Merge branch 'master' into CCSDt_interface

docs/source/quickstart/edmetfinite.py

@@ -26,7 +26,7 @@
 # The KS object needs to be converted to a HF object:
 lda = lda.to_rhf()
 
-emb_lda = vayesta.edmet.EDMET(lda, dmet_threshold=1e-12, solver="EBCCSD",
+emb_lda = vayesta.edmet.EDMET(lda, bath_options=dict(dmet_threshold=1e-12), solver="CCSD-S-1-1",
                               oneshot=True, make_dd_moments=False)
 
 with emb_lda.iao_fragmentation() as f:
@@ -42,7 +42,7 @@
 # The KS opbject needs to be converted to a HF object:
 b3lyp = b3lyp.to_rhf()
 
-emb_b3lyp = vayesta.edmet.EDMET(b3lyp, dmet_threshold=1e-12, solver="EBCCSD", oneshot=True, make_dd_moments=False)
+emb_b3lyp = vayesta.edmet.EDMET(b3lyp, bath_options=dict(dmet_threshold=1e-12), solver="CCSD-S-1-1", oneshot=True, make_dd_moments=False)
 with emb_b3lyp.iao_fragmentation() as f:
     f.add_all_atomic_fragments()
 emb_b3lyp.kernel()
@@ -51,7 +51,7 @@
 hf = pyscf.scf.RHF(mol).density_fit()
 hf.kernel()
 
-emb_hf = vayesta.edmet.EDMET(hf, dmet_threshold=1e-12, solver="EBCCSD",
+emb_hf = vayesta.edmet.EDMET(hf, bath_options=dict(dmet_threshold=1e-12), solver="CCSD-S-1-1",
                              oneshot=True, make_dd_moments=False)
 
 with emb_hf.iao_fragmentation() as f:
@@ -65,7 +65,7 @@
 print("E(LDA)=              %+16.8f Ha" % lda.e_tot)
 print("E(Emb. CCSD @LDA)=   %+16.8f Ha" % emb_lda.e_tot)
 print("E(HF @LDA)=          %+16.8f Ha" % emb_lda.e_mf)
-print("E(CCSD)=             %+16.8f Ha" % cc.)
+print("E(CCSD)=             %+16.8f Ha" % cc.e_tot)
 print("E(B3LYP)=            %+16.8f Ha" % b3lyp.e_tot)
 print("E(HF)=               %+16.8f Ha" % hf.e_tot)
 print("E(HF @B3LYP)=        %+16.8f Ha" % emb_b3lyp.e_mf)

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/bath/dmet.py

@@ -47,7 +47,8 @@ def kernel(self):
         # --- Separate occupied and virtual in cluster
         cluster = [self.c_frag, c_dmet]
         self.base._check_orthonormal(*cluster, mo_name='cluster MO')
-        c_cluster_occ, c_cluster_vir = self.fragment.diagonalize_cluster_dm(*cluster, tol=2*self.dmet_threshold)
+        tol = self.base.opts.bath_options['occupation_tolerance']
+        c_cluster_occ, c_cluster_vir = self.fragment.diagonalize_cluster_dm(*cluster, tol=tol)
         # Canonicalize
         c_cluster_occ = self.fragment.canonicalize_mo(c_cluster_occ)[0]
         c_cluster_vir = self.fragment.canonicalize_mo(c_cluster_vir)[0]

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 @@ class Options(OptionsBase):
     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 @@ class Results:
         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
@@ -517,8 +520,8 @@ def diagonalize_cluster_dm(self, *mo_coeff, dm1=None, norm=2, tol=1e-4):
         sc = np.dot(self.base.get_ovlp(), c_cluster)
         dm = dot(sc.T, dm1, sc)
         e, r = np.linalg.eigh(dm)
-        if tol and not np.allclose(np.fmin(abs(e), abs(e-norm)), 0, atol=tol, rtol=0):
-            raise RuntimeError("Eigenvalues of cluster-DM not all close to 0 or %d:\n%s" % (norm, e))
+        if tol and not np.allclose(np.fmin(abs(e), abs(e-norm)), 0, atol=2*tol, rtol=0):
+            self.log.warn("Eigenvalues of cluster-DM not all close to 0 or %d:\n%s" % (norm, e))
         e, r = e[::-1], r[:,::-1]
         c_cluster = np.dot(c_cluster, r)
         c_cluster = fix_orbital_sign(c_cluster)[0]
@@ -863,7 +866,7 @@ def get_orbitals(occtype):
         def check_occup(mo_coeff, expected):
             occup = self.get_mo_occupation(mo_coeff)
             # RHF
-            atol = self.opts.bath_options['dmet_threshold']
+            atol = self.opts.bath_options['occupation_tolerance']
             if np.ndim(occup[0]) == 0:
                 assert np.allclose(occup, 2*expected, rtol=0, atol=2*atol)
             else:
@@ -1036,10 +1039,25 @@ def get_solver(self, solver=None):
             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
@@ -67,7 +67,7 @@ class Options(OptionsBase):
     # --- Bath options
     bath_options: dict = OptionsBase.dict_with_defaults(
         # General
-        bathtype='dmet', canonicalize=True,
+        bathtype='dmet', canonicalize=True, occupation_tolerance=1e-8,
         # DMET bath
         dmet_threshold=1e-8,
         # R2 bath
@@ -105,15 +105,17 @@ class Options(OptionsBase):
             # EBFCI/EBCCSD
             max_boson_occ=2,
             # EBCC
-            ansatz=None, store_as_ccsd=None,
+            ansatz=None, store_as_ccsd=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 @@ def e_nuc(self):
         """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 @@ def get_eris_object(self, postscf, fock=None):
     @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_fragments(self, *args, **kwargs):
     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 @@ def get_eris_object(self, postscf, fock=None):
         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