code review suggestions

qutip/solver/heom/bofin_solvers.py

@@ -105,7 +105,6 @@ class HierarchyADOs:
     labels: list of tuples
         A list of the ADO labels within the hierarchy.
     """
-
     def __init__(self, exponents, max_depth):
         self.exponents = exponents
         self.max_depth = max_depth
@@ -569,18 +568,19 @@ class HEOMSolver(Solver):
         or bosonic baths.
 
     odd_parity : Bool
-        For fermionic baths only.  Default is "False".  If set to "True",
-        the RHS construction differs slightly (it implies the RHS is acting on
-        an initial density operator which has odd parity in terms of its
-        representation in fermionic annhilation operators). In other words, a
-        state with a well defined number of excitations has even parity, but
-        one with superpositions of the number of excitations has odd parity.
-        This is normally prohibited in 'normal' electronic systems, where
-        charge is conserved, but is useful for the construction of
-        electron spectrum (under linear response assumptions).
+        For fermionic baths only. Default is "False". "Parity" refers to the
+        parity of the initial system state used with the HEOM. An example of
+        an odd parity state is one made from applying an odd number of
+        fermionic creation operators to a physical density operator.
+        Physical systems have even parity, but allowing the generalization
+        to odd-parity states allows one to calculate useful physical quantities
+        like the system power spectrum or density of states.
+        The form of the HEOM differs depending on the parity of the initial
+        system state, so if this option is set to "True", a different RHS is
+        constructed, which can then be used with a system state of odd parity.
 
     max_depth : int
-        The maximum depth of the heirarchy (i.e. the maximum number of bath
+        The maximum depth of the hierarchy (i.e. the maximum number of bath
         exponent "excitations" to retain).
 
     options : dict, optional
@@ -767,9 +767,11 @@ def _grad_prev_fermionic(self, he_n, k):
         ]
 
         n_excite = sum(he_fermionic_n)
-        n_excite_before_m = sum(he_fermionic_n[:k])
         sign1 = (-1) ** (n_excite + 1 - self.odd_parity)
+
+        n_excite_before_m = sum(he_fermionic_n[:k])
         sign2 = (-1)**(n_excite_before_m + self.odd_parity)
+
         sigma_bar_k = k + self.ados.sigma_bar_k_offset[k]
 
         if self.ados.exponents[k].type == BathExponent.types["+"]:
@@ -812,8 +814,9 @@ def _grad_next_fermionic(self, he_n, k):
             for i, exp in zip(he_n, self.ados.exponents)
         ]
         n_excite = sum(he_fermionic_n)
-        n_excite_before_m = sum(he_fermionic_n[:k])
         sign1 = (-1) ** (n_excite + 1 - self.odd_parity)
+
+        n_excite_before_m = sum(he_fermionic_n[:k])
         sign2 = (-1)**(n_excite_before_m + self.odd_parity)
 
         if self.ados.exponents[k].type == BathExponent.types["+"]:
@@ -1244,7 +1247,6 @@ class HSolverDL(HEOMSolver):
         operator). See :meth:`BosonicBath.combine` for details.
         Keyword only. Default: True.
     """
-
     def __init__(
         self, H_sys, coup_op, coup_strength, temperature,
         N_cut, N_exp, cut_freq, *, bnd_cut_approx=False, options=None,
@@ -1293,7 +1295,6 @@ class _GatherHEOMRHS:
         nhe : int
             The number of ADOs in the hierarchy.
     """
-
     def __init__(self, f_idx, block, nhe):
         self._block_size = block
         self._n_blocks = nhe

qutip/tests/solver/heom/test_bofin_solvers.py

@@ -9,7 +9,7 @@
 
 from qutip import (
     basis, destroy, expect, liouvillian, qeye, sigmax, sigmaz,
-    tensor, Qobj, QobjEvo, fcreate, fdestroy
+    tensor, Qobj, QobjEvo, fdestroy
 )
 from qutip.core import data as _data
 from qutip.solver.heom.bofin_baths import (
@@ -232,7 +232,6 @@ class DrudeLorentzPureDephasingModel:
     """ Analytic Drude-Lorentz pure-dephasing model for testing the HEOM
         solver.
     """
-
     def __init__(self, lam, gamma, T, Nk):
         self.lam = lam
         self.gamma = gamma
@@ -290,7 +289,6 @@ class UnderdampedPureDephasingModel:
     """ Analytic Drude-Lorentz pure-dephasing model for testing the HEOM
         solver.
     """
-
     def __init__(self, lam,  gamma, w0, T, Nk):
         self.lam = lam
         self.gamma = gamma
@@ -333,7 +331,6 @@ class BosonicMode:
     """ A description of a bosonic mode for inclusion in a
         DiscreteLevelCurrentModel.
     """
-
     def __init__(self, N, Lambda, Omega, gamma_b):
         self.N = N
         self.Lambda = Lambda
@@ -359,7 +356,6 @@ class DiscreteLevelCurrentModel:
     """ Analytic discrete level current model for testing the HEOM solver
         with a fermionic bath (and optionally a bosonic mode).
     """
-
     def __init__(self, gamma, W, T, lmax, theta=2., e1=1., bosonic_mode=None):
         # single fermion
         self.e1 = e1  # energy
@@ -1158,20 +1154,17 @@ def test_parity(self):
         bath2 = LorentzianPadeBath(
             Q=d_2, gamma=2 * Gamma, w=W, mu=mu, T=1 / beta, Nk=Nk,
             tag="Lead 2")
-        resultHEOMPade = HEOMSolver(L, [bath1, bath2], depth, odd_parity=True)
-        rhoss, _ = resultHEOMPade.steady_state(use_mkl=False)
+        solver = HEOMSolver(L, [bath1, bath2], depth, odd_parity=True)
+        rhoss, _ = solver.steady_state(use_mkl=False)
         rhoss = rhoss.full()
         expected_odd = np.diag([-0.18472, 0.68472, 0.68472, -0.18472])
-        expected = np.diag([0.10623, 0.39376, 0.39376, 0.10623])
-        assert np.isclose(rhoss, expected_odd,atol=1e-5).all()
-        resultHEOMPade = HEOMSolver(L, [bath1, bath2], depth, odd_parity=False)
-        rhoss, _ = resultHEOMPade.steady_state(use_mkl=False)
-        rhoss = rhoss.full()
-        assert np.isclose(rhoss, expected,atol=1e-5).all()
-
-
-
+        np.testing.assert_allclose(rhoss, expected_odd, atol=1e-5)
 
+        solver = HEOMSolver(L, [bath1, bath2], depth, odd_parity=False)
+        rhoss, _ = solver.steady_state(use_mkl=False)
+        rhoss = rhoss.full()
+        expected = np.diag([0.10623, 0.39376, 0.39376, 0.10623])
+        np.testing.assert_allclose(rhoss, expected, atol=1e-5)
 class TestHeomsolveFunction:
     @pytest.mark.parametrize(['evo'], [
         pytest.param("qobj", id="qobj"),