diff --git a/dyson/expressions/ccsd.py b/dyson/expressions/ccsd.py
index 4af07bd..bda91a0 100644
--- a/dyson/expressions/ccsd.py
+++ b/dyson/expressions/ccsd.py
@@ -3,8 +3,7 @@
 """
 
 import numpy as np
-from pyscf import ao2mo, cc, lib
-
+from pyscf import ao2mo, cc, lib, scf, pbc
 from dyson import util
 from dyson.expressions import BaseExpression
 
@@ -20,22 +19,33 @@ def __init__(self, *args, ccsd=None, t1=None, t2=None, l1=None, l2=None, **kwarg
         BaseExpression.__init__(self, *args, **kwargs)
 
         if ccsd is None:
-            ccsd = cc.CCSD(self.mf, mo_coeff=self.mo_coeff, mo_occ=self.mo_occ)
-            ccsd.verbose = 0
-
-        # Update amplitudes if provided as kwargs
-        self.t1 = t1 if t1 is not None else ccsd.t1
-        self.t2 = t2 if t2 is not None else ccsd.t2
-        self.l1 = l1 if l1 is not None else ccsd.l1
-        self.l2 = l2 if l2 is not None else ccsd.l2
+            if isinstance(self.mf, scf.hf.RHF):
+                ccsd = cc.CCSD(self.mf, mo_coeff=self.mo_coeff, mo_occ=self.mo_occ)
+            elif isinstance(self.mf, pbc.scf.hf.RHF):
+                ccsd = pbc.cc.CCSD(self.mf, mo_coeff=self.mo_coeff, mo_occ=self.mo_occ)
+            else:
+                raise NotImplementedError("EOM-CCSD not implemented for this type of mean-field object.")
+
+            #ccsd = cc.CCSD(self.mf, mo_coeff=self.mo_coeff, mo_occ=self.mo_occ)
+            ccsd.t1 = t1
+            ccsd.t2 = t2
+            ccsd.l1 = l1
+            ccsd.l2 = l2 
 
         # Solve CCSD if amplitudes are not provided
         if ccsd.t1 is None or ccsd.t2 is None:
             ccsd.kernel()
             self.t1 = ccsd.t1
             self.t2 = ccsd.t2
+        else:
+            self.t1 = ccsd.t1
+            self.t2 = ccsd.t2
+        
         if ccsd.l1 is None or ccsd.l2 is None:
             self.l1, self.l2 = ccsd.solve_lambda()
+        else:
+            self.l1 = ccsd.l1
+            self.l2 = ccsd.l2
 
         self.eris = ccsd.ao2mo()
         self.imds = cc.eom_rccsd._IMDS(ccsd, eris=self.eris)
@@ -114,23 +124,34 @@ def __init__(self, *args, ccsd=None, t1=None, t2=None, l1=None, l2=None, **kwarg
         BaseExpression.__init__(self, *args, **kwargs)
 
         if ccsd is None:
-            ccsd = cc.CCSD(self.mf, mo_coeff=self.mo_coeff, mo_occ=self.mo_occ)
-            ccsd.verbose = 0
-
-        # Update amplitudes if provided as kwargs
-        self.t1 = t1 if t1 is not None else ccsd.t1
-        self.t2 = t2 if t2 is not None else ccsd.t2
-        self.l1 = l1 if l1 is not None else ccsd.l1
-        self.l2 = l2 if l2 is not None else ccsd.l2
+            if isinstance(self.mf, scf.hf.RHF):
+                ccsd = cc.CCSD(self.mf, mo_coeff=self.mo_coeff, mo_occ=self.mo_occ)
+            elif isinstance(self.mf, pbc.scf.hf.RHF):
+                ccsd = pbc.cc.CCSD(self.mf, mo_coeff=self.mo_coeff, mo_occ=self.mo_occ)
+            else:
+                raise NotImplementedError("momCCSD not implemented for this type of mean-field object.")
+            #ccsd = cc.CCSD(self.mf, mo_coeff=self.mo_coeff, mo_occ=self.mo_occ)
+            # Use provided amplitudes if available
+            ccsd.t1 = t1
+            ccsd.t2 = t2
+            ccsd.l1 = l1
+            ccsd.l2 = l2 
 
         # Solve CCSD if amplitudes are not provided
         if ccsd.t1 is None or ccsd.t2 is None:
             ccsd.kernel()
             self.t1 = ccsd.t1
             self.t2 = ccsd.t2
+        else:
+            self.t1 = ccsd.t1
+            self.t2 = ccsd.t2
+        
         if ccsd.l1 is None or ccsd.l2 is None:
             self.l1, self.l2 = ccsd.solve_lambda()
-
+        else:
+            self.l1 = ccsd.l1
+            self.l2 = ccsd.l2
+            
         self.eris = ccsd.ao2mo()
         self.imds = cc.eom_rccsd._IMDS(ccsd, eris=self.eris)
         self.imds.make_ea()
diff --git a/dyson/expressions/expression.py b/dyson/expressions/expression.py
index 9b05f3f..897368d 100644
--- a/dyson/expressions/expression.py
+++ b/dyson/expressions/expression.py
@@ -152,9 +152,10 @@ def build_gf_moments(self, nmom, store_vectors=True, left=False):
                     u = apply_hamiltonian(u)
 
         if left:
-            t = t.transpose(1, 2).conj()
+            t = t.transpose(0,2,1).conj()
 
         return t
+    
 
     def build_gf_chebyshev_moments(self, nmom, store_vectors=True, left=False, scaling=None):
         """Build moments of the Green's function using Chebyshev polynomials.
diff --git a/dyson/lehmann.py b/dyson/lehmann.py
index d20685d..2c71c1c 100644
--- a/dyson/lehmann.py
+++ b/dyson/lehmann.py
@@ -108,7 +108,7 @@ def moment(self, order):
 
         couplings_l, couplings_r = self._unpack_couplings()
 
-        moment = np.einsum(
+        moment = lib.einsum(
             "pk,qk,nk->npq",
             couplings_l,
             couplings_r.conj(),
@@ -447,21 +447,25 @@ def as_perturbed_mo_energy(self):
 
         return mo_energy
 
-    def on_grid(self, grid, eta=1e-1, ordering="time-ordered"):
+    def on_grid(self, grid, eta=1e-1, ordering="time-ordered", axis='real'):
         """
-        Return the Lehmann representation realised on a real frequency
+        Return the Lehmann representation realised on a frequency
         grid.
 
         Parameters
         ----------
         grid : numpy.ndarray
-            Array of real frequency points.
+            Array of frequency points.
         eta : float, optional
             Broadening parameter.  Default value is `1e-1`.
+            Only relevant for real axis.
         ordering : str, optional
             Time ordering.  Can be one of `{"time-ordered",
             "advanced", "retarded"}`.  Default value is
-            `"time-ordered"`.
+            `"time-ordered"`. 
+        axis : str, optional
+            Frequency axis.  Can be one of `{"real", "imag"}`.  Default
+            value is `"real"`.
 
         Returns
         -------
@@ -480,7 +484,12 @@ def on_grid(self, grid, eta=1e-1, ordering="time-ordered"):
 
         couplings_l, couplings_r = self._unpack_couplings()
 
-        denom = 1.0 / lib.direct_sum("w+k-k->wk", grid, signs * 1.0j * eta, self.energies)
+        if axis == "real":
+            denom = 1.0 / lib.direct_sum("w+k-k->wk", grid, signs * 1.0j * eta, self.energies)
+        elif axis == "imag":
+            denom = 1.0 / lib.direct_sum("w-k->wk", 1j*grid, self.energies)
+        else:
+            raise ValueError("axis = {}".format(axis))
         f = lib.einsum("pk,qk,wk->wpq", couplings_l, couplings_r.conj(), denom)
 
         return f
diff --git a/dyson/solvers/__init__.py b/dyson/solvers/__init__.py
index 916405d..3807718 100644
--- a/dyson/solvers/__init__.py
+++ b/dyson/solvers/__init__.py
@@ -6,6 +6,6 @@
 from dyson.solvers.mblgf import MBLGF, MixedMBLGF
 from dyson.solvers.kpmgf import KPMGF
 from dyson.solvers.cpgf import CPGF
-from dyson.solvers.chempot import AufbauPrinciple, AuxiliaryShift
+from dyson.solvers.chempot import AufbauPrinciple, AufbauPrincipleBisect, AuxiliaryShift
 from dyson.solvers.density import DensityRelaxation
 from dyson.solvers.self_consistent import SelfConsistent
diff --git a/dyson/solvers/chempot.py b/dyson/solvers/chempot.py
index 061f4b4..c09da5a 100644
--- a/dyson/solvers/chempot.py
+++ b/dyson/solvers/chempot.py
@@ -119,6 +119,57 @@ def get_self_energy(self):
     def get_greens_function(self):
         return self.gf.copy(chempot=self.chempot, deep=False)
 
+class AufbauPrincipleBisect(AufbauPrinciple):
+
+    def _kernel(self):
+        energies = self.gf.energies
+        couplings_l, couplings_r = self.gf._unpack_couplings()
+        weights = self.gf.weights()
+        low, high = 0, self.gf.naux
+        mid = high//2
+        for iter in range(100):           
+            sum = self.occupancy * weights[:mid].sum()
+            self.log.debug("Number of electrons [0:%d] = %.6f", iter + 1, sum) 
+            if sum < self.nelec:
+                low = mid
+                mid = mid + (high-low)//2
+            else:
+                high = mid
+                mid = mid - (high-low)//2
+
+            if low == mid or high == mid:
+                break
+
+        n_low, n_high = self.occupancy * weights[:low].sum(), self.occupancy * weights[:high].sum()
+        # assert n_low < self.nelec
+        # assert n_high > self.nelec
+        
+        if abs(n_low - self.nelec) < abs(n_high - self.nelec):
+            homo = low - 1 
+            error = self.nelec - n_low
+        else:
+            homo = high - 1
+            error = self.nelec - n_high
+
+        try:
+            lumo = homo + 1
+            chempot = 0.5 * (energies[homo] + energies[lumo])
+        except:
+            raise ValueError("Failed to find Fermi energy.")
+        self.log.info('HOMO LUMO %s %s'%(homo, lumo))
+        self.log.info("HOMO = %.6f", energies[homo])
+        self.log.info("LUMO = %.6f", energies[lumo])
+        self.log.info("Chemical potential = %.6f", chempot)
+        self.log.info("Error in nelec = %.3g", error)
+
+        self.converged = True
+        self.homo = energies[homo]
+        self.lumo = energies[lumo]
+        self.chempot = chempot
+        self.error = error
+
+        return chempot, error
+
 
 class AuxiliaryShift(BaseSolver):
     """
diff --git a/tests/solvers/test_aufbau.py b/tests/solvers/test_aufbau.py
index c5f43d9..0a7a0f2 100644
--- a/tests/solvers/test_aufbau.py
+++ b/tests/solvers/test_aufbau.py
@@ -8,7 +8,7 @@
 from pyscf import gto, scf, agf2
 import numpy as np
 
-from dyson import util, AufbauPrinciple, NullLogger
+from dyson import util, AufbauPrinciple, AufbauPrincipleBisect, NullLogger, MBLGF
 from dyson.lehmann import Lehmann
 
 
@@ -58,5 +58,25 @@ def test_agf2(self):
         self.assertAlmostEqual(solver.error, 0.017171058925, 7)
 
 
+@pytest.mark.regression
+class AufbauPrincipleBisect_Tests(unittest.TestCase):
+    def test_wrt_AufbauPrinciple(self):
+        for i in range(10):
+            n = 100
+            moms = np.random.random((16, n, n))
+            moms = moms + moms.transpose(0, 2, 1)
+            mblgf = MBLGF(moms)
+            mblgf.kernel()
+            gf = mblgf.get_greens_function()
+            nelec = 25
+
+            solver = AufbauPrinciple(gf, nelec, occupancy=2)
+            solver.kernel()
+
+            solver_bisect = AufbauPrincipleBisect(gf, nelec, occupancy=2)
+            solver_bisect.kernel()
+
+            assert np.allclose(solver.chempot, solver_bisect.chempot)
+
 if __name__ == "__main__":
     unittest.main()