MesoHOPS 1.4.3

1. Fixed a bug in which integer overflow was setting flux up along modes with depth 15 to 0, terminating adaptive auxiliary basis construction.
2. Added a test that confirms the full auxiliary basis is eventually constructed when the derivative error bound is sufficiently small.
As a result, the maximum hierarchy depth of 255 is now compatible with adaptivity.

Author: CaptainExasperated

mesohops/dynamics/hops_basis.py

@@ -720,13 +720,13 @@ def M2_mode_from_state(self):
     @property
     def K2_aux_by_mode(self):
         # Construct the array values of k[n] in the space of [mode, aux]
-        K2_aux_by_mode = np.zeros([self.n_hmodes, self.n_hier], dtype=np.uint8)
+        K2_aux_by_mode = np.zeros([self.n_hmodes, self.n_hier], dtype=np.uint16)
         for aux in self.hierarchy.auxiliary_list:
             array_index = np.array([list(self.list_absindex_mode).index(mode)
                                     for (mode, value) in aux.tuple_aux_vec
                                     if mode in self.list_absindex_mode],
                                    dtype=int)
-            array_values = [np.uint8(value) for (mode, value) in aux.tuple_aux_vec
+            array_values = [np.uint16(value) for (mode, value) in aux.tuple_aux_vec
                             if mode in self.list_absindex_mode]
             K2_aux_by_mode[array_index, aux._index] = array_values
         return K2_aux_by_mode

mesohops/testing/test_hops_basis.py

@@ -1468,3 +1468,62 @@ def test_M2_mode_from_state():
     )
 
     assert np.allclose(M2_mode_from_state_test, M2_mode_from_state_known)
+
+
+def test_hier_max_adaptive():
+    """
+    Tests that auxiliaries are constructed in full for calculations with small
+    delta_A up to the maximum hierarchy depth of 255 using a monomer model.
+    """
+    nsite = 1
+    e_lambda = 300
+    gamma = 60
+    temp = 300
+    hs = np.zeros([1, 1], np.complex128)
+    (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0, temp)
+    gw_sysbath = [[g_0, w_0]]
+    lop_list = np.ones([1, 1, 1])
+
+    sys_param = {
+        "HAMILTONIAN": hs,
+        "GW_SYSBATH": gw_sysbath,
+        "L_HIER": lop_list,
+        "L_NOISE1": lop_list,
+        "ALPHA_NOISE1": bcf_exp,
+        "PARAM_NOISE1": gw_sysbath,
+    }
+
+    eom_param = {"EQUATION_OF_MOTION": "NORMALIZED NONLINEAR"}
+
+    noise_param = {
+        "SEED": None,
+        "MODEL": "ZERO",
+        "TLEN": 3000.0,
+        "TAU": 0.5,
+    }
+
+    integrator_param = {
+        "INTEGRATOR": "RUNGE_KUTTA",
+        'EARLY_ADAPTIVE_INTEGRATOR': 'INCH_WORM',
+        'EARLY_INTEGRATOR_STEPS': 5,
+        'INCHWORM_CAP': 5,
+        'STATIC_BASIS': None
+    }
+
+    psi_0 = np.array([1.0], dtype=np.complex128)
+    psi_0 = psi_0 / np.linalg.norm(psi_0)
+
+    hops_ad = HOPS(
+        sys_param,
+        noise_param=noise_param,
+        hierarchy_param={"MAXHIER": 255},
+        eom_param=eom_param,
+        integration_param=integrator_param,
+    )
+    # Set delta_a to floating point minimum and run until all auxiliaries can be populated
+    hops_ad.make_adaptive(1e-323, 0, 1)
+    hops_ad.initialize(psi_0)
+    assert hops_ad.basis.n_hier < 256
+
+    hops_ad.propagate(300, 1.0)
+    assert hops_ad.basis.n_hier == 256