Merge pull request BoothGroup#9 from BoothGroup/pbc

k-space TDA

.github/workflows/ci.yaml

@@ -9,7 +9,7 @@ on:
 
 jobs:
   build:
-    name: python "3.10" on ubuntu-latest
+    name: python "3.11" on ubuntu-latest
     runs-on: ubuntu-latest
 
     strategy:
@@ -20,7 +20,7 @@ jobs:
       - name: Set up python
         uses: actions/setup-python@v2
         with:
-          python-version: "3.10"
+          python-version: "3.11"
       - name: Upgrade pip
         run: |
           python -m pip install --upgrade pip

examples/03-qsgw.py

@@ -34,3 +34,12 @@
 gw = qsGW(mf)
 gw.solver_options = dict(optimise_chempot=True, fock_loop=True)
 gw.kernel(nmom_max=3)
+
+# qsGW finds a self-consistent Green's function in the space defined
+# by the moment expansion, and via quasiparticle self-consistency
+# obtains a (different) set of single particle Koopmans states. Both
+# of these can be accessed via the `gw` object:
+print("QP energies:")
+print(gw.qp_energy)
+print("GF energies:")
+print(gw.gf.energy)

examples/11-srg_qsgw_s_parameter.py

@@ -40,7 +40,7 @@
 
 # GW
 gw = GW(mf)
-_, gf, se = gw.kernel(nmom_max)
+_, gf, se, _ = gw.kernel(nmom_max)
 gf.remove_uncoupled(tol=0.8)
 if which == "ip":
     gw_eta = -gf.get_occupied().energy.max() * HARTREE2EV
@@ -50,12 +50,12 @@
 # qsGW
 gw = qsGW(mf)
 gw.eta = 0.05
-_, gf, se = gw.kernel(nmom_max)
+_, gf, se, _ = gw.kernel(nmom_max)
 gf.remove_uncoupled(tol=0.8)
 if which == "ip":
-    qsgw_eta = -gf.get_occupied().energy.max() * HARTREE2EV
+    qsgw_eta = -np.max(gw.qp_energy[mf.mo_occ > 0]) * HARTREE2EV
 else:
-    qsgw_eta = gf.get_virtual().energy.min() * HARTREE2EV
+    qsgw_eta = np.min(gw.qp_energy[mf.mo_occ == 0]) * HARTREE2EV
 
 # SRG-qsGW
 s_params = sorted(list(data.keys()))[::-1]
@@ -72,16 +72,16 @@
     gw.diis_space = 10
     gw.conv_tol = 1e-5
     gw.conv_tol_moms = 1
-    conv, gf, se = gw.kernel(nmom_max, moments=moments)
+    conv, gf, se, _ = gw.kernel(nmom_max, moments=moments)
     moments = (
             se.get_occupied().moment(range(nmom_max+1)),
             se.get_virtual().moment(range(nmom_max+1)),
     )
     gf.remove_uncoupled(tol=0.8)
     if which == "ip":
-        qsgw_srg.append(-gf.get_occupied().energy.max() * HARTREE2EV)
+        qsgw_srg.append(-np.max(gw.qp_energy[mf.mo_occ > 0]) * HARTREE2EV)
     else:
-        qsgw_srg.append(gf.get_virtual().energy.min() * HARTREE2EV)
+        qsgw_srg.append(np.min(gw.qp_energy[mf.mo_occ == 0]) * HARTREE2EV)
 
 qsgw_srg = np.array(qsgw_srg)
 

examples/30-ktda.py

@@ -0,0 +1,44 @@
+"""Example of running k-space GW calculations with dTDA screening.
+"""
+
+import numpy as np
+from pyscf.pbc import gto, dft
+from momentGW.pbc.gw import KGW
+from momentGW.pbc.qsgw import qsKGW
+from momentGW.pbc.evgw import evKGW
+from momentGW.pbc.scgw import scKGW
+
+cell = gto.Cell()
+cell.atom = "He 0 0 0; He 1 1 1"
+cell.a = np.eye(3) * 3
+cell.basis = "6-31g"
+cell.verbose = 5
+cell.build()
+
+kpts = cell.make_kpts([2, 2, 2])
+
+mf = dft.KRKS(cell, kpts)
+mf = mf.density_fit()
+mf.xc = "hf"
+mf.kernel()
+
+# KGW
+gw = KGW(mf)
+gw.polarizability = "dtda"
+gw.kernel(nmom_max=5)
+
+# qsKGW
+gw = qsKGW(mf)
+gw.polarizability = "dtda"
+gw.srg = 100
+gw.kernel(nmom_max=1)
+
+# evKGW
+gw = evKGW(mf)
+gw.polarizability = "dtda"
+gw.kernel(nmom_max=1)
+
+# scKGW
+gw = scKGW(mf)
+gw.polarizability = "dtda"
+gw.kernel(nmom_max=1)

examples/31-kpt_interpolation.py

@@ -0,0 +1,55 @@
+"""Example of interpolation of a GW calculation onto a new k-point mesh.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from pyscf.pbc import gto, dft
+from momentGW.pbc.gw import KGW
+
+cell = gto.Cell()
+cell.atom = "H 0 0 0; H 0 0 1"
+cell.a = np.array([[25, 0, 0], [0, 25, 0], [0, 0, 2]])
+cell.basis = "sto6g"
+cell.max_memory = 1e10
+cell.verbose = 5
+cell.build()
+
+kmesh1 = [1, 1, 3]
+kmesh2 = [1, 1, 9]
+kpts1 = cell.make_kpts(kmesh1)
+kpts2 = cell.make_kpts(kmesh2)
+
+mf1 = dft.KRKS(cell, kpts1, xc="hf")
+mf1 = mf1.density_fit(auxbasis="weigend")
+mf1.exxdiv = None
+mf1.conv_tol = 1e-10
+mf1.kernel()
+
+mf2 = dft.KRKS(cell, kpts2, xc="hf")
+mf2 = mf2.density_fit(mf1.with_df.auxbasis)
+mf2.exxdiv = mf1.exxdiv
+mf2.conv_tol = mf1.conv_tol
+mf2.kernel()
+
+gw1 = KGW(mf1)
+gw1.polarizability = "dtda"
+gw1.kernel(5)
+
+gw2 = gw1.interpolate(mf2, 5)
+
+e1 = gw1.qp_energy
+e2 = gw2.qp_energy
+
+
+def get_xy(kpts, e):
+    kpts = kpts.wrap_around(kpts._kpts)[:, 2]
+    arg = np.argsort(kpts)
+    return kpts[arg], np.array(e)[arg]
+
+plt.figure()
+plt.plot(*get_xy(gw1.kpts, e1), "C0.-", label="original")
+plt.plot(*get_xy(gw2.kpts, e2), "C2.-", label="interpolated")
+handles, labels = plt.gca().get_legend_handles_labels()
+by_label = dict(zip(labels, handles))
+plt.legend(by_label.values(), by_label.keys())
+plt.show()

momentGW/__init__.py

@@ -52,3 +52,7 @@
 from momentGW.evgw import evGW
 from momentGW.scgw import scGW
 from momentGW.qsgw import qsGW
+from momentGW.pbc.gw import KGW
+from momentGW.pbc.evgw import evKGW
+from momentGW.pbc.qsgw import qsKGW
+from momentGW.pbc.scgw import scKGW

momentGW/base.py

@@ -2,6 +2,8 @@
 Base class for moment-constrained GW solvers.
 """
 
+import warnings
+
 import numpy as np
 from pyscf import lib
 from pyscf.agf2 import mpi_helper
@@ -82,6 +84,9 @@ def __init__(self, mf, **kwargs):
         self.stdout = self.mol.stdout
         self.max_memory = 1e10
 
+        if kwargs.pop("vhf_df", None) is not None:
+            warnings.warn("Keyword argument vhf_df is deprecated.", DeprecationWarning)
+
         for key, val in kwargs.items():
             if not hasattr(self, key):
                 raise AttributeError("%s has no attribute %s", self.name, key)
@@ -97,6 +102,7 @@ def __init__(self, mf, **kwargs):
         self.converged = None
         self.se = None
         self.gf = None
+        self._qp_energy = None
 
         self._keys = set(self.__dict__.keys()).union(self._opts)
 
@@ -118,10 +124,60 @@ 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."""
 
+        if t_prev is None:
+            t_prev = np.zeros_like(t)
+
         error = 0
         for a, b in zip(t, t_prev):
             a = a / max(np.max(np.abs(a)), 1)
@@ -141,6 +197,63 @@ def _gf_to_occ(gf):
 
         return occ
 
+    @staticmethod
+    def _gf_to_energy(gf):
+        """Return the `energy` attribute of a `gf`. Allows hooking in
+        `pbc` methods to retain syntax.
+        """
+        return gf.energy
+
+    @staticmethod
+    def _gf_to_coupling(gf):
+        """Return the `coupling` attribute of a `gf`. Allows hooking in
+        `pbc` methods to retain syntax.
+        """
+        return gf.coupling
+
+    def _gf_to_mo_energy(self, gf):
+        """Find the poles of a GF which best overlap with the MOs.
+
+        Parameters
+        ----------
+        gf : GreensFunction
+            Green's function object.
+
+        Returns
+        -------
+        mo_energy : ndarray
+            Updated MO energies.
+        """
+
+        check = set()
+        mo_energy = np.zeros_like(self.mo_energy)
+
+        for i in range(self.nmo):
+            arg = np.argmax(gf.coupling[i] ** 2)
+            mo_energy[i] = gf.energy[arg]
+            check.add(arg)
+
+        if len(check) != self.nmo:
+            logger.warn(self, "Inconsistent quasiparticle weights!")
+
+        return mo_energy
+
+    @property
+    def qp_energy(self):
+        """
+        Return the quasiparticle energies. For most GW methods, this
+        simply consists of the poles of the `self.gf` that best
+        overlap with the MOs, in order. In some methods such as qsGW,
+        these two quantities are not the same.
+        """
+
+        if self._qp_energy is not None:
+            return self._qp_energy
+
+        qp_energy = self._gf_to_mo_energy(self.gf)
+
+        return qp_energy
+
     @property
     def mol(self):
         return self._scf.mol
@@ -151,6 +264,17 @@ def with_df(self):
             raise ValueError("GW solvers require density fitting.")
         return self._scf.with_df
 
+    @property
+    def has_fock_loop(self):
+        """
+        Returns a boolean indicating whether the solver requires a Fock
+        loop. For most GW methods, this is simply `self.fock_loop`. In
+        some methods such as qsGW, a Fock loop is required with or
+        without `self.fock_loop` for the quasiparticle self-consistency,
+        with this property acting as a hook to indicate this.
+        """
+        return self.fock_loop
+
     get_nmo = get_nmo
     get_nocc = get_nocc
     get_frozen_mask = get_frozen_mask