MesoHops v1.3

This commit is a major update over the previous version of MesoHOPS including:
1. The implementation of an adaptive algorithm that works for arbitrary diagonal L-operators. In practice, this means the adHOPS can now be used for multi-particle dynamics!
2. A substantial speed improvement (3-8x faster for large systems)
3. A modest further improvement in memory when managing a large number of auxiliary wave functions
4. Some behind the scenes simplifications in the code design that will allow enable future generalizations
5. The noise class has been condensed and now has some additional options for how the uncorrelated noise is initially generated

Starting with this commit, we are implementing levels of testing. For most purposes, pytest -level 1 is sufficient to ensure the code is behaving as expected. pytest -level 2 invokes a much slower set of tests that ensures the random number generators are performing within expected tolerances.

Author: digbennett

mesohops/dynamics/basis_functions_adaptive.py

@@ -181,7 +181,7 @@ def error_flux_up(Φ, n_state, n_hier, n_hmodes, list_w, K2_aux_bymode,
     if (type == "H"):
         # \sum_s |L_m[s,s] \psi_k[s]|^2 = \sum_s |M2[m,s] * \psi_k[s]|^2
         # \sum_s |M2[m,s]|^2 * |\psi_k[s]|^2
-        P2_pop_modes = (np.abs(M2_mode_from_state)**2 @ (np.abs(C2_phi) ** 2))
+        P2_pop_modes = (np.abs(M2_mode_from_state).power(2) @ (np.abs(C2_phi) ** 2))
         E2_error = np.abs(W2_bymode) ** 2 * (1 + K2_aux_bymode) ** 2 * P2_pop_modes / hbar ** 2
         if F2_filter is not None:
             E2_error = F2_filter*E2_error
@@ -191,7 +191,7 @@ def error_flux_up(Φ, n_state, n_hier, n_hmodes, list_w, K2_aux_bymode,
         E2_error = np.abs(np.abs(W2_bymode) * (1 + K2_aux_bymode)) ** 2
         if F2_filter is not None:
             E2_error *= F2_filter
-        E2_error = np.transpose(M2_mode_from_state) @ E2_error
+        E2_error = np.transpose(M2_mode_from_state.power(2)) @ E2_error
         E2_error = P2_pop_state * E2_error / hbar ** 2
 
     else:
@@ -263,7 +263,7 @@ def error_flux_down(Φ, n_state, n_hier, n_hmodes, list_g, list_w, M2_mode_from_
         # Hierarchy Type Downward Flux
         # ============================
         D2_mode_from_state = np.zeros([n_hmodes,n_state])
-        D2_mode_from_state[:,:] = M2_mode_from_state - E1_Lm.reshape([len(E1_Lm),1])
+        D2_mode_from_state[:,:] = M2_mode_from_state.toarray() - E1_Lm.reshape([len(E1_Lm),1])
         E2_flux_down_error = (
                 np.real(
                     (np.abs(G2_bymode / W2_bymode) ** 2)
@@ -278,7 +278,7 @@ def error_flux_down(Φ, n_state, n_hier, n_hmodes, list_g, list_w, M2_mode_from_
         # State Type Downward Flux
         # ========================
         D2_state_from_mode = np.zeros([n_state,n_hmodes])
-        D2_state_from_mode[:,:] = np.transpose(M2_mode_from_state) - E1_Lm
+        D2_state_from_mode[:,:] = np.transpose(M2_mode_from_state).toarray() - E1_Lm
         E2_flux_down_error = np.abs(G2_bymode / W2_bymode)**2
         if F2_filter is not None:
             E2_flux_down_error *= F2_filter

mesohops/dynamics/eom_functions.py

@@ -3,7 +3,7 @@
 
 
 __title__ = "EOM Functions"
-__author__ = "D. I. G. Bennett"
+__author__ = "D. I. G. Bennett, J. K. Lynd"
 __version__ = "1.2"
 
 
@@ -191,3 +191,99 @@ def calc_norm_corr(
         delta += (np.conj(phi_0) @ phi_1) * list_avg_L2[l_ind]
 
     return np.real(delta)
+
+def calc_LT_corr(
+    list_LT_coeff, list_L2, list_avg_L2, physical = False
+):
+    """
+    This function computes the low-temperature correction factor associated with each
+    member of the hierarchy in the nonlinear equation of motion. The factor is given by
+    the sum over the low-temperature correction coefficients and associated
+    L-operators c_n and L_n:
+    \sum_n (2<L_n>Re[c_n] - L_nc_n)L_n
+    Where c_n is the nth low-temperature correction factor, and L_n is the nth
+    L-operator associated with that factor.
+
+    PARAMETERS
+    ----------
+    1. list_LT_coeff : list
+                       relative list of low-temperature coefficients
+    2. list_L2 : list
+                 relative list of active L operators
+    3. list_avg_L2 : list
+                     relative list of the expectation values of the active L operators
+    4. physical : boolean
+                  Determines if this correction is to be applied to the physical
+                  wavefunction or elsewhere
+    RETURNS
+    -------
+    1. C2_corr : np.array
+                 the low-temperature correction self-derivative term
+    """
+    # Adds the terminated flux from above to the physical wavefunction
+    if physical:
+        return np.sum(np.array([(list_LT_coeff[n]*list_avg_L2[n]*np.eye(
+        list_L2[n].shape[0]) - list_LT_coeff[n]*list_L2[n])@list_L2[n] for n in
+                            range(len(list_LT_coeff))]), axis=0)
+    # Adds the delta-approximated noise memory drift to all members of the hierarchy
+    else:
+        return np.sum(np.array([(np.conj(list_LT_coeff[n])*list_avg_L2[n]*np.eye(
+        list_L2[n].shape[0]))@list_L2[n] for n in range(len(list_LT_coeff))]), axis=0)
+
+def calc_LT_corr_to_norm_corr(
+    list_LT_coeff, list_avg_L2, list_avg_L2_sq
+):
+    """
+    This function computes the low-temperature correction to the normalization factor in
+    the normalized nonlinear equation of motion. The correction is given by the sum
+    over the low-temperature correction coefficients and associated L-operators c_n
+    and L_n:
+    \sum_n Re[c_n](2<L_n>^2 - <L_n^2>)
+    Where c_n is the nth low-temperature correction factor, and L_n is the nth
+    L-operator associated with that factor.
+
+    PARAMETERS
+    ----------
+    1. list_LT_coeff : list
+                       relative list of low-temperature coefficients
+    2. list_avg_L2 : list
+                     relative list of the expectation values of the L operators
+    3. list_avg_L2_sq : list
+                        relative list of the expectation values of the squared L
+                       operators
+    RETURNS
+    -------
+    1. delta_corr : float
+                    the low-temperature correction to the normalization correction
+                    factor
+    """
+    return np.sum(np.array([
+    np.real(list_LT_coeff[j])*(2*list_avg_L2[j]**2 - list_avg_L2_sq[j])
+        for j in range(len(list_LT_coeff))]))
+
+def calc_LT_corr_linear(
+    list_LT_coeff, list_L2
+):
+    """
+    This function computes the low-temperature correction factor associated with each
+    member of the hierarchy in the linear equation of motion. The factor is given by
+    the sum over the low-temperature correction coefficients and associated
+    L-operators c_n and L_n:
+    -\sum_n c_nL_n^2
+    Where c_n is the nth low-temperature correction factor, and L_n is the nth
+    L-operator associated with that factor.
+    NOTE: this correction should only be applied to the physical wavefunction.
+
+    PARAMETERS
+    ----------
+    1. list_LT_coeff : list
+                       relative list of low-temperature coefficients
+    2. list_L2 : list
+                 relative list of L operators
+    RETURNS
+    -------
+    1. C2_corr : np.array
+                 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)