cRPA Screening and external RPA corrections

examples/ewf/molecules/08-screening-rpa-corrections.py

@@ -0,0 +1,49 @@
+import pyscf.gto
+import pyscf.scf
+import pyscf.cc
+import vayesta.ewf
+
+
+mol = pyscf.gto.Mole()
+mol.atom = """
+O  0.0000   0.0000   0.1173
+H  0.0000   0.7572  -0.4692
+H  0.0000  -0.7572  -0.4692
+"""
+mol.basis = 'cc-pVTZ'
+mol.output = 'pyscf.out'
+mol.build()
+
+# Hartree-Fock
+mf = pyscf.scf.RHF(mol).density_fit()
+mf.kernel()
+
+# Reference full system CCSD:
+cc = pyscf.cc.CCSD(mf)
+cc.kernel()
+
+eta = 1e-6
+
+# Embedded CCSD calculation with bare interactions and no energy correction.
+emb_bare = vayesta.ewf.EWF(mf, bath_options=dict(threshold=eta))
+emb_bare.kernel()
+
+# Embedded CCSD with mRPA screened interactions and RPA cumulant approximation for nonlocal correlations.
+emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=eta), screening="mrpa", ext_rpa_correction="cumulant")
+emb.kernel()
+e_nonlocal_cumulant = emb.e_nonlocal
+# Embedded CCSD with mRPA screened interactions and delta RPA approximation for nonlocal correlations.
+emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=eta), screening="mrpa", ext_rpa_correction="erpa")
+emb.kernel()
+e_nonlocal_erpa = emb.e_nonlocal
+
+# Note that mRPA screening and external corrections often cancel with each other in the case of the energy.
+print("E(CCSD)=                              %+16.8f Ha" % cc.e_tot)
+print("E(RPA)=                               %+16.8f Ha  (error= %+.8f Ha)" % (emb.e_mf + emb.e_rpa,
+                                                                               emb.e_mf + emb.e_rpa - cc.e_tot))
+print("E(Emb. CCSD)=                         %+16.8f Ha  (error= %+.8f Ha)" % (emb_bare.e_tot, emb_bare.e_tot-cc.e_tot))
+print("E(Emb. Screened CCSD)=                %+16.8f Ha  (error= %+.8f Ha)" % (emb.e_tot, emb.e_tot-cc.e_tot))
+print("E(Emb. Screened CCSD + \Delta E_k)=   %+16.8f Ha  (error= %+.8f Ha)" % (emb.e_tot+e_nonlocal_cumulant,
+                                                                               emb.e_tot+e_nonlocal_cumulant-cc.e_tot))
+print("E(Emb. Screened CCSD + \Delta RPA)=   %+16.8f Ha  (error= %+.8f Ha)" % (emb.e_tot+e_nonlocal_erpa,
+                                                                               emb.e_tot+e_nonlocal_erpa-cc.e_tot))

vayesta/core/qemb/fragment.py

@@ -9,7 +9,7 @@
 import pyscf.lo
 # --- Internal
 from vayesta.core.util import (OptionsBase, cache, deprecated, dot, einsum, energy_string, fix_orbital_sign, hstack,
-                               log_time, time_string, timer)
+                               log_time, time_string, timer, AbstractMethodError)
 from vayesta.core import spinalg
 from vayesta.core.types import Cluster
 from vayesta.core.symmetry import SymmetryIdentity
@@ -28,6 +28,7 @@
 from vayesta.misc.cubefile import CubeFile
 from vayesta.mpi import mpi
 from vayesta.solver import get_solver_class, check_solver_config, ClusterHamiltonian
+from vayesta.core.screening.screening_crpa import get_frag_W
 
 # Get MPI rank of fragment
 get_fragment_mpi_rank = lambda *args : args[0].mpi_rank
@@ -45,6 +46,7 @@
     store_eris: bool = None     # If True, ERIs will be cached in Fragment.hamil
     dm_with_frozen: bool = None # TODO: is still used?
     screening: typing.Optional[str] = None
+    match_cluster_fock: bool = None
     # Fragment specific
     # -----------------
     # Auxiliary fragments are treated before non-auxiliary fragments, but do not contribute to expectation values
@@ -73,6 +75,7 @@
         converged: bool = None      # True, if solver reached convergence criterion or no convergence required (eg. MP2 solver)
         # --- Energies
         e_corr: float = None        # Fragment correlation energy contribution
+        e_corr_rpa: float = None    # Fragment correlation energy contribution at the level of RPA.
         # --- Wave-function
         wf: WaveFunction = None     # WaveFunction object (MP2, CCSD,...)
         pwf: WaveFunction = None    # Fragment-projected wave function
@@ -1036,10 +1039,25 @@
             cluster_solver.hamil.add_screening(self._seris_ov)
         return cluster_solver
 
+    def get_local_rpa_correction(self, hamil=None):
+        e_loc_rpa = None
+        if self.base.opts.ext_rpa_correction:
+            hamil = hamil or self.hamil
+            cumulant = self.base.opts.ext_rpa_correction == "cumulant"
+            if self.base.opts.ext_rpa_correction not in ["erpa", "cumulant"]:
+                raise ValueError("Unknown external rpa correction %s specified.")
+            e_loc_rpa = hamil.calc_loc_erpa(*self._seris_ov[1:], cumulant)
+        return e_loc_rpa
+
     def get_frag_hamil(self):
+        if self.opts.screening is not None:
+            if "crpa_full" in self.opts.screening:
+                self.bos_freqs, self.couplings = get_frag_W(self.mf, self, pcoupling=("pcoupled" in self.opts.screening),
+                                                            only_ov_screened=("ov" in self.opts.screening),
+                                                            log=self.log)
         # This detects based on fragment what kind of Hamiltonian is appropriate (restricted and/or EB).
         return ClusterHamiltonian(self, self.mf, self.log, screening=self.opts.screening,
-                                  cache_eris=self.opts.store_eris)
+                                  cache_eris=self.opts.store_eris, match_fock=self.opts.match_cluster_fock)
 
     def get_solver_options(self, *args, **kwargs):
         raise AbstractMethodError

vayesta/core/qemb/qemb.py

@@ -25,7 +25,7 @@
 from vayesta.core.ao2mo import postscf_ao2mo
 from vayesta.core.ao2mo import postscf_kao2gmo
 from vayesta.core.scmf import PDMET, Brueckner
-from vayesta.core.qemb.scrcoulomb import build_screened_eris
+from vayesta.core.screening.screening_moment import build_screened_eris
 from vayesta.mpi import mpi
 from vayesta.core.qemb.register import FragmentRegister
 from vayesta.rpa import ssRIRPA
@@ -105,15 +105,17 @@
             # EBFCI/EBCCSD
             max_boson_occ=2,
             # EBCC
-            ansatz=None,
+            ansatz=None, fermion_wf=False,
             # Dump
             dumpfile='clusters.h5',
             # MP2
             compress_cderi=False)
     # --- Other
     symmetry_tol: float = 1e-6              # Tolerance (in Bohr) for atomic positions
     symmetry_mf_tol: float = 1e-5           # Tolerance for mean-field solution
-    screening: Optional[str] = None
+    screening: Optional[str] = None         # What form of screening to use in clusters.
+    ext_rpa_correction: Optional[str] = None
+    match_cluster_fock: bool = False
 
 class Embedding:
     """Base class for quantum embedding methods.
@@ -518,6 +520,12 @@
         """Nuclear-repulsion energy per unit cell (not folded supercell)."""
         return self.mol.energy_nuc()/self.ncells
 
+    @property
+    def e_nonlocal(self):
+        if self.opts.ext_rpa_correction is None:
+            return 0.0
+        return self.e_rpa - sum([x.results.e_corr_rpa * x.symmetry_factor for x in self.get_fragments(sym_parent=None)])
+
     # Embedding properties
 
     @property
@@ -755,13 +763,33 @@
     @log_method()
     @with_doc(build_screened_eris)
     def build_screened_eris(self, *args, **kwargs):
-        scrfrags = [x.opts.screening for x in self.fragments if x.opts.screening is not None]
-        if len(scrfrags) > 0:
-            # Calculate total dRPA energy in N^4 time; this is cheaper than screening calculations.
-            rpa = ssRIRPA(self.mf, log=self.log)
-            self.e_nonlocal, energy_error = rpa.kernel_energy(correction='linear')
-            if scrfrags.count("mrpa") > 0:
-                build_screened_eris(self, *args, **kwargs)
+        nmomscr = len([x.opts.screening for x in self.fragments if x.opts.screening == "mrpa"])
+        if self.opts.ext_rpa_correction:
+            cumulant = self.opts.ext_rpa_correction == "cumulant"
+            if nmomscr < self.nfrag:
+                raise NotImplementedError("External dRPA correction currently requires all fragments use mrpa screening.")
+
+            if self.opts.ext_rpa_correction not in ["erpa", "cumulant"]:
+                raise ValueError("Unknown external rpa correction %s specified.")
+            l_ = self.get_cderi(self.mf.mo_coeff)[0]
+            rpa = ssRIRPA(self.mf, log=self.log, Lpq=l_)
+            if cumulant:
+                l_ = l_[:, :self.nocc, self.nocc:].reshape((l_.shape[0], -1))
+                l_ = np.concatenate([l_, l_], axis=1)
+
+                m0 = rpa.kernel_moms(0, target_rot=l_)[0][0]
+                # Deduct effective mean-field contribution and project the RHS and we're done.
+                self.e_rpa = 0.5 * einsum("pq,pq->", m0 - l_, l_)
+            else:
+                # Calculate total dRPA energy in N^4 time; this is cheaper than screening calculations.
+                self.e_rpa, energy_error = rpa.kernel_energy(correction='linear')
+            self.log.info("Set total RPA correlation energy contribution as %s", energy_string(self.e_rpa))
+        if nmomscr > 0:
+            self.log.info("")
+            self.log.info("SCREENED INTERACTION SETUP")
+            self.log.info("==========================")
+            with log_time(self.log.timing, "Time for screened interation setup: %s"):
+                build_screened_eris(self, *args, store_m0=self.opts.ext_rpa_correction, **kwargs)
 
     # Symmetry between fragments
     # --------------------------
@@ -1561,6 +1589,7 @@
     def _reset(self):
         self.e_corr = None
         self.converged = False
+        self.e_rpa = None
 
     def reset(self, *args, **kwargs):
         self._reset()

vayesta/core/qemb/uqemb.py

@@ -196,6 +196,13 @@
         eris = postscf_ao2mo(postscf, fock=fock, mo_energy=mo_energy, e_hf=e_hf)
         return eris
 
+    @log_method()
+    @with_doc(Embedding.build_screened_eris)
+    def build_screened_eris(self, *args, **kwargs):
+        if self.opts.ext_rpa_correction is not None:
+            raise NotImplementedError("External RPA correlation energy only implemented for restricted references!")
+        super().build_screened_eris(*args, **kwargs)
+
     def update_mf(self, mo_coeff, mo_energy=None, veff=None):
         """Update underlying mean-field object."""
         # Chech orthonormal MOs