diff --git a/src/simsopt/geo/boozersurface.py b/src/simsopt/geo/boozersurface.py
index 32a447023..86c29b1fb 100644
--- a/src/simsopt/geo/boozersurface.py
+++ b/src/simsopt/geo/boozersurface.py
@@ -85,8 +85,6 @@ def __init__(self, biotsavart, surface, label, targetlabel, constraint_weight=No
         # set the default options now
         if 'verbose' not in options:
             options['verbose'] = True
-        if 'G' not in options:
-            options['G'] = None
         
         # default solver options for the BoozerExact and BoozerLS solvers
         if self.boozer_type == 'exact':
@@ -120,7 +118,7 @@ def run_code(self, iota, G=None):
         Args:
             boozer_type: either 'exact' or 'ls', to indicate whether a BoozerExact or a BoozerLS surface is to be computed.
             iota: guess for value of rotational transform on the surface,
-            G: guess for value of G on surface, defaults to None.
+            G: guess for value of G on surface, defaults to None.  Note that if None is used, then the coil currents must be fixed.
             options: dictionary of solver options. If keyword is not specified, then a default 
                     value is used.  Possible keywords are:
                     `verbose`: display convergence information
@@ -134,6 +132,10 @@ def run_code(self, iota, G=None):
         if not self.need_to_run_code:
             return
         
+        # for coil optimizations, the gradient calculations of the objective assume that the coil currents are fixed.
+        if G is None:
+            assert np.all([c.current.is_fixed(0) for c in self.biotsavart.coils])
+
         # BoozerExact default solver
         if self.boozer_type == 'exact':
             res = self.solve_residual_equation_exactly_newton(iota=iota, G=G, tol=self.options['newton_tol'], maxiter=self.options['newton_maxiter'], verbose=self.options['verbose'])
diff --git a/tests/geo/surface_test_helpers.py b/tests/geo/surface_test_helpers.py
index 3e9490ba6..12d9fcfdd 100644
--- a/tests/geo/surface_test_helpers.py
+++ b/tests/geo/surface_test_helpers.py
@@ -90,6 +90,11 @@ def get_boozer_surface(label="Volume", nphi=None, ntheta=None, boozer_type='exac
     current_sum = sum(abs(c.current.get_value()) for c in coils)
     G0 = 2. * np.pi * current_sum * (4 * np.pi * 10**(-7) / (2 * np.pi)) if optimize_G else None
 
+    # currents need to be fixed if optimize_G is None
+    if optimize_G is False:
+        for c in bs.coils:
+            c.current.fix_all()
+
     ## RESOLUTION DETAILS OF SURFACE ON WHICH WE OPTIMIZE FOR QA
     mpol = 6 if boozer_type == 'exact' else 3
     ntor = 6 if boozer_type == 'exact' else 3