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

@@ -105,6 +105,7 @@ 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
@@ -347,6 +348,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
@@ -566,6 +568,17 @@ 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".  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).
+
     max_depth : int
         The maximum depth of the heirarchy (i.e. the maximum number of bath
         exponent "excitations" to retain).
@@ -602,9 +615,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")
 
@@ -753,11 +767,9 @@ def _grad_prev_fermionic(self, he_n, k):
         ]
 
         n_excite = sum(he_fermionic_n)
-        sign1 = (-1) ** (n_excite + 1)
-
         n_excite_before_m = sum(he_fermionic_n[:k])
-        sign2 = (-1) ** (n_excite_before_m)
-
+        sign1 = (-1) ** (n_excite + 1 - self.odd_parity)
+        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["+"]:
@@ -800,10 +812,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)
-        sign1 = (-1) ** (n_excite + 1)
-
         n_excite_before_m = sum(he_fermionic_n[:k])
-        sign2 = (-1) ** (n_excite_before_m)
+        sign1 = (-1) ** (n_excite + 1 - self.odd_parity)
+        sign2 = (-1)**(n_excite_before_m + self.odd_parity)
 
         if self.ados.exponents[k].type == BathExponent.types["+"]:
             if sign1 == -1:
@@ -1233,6 +1244,7 @@ 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,
@@ -1281,6 +1293,7 @@ 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
+    tensor, Qobj, QobjEvo, fcreate, 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"
     )
 
 
@@ -232,6 +232,7 @@ 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
@@ -289,6 +290,7 @@ 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
@@ -331,6 +333,7 @@ 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
@@ -356,6 +359,7 @@ 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
@@ -1127,6 +1131,46 @@ 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")
+        resultHEOMPade = HEOMSolver(L, [bath1, bath2], depth, odd_parity=True)
+        rhoss, _ = resultHEOMPade.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()
+
+
+
 
 class TestHeomsolveFunction:
     @pytest.mark.parametrize(['evo'], [