Improve inheritence of kernels

momentGW/base.py

@@ -124,6 +124,53 @@ def build_se_moments(self, *args, **kwargs):
     def solve_dyson(self, *args, **kwargs):
         raise NotImplementedError
 
+    def _kernel(self, *args, **kwargs):
+        raise NotImplementedError
+
+    def kernel(
+        self,
+        nmom_max,
+        mo_energy=None,
+        mo_coeff=None,
+        moments=None,
+        integrals=None,
+    ):
+        if mo_coeff is None:
+            mo_coeff = self.mo_coeff
+        if mo_energy is None:
+            mo_energy = self.mo_energy
+
+        cput0 = (logger.process_clock(), logger.perf_counter())
+        self.dump_flags()
+        logger.info(self, "nmom_max = %d", nmom_max)
+
+        self.converged, self.gf, self.se, self._qp_energy = self._kernel(
+            nmom_max,
+            mo_energy,
+            mo_coeff,
+            integrals=integrals,
+        )
+
+        gf_occ = self.gf.get_occupied()
+        gf_occ.remove_uncoupled(tol=1e-1)
+        for n in range(min(5, gf_occ.naux)):
+            en = -gf_occ.energy[-(n + 1)]
+            vn = gf_occ.coupling[:, -(n + 1)]
+            qpwt = np.linalg.norm(vn) ** 2
+            logger.note(self, "IP energy level %d E = %.16g  QP weight = %0.6g", n, en, qpwt)
+
+        gf_vir = self.gf.get_virtual()
+        gf_vir.remove_uncoupled(tol=1e-1)
+        for n in range(min(5, gf_vir.naux)):
+            en = gf_vir.energy[n]
+            vn = gf_vir.coupling[:, n]
+            qpwt = np.linalg.norm(vn) ** 2
+            logger.note(self, "EA energy level %d E = %.16g  QP weight = %0.6g", n, en, qpwt)
+
+        logger.timer(self, self.name, *cput0)
+
+        return self.converged, self.gf, self.se, self.qp_energy
+
     @staticmethod
     def _moment_error(t, t_prev):
         """Compute scaled error between moments."""

momentGW/evgw.py

@@ -50,6 +50,9 @@ def kernel(
         Green's function object
     se : pyscf.agf2.SelfEnergy
         Self-energy object
+    qp_energy : numpy.ndarray
+        Quasiparticle energies. Always None for evGW, returned for
+        compatibility with other evGW methods.
     """
 
     logger.warn(gw, "evGW is untested!")
@@ -118,7 +121,7 @@ def kernel(
         if conv:
             break
 
-    return conv, gf, se
+    return conv, gf, se, None
 
 
 class evGW(GW):
@@ -178,6 +181,8 @@ class evGW(GW):
     def name(self):
         return "evG%sW%s" % ("0" if self.g0 else "", "0" if self.w0 else "")
 
+    _kernel = kernel
+
     def check_convergence(self, mo_energy, mo_energy_prev, th, th_prev, tp, tp_prev):
         """Check for convergence, and print a summary of changes."""
 
@@ -201,53 +206,3 @@ def check_convergence(self, mo_energy, mo_energy_prev, th, th_prev, tp, tp_prev)
                 max(error_th, error_tp) < self.conv_tol_moms,
             )
         )
-
-    def kernel(
-        self,
-        nmom_max,
-        mo_energy=None,
-        mo_coeff=None,
-        moments=None,
-        integrals=None,
-    ):
-        if mo_coeff is None:
-            mo_coeff = self.mo_coeff
-        if mo_energy is None:
-            mo_energy = self.mo_energy
-
-        cput0 = (logger.process_clock(), logger.perf_counter())
-        self.dump_flags()
-        logger.info(self, "nmom_max = %d", nmom_max)
-
-        self.converged, self.gf, self.se = kernel(
-            self,
-            nmom_max,
-            mo_energy,
-            mo_coeff,
-            integrals=integrals,
-        )
-
-        gf_occ = self.gf.get_occupied()
-        gf_occ.remove_uncoupled(tol=1e-1)
-        for n in range(min(5, gf_occ.naux)):
-            en = -gf_occ.energy[-(n + 1)]
-            vn = gf_occ.coupling[:, -(n + 1)]
-            qpwt = np.linalg.norm(vn) ** 2
-            logger.note(self, "IP energy level %d E = %.16g  QP weight = %0.6g", n, en, qpwt)
-
-        gf_vir = self.gf.get_virtual()
-        gf_vir.remove_uncoupled(tol=1e-1)
-        for n in range(min(5, gf_vir.naux)):
-            en = gf_vir.energy[n]
-            vn = gf_vir.coupling[:, n]
-            qpwt = np.linalg.norm(vn) ** 2
-            logger.note(self, "EA energy level %d E = %.16g  QP weight = %0.6g", n, en, qpwt)
-
-        if self.converged:
-            logger.note(self, "%s converged", self.name)
-        else:
-            logger.note(self, "%s failed to converge", self.name)
-
-        logger.timer(self, self.name, *cput0)
-
-        return self.converged, self.gf, self.se

momentGW/gw.py

@@ -53,12 +53,15 @@ def kernel(
     Returns
     -------
     conv : bool
-        Convergence flag. Always True for AGW, returned for
+        Convergence flag. Always True for GW, returned for
         compatibility with other GW methods.
     gf : pyscf.agf2.GreensFunction
         Green's function object
     se : pyscf.agf2.SelfEnergy
         Self-energy object
+    qp_energy : numpy.ndarray
+        Quasiparticle energies. Always None for GW, returned for
+        compatibility with other GW methods.
     """
 
     if integrals is None:
@@ -85,7 +88,7 @@ def kernel(
     gf, se = gw.solve_dyson(th, tp, se_static, integrals=integrals)
     conv = True
 
-    return conv, gf, se
+    return conv, gf, se, None
 
 
 class GW(BaseGW):
@@ -98,6 +101,8 @@ class GW(BaseGW):
     def name(self):
         return "G0W0"
 
+    _kernel = kernel
+
     def build_se_static(self, integrals, mo_coeff=None, mo_energy=None):
         """Build the static part of the self-energy, including the
         Fock matrix.
@@ -301,48 +306,3 @@ def moment_error(self, se_moments_hole, se_moments_part, se):
         )
 
         return eh, ep
-
-    def kernel(
-        self,
-        nmom_max,
-        mo_energy=None,
-        mo_coeff=None,
-        moments=None,
-        integrals=None,
-    ):
-        if mo_coeff is None:
-            mo_coeff = self.mo_coeff
-        if mo_energy is None:
-            mo_energy = self.mo_energy
-
-        cput0 = (logger.process_clock(), logger.perf_counter())
-        self.dump_flags()
-        logger.info(self, "nmom_max = %d", nmom_max)
-
-        self.converged, self.gf, self.se = kernel(
-            self,
-            nmom_max,
-            mo_energy,
-            mo_coeff,
-            integrals=integrals,
-        )
-
-        gf_occ = self.gf.get_occupied()
-        gf_occ.remove_uncoupled(tol=1e-1)
-        for n in range(min(5, gf_occ.naux)):
-            en = -gf_occ.energy[-(n + 1)]
-            vn = gf_occ.coupling[:, -(n + 1)]
-            qpwt = np.linalg.norm(vn) ** 2
-            logger.note(self, "IP energy level %d E = %.16g  QP weight = %0.6g", n, en, qpwt)
-
-        gf_vir = self.gf.get_virtual()
-        gf_vir.remove_uncoupled(tol=1e-1)
-        for n in range(min(5, gf_vir.naux)):
-            en = gf_vir.energy[n]
-            vn = gf_vir.coupling[:, n]
-            qpwt = np.linalg.norm(vn) ** 2
-            logger.note(self, "EA energy level %d E = %.16g  QP weight = %0.6g", n, en, qpwt)
-
-        logger.timer(self, self.name, *cput0)
-
-        return self.converged, self.gf, self.se

momentGW/pbc/base.py

@@ -81,6 +81,50 @@ def __init__(self, mf, **kwargs):
 
         self._keys = set(self.__dict__.keys()).union(self._opts)
 
+    def kernel(
+        self,
+        nmom_max,
+        mo_energy=None,
+        mo_coeff=None,
+        moments=None,
+        integrals=None,
+    ):
+        if mo_coeff is None:
+            mo_coeff = self.mo_coeff
+        if mo_energy is None:
+            mo_energy = self.mo_energy
+
+        cput0 = (logger.process_clock(), logger.perf_counter())
+        self.dump_flags()
+        logger.info(self, "nmom_max = %d", nmom_max)
+
+        self.converged, self.gf, self.se, self._qp_energy = self._kernel(
+            nmom_max,
+            mo_energy,
+            mo_coeff,
+            integrals=integrals,
+        )
+
+        gf_occ = self.gf[0].get_occupied()
+        gf_occ.remove_uncoupled(tol=1e-1)
+        for n in range(min(5, gf_occ.naux)):
+            en = -gf_occ.energy[-(n + 1)]
+            vn = gf_occ.coupling[:, -(n + 1)]
+            qpwt = np.linalg.norm(vn) ** 2
+            logger.note(self, "IP energy level (Γ) %d E = %.16g  QP weight = %0.6g", n, en, qpwt)
+
+        gf_vir = self.gf[0].get_virtual()
+        gf_vir.remove_uncoupled(tol=1e-1)
+        for n in range(min(5, gf_vir.naux)):
+            en = gf_vir.energy[n]
+            vn = gf_vir.coupling[:, n]
+            qpwt = np.linalg.norm(vn) ** 2
+            logger.note(self, "EA energy level (Γ) %d E = %.16g  QP weight = %0.6g", n, en, qpwt)
+
+        logger.timer(self, self.name, *cput0)
+
+        return self.converged, self.gf, self.se, self.qp_energy
+
     @staticmethod
     def _gf_to_occ(gf):
         return tuple(BaseGW._gf_to_occ(g) for g in gf)

momentGW/pbc/evgw.py

@@ -12,7 +12,7 @@
 from pyscf.pbc import dft, gto
 from pyscf.pbc.tools import k2gamma
 
-from momentGW.evgw import evGW, kernel
+from momentGW.evgw import evGW
 from momentGW.pbc.gw import KGW
 
 
@@ -51,48 +51,3 @@ def check_convergence(self, mo_energy, mo_energy_prev, th, th_prev, tp, tp_prev)
                 max(error_th, error_tp) < self.conv_tol_moms,
             )
         )
-
-    def kernel(
-        self,
-        nmom_max,
-        mo_energy=None,
-        mo_coeff=None,
-        moments=None,
-        integrals=None,
-    ):
-        if mo_coeff is None:
-            mo_coeff = self.mo_coeff
-        if mo_energy is None:
-            mo_energy = self.mo_energy
-
-        cput0 = (logger.process_clock(), logger.perf_counter())
-        self.dump_flags()
-        logger.info(self, "nmom_max = %d", nmom_max)
-
-        self.converged, self.gf, self.se = kernel(
-            self,
-            nmom_max,
-            mo_energy,
-            mo_coeff,
-            integrals=integrals,
-        )
-
-        gf_occ = self.gf[0].get_occupied()
-        gf_occ.remove_uncoupled(tol=1e-1)
-        for n in range(min(5, gf_occ.naux)):
-            en = -gf_occ.energy[-(n + 1)]
-            vn = gf_occ.coupling[:, -(n + 1)]
-            qpwt = np.linalg.norm(vn) ** 2
-            logger.note(self, "IP energy level (Γ) %d E = %.16g  QP weight = %0.6g", n, en, qpwt)
-
-        gf_vir = self.gf[0].get_virtual()
-        gf_vir.remove_uncoupled(tol=1e-1)
-        for n in range(min(5, gf_vir.naux)):
-            en = gf_vir.energy[n]
-            vn = gf_vir.coupling[:, n]
-            qpwt = np.linalg.norm(vn) ** 2
-            logger.note(self, "EA energy level (Γ) %d E = %.16g  QP weight = %0.6g", n, en, qpwt)
-
-        logger.timer(self, self.name, *cput0)
-
-        return self.converged, self.gf, self.se

momentGW/pbc/gw.py

@@ -10,7 +10,7 @@
 from pyscf.lib import logger
 from pyscf.pbc import scf
 
-from momentGW.gw import GW, kernel
+from momentGW.gw import GW
 from momentGW.pbc.base import BaseKGW
 from momentGW.pbc.fock import fock_loop, minimize_chempot, search_chempot
 from momentGW.pbc.ints import KIntegrals
@@ -207,48 +207,3 @@ def make_rdm1(self, gf=None):
             gf = [GreensFunction(self.mo_energy, np.eye(self.nmo))]
 
         return np.array([g.make_rdm1() for g in gf])
-
-    def kernel(
-        self,
-        nmom_max,
-        mo_energy=None,
-        mo_coeff=None,
-        moments=None,
-        integrals=None,
-    ):
-        if mo_coeff is None:
-            mo_coeff = self.mo_coeff
-        if mo_energy is None:
-            mo_energy = self.mo_energy
-
-        cput0 = (logger.process_clock(), logger.perf_counter())
-        self.dump_flags()
-        logger.info(self, "nmom_max = %d", nmom_max)
-
-        self.converged, self.gf, self.se = kernel(
-            self,
-            nmom_max,
-            mo_energy,
-            mo_coeff,
-            integrals=integrals,
-        )
-
-        gf_occ = self.gf[0].get_occupied()
-        gf_occ.remove_uncoupled(tol=1e-1)
-        for n in range(min(5, gf_occ.naux)):
-            en = -gf_occ.energy[-(n + 1)]
-            vn = gf_occ.coupling[:, -(n + 1)]
-            qpwt = np.linalg.norm(vn) ** 2
-            logger.note(self, "IP energy level (Γ) %d E = %.16g  QP weight = %0.6g", n, en, qpwt)
-
-        gf_vir = self.gf[0].get_virtual()
-        gf_vir.remove_uncoupled(tol=1e-1)
-        for n in range(min(5, gf_vir.naux)):
-            en = gf_vir.energy[n]
-            vn = gf_vir.coupling[:, n]
-            qpwt = np.linalg.norm(vn) ** 2
-            logger.note(self, "EA energy level (Γ) %d E = %.16g  QP weight = %0.6g", n, en, qpwt)
-
-        logger.timer(self, self.name, *cput0)
-
-        return self.converged, self.gf, self.se