Head and wings corrections for k-dTDA

examples/20-finite_size_corrections.py

@@ -0,0 +1,46 @@
+"""
+Examples of finite size corrections for `momentGW` calculations for periodic solids.
+"""
+
+import numpy as np
+from pyscf.pbc import dft, gto
+
+from momentGW import KGW, KUGW
+
+# Define a unit cell
+cell = gto.Cell()
+cell.atom = "He 0 0 0; He 1 1 1"
+cell.basis = "6-31g"
+cell.a = np.eye(3) * 3
+cell.verbose = 5
+cell.build()
+kpts = cell.make_kpts([3, 1, 1])
+
+# Run a DFT calculation
+mf = dft.KRKS(cell, kpts)
+mf = mf.density_fit()
+mf.xc = "b3lyp"
+mf.kernel()
+
+# For any `pbc` calculation a finite size correction can be added for
+# the divergence associated with G = 0, q=0 for GDF Coulomb integrals.
+# This is done by setting the `fsc` attribute to any combination of
+# "H", "W", and/or "B" for the Head, Wings and Body portion of the
+# correction. The Ewald correction is added in all cases. The default
+# is `None` which disables the correction.
+
+# RHF reference
+gw = KGW(mf)
+gw.polarizability = "dRPA"
+gw.fsc = "HWB"
+gw.kernel(nmom_max=3)
+
+gw = KGW(mf)
+gw.polarizability = "dTDA"
+gw.fsc = "HW"
+gw.kernel(nmom_max=3)
+
+gw = KGW(mf)
+gw.polarizability = "dRPA"
+gw.fsc = "H"
+gw.kernel(nmom_max=3)

momentGW/gw.py

@@ -144,9 +144,37 @@ def name(self):
         polarizability = self.polarizability.upper().replace("DTDA", "dTDA").replace("DRPA", "dRPA")
         return f"{polarizability}-G0W0"
 
+    def get_veff(self, integrals, dm=None, **kwargs):
+        """Get the effective potential.
+
+        Parameters
+        ----------
+        integrals : Integrals
+            Integrals object.
+        dm : numpy.ndarray, optional
+            Density matrix. If `None`, determine using `self.make_rdm1`.
+            Default value is `None`.
+        **kwargs : dict, optional
+            Additional keyword arguments passed to the integrals object.
+
+        Returns
+        -------
+        veff : numpy.ndarray
+            Effective potential.
+        """
+
+        # Get the density matrix
+        if dm is None:
+            dm = self.make_rdm1()
+
+        # Get the effective potential
+        veff = integrals.get_veff(dm, **kwargs)
+
+        return veff
+
     @logging.with_timer("Static self-energy")
     @logging.with_status("Building static self-energy")
-    def build_se_static(self, integrals):
+    def build_se_static(self, integrals, force_build=False):
         """
         Build the static part of the self-energy, including the Fock
         matrix.
@@ -155,6 +183,10 @@ def build_se_static(self, integrals):
         ----------
         integrals : Integrals
             Integrals object.
+        force_build : bool, optional
+            If `True`, force the static self-energy to be built, even if
+            the underlying mean-field object is Hartree--Fock. Default
+            value is `False`.
 
         Returns
         -------
@@ -168,15 +200,15 @@ def build_se_static(self, integrals):
         dm = self._scf.make_rdm1(mo_coeff=self._mo_coeff)
 
         # Get the contribution from the exchange-correlation potential
-        if getattr(self._scf, "xc", "hf") == "hf":
+        if getattr(self._scf, "xc", "hf") == "hf" and not force_build:
             se_static = np.zeros_like(dm)
             se_static = se_static[..., mask, :][..., :, mask]
         else:
             with util.SilentSCF(self._scf):
                 veff = self._scf.get_veff(None, dm)[..., mask, :][..., :, mask]
                 vj = self._scf.get_j(None, dm)[..., mask, :][..., :, mask]
 
-            vhf = integrals.get_veff(dm, j=vj, basis="ao")
+            vhf = self.get_veff(integrals, dm=dm, j=vj, basis="ao")
             se_static = vhf - veff
             se_static = util.einsum(
                 "...pq,...pi,...qj->...ij", se_static, np.conj(self.mo_coeff), self.mo_coeff

momentGW/pbc/base.py

@@ -53,14 +53,14 @@ class BaseKGW(BaseGW):
     thc_opts : dict, optional
         Dictionary of options to be used for THC calculations. Current
         implementation requires a filepath to import the THC integrals.
-    fc : bool, optional
+    fsc : bool, optional
         If `True`, apply finite size corrections. Default value is
         `False`.
     """
 
     _defaults = OrderedDict(
         **BaseGW._defaults,
-        fc=False,
+        fsc=None,
     )
     _defaults["compression"] = None
 
@@ -71,7 +71,7 @@ def __init__(self, mf, **kwargs):
         super().__init__(mf, **kwargs)
 
         # Options
-        self.fc = False
+        self.fsc = None
 
         # Attributes
         self._kpts = KPoints(self.cell, getattr(mf, "kpts", np.zeros((1, 3))))

momentGW/pbc/gw.py

@@ -56,7 +56,7 @@ class KGW(BaseKGW, GW):
     thc_opts : dict, optional
         Dictionary of options to be used for THC calculations. Current
         implementation requires a filepath to import the THC integrals.
-    fc : bool, optional
+    fsc : bool, optional
         If `True`, apply finite size corrections. Default value is
         `False`.
     """
@@ -69,9 +69,44 @@ def name(self):
         polarizability = self.polarizability.upper().replace("DTDA", "dTDA").replace("DRPA", "dRPA")
         return f"{polarizability}-KG0W0"
 
+    def get_veff(self, integrals, dm=None, **kwargs):
+        """Get the effective potential.
+
+        Parameters
+        ----------
+        integrals : KIntegrals
+            Integrals object.
+        dm : numpy.ndarray, optional
+            Density matrix at each k-point. If `None`, determine using
+            `self.make_rdm1`. Default value is `None`.
+        **kwargs : dict, optional
+            Additional keyword arguments passed to the integrals object.
+
+        Returns
+        -------
+        veff : numpy.ndarray
+            Effective potential at each k-point.
+        """
+
+        # Get the density matrix
+        if dm is None:
+            dm = self.make_rdm1()
+
+        # Set the default options
+        if "ewald" not in kwargs:
+            if self.fsc is not None:
+                kwargs = {**kwargs, "ewald": True}
+            else:
+                kwargs = {**kwargs, "ewald": False}
+
+        # Get the effective potential
+        veff = integrals.get_veff(dm, **kwargs)
+
+        return veff
+
     @logging.with_timer("Static self-energy")
     @logging.with_status("Building static self-energy")
-    def build_se_static(self, integrals):
+    def build_se_static(self, integrals, **kwargs):
         """
         Build the static part of the self-energy, including the Fock
         matrix.
@@ -80,14 +115,29 @@ def build_se_static(self, integrals):
         ----------
         integrals : KIntegrals
             Integrals object.
+        **kwargs : dict, optional
+            Additional keyword arguments.
 
         Returns
         -------
         se_static : numpy.ndarray
             Static part of the self-energy at each k-point. If
             `self.diagonal_se`, non-diagonal elements are set to zero.
         """
-        return super().build_se_static(integrals)
+        if self.fsc is not None:
+            if len(list(self.fsc)) > 3:
+                raise ValueError(
+                    "Finite size corrections require as an input a combination of H, W and B "
+                    "for the different finite size corrections (H - Head, W - Wing, B - Body)")
+            for i, letter in enumerate(list(self.fsc)):
+                if letter not in ["H", "W", "B"]:
+                    raise ValueError(
+                        "Finite size corrections require as an input a combination of H, W and B "
+                        "for the different finite size corrections (H - Head, W - Wing, B - Body)")
+            kwargs = {**kwargs, "force_build": True}
+        else:
+            kwargs = {**kwargs, "force_build": False}
+        return super().build_se_static(integrals, **kwargs)
 
     def build_se_moments(self, nmom_max, integrals, **kwargs):
         """Build the moments of the self-energy.
@@ -112,10 +162,10 @@ def build_se_moments(self, nmom_max, integrals, **kwargs):
         """
 
         if self.polarizability.lower() == "dtda":
-            tda = dTDA(self, nmom_max, integrals, **kwargs)
+            tda = dTDA(self, nmom_max, integrals, fsc=self.fsc, **kwargs)
             return tda.kernel()
         if self.polarizability.lower() == "drpa":
-            rpa = dRPA(self, nmom_max, integrals, **kwargs)
+            rpa = dRPA(self, nmom_max, integrals, fsc=self.fsc, **kwargs)
             return rpa.kernel()
         elif self.polarizability.lower() == "thc-dtda":
             tda = thc.dTDA(self, nmom_max, integrals, **kwargs)
@@ -154,6 +204,8 @@ def ao2mo(self, transform=True):
                 compression=self.compression,
                 compression_tol=self.compression_tol,
                 store_full=self.fock_loop,
+                mo_energy_w=self.mo_energy,
+                fsc=self.fsc,
                 input_path=self.thc_opts["file_path"],
             )
 
@@ -334,7 +386,7 @@ def make_rdm1(self, gf=None):
 
     @logging.with_timer("Energy")
     @logging.with_status("Calculating energy")
-    def energy_hf(self, gf=None, integrals=None):
+    def energy_hf(self, gf=None, integrals=None, **kwargs):
         """Calculate the one-body (Hartree--Fock) energy.
 
         Parameters
@@ -367,7 +419,7 @@ def energy_hf(self, gf=None, integrals=None):
                 "kpq,kpi,kqj->kij", self._scf.get_hcore(), self.mo_coeff.conj(), self.mo_coeff
             )
         rdm1 = self.make_rdm1()
-        fock = integrals.get_fock(rdm1, h1e)
+        fock = integrals.get_fock(rdm1, h1e, **kwargs)
 
         # Calculate the Hartree--Fock energy at each k-point
         e_1b = sum(energy.hartree_fock(rdm1[k], fock[k], h1e[k]) for k in self.kpts.loop(1))