Support periodic supercell CCSD and add bisection based Aufbau

dyson/expressions/ccsd.py

@@ -3,7 +3,7 @@
 """
 
 import numpy as np
-from pyscf import ao2mo, cc, lib
+from pyscf import ao2mo, cc, lib, pbc, scf
 
 from dyson import util
 from dyson.expressions import BaseExpression
@@ -20,22 +20,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(
+                    "EOM-CCSD not implemented for this type of mean-field object."
+                )
+
+            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,22 +126,35 @@ 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)

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(),

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

dyson/solvers/chempot.py

@@ -120,6 +120,71 @@ def get_greens_function(self):
         return self.gf.copy(chempot=self.chempot, deep=False)
 
 
+class AufbauPrincipleBisect(AufbauPrinciple):
+    """
+    Fill a series of orbitals according to the Aufbau principle using a bisection algorithim.
+
+    Parameters
+    ----------
+    *args : tuple
+        Input arguments.  Either `(gf, nelec)` or `(fock, se, nelec)`
+        where `gf` is the Lehmann representation of the Green's
+        function, `fock` is the Fock matrix, `se` is the Lehmann
+        representation of the self-energy and `nelec` is the number
+        of electrons in the physical space.
+    occupancy : int, optional
+        Occupancy of each state, i.e. `2` for a restricted reference
+        and `1` for other references.  Default value is `2`.
+    """
+
+    def _kernel(self):
+        energies = self.gf.energies
+        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()
+
+        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):
     """
     Shift the self-energy auxiliaries with respect to the Green's

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,26 @@ 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, log=NullLogger())
+            solver.kernel()
+
+            solver_bisect = AufbauPrincipleBisect(gf, nelec, occupancy=2, log=NullLogger())
+            solver_bisect.kernel()
+
+            assert np.allclose(solver.chempot, solver_bisect.chempot)
+
+
 if __name__ == "__main__":
     unittest.main()