diff --git a/f90wrap_helpers/classes.py b/f90wrap_helpers/classes.py
index 5cd2e16..485464d 100644
--- a/f90wrap_helpers/classes.py
+++ b/f90wrap_helpers/classes.py
@@ -42,10 +42,14 @@ def set_linear_power_spectrum(self, k, z, pofk_lin):
         self.k_lin = np.asfortranarray(k, dtype=np.float64)
 
         pyhmcode.cosmology_functions.init_external_linear_power_tables(self, self.k_lin, self.a_lin, pk_lin)
-        pyhmcode.cosmology_functions.init_cosmology(self)
+        self._update()
 
         self._has_linear_power_set = True
 
+    def _update(self):
+        pyhmcode.cosmology_functions.init_cosmology(self)
+
+
 class Halomodel(pyhmcode.hmx.halomod):
     def __new__(cls, mode, verbose=False):
         # Need to use __new__ to intercept creation of the instance and use
@@ -70,8 +74,6 @@ def As(self, As):
         _pyhmcode.f90wrap_halomod__set__as(self._handle, As)
 
 
-
-
 def calculate_nonlinear_power_spectrum(cosmology, halomodel, fields=None,
                                        return_halo_terms=False, verbose=False):
     """Compute power spectra using HMCode or HMx.
@@ -102,6 +104,9 @@ def calculate_nonlinear_power_spectrum(cosmology, halomodel, fields=None,
     if not cosmology._has_linear_power_set:
         raise RuntimeError("Cosmology has no linear power spectrum set.")
 
+    # Update Fortran cosmology internal state 
+    cosmology._update()
+
     if halomodel.hmcode_version in [HMcode2015, HMcode2016, HMcode2020, HMcode2020_feedback]:
         fields = np.array([field_dmonly])
     elif fields is None: