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

@@ -566,6 +566,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 +613,9 @@ 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()
-
+        self.odd_parity = odd_parity
         if not isinstance(H, (Qobj, QobjEvo)):
             raise TypeError("The Hamiltonian (H) must be a Qobj or QobjEvo")
 
@@ -745,6 +756,16 @@ def _grad_prev_bosonic(self, he_n, k):
             )
         return op
 
+    def _get_signs(self, n_excite, n_excite_before_m):
+        """gets the correct signs according to the parity"""
+        if self.odd_parity:
+            sign1 = (-1) ** (n_excite)
+            sign2 = (-1) ** (n_excite_before_m + 1)
+        else:
+            sign1 = (-1) ** (n_excite + 1)
+            sign2 = (-1) ** (n_excite_before_m)
+        return sign1, sign2
+
     def _grad_prev_fermionic(self, he_n, k):
         ck = self.ados.ck
         he_fermionic_n = [
@@ -753,11 +774,10 @@ 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, sign2 = self._get_signs(
+            n_excite=n_excite, n_excite_before_m=n_excite_before_m
+            )
         sigma_bar_k = k + self.ados.sigma_bar_k_offset[k]
 
         if self.ados.exponents[k].type == BathExponent.types["+"]:
@@ -800,11 +820,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, sign2 = self._get_signs(
+            n_excite=n_excite, n_excite_before_m=n_excite_before_m)
         if self.ados.exponents[k].type == BathExponent.types["+"]:
             if sign1 == -1:
                 op = _data.mul(self._s_pre_minus_post_Q[k], -1j * sign2)