HEOM Parity support

doc/changes/2261.feature

@@ -0,0 +1 @@
+It adds odd parity support to HEOM's fermionic solver

qutip/solver/heom/bofin_solvers.py

@@ -347,6 +347,7 @@ class HierarchyADOsState:
     See :class:`HierarchyADOs` for a full list of the available attributes
     and methods.
     """
+
     def __init__(self, rho, ados, ado_state):
         self.rho = rho
         self._ado_state = ado_state
@@ -574,8 +575,20 @@ class HEOMSolver(Solver):
         If multiple baths are given, they must all be either fermionic
         or bosonic baths.
 
+    odd_parity : Bool
+        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
@@ -610,9 +623,10 @@ class HEOMSolver(Solver):
         "state_data_type": "dense",
     }
 
-    def __init__(self, H, bath, max_depth, *, options=None):
+    def __init__(self, H, bath, max_depth, *, odd_parity=False, options=None):
         _time_start = time()
-
+        # we call bool here because odd_parity will be used in arithmetic
+        self.odd_parity = bool(odd_parity)
         if not isinstance(H, (Qobj, QobjEvo)):
             raise TypeError("The Hamiltonian (H) must be a Qobj or QobjEvo")
 
@@ -761,10 +775,10 @@ def _grad_prev_fermionic(self, he_n, k):
         ]
 
         n_excite = sum(he_fermionic_n)
-        sign1 = (-1) ** (n_excite + 1)
+        sign1 = (-1) ** (n_excite + 1 - self.odd_parity)
 
         n_excite_before_m = sum(he_fermionic_n[:k])
-        sign2 = (-1) ** (n_excite_before_m)
+        sign2 = (-1)**(n_excite_before_m + self.odd_parity)
 
         sigma_bar_k = k + self.ados.sigma_bar_k_offset[k]
 
@@ -808,10 +822,10 @@ def _grad_next_fermionic(self, he_n, k):
             for i, exp in zip(he_n, self.ados.exponents)
         ]
         n_excite = sum(he_fermionic_n)
-        sign1 = (-1) ** (n_excite + 1)
+        sign1 = (-1) ** (n_excite + 1 - self.odd_parity)
 
         n_excite_before_m = sum(he_fermionic_n[:k])
-        sign2 = (-1) ** (n_excite_before_m)
+        sign2 = (-1)**(n_excite_before_m + self.odd_parity)
 
         if self.ados.exponents[k].type == BathExponent.types["+"]:
             if sign1 == -1:

qutip/tests/solver/heom/test_bofin_solvers.py

@@ -9,7 +9,7 @@
 
 from qutip import (
     basis, destroy, expect, liouvillian, qeye, sigmax, sigmaz, sigmay,
-    tensor, Qobj, QobjEvo, fidelity
+    tensor, Qobj, QobjEvo, fidelity, fdestroy
 )
 from qutip.core import data as _data
 from qutip.solver.heom.bofin_baths import (
@@ -54,8 +54,8 @@ def assert_raises_steady_state_time_dependent(hsolver):
     with pytest.raises(ValueError) as err:
         hsolver.steady_state()
     assert str(err.value) == (
-       "A steady state cannot be determined for a time-dependent"
-       " system"
+        "A steady state cannot be determined for a time-dependent"
+        " system"
     )
 
 
@@ -1145,7 +1145,44 @@ def test_solving_with_step(self):
 
         assert states[-1] == ado_state.extract(0)
 
-
+    def test_parity(self):
+        depth = 2
+        Nk = 2
+        # system: two fermions
+        N = 2
+        d_1 = fdestroy(N, 0)
+        d_2 = fdestroy(N, 1)
+        # bath params:
+        mu = 0.  # chemical potential
+        Gamma = 1  # coupling strenght
+        W = 2.5  # bath width
+        # system params:
+        # coulomb repulsion
+        U = 3 * np.pi * Gamma
+        # impurity energy
+        w0 = - U / 2.
+        beta = 1 / (0.2 * Gamma)  # Inverse Temperature
+        H = w0 * (d_1.dag() * d_1 + d_2.dag()
+                  * d_2) + U * d_1.dag() * d_1 * d_2.dag() * d_2
+
+        L = liouvillian(H)
+        bath1 = LorentzianPadeBath(
+            Q=d_1, gamma=2 * Gamma, w=W, mu=mu, T=1 / beta, Nk=Nk,
+            tag="Lead 1")
+        bath2 = LorentzianPadeBath(
+            Q=d_2, gamma=2 * Gamma, w=W, mu=mu, T=1 / beta, Nk=Nk,
+            tag="Lead 2")
+        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])
+        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"),