MesoHOPS v1.3.1

This update provides a number of quality-of-life improvements (e.g. a larger variety of ways to define the noise trajectory) in the code base and fixes one problem with memory scaling introduced by CPU efficient (but linear scaling memory) algorithm for indexing auxiliary wave functions. The code is also more robust against some errors that were arising in edge-cases.

Behind the scenes, a number of code updates have provided a more consistent set of naming conventions.

Author: Doran Raccah

mesohops/dynamics/basis_functions_adaptive.py

@@ -196,7 +196,7 @@ def error_flux_up(Φ, n_state, n_hier, n_hmodes, list_w, K2_aux_bymode,
 
     else:
         E2_error = 0
-        raise UnsupportedRequest(type, f"in error_flux_up received type {type}")
+        raise UnsupportedRequest("type {}".format(type), "error_flux_up")
 
     return E2_error
 
@@ -288,7 +288,7 @@ def error_flux_down(Φ, n_state, n_hier, n_hmodes, list_g, list_w, M2_mode_from_
         )
     else:
         E2_flux_down_error = 0
-        raise UnsupportedRequest(type, f"in error_flux_down received type {type}")
+        raise UnsupportedRequest("type {}".format(type), "error_flux_down")
 
     return E2_flux_down_error
 

mesohops/dynamics/bath_corr_functions.py

@@ -72,6 +72,59 @@ def bcf_convert_sdl_to_exp(lambda_sdl, gamma_sdl, omega_sdl, temp):
 
     return (g_exp, w_exp)
 
+def bcf_convert_dl_ud_to_exp(lambda_dl, gamma_dl, omega_dl, temp):
+    """
+    Converts underdamped Drude-Lorentz spectral density parameters to the exponential
+    equivalent. Assumes that omega_dl (the underdamped frequency) is larger than
+    gamma_dl (the reorganization timescale). Does not account for Matsubara modes.
+
+    Parameters
+    ----------
+    1. lambda_sdl : float [unit: cm^-1]
+                    Reorganization energy.
+
+    2. gamma_sdl : float [unit: cm^-1]
+                   Reorganization time scale.
+
+    3. omega_sdl : float [unit: cm^-1]
+                   Vibrational frequency.
+
+    4. temp : float [unit: K]
+              Temperature.
+
+    Returns
+    -------
+    1. list_modes: list(complex)
+                   List of the exponential modes that comprise the correlation
+                   function, alternating gs and ws (complex, [cm^-2] and [cm^-1],
+                   representing the constant prefactor and exponential decay rate,
+                   respectively)
+    """
+    beta = 1 / (kB * temp)
+    xi = np.sqrt(omega_dl**2 - gamma_dl**2)
+    prefactor_base = (lambda_dl * omega_dl**2)/(2 * xi)
+    w_1 = xi + 1j*gamma_dl
+    w_2 = -1*xi + 1j*gamma_dl
+    g_1 = prefactor_base
+    g_2 = -prefactor_base
+    coth_w_1 = 1 / np.tanh(w_1*beta/2)
+    coth_w_2 = 1 / np.tanh(w_2*beta/2)
+    g_1 += (coth_w_1 - np.conj(coth_w_2)) * prefactor_base
+    g_2 += (-coth_w_2 + np.conj(coth_w_1)) * prefactor_base
+
+    # # Test to prove that this expression is the same as the high-temperature
+    # # approximation of equation S52 from Bennet et al., Supplementary Information:
+    # # Mechanistic regimes of vibronic transport in a heterodimer and the design
+    # # principle of incoherent vibronic transport in phycobiliproteins, J. Phys. Chem.
+    # # Lett., 2018, https://doi.org/10.1021/acs.jpclett.8b00844:
+    # t_axis = np.arange(0, 0.21, 0.01)
+    # exp_form = g_1*np.exp(1j*w_1*t_axis) + g_2*np.exp(1j*w_2*t_axis)
+    # analytic_form = (lambda_dl*(gamma_dl**2 + xi**2)/xi) * np.exp(-1*gamma_dl*t_axis)\
+    #                 * (2*(np.sin(beta*gamma_dl)*np.sin(xi*t_axis) + np.sinh(
+    #     beta*xi)*np.cos(xi*t_axis))/(np.cosh(beta*xi)-np.cos(beta*gamma_dl)) +
+    #                    1j*np.sin(xi*t_axis))
+    return [g_1, -1j*w_1, g_2, -1j*w_2]
+
 
 def bcf_convert_dl_to_exp_with_Matsubara(lambda_dl, gamma_dl, temp, k_matsubara):
     """
@@ -189,4 +242,4 @@ def E_tilde_k(w):
     w_exp_conj = np.conj(w_exp)
 
     return [g_exp - 1j*g_exp_im, w_exp, g_exp_conj - 1j*g_exp_im_conj, w_exp_conj] + \
-           list_mats_modes
+           list_mats_modes

mesohops/dynamics/eom_functions.py

@@ -9,7 +9,7 @@
 
 def operator_expectation(oper, vec, flag_gcorr = False):
     """
-    Calculates the expectation value of a sparse operator
+    Calculates the expectation value of a sparse operator.
 
     Parameters
     ----------
@@ -63,10 +63,10 @@ def compress_zmem(z_mem, list_index_L2_by_mode, list_absindex_mode):
     return dz_hat
 
 
-def calc_delta_zmem(z_mem, list_avg_L2, list_g, list_w, list_index_L2_by_mode,
-                    list_absindex_mode, list_index_L2_active):
+def calc_delta_zmem(z_mem, list_avg_L2, list_g, list_w, list_absindex_L2_by_mode,
+                    list_absindex_mode, list_absindex_L2_active):
     """
-    This updates the memory term. The form of the equation depends on expanding the
+    Updates the memory term. The form of the equation depends on expanding the
     memory integral assuming an exponential expansion.
 
     NOTE: THIS ASSUMES THE NOISE HAS EXPONENTIAL FORM.
@@ -77,9 +77,9 @@ def calc_delta_zmem(z_mem, list_avg_L2, list_g, list_w, list_index_L2_by_mode,
           list_avg_L2 : relative
           list_g : absolute
           list_w : absolute
-          list_index_L2_by_mode : relative
+          list_absindex_L2_by_mode : absolute
           list_absindex_mode : mapping from relative-->absolute
-          list_index_L2_active : relative
+          list_absindex_L2_active : absolute
           
     Parameters
     ----------
@@ -90,23 +90,24 @@ def calc_delta_zmem(z_mem, list_avg_L2, list_g, list_w, list_index_L2_by_mode,
                      Relative list of the expectation values of the L operators.
 
     3. list_g : list(complex)
-                List of pre exponential factors for bath correlation functions [
-                absolute].
+                List of pre exponential factors for bath correlation functions [unit:
+                cm^-2].
 
     4. list_w :  list(complex)
-                 List of exponents for bath correlation functions (w = γ+iΩ) [absolute].
+                 List of exponents for bath correlation functions (w = γ+iΩ) [unit:
+                 cm^-1].
 
-    5. list_index_L2_by_mode : list(int)
-                               List of length equal to the number of 'modes' in the
-                               current hierarchy basis and each entry is an index for
-                               the absolute list_L2.
+    5. list_absindex_L2_by_mode : list(int)
+                                  List of indices for the absolute list of L-operators
+                                  to match L-operators to the associated absolute mode
+                                  index.
 
     6. list_absindex_mode : list(int)
                             List of the absolute indices  of the modes in current basis.
 
-    7. list_index_L2_active : list(int)
-                              List of relative indices of L-operators that have any
-                              non-zero values.
+    7. list_absindex_L2_active : list(int)
+                                 List of absolute indices of L-operators that have any 
+                                 non-zero values.
 
     Returns
     -------
@@ -121,11 +122,12 @@ def calc_delta_zmem(z_mem, list_avg_L2, list_g, list_w, list_index_L2_by_mode,
     )
 
     # Loop over modes in the current basis
-    for (relindex_mode, absindex_mode) in enumerate(list_absindex_mode):
-        index_L2 = list_index_L2_by_mode[relindex_mode]
-        if index_L2 in list_index_L2_active:
-            l_avg = list_avg_L2[list_index_L2_active.index(index_L2)]
-        else:
+    for absindex_mode in list_absindex_mode:
+        absindex_L2 = list_absindex_L2_by_mode[absindex_mode]
+        try:
+            relindex_L2 = list(list_absindex_L2_active).index(absindex_L2)
+            l_avg = list_avg_L2[relindex_L2]
+        except:
             l_avg = 0
         delta_z_mem[absindex_mode] = (
             l_avg * np.conj(list_g[absindex_mode])
@@ -286,4 +288,4 @@ def calc_LT_corr_linear(
                  the low-temperature correction self-derivative term
     """
     return -1*np.sum(np.array([(list_LT_coeff[n])*list_L2[n]@list_L2[n] for n in
-                            range(len(list_LT_coeff))]),axis=0)
+                            range(len(list_LT_coeff))]),axis=0)

mesohops/dynamics/eom_hops_ksuper.py

@@ -284,12 +284,12 @@ def _add_crossterms_stable_K(
         for l_mode_abs in list_new_mode:
             #For each new mode, we look through our dictionary of mode connections between stable auxiliaries and add the corresponding entry to K
             #super operators.
-            if(len(hierarchy.stable_aux_hash_conn_by_mode[l_mode_abs]) > 0):
-                list_hashes, list_aux_conn_hash_and_value = list(zip(*hierarchy.stable_aux_hash_conn_by_mode[
+            if(len(hierarchy.stable_aux_id_conn_by_mode[l_mode_abs]) > 0):
+                list_ids, list_aux_conn_id_and_value = list(zip(*hierarchy.stable_aux_id_conn_by_mode[
                     l_mode_abs].items()))
-                list_aux_indices = [hierarchy._aux_index(hierarchy._aux_by_hash[hash_]) for hash_ in list_hashes]
-                list_hashes_p1, list_aux_mode_values = list(zip(*list_aux_conn_hash_and_value))
-                list_aux_indices_p1 = [hierarchy._aux_index(hierarchy._aux_by_hash[hash_]) for hash_ in list_hashes_p1]
+                list_aux_indices = [hierarchy._aux_index(hierarchy.dict_aux_by_id[id_]) for id_ in list_ids]
+                list_ids_p1, list_aux_mode_values = list(zip(*list_aux_conn_id_and_value))
+                list_aux_indices_p1 = [hierarchy._aux_index(hierarchy.dict_aux_by_id[id_]) for id_ in list_ids_p1]
                 l_mod = list(mode.list_absindex_mode).index(l_mode_abs)
                 i_lop = mode.list_index_L2_by_hmode[l_mod]
                 Kp1_data.extend(