diff --git a/examples/ewf/molecules/08-screening-rpa-corrections.py b/examples/ewf/molecules/08-screening-rpa-corrections.py new file mode 100644 index 000000000..17375d4ef --- /dev/null +++ b/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)) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index a6a152de5..72deaf53d 100644 --- a/vayesta/core/qemb/fragment.py +++ b/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 @@ -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 diff --git a/vayesta/core/qemb/qemb.py b/vayesta/core/qemb/qemb.py index d28b89d13..d0a1dd8c3 100644 --- a/vayesta/core/qemb/qemb.py +++ b/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,7 +105,7 @@ class Options(OptionsBase): # EBFCI/EBCCSD max_boson_occ=2, # EBCC - ansatz=None, + ansatz=None, fermion_wf=False, # Dump dumpfile='clusters.h5', # MP2 @@ -113,7 +113,9 @@ class Options(OptionsBase): # --- 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() diff --git a/vayesta/core/qemb/uqemb.py b/vayesta/core/qemb/uqemb.py index fba6e8244..46cce67f9 100644 --- a/vayesta/core/qemb/uqemb.py +++ b/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 diff --git a/vayesta/core/screening/__init__.py b/vayesta/core/screening/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/vayesta/core/screening/screening_crpa.py b/vayesta/core/screening/screening_crpa.py new file mode 100644 index 000000000..0de9fd65d --- /dev/null +++ b/vayesta/core/screening/screening_crpa.py @@ -0,0 +1,223 @@ +import scipy.linalg + +from vayesta.rpa import ssRPA +from .screening_moment import _get_target_rot +import copy +from vayesta.core.util import * +import numpy as np +from pyscf import lib +from pyscf.lib import logger +from pyscf.ao2mo import _ao2mo +from pyscf import __config__ + + +class cRPAError(RuntimeError): + pass + +def get_frag_W(mf, fragment, pcoupling=True, only_ov_screened=False, log=None): + """Generates screened coulomb interaction due to screening at the level of cRPA. + Note that this currently scales as O(N_frag N^6), so is not practical without further refinement. + + Parameters + ---------- + mf : pyscf.scf object + Mean-field instance. + fragment : vayesta.qemb.Fragment subclass + Fragments for the calculation, used to define local interaction space. + log : logging.Logger, optional + Logger object. If None, the logger of the `emb` object is used. Default: None. + + Returns + ------- + freqs : np.array + Effective bosonic frequencies for cRPA screening. + couplings : np.array + Effective bosonic couplings for cRPA screening. + """ + + log.info("Generating screened interaction via frequency dependent cRPA.") + try: + l_a, l_b, crpa = set_up_W_crpa(mf, fragment, pcoupling, only_ov_screened, log=log) + except cRPAError as e: + freqs = np.zeros((0,)) + nmo = mf.mo_coeff.shape[1] + couplings = (np.zeros((0, nmo, nmo)), np.zeros((0, nmo, nmo))) + else: + freqs = crpa.freqs_ss + couplings = (l_a, l_b) + log.info("cRPA resulted in %d poles", len(freqs)) + + return freqs, couplings + + +def get_frag_deltaW(mf, fragment, pcoupling=True, only_ov_screened=False, log=None): + """Generates change in coulomb interaction due to screening at the level of static limit of cRPA. + Note that this currently scales as O(N_frag N^6), so is not practical without further refinement. + + Parameters + ---------- + mf : pyscf.scf object + Mean-field instance. + fragment : vayesta.qemb.Fragment subclass + Fragments for the calculation, used to define local interaction space. + log : logging.Logger, optional + Logger object. If None, the logger of the `emb` object is used. Default: None. + + Returns + ------- + deltaW : np.array + Change in cluster local coulomb interaction due to cRPA screening. + """ + + log.info("Generating screened interaction via static limit of cRPA.") + log.warning("This is poorly defined for non-CAS fragmentations.") + log.warning("This implementation is expensive, with O(N^6) computational cost per cluster.") + try: + l_a, l_b, crpa = set_up_W_crpa(mf, fragment, pcoupling, only_ov_screened=only_ov_screened, log=log) + except cRPAError as e: + nmo = mf.mo_coeff.shape[1] + delta_w = tuple([np.zeros([nmo] * 4)] * 4) + crpa = None + else: + # Have a factor of -2 due to negative value of RPA dd response, and summation of + # the excitation and deexcitation branches of the dd response. + static_fac = - 1.0 * (crpa.freqs_ss ** (-1)) + + delta_w = ( + einsum("npq,n,nrs->pqrs", l_a, static_fac, l_a) + einsum("nqp,n,nsr->pqrs", l_a, static_fac, l_a), + einsum("npq,n,nrs->pqrs", l_a, static_fac, l_b) + einsum("nqp,n,nsr->pqrs", l_a, static_fac, l_b), + einsum("npq,n,nrs->pqrs", l_b, static_fac, l_b) + einsum("nqp,n,nsr->pqrs", l_b, static_fac, l_b), + ) + return delta_w, crpa + +def set_up_W_crpa(mf, fragment, pcoupling=True, only_ov_screened=False, log=None): + is_rhf = np.ndim(mf.mo_coeff[1]) == 1 + if not hasattr(mf, "with_df"): + raise NotImplementedError("Screened interactions require density-fitting.") + crpa, rot_loc, rot_crpa = get_crpa(mf, fragment, log) + # Now need to calculate interactions. + nmo = mf.mo_coeff.shape[1] + nocc = sum(mf.mo_occ.T > 0) + c = mf.mo_coeff + if is_rhf: + nocc = (nocc, nocc) + c = (c, c) + + # First calculate alpha contribution. + lov_a = ao2mo(mf, mo_coeff=c[0], ijslice=(0, nocc[0], nocc[0], nmo)).reshape((-1, crpa.ova)) + lov_a = dot(lov_a, crpa.ov_rot[0].T) + l_aux = dot(lov_a, crpa.XpY_ss[0]) + del lov_a + lov_b = ao2mo(mf, mo_coeff=c[1], ijslice=(0, nocc[1], nocc[1], nmo)).reshape((-1, crpa.ovb)) + lov_b = dot(lov_b, crpa.ov_rot[1].T) + + # This is a decomposition of the cRPA reducible dd response in the auxilliary basis + l_aux += dot(lov_b, crpa.XpY_ss[1]) + del lov_b + + # This is expensive, and we'd usually want to avoid doing it twice unnecessarily, but other things will be worse. + # Now calculate the coupling back to the fragment itself. + c_act = fragment.cluster.c_active + if is_rhf: + c_act = (c_act, c_act) + + # Calculating this quantity scales as O(N^3), rather than O(N^4) if we explicitly calculated the cRPA dd response + # in the auxiliary basis. + lpqa_loc = ao2mo(mf, mo_coeff=c_act[0]) + l_a = einsum("npq,nm->mpq", lpqa_loc, l_aux) + del lpqa_loc + lpqb_loc = ao2mo(mf, mo_coeff=c_act[1]) + l_b = einsum("npq,nm->mpq", lpqb_loc, l_aux) + del lpqb_loc + + if pcoupling: + # Need to calculate additional contribution to the couplings resulting from rotation of the irreducible + # polarisability. + # Generate the full-system matrix of orbital energy differences. + eps = ssRPA(mf)._gen_eps() + # First, generate epsilon couplings between cluster and crpa spaces. + eps_fb = [einsum("p,qp,rp->qr", e, l, nl) for e, l, nl in zip(eps, rot_loc, crpa.ov_rot)] + # Then generate X and Y values for this correction. + x_crpa = [(p + m)/2 for p, m in zip(crpa.XpY_ss, crpa.XmY_ss)] + y_crpa = [(p - m)/2 for p, m in zip(crpa.XpY_ss, crpa.XmY_ss)] + # Contract with epsilon values + a_fb = [dot(e, x) for x, e in zip(x_crpa, eps_fb)] + b_fb = [dot(e, y) for y, e in zip(y_crpa, eps_fb)] + no = fragment.cluster.nocc_active + if isinstance(no, int): + no = (no, no) + nv = fragment.cluster.nvir_active + if isinstance(nv, int): + nv = (nv, nv) + l_a[:, :no[0], no[0]:] += a_fb[0].T.reshape((a_fb[0].shape[-1], no[0], nv[0])) + l_b[:, :no[1], no[1]:] += a_fb[1].T.reshape((a_fb[1].shape[-1], no[1], nv[1])) + + l_a[:, no[0]:, :no[0]] += b_fb[0].T.reshape((b_fb[0].shape[-1], no[0], nv[0])).transpose(0,2,1) + l_b[:, no[1]:, :no[1]] += b_fb[1].T.reshape((b_fb[1].shape[-1], no[1], nv[1])).transpose(0,2,1) + + if only_ov_screened: + # Zero out all contributions screening oo or vv contributions. + no = fragment.cluster.nocc_active + if isinstance(no, int): + no = (no, no) + l_a[:, no[0]:, no[0]:] = 0.0 + l_a[:, :no[0], :no[0]] = 0.0 + l_b[:, no[1]:, no[1]:] = 0.0 + l_b[:, :no[1], :no[1]] = 0.0 + return l_a, l_b, crpa + + +def get_crpa(orig_mf, f, log): + + def construct_loc_rot(f): + """Constructs the rotation of the overall mean-field space into which """ + ro = f.get_overlap("cluster[occ]|mo[occ]") + rv = f.get_overlap("cluster[vir]|mo[vir]") + + if isinstance(ro, np.ndarray): + ro = (ro, ro) + if isinstance(rv, np.ndarray): + rv = (rv, rv) + + rot_ova = einsum("Ij,Ab->IAjb", ro[0], rv[0]) + rot_ova = rot_ova.reshape((rot_ova.shape[0] * rot_ova.shape[1], -1)) + + rot_ovb = einsum("Ij,Ab->IAjb", ro[1], rv[1]) + rot_ovb = rot_ovb.reshape((rot_ovb.shape[0] * rot_ovb.shape[1], -1)) + return rot_ova, rot_ovb + + rot_loc = construct_loc_rot(f) + rot_ov = scipy.linalg.null_space(rot_loc[0]).T, scipy.linalg.null_space(rot_loc[1]).T + if rot_ov[0].shape[0] == 0 and rot_ov[1].shape[0] == 0: + log.warning("cRPA space contains no excitations! Interactions will be unscreened.") + raise cRPAError("cRPA space contains no excitations!") + # RPA calculation and new_mf will contain all required information for the response. + crpa = ssRPA(orig_mf, ov_rot=rot_ov) + crpa.kernel() + return crpa, rot_loc, rot_ov + + +def ao2mo(mf, mo_coeff=None, ijslice=None): + """Get MO basis density-fitted integrals. + """ + + if mo_coeff is None: + mo_coeff = mf.mo_coeff + nmo = mo_coeff.shape[1] + naux = mf.with_df.get_naoaux() + mem_incore = (2 * nmo ** 2 * naux) * 8 / 1e6 + mem_now = lib.current_memory()[0] + + mo = np.asarray(mo_coeff, order='F') + if ijslice is None: + ijslice = (0, nmo, 0, nmo) + + finshape = (naux, ijslice[1] - ijslice[0], ijslice[3] - ijslice[2]) + + Lpq = None + if (mem_incore + mem_now < 0.99 * mf.max_memory) or mf.mol.incore_anyway: + Lpq = _ao2mo.nr_e2(mf.with_df._cderi, mo, ijslice, aosym='s2', out=Lpq) + return Lpq.reshape(*finshape) + else: + logger.warn(mf, 'Memory may not be enough!') + raise NotImplementedError diff --git a/vayesta/core/qemb/scrcoulomb.py b/vayesta/core/screening/screening_moment.py similarity index 78% rename from vayesta/core/qemb/scrcoulomb.py rename to vayesta/core/screening/screening_moment.py index dacbefbc3..5ad5844f6 100644 --- a/vayesta/core/qemb/scrcoulomb.py +++ b/vayesta/core/screening/screening_moment.py @@ -2,12 +2,12 @@ import numpy as np import scipy import scipy.linalg -from vayesta.rpa import ssRIRPA +from vayesta.rpa import ssRIRPA, ssRPA from vayesta.core.util import dot, einsum from vayesta.mpi import mpi -def build_screened_eris(emb, fragments=None, cderi_ov=None, calc_e=True, npoints=48, log=None): +def build_screened_eris(emb, fragments=None, cderi_ov=None, store_m0=True, npoints=48, log=None): """Generates renormalised coulomb interactions for use in local cluster calculations. Currently requires unrestricted system. @@ -21,8 +21,8 @@ def build_screened_eris(emb, fragments=None, cderi_ov=None, calc_e=True, npoints cderi_ov : np.array or tuple of np.array, optional. Cholesky-decomposed ERIs in the particle-hole basis of mf. If mf is unrestricted this should be a list of arrays corresponding to the different spin channels. - calc_ecorrection : bool, optional. - Whether to calculate a nonlocal energy correction at the level of RPA + store_m0 : bool, optional. + Whether to store the local zeroth moment in the fragment class for use later. npoints : int, optional Number of points for numerical integration. Default: 48. log : logging.Logger, optional @@ -45,51 +45,35 @@ def build_screened_eris(emb, fragments=None, cderi_ov=None, calc_e=True, npoints if fragments is None: fragments = emb.get_fragments(active=True, sym_parent=None, mpi_rank=mpi.rank) fragments = [f for f in fragments if f.opts.screening == 'mrpa'] - if emb.spinsym != 'unrestricted': - raise NotImplementedError("Screened interactions require a spin-unrestricted formalism.") if emb.df is None: raise NotImplementedError("Screened interactions require density-fitting.") r_occs = [f.get_overlap('mo[occ]|cluster[occ]') for f in fragments] r_virs = [f.get_overlap('mo[vir]|cluster[vir]') for f in fragments] target_rots, ovs_active = _get_target_rot(r_occs, r_virs) - rpa = ssRIRPA(emb.mf, log=log, Lpq=cderi_ov) - if calc_e: - # This scales as O(N^4) - erpa, energy_error = rpa.kernel_energy(correction='linear') - else: - erpa = None - - tr = np.concatenate(target_rots, axis=0) - if sum(sum(ovs_active)) > 0: - # Computation scales as O(N^4) - moms_interact, est_errors = rpa.kernel_moms(0, tr, npoints=npoints) - momzero_interact = moms_interact[0] - else: - momzero_interact = np.zeros_like(np.concatenate(tr, axis=0)) - - # Now need to separate into local contributions - n = 0 - local_moments = [] - for nov, rot in zip(ovs_active, target_rots): - # Computation costs O(N^2 N_clus^2) - # Get corresponding section of overall moment, then project to just local contribution. - local_moments += [dot(momzero_interact[n:n+sum(nov)], rot.T)] - n += sum(nov) + local_moments = calc_moms_RIRPA(emb.mf, target_rots, ovs_active, log, cderi_ov, npoints) + # Could generate moments using N^6 moments instead, but just for debugging. + #local_moments, erpa = calc_moms_RPA(emb.mf, target_rots, ovs_active, log, cderi_ov, calc_e, npoints) # Then construct the RPA coupling matrix A-B, given by the diagonal matrix of energy differences. no = np.array(sum(emb.mf.mo_occ.T > 0)) - norb = np.array(emb.mo_coeff).shape[1] + if no.size == 1: no = np.array([int(no), int(no)]) + norb = emb.mo_coeff[0].shape[0] nv = norb - no + mo_e = emb.mf.mo_energy + + if isinstance(mo_e[0], float): + mo_e = np.array([mo_e, mo_e]) + def get_eps_singlespin(no_, nv_, mo_energy): eps = np.zeros((no_, nv_)) eps = eps + mo_energy[no_:] eps = (eps.T - mo_energy[:no_]).T eps = eps.reshape(-1) return eps - eps = np.concatenate([get_eps_singlespin(no[0], nv[0], emb.mf.mo_energy[0]), - get_eps_singlespin(no[1], nv[1], emb.mf.mo_energy[1])]) + eps = np.concatenate([get_eps_singlespin(no[0], nv[0], mo_e[0]), + get_eps_singlespin(no[1], nv[1], mo_e[1])]) # And use this to perform inversion to calculate interaction in cluster. seris_ov = [] @@ -97,7 +81,11 @@ def get_eps_singlespin(no_, nv_, mo_energy): amb = einsum("pn,qn,n->pq", rot, rot, eps) # O(N^2 N_clus^4) # Everything from here on is independent of system size, scaling at most as O(N_clus^6) # (arrays have side length equal to number of cluster single-particle excitations). - mominv = np.linalg.inv(mom) + e, c = np.linalg.eigh(mom) + if min(e) < 1e-4: + log.warning("Small eigenvalue of local rpa moment in %s: %e", f.name, min(e)) + + mominv = einsum("pn,n,qn->pq", c, e**(-1), c) apb = dot(mominv, amb, mominv) # This is the renormalised coulomb kernel in the cluster. @@ -108,23 +96,58 @@ def get_eps_singlespin(no_, nv_, mo_energy): # indices. no = f.cluster.nocc_active nv = f.cluster.nvir_active + if isinstance(no, int): + no = (no, no) + nv = (nv, nv) + kcaa = kc[:ova, :ova].reshape((no[0], nv[0], no[0], nv[0])) kcab = kc[:ova, ova:].reshape((no[0], nv[0], no[1], nv[1])) kcbb = kc[ova:, ova:].reshape((no[1], nv[1], no[1], nv[1])) - if kcaa.shape == kcbb.shape: - # This is not that meaningful, since the alpha and beta basis themselves may be different: - log.info("Screened interations in %s: spin-symmetry= %.3e spin-dependence= %.3e", - f.name, abs(kcaa-kcbb).max(), abs((kcaa+kcbb)/2-kcab).max()) kc = (kcaa, kcab, kcbb) - f._seris_ov = kc + f._seris_ov = (kc, mom, amb) if store_m0 else (kc,) seris_ov.append(kc) - return seris_ov, erpa + +def calc_moms_RIRPA(mf, target_rots, ovs_active, log, cderi_ov, npoints): + rpa = ssRIRPA(mf, log=log, Lpq=cderi_ov) + + tr = np.concatenate(target_rots, axis=0) + if sum(sum(ovs_active)) > 0: + # Computation scales as O(N^4) + moms_interact, est_errors = rpa.kernel_moms(0, tr, npoints=npoints) + momzero_interact = moms_interact[0] + else: + momzero_interact = np.zeros_like(np.concatenate(tr, axis=0)) + + # Now need to separate into local contributions + n = 0 + local_moments = [] + for nov, rot in zip(ovs_active, target_rots): + # Computation costs O(N^2 N_clus^2) + # Get corresponding section of overall moment, then project to just local contribution. + mom = dot(momzero_interact[n:n+sum(nov)], rot.T) + # This isn't exactly symmetric due to numerical integration, so enforce here. + mom = (mom + mom.T) / 2 + local_moments += [mom] + n += sum(nov) + + return local_moments + +def calc_moms_RPA(mf, target_rots, ovs_active, log, cderi_ov, calc_e, npoints): + rpa = ssRPA(mf, log=log) + erpa = rpa.kernel() + mom0 = rpa.gen_moms(0)[0] + local_moments = [] + for rot in target_rots: + local_moments += [dot(rot, mom0, rot.T)] + return local_moments, erpa def get_screened_eris_full(eris, seris_ov, copy=True, log=None): """Build full array of screened ERIs, given the bare ERIs and screening.""" + log.info("Generating screened interaction to conserve zeroth moment of the dd response.") + def replace_ov(full, ov, spins): out = full.copy() if copy else full no1, no2 = ov.shape[0], ov.shape[2] @@ -134,12 +157,11 @@ def replace_ov(full, ov, spins): out[v1,o1,o2,v2] = ov.transpose([1, 0, 2, 3]) out[o1,v1,v2,o2] = ov.transpose([0, 1, 3, 2]) out[v1,o1,v2,o2] = ov.transpose([1, 0, 3, 2]) - if log: - maxidx = np.unravel_index(np.argmax(abs(full-out)), full.shape) - log.info("Maximally screened element of W(%2s|%2s): V= %.3e -> W= %.3e (delta= %.3e)", - 2*spins[0], 2*spins[1], full[maxidx], out[maxidx], out[maxidx]-full[maxidx]) return out + if isinstance(eris, np.ndarray): + eris = (eris, eris, eris) + seris = (replace_ov(eris[0], seris_ov[0], 'aa'), replace_ov(eris[1], seris_ov[1], 'ab'), replace_ov(eris[2], seris_ov[2], 'bb')) @@ -226,6 +248,10 @@ def get_target_rot_spat(ro, rv): for i, (r_o, r_v) in enumerate(zip(r_active_occs, r_active_virs)): + if isinstance(r_o, np.ndarray): + r_o = (r_o, r_o) + r_v = (r_v, r_v) + arot, ova = get_target_rot_spat(r_o[0], r_v[0]) brot, ovb = get_target_rot_spat(r_o[1], r_v[1]) diff --git a/vayesta/dmet/dmet.py b/vayesta/dmet/dmet.py index 911268450..131b6abd0 100644 --- a/vayesta/dmet/dmet.py +++ b/vayesta/dmet/dmet.py @@ -128,11 +128,10 @@ def kernel(self): if type(nelec_mf) == tuple: nelec_mf = sum(nelec_mf) - if self.opts.screening == 'mrpa': - for f in self.get_fragments(sym_parent=None): - f.make_bath() - f.make_cluster() - self.build_screened_eris() + for f in self.get_fragments(sym_parent=None): + f.make_bath() + f.make_cluster() + self.build_screened_eris() def electron_err(cpt, construct_bath=False): err = self.calc_electron_number_defect(cpt, nelec_mf, sym_parents, nsym, construct_bath) diff --git a/vayesta/dmet/fragment.py b/vayesta/dmet/fragment.py index ea57d83b9..cff53c87b 100644 --- a/vayesta/dmet/fragment.py +++ b/vayesta/dmet/fragment.py @@ -67,6 +67,7 @@ def kernel(self, solver=None, init_guess=None, seris_ov=None, construct_bath=Tru return None cluster_solver = self.get_solver(solver) + # Chemical potential if chempot is not None: px = self.get_fragment_projector(self.cluster.c_active) diff --git a/vayesta/edmet/edmet.py b/vayesta/edmet/edmet.py index 5e2dda8e9..6e73338ba 100644 --- a/vayesta/edmet/edmet.py +++ b/vayesta/edmet/edmet.py @@ -54,6 +54,14 @@ def eps(self): eps = eps.reshape(-1) return eps, eps + @property + def e_nonlocal(self): + return self._e_nonlocal or 0.0 + + @e_nonlocal.setter + def e_nonlocal(self, value): + self._e_nonlocal = value + def check_solver(self, solver): is_uhf = np.ndim(self.mo_coeff[1]) == 2 is_eb = True @@ -375,4 +383,8 @@ def run_exact_full_ac(self, xc_kernel=None, deg=5, calc_local=False, cluster_con else: raise NotImplementedError + def _reset(self): + super()._reset() + self._e_nonlocal = None + REDMET = EDMET diff --git a/vayesta/ewf/ewf.py b/vayesta/ewf/ewf.py index cf81c3c05..adabbd690 100644 --- a/vayesta/ewf/ewf.py +++ b/vayesta/ewf/ewf.py @@ -27,6 +27,8 @@ class Options(Embedding.Options): # --- Bath settings bath_options: dict = Embedding.Options.change_dict_defaults('bath_options', bathtype='mp2', threshold=1e-8) + solver_options: dict = Embedding.Options.change_dict_defaults('solver_options', + fermion_wf=True) #ewdmet_max_order: int = 1 # If multiple bno thresholds are to be calculated, we can project integrals and amplitudes from a previous larger cluster: project_eris: bool = False # Project ERIs from a pervious larger cluster (corresponding to larger eta), can result in a loss of accuracy especially for large basis sets! @@ -139,12 +141,7 @@ def kernel(self): self.communicate_clusters() # --- Screened Coulomb interaction - if any(x.opts.screening is not None for x in fragments): - self.log.info("") - self.log.info("SCREENING INTERACTIONS") - self.log.info("======================") - with log_time(self.log.timing, "Time for screened interations: %s"): - self.build_screened_eris() + self.build_screened_eris() # --- Loop over fragments with no symmetry parent and with own MPI rank self.log.info("") @@ -294,6 +291,7 @@ def _make_rdm2_ccsd_proj_lambda(self, *args, **kwargs): def get_e_corr(self, functional=None, **kwargs): functional = (functional or self.opts.energy_functional) + if functional == 'projected': self.log.warning("functional='projected' is deprecated; use functional='wf' instead.") functional = 'wf' diff --git a/vayesta/ewf/fragment.py b/vayesta/ewf/fragment.py index e9aacfffc..e0de257f3 100644 --- a/vayesta/ewf/fragment.py +++ b/vayesta/ewf/fragment.py @@ -190,7 +190,8 @@ def kernel(self, solver=None, init_guess=None): # Create solver object cluster_solver = self.get_solver(solver) - + # Calculate cluster energy at the level of RPA. + e_corr_rpa = self.get_local_rpa_correction(cluster_solver.hamil) # --- Chemical potential cpt_frag = self.base.opts.global_frag_chempot if self.opts.nelectron_target is not None: @@ -253,7 +254,7 @@ def kernel(self, solver=None, init_guess=None): # --- Add to results data class self._results = results = self.Results(fid=self.id, n_active=cluster.norb_active, - converged=cluster_solver.converged, wf=wf, pwf=pwf, moms=moms) + converged=cluster_solver.converged, wf=wf, pwf=pwf, moms=moms, e_corr_rpa=e_corr_rpa) self.hamil = cluster_solver.hamil diff --git a/vayesta/rpa/ssrpa.py b/vayesta/rpa/ssrpa.py index 4dbb0d152..0dd042d00 100644 --- a/vayesta/rpa/ssrpa.py +++ b/vayesta/rpa/ssrpa.py @@ -94,8 +94,13 @@ def kernel(self, xc_kernel=None, alpha=1.0): AmBrtinv = einsum("pn,n,qn->pq", c2, e2 ** (-0.5), c2) XpY = np.einsum("n,pq,qn->pn", self.freqs_ss ** (-0.5), AmBrt, c) XmY = np.einsum("n,pq,qn->pn", self.freqs_ss ** (0.5), AmBrtinv, c) - self.XpY_ss = (XpY[: self.ova], XpY[self.ova :]) - self.XmY_ss = (XmY[: self.ova], XmY[self.ova :]) + + nov = self.ova + if self.ov_rot is not None: + nov = self.ov_rot[0].shape[0] + + self.XpY_ss = (XpY[:nov], XpY[nov:]) + self.XmY_ss = (XmY[:nov], XmY[nov:]) self.freqs_sf = None self.log.timing("Time to solve RPA problem: %s", time_string(timer() - t0)) @@ -199,8 +204,7 @@ def get_val_alpha(alpha): return e - def _gen_arrays(self, xc_kernel=None, alpha=1.0): - t0 = timer() + def _gen_eps(self): # Only have diagonal components in canonical basis. eps = np.zeros((self.nocc, self.nvir)) eps = eps + self.mf.mo_energy[self.nocc :] @@ -208,18 +212,30 @@ def _gen_arrays(self, xc_kernel=None, alpha=1.0): eps = eps.reshape((self.ova,)) if self.ov_rot is not None: - eps = einsum("pn,n,qn->pq", self.ov_rot, eps, self.ov_rot) - # want to ensure diagonalise eps in our new basis; nb this also rotates the existing rotation for - # consistency. - eps, c = scipy.linalg.eigh(eps) - self.ov_rot = dot(c.T, self.ov_rot) + if self.ov_rot[0].shape[0] > 0: + epsa = einsum("pn,n,qn->pq", self.ov_rot[0], eps, self.ov_rot[0]) + epsa, ca = scipy.linalg.eigh(epsa) + self.ov_rot = (dot(ca.T, self.ov_rot[0]), self.ov_rot[1]) + if self.ov_rot[1].shape[0] > 0: + epsb = einsum("pn,n,qn->pq", self.ov_rot[1], eps, self.ov_rot[1]) + epsb, cb = scipy.linalg.eigh(epsb) + self.ov_rot = (self.ov_rot[0], dot(cb.T, self.ov_rot[1])) + else: + epsa = epsb = eps + return epsa, epsb + + + def _gen_arrays(self, xc_kernel=None, alpha=1.0): + t0 = timer() + + epsa, epsb = self._gen_eps() - AmB = np.concatenate([eps, eps]) + AmB = np.concatenate([epsa, epsb]) fullv = self.get_k() ApB = 2 * fullv * alpha if self.ov_rot is not None: - fullrot = scipy.linalg.block_diagonal(self.ov_rot, self.ov_rot) + fullrot = scipy.linalg.block_diag(self.ov_rot[0], self.ov_rot[1]) ApB = dot(fullrot, ApB, fullrot.T) # At this point AmB is just epsilon so add in. @@ -239,7 +255,7 @@ def _gen_arrays(self, xc_kernel=None, alpha=1.0): M = dot(AmBrt, ApB, AmBrt) self.log.timing("Time to build RPA arrays: %s", time_string(timer() - t0)) - return M, AmB, ApB, (eps, eps), fullv + return M, AmB, ApB, (epsa, epsb), fullv def get_k(self): eris = self.ao2mo() @@ -319,30 +335,21 @@ def get_xc_contribs(self, xc_kernel, c_o, c_v, alpha=1.0): return ApB * alpha, AmB * alpha def gen_moms(self, max_mom, xc_kernel=None): - res = np.zeros((max_mom + 1, self.ov, self.ov)) + nova = self.ova + nov = self.ov + if self.ov_rot is not None: + nova = self.ov_rot[0].shape[0] + nov = nova + self.ov_rot[1].shape[0] + res = np.zeros((max_mom+1, nov, nov)) for x in range(max_mom + 1): # Have different spin components in general; these are alpha-alpha, alpha-beta and beta-beta. - res[x, : self.ova, : self.ova] = np.einsum( - "pn,n,qn->pq", self.XpY_ss[0], self.freqs_ss**x, self.XpY_ss[0] - ) - res[x, self.ova :, self.ova :] = np.einsum( - "pn,n,qn->pq", self.XpY_ss[1], self.freqs_ss**x, self.XpY_ss[1] - ) - res[x, : self.ova, self.ova :] = np.einsum( - "pn,n,qn->pq", self.XpY_ss[0], self.freqs_ss**x, self.XpY_ss[1] - ) - res[x, self.ova :, : self.ova] = res[x][: self.ova, self.ova :].T + res[x, :nova, :nova] = np.einsum("pn,n,qn->pq", self.XpY_ss[0], self.freqs_ss ** x, self.XpY_ss[0]) + res[x, nova:, nova:] = np.einsum("pn,n,qn->pq", self.XpY_ss[1], self.freqs_ss ** x, self.XpY_ss[1]) + res[x, :nova, nova:] = np.einsum("pn,n,qn->pq", self.XpY_ss[0], self.freqs_ss ** x, self.XpY_ss[1]) + res[x, nova:, :nova] = res[x][:nova, nova:].T return res - @property - def mo_coeff(self): - return self.mf.mo_coeff - - @property - def nao(self): - return self.mf.mol.nao - def ao2mo(self, mo_coeff=None, compact=False): """Get the ERIs in MO basis""" mo_coeff = self.mo_coeff if mo_coeff is None else mo_coeff diff --git a/vayesta/rpa/ssurpa.py b/vayesta/rpa/ssurpa.py index ec6a62383..48a353e20 100644 --- a/vayesta/rpa/ssurpa.py +++ b/vayesta/rpa/ssurpa.py @@ -3,7 +3,7 @@ import scipy.linalg from timeit import default_timer as timer -from vayesta.core.util import dot, time_string +from vayesta.core.util import dot, time_string, einsum class ssURPA(ssRPA): @@ -43,9 +43,7 @@ def mo_coeff_vir(self): na, nb = self.nocc return self.mo_coeff[0][:, na:], self.mo_coeff[1][:, nb:] - def _gen_arrays(self, xc_kernel=None, alpha=1.0): - t0 = timer() - + def _gen_eps(self): nocc_a, nocc_b = self.nocc nvir_a, nvir_b = self.nvir # Only have diagonal components in canonical basis. @@ -58,6 +56,11 @@ def _gen_arrays(self, xc_kernel=None, alpha=1.0): epsb = epsb + self.mf.mo_energy[1][nocc_b:] epsb = (epsb.T - self.mf.mo_energy[1][:nocc_b]).T epsb = epsb.reshape((self.ovb,)) + return epsa, epsb + + def _gen_arrays(self, xc_kernel=None, alpha=1.0): + t0 = timer() + epsa, epsb = self._gen_eps() if self.ov_rot is not None: epsa = einsum("pn,n,qn->pq", self.ov_rot[0], epsa, self.ov_rot[0]) @@ -78,7 +81,7 @@ def _gen_arrays(self, xc_kernel=None, alpha=1.0): fullv = self.get_k() ApB = 2 * fullv * alpha if self.ov_rot is not None: - fullrot = scipy.linalg.block_diagonal(self.ov_rot[0], self.ov_rot[1]) + fullrot = scipy.linalg.block_diag(self.ov_rot[0], self.ov_rot[1]) ApB = dot(fullrot, ApB, fullrot.T) # At this point AmB is just epsilon so add in. diff --git a/vayesta/solver/ebcc.py b/vayesta/solver/ebcc.py index 93a72938e..4ca386a82 100644 --- a/vayesta/solver/ebcc.py +++ b/vayesta/solver/ebcc.py @@ -172,9 +172,12 @@ class EB_REBCC_Solver(REBCC_Solver): @dataclasses.dataclass class Options(REBCC_Solver.Options): ansatz: str = "CCSD-S-1-1" + fermion_wf: bool = False @property def is_fCCSD(self): + if self.opts.fermion_wf: + return super().is_fCCSD return False def get_nonnull_solver_opts(self): diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 2b7832230..153ddd9e8 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -4,12 +4,11 @@ import numpy as np import pyscf.lib import pyscf.scf -import scipy.linalg -from vayesta.core.qemb import scrcoulomb +from vayesta.core.util import dot, einsum, OptionsBase, break_into_lines, log_time, energy_string +from vayesta.core.screening import screening_moment, screening_crpa from vayesta.core.types import Orbitals -from vayesta.core.util import dot, einsum, OptionsBase, break_into_lines, log_time -from vayesta.rpa import ssRPA +from typing import Optional def is_ham(ham): @@ -60,6 +59,7 @@ class RClusterHamiltonian: class Options(OptionsBase): screening: Optional[str] = None cache_eris: bool = True + match_fock: bool = True @property def _scf_class(self): @@ -112,7 +112,11 @@ def get_integrals(self, bare_eris=None, with_vext=True): seris = self.get_eris_screened() return heff, seris - def get_fock(self, with_vext=True, use_seris=True, with_exxdiv=False): + def get_fock(self, with_vext=True, use_seris=None, with_exxdiv=False): + if use_seris is None: + # Actual default depends on self.opts.match_fock; to reproduce fock we will effectively modify heff + use_seris = not self.opts.match_fock + c = self.cluster.c_active fock = dot(c.T, self._fragment.base.get_fock(with_exxdiv=with_exxdiv), c) if with_vext and self.v_ext is not None: @@ -133,9 +137,9 @@ def get_fock(self, with_vext=True, use_seris=True, with_exxdiv=False): def get_heff(self, eris=None, fock=None, with_vext=True, with_exxdiv=False): if eris is None: - eris = self.get_eris_bare() + eris = self.get_eris_screened() if fock is None: - fock = self.get_fock(with_vext=False, use_seris=False, with_exxdiv=with_exxdiv) + fock = self.get_fock(with_vext=False, with_exxdiv=with_exxdiv) occ = np.s_[:self.cluster.nocc_active] v_act = 2 * einsum('iipq->pq', eris[occ, occ]) - einsum('iqpi->pq', eris[occ, :, :, occ]) h_eff = fock - v_act @@ -325,17 +329,20 @@ def pad(a, diag_val): else: # Just set up heff using standard bare eris. bare_eris = self.get_eris_bare() - heff = pad_to_match(self.get_heff(eris=bare_eris, with_vext=True), dummy_energy) if force_bare_eris: clusmf._eri = pad_to_match(bare_eris, 0.0) else: clusmf._eri = pad_to_match(self.get_eris_screened()) + heff = pad_to_match(self.get_heff(with_vext=True), dummy_energy) + + clusmf.get_hcore = lambda *args, **kwargs: heff if overwrite_fock: - clusmf.get_fock = lambda *args, **kwargs: pad_to_match(self.get_fock(with_vext=True), dummy_energy) - clusmf.get_veff = lambda *args, **kwargs: np.array(clusmf.get_fock(*args, **kwargs)) - np.array( - clusmf.get_hcore()) + clusmf.get_fock = lambda *args, **kwargs: pad_to_match( + self.get_fock(with_vext=True, use_seris=not force_bare_eris), dummy_energy) + clusmf.get_veff = lambda *args, **kwargs: np.array(clusmf.get_fock(*args, **kwargs)) - \ + np.array(clusmf.get_hcore()) return clusmf, orbs_to_freeze @@ -358,48 +365,145 @@ def get_clus_mf_info(self, ao_basis=False, with_vext=True, with_exxdiv=False): # Functionality for use with screened interactions and external corrections. - def add_screening(self, seris_intermed=None): - """Add screened interactions into the Hamiltonian.""" - if self.opts.screening == "mrpa": - assert (seris_intermed is not None) + def calc_loc_erpa(self, m0, amb, only_cumulant=False): - # Use bare coulomb interaction from hamiltonian; this could well be cached in future. - bare_eris = self.get_eris_bare() + no, nv = self.cluster.nocc_active, self.cluster.nvir_active + nov = no * nv + # Bare coulomb interaction in cluster ov-ov space. + v = self.get_eris_bare()[:no, no:, :no, no:].reshape((nov, nov)) + ro = self._fragment.get_overlap("fragment|cluster-occ") + po = dot(ro.T, ro) + + def gen_spin_components(mat): + return mat[:nov, :nov], mat[:nov, nov:], mat[nov:, nov:] + + m0_aa, m0_ab, m0_bb = gen_spin_components(m0) - self._seris = scrcoulomb.get_screened_eris_full(bare_eris, seris_intermed) + if only_cumulant: - elif self.opts.screening == "crpa": - raise NotImplementedError() + def compute_e_rrpa(proj): + def pr(m): + m = m.reshape((no, nv, no * nv)) + m = np.tensordot(proj, m, axes=[0, 0]) + return m.reshape((no * nv, no * nv)) + + erpa = 0.5 * einsum("pq,qp->", pr(m0_aa + m0_ab + m0_ab.T + m0_bb - 2 * np.eye(nov)), v) + self.log.info("Computed fragment RPA cumulant energy contribution for cluster %s as %s", + self._fragment.id, + energy_string(erpa)) + return erpa else: - raise ValueError("Unknown cluster screening protocol: %s" % self.opts.screening) + d_aa, d_ab, d_bb = gen_spin_components(amb) + # This should be zero. + assert (abs(d_ab).max() < 1e-10) + + def compute_e_rrpa(proj): + def pr(m): + m = m.reshape((no, nv, no * nv)) + m = np.tensordot(proj, m, axes=[0, 0]) + return m.reshape((no * nv, no * nv)) + + erpa = 0.5 * (einsum("pq,qp->", pr(m0_aa), d_aa) + einsum("pq,qp->", pr(m0_bb), d_bb)) + erpa += einsum("pq,qp->", pr(m0_aa + m0_ab + m0_ab.T + m0_bb), v) + erpa -= 0.5 * (pr(d_aa + v + d_bb + v).trace()) + self.log.info("Computed fragment RPA energy contribution for cluster %s as %s", self._fragment.id, + energy_string(erpa)) + return erpa + + compute_e_rrpa(np.eye(no)) + + return compute_e_rrpa(po) + + def add_screening(self, seris_intermed=None): + """ + Adds appropriate screening according to the value of self.opts.screening. + -`None`: gives bare interactions, but this function shouldn't be called in that case. + -'mrpa': moment-conserving interactions. + -'crpa': gives cRPA interactions. Including 'ov' after 'crpa' will only apply cRPA screening in the o-v channel. + Including 'pcoupling' similarly will account for the polarisability in non-canonical cluster spaces. - def calc_loc_erpa(self): + seris_intermed is only required for mRPA interactions. + """ + self._seris = self._add_screening(seris_intermed, spin_integrate=True) - clusmf = self.to_pyscf_mf(force_bare_eris=True) - clusrpa = ssRPA(clusmf) - M, AmB, ApB, eps, v = clusrpa._gen_arrays() - erpa = clusrpa.kernel() - m0 = clusrpa.gen_moms(0) + def _add_screening(self, seris_intermed=None, spin_integrate=True): - def get_product_projector(): - nocc = self.nelec - nvir = tuple([x - y for x, y in zip(self.ncas, self.nelec)]) - p_occ_frag = self.target_space_projector(self.cluster.c_active_occ) - if (not isinstance(p_occ_frag, tuple)) and np.ndim(p_occ_frag) == 2: - p_occ_frag = (p_occ_frag, p_occ_frag) + def spin_integrate_and_report(m, warn_threshold=1e-6): + spat = (m[0] + m[1] + m[2] + m[1].transpose((2, 3, 0, 1))) / 4.0 - def get_product_projector(p_o, p_v, no, nv): - return einsum("ij,ab->iajb", p_o, p_v).reshape((no * nv, no * nv)) + dev = [abs(x - spat).max() for x in m] + [abs(m[2].transpose(2, 3, 0, 1) - spat).max()] + self.log.info("Largest change in screened interactions due to spin integration: %e", max(dev)) + if max(dev) > warn_threshold: + self.log.warning("Significant change in screened interactions due to spin integration: %e", max(dev)) + return spat - pa = get_product_projector(p_occ_frag[0], np.eye(nvir[0]), nocc[0], nvir[0]) - pb = get_product_projector(p_occ_frag[1], np.eye(nvir[1]), nocc[1], nvir[1]) - return scipy.linalg.block_diag(pa, pb) + if self.opts.screening is None: + raise ValueError("Attempted to add screening to fragment with no screening protocol specified.") + if self.opts.screening == "mrpa": + assert (seris_intermed is not None) + # Use bare coulomb interaction from hamiltonian. + bare_eris = self.get_eris_bare() + seris = screening_moment.get_screened_eris_full(bare_eris, seris_intermed[0], log=self.log) + if spin_integrate: + seris = spin_integrate_and_report(seris) + elif self.opts.screening[:4] == "crpa": + bare_eris = self.get_eris_bare() + delta, crpa = screening_crpa.get_frag_deltaW(self.orig_mf, self._fragment, + pcoupling=("pcoupled" in self.opts.screening), + only_ov_screened=("ov" in self.opts.screening), + log=self.log) + if "store" in self.opts.screening: + self.log.warning("Storing cRPA object in Hamiltonian- O(N^4) memory cost!") + self.crpa = crpa + if "full" in self.opts.screening: + # Add a check just in case. + self.log.critical("Static screening of frequency-dependent interactions not supported") + self.log.critical("This statement should be impossible to reach!") + raise ValueError("Static screening of frequency-dependent interactions not supported") + else: + if spin_integrate: + delta = spin_integrate_and_report(delta) + seris = bare_eris + delta + else: + seris = tuple([x + y for x, y in zip(bare_eris, delta)]) + else: + raise ValueError("Unknown cluster screening protocol: %s" % self.opts.screening) - proj = get_product_projector() - eloc = 0.5 * einsum("pq,qr,rp->", proj, m0, ApB) - einsum("pq,qp->", proj, ApB + AmB) - return eloc + def report_screening(screened, bare, spins): + maxidx = np.unravel_index(np.argmax(abs(screened - bare)), bare.shape) + if spins is None: + wstring = "W" + else: + wstring = "W(%2s|%2s)" % (2 * spins[0], 2 * spins[1]) + self.log.info( + "Maximally screened element of %s: V= %.3e -> W= %.3e (delta= %.3e)", + wstring, bare[maxidx], screened[maxidx], screened[maxidx] - bare[maxidx]) + # self.log.info( + # " Corresponding norms%s: ||V||= %.3e, ||W||= %.3e, ||delta||= %.3e", + # " " * len(wstring), np.linalg.norm(bare), np.linalg.norm(screened), + # np.linalg.norm(screened-bare)) + + if spin_integrate: + report_screening(seris, bare_eris, None) + else: + report_screening(seris[0], bare_eris[0], "aa") + report_screening(seris[1], bare_eris[1], "ab") + report_screening(seris[2], bare_eris[2], "bb") + + def get_sym_breaking(norm_aa, norm_ab, norm_bb): + spinsym = abs(norm_aa - norm_bb) / ((norm_aa + norm_bb) / 2) + spindep = abs((norm_aa + norm_bb) / 2 - norm_ab) / ((norm_aa + norm_bb + norm_ab) / 3) + return spinsym, spindep + + bss, bsd = get_sym_breaking(*[np.linalg.norm(x) for x in bare_eris]) + sss, ssd = get_sym_breaking(*[np.linalg.norm(x) for x in seris]) + dss, dsd = get_sym_breaking(*[np.linalg.norm(x - y) for x, y in zip(bare_eris, seris)]) + + self.log.info("Proportional spin symmetry breaking in norms: V= %.3e, W= %.3e, (W-V= %.3e)", bss, sss, dss) + self.log.info("Proportional spin dependence in norms: V= %.3e, W= %.3e, (W-V= %.3e)", bsd, ssd, dsd) + return seris def assert_equal_spin_channels(self, message=""): na, nb = self.ncas @@ -613,6 +717,10 @@ def __contains__(self, item): fock = self.get_fock(with_vext=with_vext, use_seris=not force_bare, with_exxdiv=with_exxdiv) return DummyERIs(getter, valid_blocks=ValidUHFKeys(), fock=fock, nocc=self.cluster.nocc_active) + def add_screening(self, seris_intermed=None): + """Add screened interactions into the Hamiltonian.""" + self._seris = self._add_screening(seris_intermed, spin_integrate=False) + class EB_RClusterHamiltonian(RClusterHamiltonian): @dataclasses.dataclass @@ -669,13 +777,17 @@ def get_polaritonic_shifted_couplings(self): couplings = tuple([x - einsum("pq,n->npq", np.eye(x.shape[1]), temp) for x in self.unshifted_couplings]) if not np.allclose(couplings[0], couplings[1]): self.log.critical("Polaritonic shifted bosonic fermion-boson couplings break cluster spin symmetry; please" - "use an unrestricted formalism.") - raise RuntimeError() + " use an unrestricted formalism.") + raise RuntimeError("Polaritonic shifted bosonic fermion-boson couplings break cluster spin symmetry; please" + " use an unrestricted formalism.") return couplings[0] def get_eb_dm_polaritonic_shift(self, dm1): return (-einsum("n,pq->pqn", self.polaritonic_shift, dm1 / 2),) * 2 + def _add_screening(self, seris_intermed=None, spin_integrate=True): + return self.get_eris_bare() + class EB_UClusterHamiltonian(UClusterHamiltonian, EB_RClusterHamiltonian): @dataclasses.dataclass @@ -702,3 +814,6 @@ def get_polaritonic_shifted_couplings(self): def get_eb_dm_polaritonic_shift(self, dm1): return tuple([-einsum("n,pq->pqn", self.polaritonic_shift, x) for x in dm1]) + + def calc_loc_erpa(self, m0, amb): + raise NotImplementedError() diff --git a/vayesta/tests/ewf/test_rpa_correction.py b/vayesta/tests/ewf/test_rpa_correction.py new file mode 100644 index 000000000..63e9265db --- /dev/null +++ b/vayesta/tests/ewf/test_rpa_correction.py @@ -0,0 +1,47 @@ +import unittest +import numpy as np +import vayesta +import vayesta.ewf +from vayesta.core.util import cache +from vayesta.tests import testsystems +from vayesta.tests.common import TestCase + + +class Test_RPA_Corrections_Ethanol_RHF(TestCase): + + system = testsystems.ethanol_631g_df + + @property + def mf(self): + return self.system.rhf() + + def get_nl_energies(self, correction, bathtype="dmet"): + + emb = vayesta.ewf.EWF(self.mf, bath_options={"bathtype":bathtype}, solver="CCSD", + screening="mrpa", ext_rpa_correction=correction) + with emb.iao_fragmentation() as f: + f.add_all_atomic_fragments() + emb.kernel() + return emb.e_nonlocal, emb.e_rpa + + def test_rpa_correction(self): + enl, erpa = self.get_nl_energies("erpa") + self.assertAlmostEqual(erpa, -0.29138256715100397) + self.assertAlmostEqual(enl, -0.3087486792144427) + + def test_cumulant_correction(self): + enl, erpa = self.get_nl_energies("cumulant") + self.assertAlmostEqual(erpa, -0.5145262339186916) + self.assertAlmostEqual(enl, -0.31553721476368396) + +class Test_RPA_Corrections_complete(Test_RPA_Corrections_Ethanol_RHF): + """Tests with a complete bath in all clusters. This should give no nonlocal correction in any case.""" + system = testsystems.water_631g_df + + def test_rpa_correction(self): + enl, erpa = self.get_nl_energies("erpa", "full") + self.assertAlmostEqual(enl, 0.0) + + def test_cumulant_correction(self): + enl, erpa = self.get_nl_energies("cumulant", "full") + self.assertAlmostEqual(enl, 0.0) diff --git a/vayesta/tests/ewf/test_screening.py b/vayesta/tests/ewf/test_screening.py index 7c7a0401b..7b9aae7ba 100644 --- a/vayesta/tests/ewf/test_screening.py +++ b/vayesta/tests/ewf/test_screening.py @@ -10,7 +10,7 @@ class TestTwoElectron(TestCase): system = testsystems.h2_ccpvdz_df - e_ref = -1.123779303361342 + e_ref = {"mrpa":-1.123779303361342, "crpa":-1.1237769151822752} @classmethod def setUpClass(cls): @@ -23,26 +23,36 @@ def tearDownClass(cls): @classmethod @cache - def emb(cls, bno_threshold, solver): + def emb(cls, bno_threshold, solver, screening): emb = vayesta.ewf.EWF(cls.mf, bath_options=dict(threshold=bno_threshold), solver=solver, - screening='mrpa', solver_options=dict(conv_tol=1e-12)) + screening=screening, solver_options=dict(conv_tol=1e-12)) emb.kernel() return emb - def test_ccsd(self): - emb = self.emb(np.inf, 'CCSD') + def test_ccsd_mrpa(self): + emb = self.emb(np.inf, 'CCSD', 'mrpa') emb.kernel() - self.assertAllclose(emb.e_tot, self.e_ref) + self.assertAllclose(emb.e_tot, self.e_ref['mrpa']) - def test_fci(self): - emb = self.emb(np.inf, 'FCI') + def test_fci_mrpa(self): + emb = self.emb(np.inf, 'FCI', 'mrpa') emb.kernel() - self.assertAllclose(emb.e_tot, self.e_ref) + self.assertAllclose(emb.e_tot, self.e_ref['mrpa']) + + def test_ccsd_crpa(self): + emb = self.emb(np.inf, 'CCSD', 'crpa') + emb.kernel() + self.assertAllclose(emb.e_tot, self.e_ref['crpa']) + + def test_fci_crpa(self): + emb = self.emb(np.inf, 'FCI', 'crpa') + emb.kernel() + self.assertAllclose(emb.e_tot, self.e_ref['crpa']) class TestTwoHole(TestTwoElectron): system = testsystems.f2_sto6g_df - e_ref = -197.84155758368854 + e_ref = {"mrpa":-197.84155758368854, "crpa":-197.83928243962046} if __name__ == '__main__': print('Running %s' % __file__) diff --git a/vayesta/tests/testsystems.py b/vayesta/tests/testsystems.py index 1c2c07617..34b6f14af 100644 --- a/vayesta/tests/testsystems.py +++ b/vayesta/tests/testsystems.py @@ -363,6 +363,7 @@ def uhf(self): water_ccpvdz_df = TestMolecule(atom=molecules.water(), basis="cc-pvdz", auxbasis="cc-pvdz-jkfit") ethanol_ccpvdz = TestMolecule(atom=molecules.ethanol(), basis="cc-pvdz") +ethanol_631g_df = TestMolecule(atom=molecules.ethanol(), basis="6-31G", auxbasis="6-31G") lih_ccpvdz = TestMolecule(atom="Li 0 0 0; H 0 0 1.4", basis="cc-pvdz") lih_631g = TestMolecule(atom="Li 0 0 0; H 0 0 1.4", basis="6-31g")