From 9a3bfd0c4d8ae8f394a0da21931c4058152745df Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Wed, 14 Sep 2022 00:04:41 +0100 Subject: [PATCH 01/33] Started adding cRPA screened coulomb interaction. --- vayesta/core/qemb/qemb.py | 2 +- vayesta/core/screening/__init__.py | 0 vayesta/core/screening/screening_crpa.py | 187 ++++++++++++++++++ .../screening_moment.py} | 4 +- vayesta/solver/hamiltonian.py | 4 +- 5 files changed, 192 insertions(+), 5 deletions(-) create mode 100644 vayesta/core/screening/__init__.py create mode 100644 vayesta/core/screening/screening_crpa.py rename vayesta/core/{qemb/scrcoulomb.py => screening/screening_moment.py} (99%) diff --git a/vayesta/core/qemb/qemb.py b/vayesta/core/qemb/qemb.py index d28b89d13..e199042b1 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 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..5d8114643 --- /dev/null +++ b/vayesta/core/screening/screening_crpa.py @@ -0,0 +1,187 @@ +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__ + + +def get_crpa_deltaW(emb, fragments=None, calc_delta_e=True, log=None): + """Generates change in 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 + ---------- + emb : Embedding + Embedding instance. + fragments : list of vayesta.qemb.Fragment subclasses, optional + List of fragments for the calculation, used to define local interaction spaces. + If None, `emb.get_fragments(sym_parent=None)` is used. Default: None. + calc_delta_e : bool, optional. + Whether to calculate a nonlocal energy correction at the level of RPA + log : logging.Logger, optional + Logger object. If None, the logger of the `emb` object is used. Default: None. + + Returns + ------- + seris_ov : list of tuples of np.array + List of spin-dependent screened (ov|ov), for each fragment provided. + delta_e: float + Delta RPA correction computed as difference between full system RPA energy and + cluster correlation energies; currently only functional in CAS fragmentations. + """ + if fragments is None: + fragments = emb.get_fragments(sym_parent=None) + 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.") + + wc = [] + + for f in fragments: + # Compute cRPA response. + new_mf, crpa = get_crpa_chi(emb.mf, f) + + # Apply static approximation to interaction. + static_factor = crpa.freqs_ss ** (-1) + # Now need to calculate interactions. + + nmo = new_mf.mo_coeff.shape[1] + nocc = sum(new_mf.mo_occ.T > 0) + + Lov_aenv = ao2mo(new_mf, mo_coeff=new_mf.mo_coeff[0], ijslice=(0, nocc[0], nocc[0], nmo[0])).reshape( + (-1, crpa.ova)) + Lov_benv = ao2mo(new_mf, mo_coeff=new_mf.mo_coeff[1], ijslice=(0, nocc[1], nocc[1], nmo[1])).reshape( + (-1, crpa.ovb)) + + # This is the coefficient of the cRPA reducible dd response in the auxilliary basis + L_aux = dot(Lov_aenv, crpa.XpY_ss[0]) + dot(Lov_benv, crpa.XpY_ss[1]) + del Lov_aenv, Lov_benv + # This is the static approximation for the screened coulomb interaction in the auxiliary basis. + chi_aux = einsum("nx,x,mx->nm", L_aux, crpa.freqs_ss ** (-1), L_aux) + # 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. + Lpqa_loc = ao2mo(new_mf, mo_coeff=f.cluster.c_active[0]) + Lpqb_loc = ao2mo(new_mf, mo_coeff=f.cluster.c_active[1]) + + deltaW = ( + einsum("npq,nm,mrs->pqrs", Lpqa_loc, chi_aux, Lpqa_loc), + einsum("npq,nm,mrs->pqrs", Lpqa_loc, chi_aux, Lpqb_loc), + einsum("npq,nm,mrs->pqrs", Lpqb_loc, chi_aux, Lpqb_loc), + ) + wc += [deltaW] + return wc + + +def get_frag_deltaW(mf, fragment, log=None): + new_mf, crpa = get_crpa_chi(mf, fragment) + + # Apply static approximation to interaction. + static_factor = crpa.freqs_ss ** (-1) + # Now need to calculate interactions. + + nmo = new_mf.mo_coeff.shape[1] + nocc = sum(new_mf.mo_occ.T > 0) + + Lov_aenv = ao2mo(new_mf, mo_coeff=new_mf.mo_coeff[0], ijslice=(0, nocc[0], nocc[0], nmo[0])).reshape((-1, crpa.ova)) + Lov_benv = ao2mo(new_mf, mo_coeff=new_mf.mo_coeff[1], ijslice=(0, nocc[1], nocc[1], nmo[1])).reshape((-1, crpa.ovb)) + + # This is the coefficient of the cRPA reducible dd response in the auxilliary basis + L_aux = dot(Lov_aenv, crpa.XpY_ss[0]) + dot(Lov_benv, crpa.XpY_ss[1]) + del Lov_aenv, Lov_benv + # This is the static approximation for the screened coulomb interaction in the auxiliary basis. + chi_aux = einsum("nx,x,mx->nm", L_aux, crpa.freqs_ss ** (-1), L_aux) + # 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. + Lpqa_loc = ao2mo(new_mf, mo_coeff=fragment.cluster.c_active[0]) + Lpqb_loc = ao2mo(new_mf, mo_coeff=fragment.cluster.c_active[1]) + + deltaW = ( + einsum("npq,nm,mrs->pqrs", Lpqa_loc, chi_aux, Lpqa_loc), + einsum("npq,nm,mrs->pqrs", Lpqa_loc, chi_aux, Lpqb_loc), + einsum("npq,nm,mrs->pqrs", Lpqb_loc, chi_aux, Lpqb_loc), + ) + return deltaW + + +def get_crpa_chi(orig_mf, f): + def construct_crpa_mf(orig_mf, f): + """Construct mean-field object upon which an RPA calculation returns the cRPA response for a cluster.""" + s = orig_mf.get_ovlp() + + def get_canon_env(c, orig_e): + # First define rotation of initial orbitals. + r = dot(orig_mf.mo_coeff.T, s, c) + # Then get environmental orbital energies + eenv = einsum("n,np,nq->pq", orig_e, r, r) + # Diagonalise to get canonicalised equivalent. + mo_energy, c_eig = np.linalg.eigh(eenv) + return mo_energy, dot(c, c_eig) + + # Do seperately to ensure no changing of mean-field. + if np.ndim(orig_mf.mo_coeff[0]) == 1: + eo, co = get_canon_env(f.cluster.c_frozen_occ, f.mo_energy) + ev, cv = get_canon_env(f.cluster.c_frozen_vir, f.mo_energy) + e = np.concatenate((eo, ev)) + c = np.concatenate((co, cv), axis=1) + occ = np.zeros_like(e) + occ[:f.cluster.nocc_frozen] = 2.0 + else: + eoa, coa = get_canon_env(f.cluster.c_frozen_occ[0], f.mo_energy[0]) + eob, cob = get_canon_env(f.cluster.c_frozen_occ[1], f.mo_energy[1]) + eva, cva = get_canon_env(f.cluster.c_frozen_vir[0], f.mo_energy[0]) + evb, cvb = get_canon_env(f.cluster.c_frozen_vir[1], f.mo_energy[1]) + + ea = np.concatenate([eoa, eva]) + eb = np.concatenate([eob, evb]) + e = np.array((ea, eb)) + + ca = np.concatenate([coa, cva], axis=1) + cb = np.concatenate([cob, cvb], axis=1) + c = np.array((ca, cb)) + occ = np.zeros_like(e) + occ[0, :f.cluster.nocc_frozen[0]] = 1.0 + occ[1, :f.cluster.nocc_frozen[1]] = 1.0 + + new_mf = copy.copy(orig_mf) + new_mf.mo_coeff = c + new_mf.mo_energy = e + new_mf.mo_occ = occ + return new_mf + + new_mf = construct_crpa_mf(orig_mf, f) + # RPA calculation and new_mf will contain all required information for the response. + crpa = ssRPA(new_mf) + crpa.kernel() + + return new_mf, crpa + + +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 99% rename from vayesta/core/qemb/scrcoulomb.py rename to vayesta/core/screening/screening_moment.py index dacbefbc3..6973ba430 100644 --- a/vayesta/core/qemb/scrcoulomb.py +++ b/vayesta/core/screening/screening_moment.py @@ -21,7 +21,7 @@ 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. + calc_delta_e : bool, optional. Whether to calculate a nonlocal energy correction at the level of RPA npoints : int, optional Number of points for numerical integration. Default: 48. @@ -117,7 +117,7 @@ def get_eps_singlespin(no_, nv_, mo_energy): 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) seris_ov.append(kc) return seris_ov, erpa diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 2b7832230..265af5e36 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -6,7 +6,7 @@ import pyscf.scf import scipy.linalg -from vayesta.core.qemb import scrcoulomb +from vayesta.core.screening import screening_moment from vayesta.core.types import Orbitals from vayesta.core.util import dot, einsum, OptionsBase, break_into_lines, log_time from vayesta.rpa import ssRPA @@ -366,7 +366,7 @@ def add_screening(self, seris_intermed=None): # Use bare coulomb interaction from hamiltonian; this could well be cached in future. bare_eris = self.get_eris_bare() - self._seris = scrcoulomb.get_screened_eris_full(bare_eris, seris_intermed) + self._seris = screening_moment.get_screened_eris_full(bare_eris, seris_intermed) elif self.opts.screening == "crpa": raise NotImplementedError() From 8abee147754481cad386726222bc00f9ff4f6af0 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Thu, 22 Sep 2022 17:58:03 +0100 Subject: [PATCH 02/33] Working cRPA calculations with arbitrary excitation space. --- vayesta/core/screening/screening_crpa.py | 153 +++++++---------------- vayesta/rpa/ssrpa.py | 18 +-- vayesta/rpa/ssurpa.py | 2 +- vayesta/solver/hamiltonian.py | 21 +++- 4 files changed, 77 insertions(+), 117 deletions(-) diff --git a/vayesta/core/screening/screening_crpa.py b/vayesta/core/screening/screening_crpa.py index 5d8114643..403e831ff 100644 --- a/vayesta/core/screening/screening_crpa.py +++ b/vayesta/core/screening/screening_crpa.py @@ -1,3 +1,5 @@ +import scipy.linalg + from vayesta.rpa import ssRPA from .screening_moment import _get_target_rot import copy @@ -9,83 +11,48 @@ from pyscf import __config__ -def get_crpa_deltaW(emb, fragments=None, calc_delta_e=True, log=None): +def get_crpa_deltaW(emb, fragments=None, log=None): + + if fragments is None: + fragments = emb.get_fragments(sym_parent=None) + 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.") + + +def get_frag_deltaW(mf, fragment, log=None): """Generates change in 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 ---------- - emb : Embedding - Embedding instance. - fragments : list of vayesta.qemb.Fragment subclasses, optional - List of fragments for the calculation, used to define local interaction spaces. - If None, `emb.get_fragments(sym_parent=None)` is used. Default: None. - calc_delta_e : bool, optional. - Whether to calculate a nonlocal energy correction at the level of RPA + 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 ------- - seris_ov : list of tuples of np.array - List of spin-dependent screened (ov|ov), for each fragment provided. - delta_e: float - Delta RPA correction computed as difference between full system RPA energy and - cluster correlation energies; currently only functional in CAS fragmentations. + deltaW : np.array + Change in cluster local coulomb interaction due to cRPA screening. """ - if fragments is None: - fragments = emb.get_fragments(sym_parent=None) - 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.") - - wc = [] - - for f in fragments: - # Compute cRPA response. - new_mf, crpa = get_crpa_chi(emb.mf, f) - - # Apply static approximation to interaction. - static_factor = crpa.freqs_ss ** (-1) - # Now need to calculate interactions. - nmo = new_mf.mo_coeff.shape[1] - nocc = sum(new_mf.mo_occ.T > 0) + crpa = get_crpa(mf, fragment) - Lov_aenv = ao2mo(new_mf, mo_coeff=new_mf.mo_coeff[0], ijslice=(0, nocc[0], nocc[0], nmo[0])).reshape( - (-1, crpa.ova)) - Lov_benv = ao2mo(new_mf, mo_coeff=new_mf.mo_coeff[1], ijslice=(0, nocc[1], nocc[1], nmo[1])).reshape( - (-1, crpa.ovb)) + # Apply static approximation to interaction. + static_factor = crpa.freqs_ss ** (-1) + # Now need to calculate interactions. - # This is the coefficient of the cRPA reducible dd response in the auxilliary basis - L_aux = dot(Lov_aenv, crpa.XpY_ss[0]) + dot(Lov_benv, crpa.XpY_ss[1]) - del Lov_aenv, Lov_benv - # This is the static approximation for the screened coulomb interaction in the auxiliary basis. - chi_aux = einsum("nx,x,mx->nm", L_aux, crpa.freqs_ss ** (-1), L_aux) - # 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. - Lpqa_loc = ao2mo(new_mf, mo_coeff=f.cluster.c_active[0]) - Lpqb_loc = ao2mo(new_mf, mo_coeff=f.cluster.c_active[1]) + nmo = mf.mo_coeff.shape[1] + nocc = sum(mf.mo_occ.T > 0) - deltaW = ( - einsum("npq,nm,mrs->pqrs", Lpqa_loc, chi_aux, Lpqa_loc), - einsum("npq,nm,mrs->pqrs", Lpqa_loc, chi_aux, Lpqb_loc), - einsum("npq,nm,mrs->pqrs", Lpqb_loc, chi_aux, Lpqb_loc), - ) - wc += [deltaW] - return wc -def get_frag_deltaW(mf, fragment, log=None): - new_mf, crpa = get_crpa_chi(mf, fragment) - # Apply static approximation to interaction. - static_factor = crpa.freqs_ss ** (-1) - # Now need to calculate interactions. - nmo = new_mf.mo_coeff.shape[1] - nocc = sum(new_mf.mo_occ.T > 0) Lov_aenv = ao2mo(new_mf, mo_coeff=new_mf.mo_coeff[0], ijslice=(0, nocc[0], nocc[0], nmo[0])).reshape((-1, crpa.ova)) Lov_benv = ao2mo(new_mf, mo_coeff=new_mf.mo_coeff[1], ijslice=(0, nocc[1], nocc[1], nmo[1])).reshape((-1, crpa.ovb)) @@ -108,57 +75,31 @@ def get_frag_deltaW(mf, fragment, log=None): return deltaW -def get_crpa_chi(orig_mf, f): - def construct_crpa_mf(orig_mf, f): - """Construct mean-field object upon which an RPA calculation returns the cRPA response for a cluster.""" - s = orig_mf.get_ovlp() - - def get_canon_env(c, orig_e): - # First define rotation of initial orbitals. - r = dot(orig_mf.mo_coeff.T, s, c) - # Then get environmental orbital energies - eenv = einsum("n,np,nq->pq", orig_e, r, r) - # Diagonalise to get canonicalised equivalent. - mo_energy, c_eig = np.linalg.eigh(eenv) - return mo_energy, dot(c, c_eig) - - # Do seperately to ensure no changing of mean-field. - if np.ndim(orig_mf.mo_coeff[0]) == 1: - eo, co = get_canon_env(f.cluster.c_frozen_occ, f.mo_energy) - ev, cv = get_canon_env(f.cluster.c_frozen_vir, f.mo_energy) - e = np.concatenate((eo, ev)) - c = np.concatenate((co, cv), axis=1) - occ = np.zeros_like(e) - occ[:f.cluster.nocc_frozen] = 2.0 - else: - eoa, coa = get_canon_env(f.cluster.c_frozen_occ[0], f.mo_energy[0]) - eob, cob = get_canon_env(f.cluster.c_frozen_occ[1], f.mo_energy[1]) - eva, cva = get_canon_env(f.cluster.c_frozen_vir[0], f.mo_energy[0]) - evb, cvb = get_canon_env(f.cluster.c_frozen_vir[1], f.mo_energy[1]) - - ea = np.concatenate([eoa, eva]) - eb = np.concatenate([eob, evb]) - e = np.array((ea, eb)) - - ca = np.concatenate([coa, cva], axis=1) - cb = np.concatenate([cob, cvb], axis=1) - c = np.array((ca, cb)) - occ = np.zeros_like(e) - occ[0, :f.cluster.nocc_frozen[0]] = 1.0 - occ[1, :f.cluster.nocc_frozen[1]] = 1.0 - - new_mf = copy.copy(orig_mf) - new_mf.mo_coeff = c - new_mf.mo_energy = e - new_mf.mo_occ = occ - return new_mf - - new_mf = construct_crpa_mf(orig_mf, f) +def get_crpa(orig_mf, f): + + def construct_crpa_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 scipy.linalg.null_space(rot_ova).T, scipy.linalg.null_space(rot_ovb).T + + rot_ov = construct_crpa_rot(f) # RPA calculation and new_mf will contain all required information for the response. - crpa = ssRPA(new_mf) + crpa = ssRPA(orig_mf, ov_rot=rot_ov) crpa.kernel() - return new_mf, crpa + return crpa def ao2mo(mf, mo_coeff=None, ijslice=None): diff --git a/vayesta/rpa/ssrpa.py b/vayesta/rpa/ssrpa.py index 4dbb0d152..5db6f00c6 100644 --- a/vayesta/rpa/ssrpa.py +++ b/vayesta/rpa/ssrpa.py @@ -208,18 +208,20 @@ 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) + epsa = einsum("pn,n,qn->pq", self.ov_rot[0], eps, self.ov_rot[0]) + epsa, ca = scipy.linalg.eigh(epsa) + epsb = einsum("pn,n,qn->pq", self.ov_rot[1], eps, self.ov_rot[1]) + epsb, cb = scipy.linalg.eigh(epsb) + self.ov_rot = (dot(ca.T, self.ov_rot[0]), dot(cb.T, self.ov_rot[1])) + else: + epsa = epsb = 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 +241,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() diff --git a/vayesta/rpa/ssurpa.py b/vayesta/rpa/ssurpa.py index ec6a62383..a8221ef8a 100644 --- a/vayesta/rpa/ssurpa.py +++ b/vayesta/rpa/ssurpa.py @@ -78,7 +78,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/hamiltonian.py b/vayesta/solver/hamiltonian.py index 265af5e36..91b5aa1e3 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -6,9 +6,10 @@ import pyscf.scf import scipy.linalg -from vayesta.core.screening import screening_moment -from vayesta.core.types import Orbitals from vayesta.core.util import dot, einsum, OptionsBase, break_into_lines, log_time +from typing import Optional +from vayesta.core.types import Orbitals +from vayesta.core.screening import screening_moment, screening_crpa from vayesta.rpa import ssRPA @@ -401,6 +402,22 @@ def get_product_projector(p_o, p_v, no, nv): eloc = 0.5 * einsum("pq,qr,rp->", proj, m0, ApB) - einsum("pq,qp->", proj, ApB + AmB) return eloc + def add_screening(self, seris_intermed=None): + """Add screened interactions into the Hamiltonian.""" + if self.opts.screening == "mrpa": + assert(seris_intermed is not None) + + # Use bare coulomb interaction from hamiltonian; this could well be cached in future. + bare_eris = self.get_eris_bare() + + self._seris = screening_moment.get_screened_eris_full(bare_eris, seris_intermed) + + elif self.opts.screening == "crpa": + raise NotImplementedError() + + else: + raise ValueError("Unknown cluster screening protocol: %s" % self.opts.screening) + def assert_equal_spin_channels(self, message=""): na, nb = self.ncas if na != nb: From 0a6369194812e0936295446954f76859b9f9ebfd Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Mon, 26 Sep 2022 15:28:37 +0100 Subject: [PATCH 03/33] Update to cRPA functionality. --- vayesta/core/qemb/qemb.py | 4 +- vayesta/core/screening/screening_crpa.py | 76 +++++++++++++----------- vayesta/rpa/ssrpa.py | 4 +- vayesta/solver/hamiltonian.py | 41 +++++++++++-- 4 files changed, 80 insertions(+), 45 deletions(-) diff --git a/vayesta/core/qemb/qemb.py b/vayesta/core/qemb/qemb.py index e199042b1..b84d53068 100644 --- a/vayesta/core/qemb/qemb.py +++ b/vayesta/core/qemb/qemb.py @@ -760,8 +760,8 @@ def build_screened_eris(self, *args, **kwargs): # 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) + if scrfrags.count("mrpa") > 0: + build_screened_eris(self, *args, **kwargs) # Symmetry between fragments # -------------------------- diff --git a/vayesta/core/screening/screening_crpa.py b/vayesta/core/screening/screening_crpa.py index 403e831ff..aba80a110 100644 --- a/vayesta/core/screening/screening_crpa.py +++ b/vayesta/core/screening/screening_crpa.py @@ -11,16 +11,6 @@ from pyscf import __config__ -def get_crpa_deltaW(emb, fragments=None, log=None): - - if fragments is None: - fragments = emb.get_fragments(sym_parent=None) - 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.") - - def get_frag_deltaW(mf, fragment, log=None): """Generates change in 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. @@ -40,39 +30,55 @@ def get_frag_deltaW(mf, fragment, log=None): Change in cluster local coulomb interaction due to cRPA screening. """ - crpa = get_crpa(mf, fragment) + is_rhf = np.ndim(mf.mo_coeff[1]) == 1 + if not hasattr(mf, "with_df"): + raise NotImplementedError("Screened interactions require density-fitting.") - # Apply static approximation to interaction. - static_factor = crpa.freqs_ss ** (-1) + crpa = get_crpa(mf, fragment) # 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[0])).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[1])).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 - - - - - - Lov_aenv = ao2mo(new_mf, mo_coeff=new_mf.mo_coeff[0], ijslice=(0, nocc[0], nocc[0], nmo[0])).reshape((-1, crpa.ova)) - Lov_benv = ao2mo(new_mf, mo_coeff=new_mf.mo_coeff[1], ijslice=(0, nocc[1], nocc[1], nmo[1])).reshape((-1, crpa.ovb)) - - # This is the coefficient of the cRPA reducible dd response in the auxilliary basis - L_aux = dot(Lov_aenv, crpa.XpY_ss[0]) + dot(Lov_benv, crpa.XpY_ss[1]) - del Lov_aenv, Lov_benv - # This is the static approximation for the screened coulomb interaction in the auxiliary basis. - chi_aux = einsum("nx,x,mx->nm", L_aux, crpa.freqs_ss ** (-1), L_aux) # 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. - Lpqa_loc = ao2mo(new_mf, mo_coeff=fragment.cluster.c_active[0]) - Lpqb_loc = ao2mo(new_mf, mo_coeff=fragment.cluster.c_active[1]) - - deltaW = ( - einsum("npq,nm,mrs->pqrs", Lpqa_loc, chi_aux, Lpqa_loc), - einsum("npq,nm,mrs->pqrs", Lpqa_loc, chi_aux, Lpqb_loc), - einsum("npq,nm,mrs->pqrs", Lpqb_loc, chi_aux, Lpqb_loc), + c_act = fragment.cluster.c_active + if is_rhf: + c_act = (c_act, c_act) + + lpqa_loc = ao2mo(mf, mo_coeff=c_act[0]) + lpqb_loc = ao2mo(mf, mo_coeff=c_act[1]) + + # 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. + l_a = einsum("npq,nm->mpq", lpqa_loc, l_aux) + del lov_a + l_b = einsum("npq,nm->mpq", lpqb_loc, l_aux) + del lov_b + + static_fac = crpa.freqs_ss ** (-1) + + delta_w = ( + einsum("npq,n,nrs->pqrs", l_a, static_fac, l_a), + einsum("npq,n,nrs->pqrs", l_a, static_fac, l_b), + einsum("npq,n,nrs->pqrs", l_b, static_fac, l_b), ) - return deltaW + return delta_w def get_crpa(orig_mf, f): diff --git a/vayesta/rpa/ssrpa.py b/vayesta/rpa/ssrpa.py index 5db6f00c6..ea6b1aa9b 100644 --- a/vayesta/rpa/ssrpa.py +++ b/vayesta/rpa/ssrpa.py @@ -94,8 +94,8 @@ 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 :]) + self.XpY_ss = (XpY[:self.ov_rot[0].shape[0]], XpY[self.ov_rot[0].shape[0]:]) + self.XmY_ss = (XmY[:self.ov_rot[0].shape[0]], XmY[self.ov_rot[0].shape[0]:]) self.freqs_sf = None self.log.timing("Time to solve RPA problem: %s", time_string(timer() - t0)) diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 91b5aa1e3..ffa3f72f9 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -404,17 +404,42 @@ def get_product_projector(p_o, p_v, no, nv): def add_screening(self, seris_intermed=None): """Add screened interactions into the Hamiltonian.""" + seris = self._add_screening(seris_intermed) + # Need to spin-integrate resultant interaction in restricted formalism. + seris_spat = (seris[0] + seris[2] + seris[1] + seris[1].transpose((2, 3, 0, 1))) / 4.0 + + dev = [abs(x - seris_spat).max() for x in seris] + [abs(seris[2].transpose(2, 3, 0, 1) - seris_spat).max()] + self.log.info("Largest change in screened interactions due to spin integration: %e", max(dev)) + self._seris = seris_spat + + def _add_screening(self, seris_intermed=None, spin_integrate=True): + + def spin_integrate_and_report(m): + spat = (m[0] + m[1] + m[2] + seris[1].transpose((2, 3, 0, 1))) / 4.0 + + 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)) + return spat + + 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; this could well be cached in future. bare_eris = self.get_eris_bare() - - self._seris = screening_moment.get_screened_eris_full(bare_eris, seris_intermed) - + seris = screening_moment.get_screened_eris_full(bare_eris, seris_intermed) + if spin_integrate: + seris = spin_integrate_and_report(seris) + return seris elif self.opts.screening == "crpa": - raise NotImplementedError() - + bare_eris = self.get_eris_bare() + delta = screening_crpa.get_frag_deltaW(self.mf, self._fragment, self.log) + 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)]) + return seris else: raise ValueError("Unknown cluster screening protocol: %s" % self.opts.screening) @@ -630,6 +655,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) + class EB_RClusterHamiltonian(RClusterHamiltonian): @dataclasses.dataclass From 3091808d4d773d95ec9dfd560de5110e4b68c091 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Tue, 27 Sep 2022 17:59:11 +0100 Subject: [PATCH 04/33] Functioning cRPA screening protocol with and without resticted interactions. --- vayesta/core/qemb/qemb.py | 6 +++++- vayesta/core/screening/screening_crpa.py | 17 +++++++++-------- vayesta/core/screening/screening_moment.py | 2 ++ vayesta/dmet/dmet.py | 9 ++++----- vayesta/dmet/fragment.py | 1 + vayesta/ewf/ewf.py | 7 +------ vayesta/solver/hamiltonian.py | 14 ++++---------- 7 files changed, 26 insertions(+), 30 deletions(-) diff --git a/vayesta/core/qemb/qemb.py b/vayesta/core/qemb/qemb.py index b84d53068..8580bf0ea 100644 --- a/vayesta/core/qemb/qemb.py +++ b/vayesta/core/qemb/qemb.py @@ -761,7 +761,11 @@ def build_screened_eris(self, *args, **kwargs): 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) + 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, **kwargs) # Symmetry between fragments # -------------------------- diff --git a/vayesta/core/screening/screening_crpa.py b/vayesta/core/screening/screening_crpa.py index aba80a110..aa65b6263 100644 --- a/vayesta/core/screening/screening_crpa.py +++ b/vayesta/core/screening/screening_crpa.py @@ -30,6 +30,8 @@ def get_frag_deltaW(mf, fragment, log=None): Change in cluster local coulomb interaction due to cRPA screening. """ + log.info("Generating screened interaction via static limit of cRPA.") + is_rhf = np.ndim(mf.mo_coeff[1]) == 1 if not hasattr(mf, "with_df"): raise NotImplementedError("Screened interactions require density-fitting.") @@ -42,13 +44,13 @@ def get_frag_deltaW(mf, fragment, log=None): 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[0])).reshape((-1, crpa.ova)) + 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[1])).reshape((-1, crpa.ovb)) + 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 @@ -61,15 +63,14 @@ def get_frag_deltaW(mf, fragment, log=None): if is_rhf: c_act = (c_act, c_act) - lpqa_loc = ao2mo(mf, mo_coeff=c_act[0]) - lpqb_loc = ao2mo(mf, mo_coeff=c_act[1]) - # 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 lov_a + del lpqa_loc + lpqb_loc = ao2mo(mf, mo_coeff=c_act[1]) l_b = einsum("npq,nm->mpq", lpqb_loc, l_aux) - del lov_b + del lpqb_loc static_fac = crpa.freqs_ss ** (-1) diff --git a/vayesta/core/screening/screening_moment.py b/vayesta/core/screening/screening_moment.py index 6973ba430..c64992db2 100644 --- a/vayesta/core/screening/screening_moment.py +++ b/vayesta/core/screening/screening_moment.py @@ -125,6 +125,8 @@ def get_eps_singlespin(no_, nv_, mo_energy): 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] 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/ewf/ewf.py b/vayesta/ewf/ewf.py index cf81c3c05..3802e90f3 100644 --- a/vayesta/ewf/ewf.py +++ b/vayesta/ewf/ewf.py @@ -139,12 +139,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("") diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index ffa3f72f9..0856d0840 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -404,18 +404,12 @@ def get_product_projector(p_o, p_v, no, nv): def add_screening(self, seris_intermed=None): """Add screened interactions into the Hamiltonian.""" - seris = self._add_screening(seris_intermed) - # Need to spin-integrate resultant interaction in restricted formalism. - seris_spat = (seris[0] + seris[2] + seris[1] + seris[1].transpose((2, 3, 0, 1))) / 4.0 - - dev = [abs(x - seris_spat).max() for x in seris] + [abs(seris[2].transpose(2, 3, 0, 1) - seris_spat).max()] - self.log.info("Largest change in screened interactions due to spin integration: %e", max(dev)) - self._seris = seris_spat + self._seris = self._add_screening(seris_intermed, spin_integrate=True) def _add_screening(self, seris_intermed=None, spin_integrate=True): def spin_integrate_and_report(m): - spat = (m[0] + m[1] + m[2] + seris[1].transpose((2, 3, 0, 1))) / 4.0 + spat = (m[0] + m[1] + m[2] + m[1].transpose((2, 3, 0, 1))) / 4.0 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)) @@ -427,7 +421,7 @@ def spin_integrate_and_report(m): assert(seris_intermed is not None) # Use bare coulomb interaction from hamiltonian; this could well be cached in future. bare_eris = self.get_eris_bare() - seris = screening_moment.get_screened_eris_full(bare_eris, seris_intermed) + seris = screening_moment.get_screened_eris_full(bare_eris, seris_intermed[0], log=self.log) if spin_integrate: seris = spin_integrate_and_report(seris) return seris @@ -657,7 +651,7 @@ def __contains__(self, item): def add_screening(self, seris_intermed=None): """Add screened interactions into the Hamiltonian.""" - self._seris = self._add_screening(seris_intermed) + self._seris = self._add_screening(seris_intermed, spin_integrate=False) class EB_RClusterHamiltonian(RClusterHamiltonian): From 04e3a2e8b4e5b3491a3c8b5b82cd0d10020c4550 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Mon, 3 Oct 2022 15:48:37 +0100 Subject: [PATCH 05/33] Added functionality to use dynamical screened coulomb interaction resulting from cRPA in clusters. --- vayesta/core/qemb/fragment.py | 5 +++ vayesta/core/screening/screening_crpa.py | 53 +++++++++++++++++++----- vayesta/solver/hamiltonian.py | 21 ++++++---- 3 files changed, 62 insertions(+), 17 deletions(-) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index a6a152de5..6d855caaa 100644 --- a/vayesta/core/qemb/fragment.py +++ b/vayesta/core/qemb/fragment.py @@ -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 @@ -1026,9 +1027,11 @@ def check_solver(self, solver): check_solver_config(is_uhf, is_eb, solver, self.log) def get_solver(self, solver=None): + if solver is None: solver = self.solver cl_ham = self.hamil + solver_cls = get_solver_class(cl_ham, solver) solver_opts = self.get_solver_options(solver) cluster_solver = solver_cls(cl_ham, **solver_opts) @@ -1038,6 +1041,8 @@ def get_solver(self, solver=None): def get_frag_hamil(self): # This detects based on fragment what kind of Hamiltonian is appropriate (restricted and/or EB). + if "crpa_full" in self.opts.screening: + self.bos_freqs, self.couplings = get_frag_W(self.mf, self, self.log) return ClusterHamiltonian(self, self.mf, self.log, screening=self.opts.screening, cache_eris=self.opts.store_eris) diff --git a/vayesta/core/screening/screening_crpa.py b/vayesta/core/screening/screening_crpa.py index aa65b6263..8d43cc00c 100644 --- a/vayesta/core/screening/screening_crpa.py +++ b/vayesta/core/screening/screening_crpa.py @@ -10,9 +10,38 @@ from pyscf.ao2mo import _ao2mo from pyscf import __config__ +def get_frag_W(mf, fragment, 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.") + + l_a, l_b, crpa = set_up_W_crpa(mf, fragment, log) + freqs = crpa.freqs_ss + log.info("cRPA resulted in %d poles", len(freqs)) + couplings = (l_a, l_b) + return freqs, couplings + def get_frag_deltaW(mf, fragment, log=None): - """Generates change in coulomb interaction due to screening at the level of cRPA. + """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 @@ -32,6 +61,18 @@ def get_frag_deltaW(mf, fragment, log=None): log.info("Generating screened interaction via static limit of cRPA.") + l_a, l_b, crpa = set_up_W_crpa(mf, fragment, log) + + static_fac = crpa.freqs_ss ** (-1) + + delta_w = ( + einsum("npq,n,nrs->pqrs", l_a, static_fac, l_a), + einsum("npq,n,nrs->pqrs", l_a, static_fac, l_b), + einsum("npq,n,nrs->pqrs", l_b, static_fac, l_b), + ) + return delta_w, crpa + +def set_up_W_crpa(mf, fragment, log=None): is_rhf = np.ndim(mf.mo_coeff[1]) == 1 if not hasattr(mf, "with_df"): raise NotImplementedError("Screened interactions require density-fitting.") @@ -71,15 +112,7 @@ def get_frag_deltaW(mf, fragment, log=None): lpqb_loc = ao2mo(mf, mo_coeff=c_act[1]) l_b = einsum("npq,nm->mpq", lpqb_loc, l_aux) del lpqb_loc - - static_fac = crpa.freqs_ss ** (-1) - - delta_w = ( - einsum("npq,n,nrs->pqrs", l_a, static_fac, l_a), - einsum("npq,n,nrs->pqrs", l_a, static_fac, l_b), - einsum("npq,n,nrs->pqrs", l_b, static_fac, l_b), - ) - return delta_w + return l_a, l_b, crpa def get_crpa(orig_mf, f): diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 0856d0840..408f4c040 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -425,15 +425,20 @@ def spin_integrate_and_report(m): if spin_integrate: seris = spin_integrate_and_report(seris) return seris - elif self.opts.screening == "crpa": + elif self.opts.screening[:4] == "crpa": bare_eris = self.get_eris_bare() - delta = screening_crpa.get_frag_deltaW(self.mf, self._fragment, self.log) - if spin_integrate: - delta = spin_integrate_and_report(delta) - seris = bare_eris + delta + delta, crpa = screening_crpa.get_frag_deltaW(self.mf, self._fragment, self.log) + if "store" in self.opts.screening: + self.crpa = crpa + if "full" in self.opts.screening: + pass else: - seris = tuple([x + y for x, y in zip(bare_eris, delta)]) - return seris + 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)]) + return seris else: raise ValueError("Unknown cluster screening protocol: %s" % self.opts.screening) @@ -716,6 +721,8 @@ def get_polaritonic_shifted_couplings(self): 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 From 80eac9d81199e9a759b4dadd2353d1124b994da6 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Tue, 4 Oct 2022 14:07:47 +0100 Subject: [PATCH 06/33] Slight refactor of static cRPA calculation. --- vayesta/core/screening/screening_crpa.py | 5 +++-- vayesta/solver/hamiltonian.py | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/vayesta/core/screening/screening_crpa.py b/vayesta/core/screening/screening_crpa.py index 8d43cc00c..fce63dac8 100644 --- a/vayesta/core/screening/screening_crpa.py +++ b/vayesta/core/screening/screening_crpa.py @@ -62,8 +62,9 @@ def get_frag_deltaW(mf, fragment, log=None): log.info("Generating screened interaction via static limit of cRPA.") l_a, l_b, crpa = set_up_W_crpa(mf, fragment, log) - - static_fac = crpa.freqs_ss ** (-1) + # 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 = - 2.0 * (crpa.freqs_ss ** (-1)) delta_w = ( einsum("npq,n,nrs->pqrs", l_a, static_fac, l_a), diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 408f4c040..7b250de7b 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -437,7 +437,7 @@ def spin_integrate_and_report(m): delta = spin_integrate_and_report(delta) seris = bare_eris + delta else: - seris = tuple([x - y for x, y in zip(bare_eris, delta)]) + seris = tuple([x + y for x, y in zip(bare_eris, delta)]) return seris else: raise ValueError("Unknown cluster screening protocol: %s" % self.opts.screening) From 0db7205ac0d61da4362635bd28d0c7a5243be683 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Tue, 4 Oct 2022 18:08:12 +0100 Subject: [PATCH 07/33] Fixed bug in ssRPA, and added option to use N^6 RPA in mRPA screening. --- vayesta/core/screening/screening_moment.py | 74 ++++++++++++++-------- vayesta/rpa/ssrpa.py | 9 ++- 2 files changed, 56 insertions(+), 27 deletions(-) diff --git a/vayesta/core/screening/screening_moment.py b/vayesta/core/screening/screening_moment.py index c64992db2..0b5bc4aea 100644 --- a/vayesta/core/screening/screening_moment.py +++ b/vayesta/core/screening/screening_moment.py @@ -2,7 +2,7 @@ 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 @@ -53,29 +53,9 @@ def build_screened_eris(emb, fragments=None, cderi_ov=None, calc_e=True, npoints 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, erpa = calc_moms_RIRPA(emb.mf, target_rots, ovs_active, log, cderi_ov, calc_e, 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)) @@ -97,7 +77,12 @@ 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) + mominv2 = np.linalg.inv(mom) + e, c = np.linalg.eigh(mom) + log.info("Minimal eigenvalue of local zeroth moment: %e", min(e)) + + mominv = einsum("pn,n,qn->pq", c, e**(-1), c) + log.info("Deviation in different inversion approaches: %e", abs(mominv - mominv2)) apb = dot(mominv, amb, mominv) # This is the renormalised coulomb kernel in the cluster. @@ -122,6 +107,45 @@ def get_eps_singlespin(no_, nv_, mo_energy): return seris_ov, erpa +def calc_moms_RIRPA(mf, target_rots, ovs_active, log, cderi_ov, calc_e, npoints): + rpa = ssRIRPA(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. + 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, erpa + +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.""" diff --git a/vayesta/rpa/ssrpa.py b/vayesta/rpa/ssrpa.py index ea6b1aa9b..076cb2b23 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.ov_rot[0].shape[0]], XpY[self.ov_rot[0].shape[0]:]) - self.XmY_ss = (XmY[:self.ov_rot[0].shape[0]], XmY[self.ov_rot[0].shape[0]:]) + + 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)) From 29387e121b1a4d954982e371446d17f3c04558dd Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Thu, 6 Oct 2022 23:30:11 +0100 Subject: [PATCH 08/33] Added standardised reporting of spin symmetry breaking and dependence for screened interactions, in orbitally invariant manner. --- vayesta/core/screening/screening_moment.py | 13 ++------ vayesta/solver/hamiltonian.py | 38 ++++++++++++++++++++-- 2 files changed, 37 insertions(+), 14 deletions(-) diff --git a/vayesta/core/screening/screening_moment.py b/vayesta/core/screening/screening_moment.py index 0b5bc4aea..290fc50e6 100644 --- a/vayesta/core/screening/screening_moment.py +++ b/vayesta/core/screening/screening_moment.py @@ -77,12 +77,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). - mominv2 = np.linalg.inv(mom) e, c = np.linalg.eigh(mom) - log.info("Minimal eigenvalue of local zeroth moment: %e", min(e)) + 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) - log.info("Deviation in different inversion approaches: %e", abs(mominv - mominv2)) apb = dot(mominv, amb, mominv) # This is the renormalised coulomb kernel in the cluster. @@ -97,10 +96,6 @@ def get_eps_singlespin(no_, nv_, mo_energy): 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, mom, amb) seris_ov.append(kc) @@ -160,10 +155,6 @@ 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 seris = (replace_ov(eris[0], seris_ov[0], 'aa'), diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 7b250de7b..f4609c979 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -414,7 +414,7 @@ def spin_integrate_and_report(m): 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)) return spat - + seris = None if self.opts.screening is None: raise ValueError("Attempted to add screening to fragment with no screening protocol specified.") if self.opts.screening == "mrpa": @@ -424,7 +424,6 @@ def spin_integrate_and_report(m): seris = screening_moment.get_screened_eris_full(bare_eris, seris_intermed[0], log=self.log) if spin_integrate: seris = spin_integrate_and_report(seris) - return seris elif self.opts.screening[:4] == "crpa": bare_eris = self.get_eris_bare() delta, crpa = screening_crpa.get_frag_deltaW(self.mf, self._fragment, self.log) @@ -438,10 +437,43 @@ def spin_integrate_and_report(m): seris = bare_eris + delta else: seris = tuple([x + y for x, y in zip(bare_eris, delta)]) - return seris else: raise ValueError("Unknown cluster screening protocol: %s" % self.opts.screening) + 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 if na != nb: From 694ed4a1bfad025189288d8b379dd31a066749c5 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Fri, 7 Oct 2022 16:05:05 +0100 Subject: [PATCH 09/33] Ensured crpa_full screening returns a fermionic wavefunction for ewf. --- vayesta/core/qemb/qemb.py | 2 +- vayesta/ewf/ewf.py | 2 ++ vayesta/solver/ebcc.py | 3 +++ 3 files changed, 6 insertions(+), 1 deletion(-) diff --git a/vayesta/core/qemb/qemb.py b/vayesta/core/qemb/qemb.py index 8580bf0ea..4424c0713 100644 --- a/vayesta/core/qemb/qemb.py +++ b/vayesta/core/qemb/qemb.py @@ -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 diff --git a/vayesta/ewf/ewf.py b/vayesta/ewf/ewf.py index 3802e90f3..11fee6151 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! 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): From 2dbdfd60361955116652810af38647f5f8b30631 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Sat, 8 Oct 2022 13:20:03 +0100 Subject: [PATCH 10/33] Added ability to use spin-integrated mRPA screened interactions. --- vayesta/core/screening/screening_moment.py | 25 +++++++++++++++++----- vayesta/solver/hamiltonian.py | 7 +++++- 2 files changed, 26 insertions(+), 6 deletions(-) diff --git a/vayesta/core/screening/screening_moment.py b/vayesta/core/screening/screening_moment.py index 290fc50e6..d5063b311 100644 --- a/vayesta/core/screening/screening_moment.py +++ b/vayesta/core/screening/screening_moment.py @@ -45,8 +45,6 @@ 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] @@ -59,17 +57,23 @@ def build_screened_eris(emb, fragments=None, cderi_ov=None, calc_e=True, 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 = [] @@ -92,6 +96,10 @@ 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])) @@ -157,6 +165,9 @@ def replace_ov(full, ov, spins): out[v1,o1,v2,o2] = ov.transpose([1, 0, 3, 2]) 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')) @@ -243,6 +254,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/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index f4609c979..abfcec328 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -408,11 +408,16 @@ def add_screening(self, seris_intermed=None): def _add_screening(self, seris_intermed=None, spin_integrate=True): - def spin_integrate_and_report(m): + 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 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)) + else: + self.log.info("Largest change in screened interactions due to spin integration: %e", max(dev)) + return spat seris = None if self.opts.screening is None: From 1340cb9a00745d720822c9439f8e23185b628f0c Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Sat, 8 Oct 2022 17:50:00 +0100 Subject: [PATCH 11/33] Fix bug in unscreened interactions. --- vayesta/core/qemb/fragment.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index 6d855caaa..0c3a1aea7 100644 --- a/vayesta/core/qemb/fragment.py +++ b/vayesta/core/qemb/fragment.py @@ -1027,11 +1027,9 @@ def check_solver(self, solver): check_solver_config(is_uhf, is_eb, solver, self.log) def get_solver(self, solver=None): - if solver is None: solver = self.solver cl_ham = self.hamil - solver_cls = get_solver_class(cl_ham, solver) solver_opts = self.get_solver_options(solver) cluster_solver = solver_cls(cl_ham, **solver_opts) @@ -1040,9 +1038,10 @@ def get_solver(self, solver=None): return cluster_solver 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, self.log) # This detects based on fragment what kind of Hamiltonian is appropriate (restricted and/or EB). - if "crpa_full" in self.opts.screening: - self.bos_freqs, self.couplings = get_frag_W(self.mf, self, self.log) return ClusterHamiltonian(self, self.mf, self.log, screening=self.opts.screening, cache_eris=self.opts.store_eris) From 921954b33c1dddcfcd30693d7d7451432d1452c7 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Sun, 9 Oct 2022 21:33:56 +0100 Subject: [PATCH 12/33] Ensured proper treatment of cRPA when no excitations are in the cRPA space. --- vayesta/core/screening/screening_crpa.py | 54 +++++++++++++++--------- vayesta/rpa/ssrpa.py | 13 +++--- 2 files changed, 43 insertions(+), 24 deletions(-) diff --git a/vayesta/core/screening/screening_crpa.py b/vayesta/core/screening/screening_crpa.py index fce63dac8..6499ccd9d 100644 --- a/vayesta/core/screening/screening_crpa.py +++ b/vayesta/core/screening/screening_crpa.py @@ -10,6 +10,10 @@ from pyscf.ao2mo import _ao2mo from pyscf import __config__ + +class cRPAError(RuntimeError): + pass + def get_frag_W(mf, fragment, 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. @@ -32,11 +36,17 @@ def get_frag_W(mf, fragment, log=None): """ log.info("Generating screened interaction via frequency dependent cRPA.") - - l_a, l_b, crpa = set_up_W_crpa(mf, fragment, log) - freqs = crpa.freqs_ss + try: + l_a, l_b, crpa = set_up_W_crpa(mf, fragment, 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)) - couplings = (l_a, l_b) + return freqs, couplings @@ -60,25 +70,29 @@ def get_frag_deltaW(mf, fragment, log=None): """ log.info("Generating screened interaction via static limit of cRPA.") - - l_a, l_b, crpa = set_up_W_crpa(mf, fragment, log) - # 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 = - 2.0 * (crpa.freqs_ss ** (-1)) - - delta_w = ( - einsum("npq,n,nrs->pqrs", l_a, static_fac, l_a), - einsum("npq,n,nrs->pqrs", l_a, static_fac, l_b), - einsum("npq,n,nrs->pqrs", l_b, static_fac, l_b), - ) + try: + l_a, l_b, crpa = set_up_W_crpa(mf, fragment, 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 = - 2.0 * (crpa.freqs_ss ** (-1)) + + delta_w = ( + einsum("npq,n,nrs->pqrs", l_a, static_fac, l_a), + einsum("npq,n,nrs->pqrs", l_a, static_fac, l_b), + einsum("npq,n,nrs->pqrs", l_b, static_fac, l_b), + ) return delta_w, crpa def set_up_W_crpa(mf, fragment, 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 = get_crpa(mf, fragment) + crpa = get_crpa(mf, fragment, log) # Now need to calculate interactions. nmo = mf.mo_coeff.shape[1] nocc = sum(mf.mo_occ.T > 0) @@ -116,7 +130,7 @@ def set_up_W_crpa(mf, fragment, log=None): return l_a, l_b, crpa -def get_crpa(orig_mf, f): +def get_crpa(orig_mf, f, log): def construct_crpa_rot(f): """Constructs the rotation of the overall mean-field space into which """ @@ -136,10 +150,12 @@ def construct_crpa_rot(f): return scipy.linalg.null_space(rot_ova).T, scipy.linalg.null_space(rot_ovb).T rot_ov = construct_crpa_rot(f) + 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 diff --git a/vayesta/rpa/ssrpa.py b/vayesta/rpa/ssrpa.py index 076cb2b23..97f38c091 100644 --- a/vayesta/rpa/ssrpa.py +++ b/vayesta/rpa/ssrpa.py @@ -213,11 +213,14 @@ def _gen_arrays(self, xc_kernel=None, alpha=1.0): eps = eps.reshape((self.ova,)) if self.ov_rot is not None: - epsa = einsum("pn,n,qn->pq", self.ov_rot[0], eps, self.ov_rot[0]) - epsa, ca = scipy.linalg.eigh(epsa) - epsb = einsum("pn,n,qn->pq", self.ov_rot[1], eps, self.ov_rot[1]) - epsb, cb = scipy.linalg.eigh(epsb) - self.ov_rot = (dot(ca.T, self.ov_rot[0]), dot(cb.T, self.ov_rot[1])) + 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 From c33fe3d9d10cae0eb95eda4b00c3eaa4ea6e3b7c Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Tue, 11 Oct 2022 18:18:15 +0100 Subject: [PATCH 13/33] Added warnings for storage of cRPA object within static cRPA calculation. --- vayesta/solver/hamiltonian.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index abfcec328..fd1a3cac6 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -334,9 +334,10 @@ def pad(a, diag_val): 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 @@ -371,7 +372,6 @@ def add_screening(self, seris_intermed=None): elif self.opts.screening == "crpa": raise NotImplementedError() - else: raise ValueError("Unknown cluster screening protocol: %s" % self.opts.screening) @@ -433,9 +433,13 @@ def spin_integrate_and_report(m, warn_threshold=1e-6): bare_eris = self.get_eris_bare() delta, crpa = screening_crpa.get_frag_deltaW(self.mf, self._fragment, 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: - pass + # 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) From ad6cd71d055886b2754bcc0df288ef3aee22c588 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Tue, 11 Oct 2022 18:27:09 +0100 Subject: [PATCH 14/33] Update RPA moment calculation to support cRPA. --- vayesta/rpa/ssrpa.py | 29 ++++++++++------------------- 1 file changed, 10 insertions(+), 19 deletions(-) diff --git a/vayesta/rpa/ssrpa.py b/vayesta/rpa/ssrpa.py index 97f38c091..829c7146d 100644 --- a/vayesta/rpa/ssrpa.py +++ b/vayesta/rpa/ssrpa.py @@ -329,30 +329,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 From 2350016c45f8901d0e0680110cc68931eb3d961f Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Wed, 12 Oct 2022 16:56:08 +0100 Subject: [PATCH 15/33] Added polarisation coupling calculation into cRPA calculation. --- vayesta/core/qemb/fragment.py | 3 +- vayesta/core/screening/screening_crpa.py | 55 ++++++++++++++++++------ vayesta/rpa/ssrpa.py | 10 ++++- vayesta/rpa/ssurpa.py | 9 ++-- vayesta/solver/hamiltonian.py | 8 ++-- 5 files changed, 62 insertions(+), 23 deletions(-) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index 0c3a1aea7..c0b0cb807 100644 --- a/vayesta/core/qemb/fragment.py +++ b/vayesta/core/qemb/fragment.py @@ -1040,7 +1040,8 @@ def get_solver(self, solver=None): 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, self.log) + self.bos_freqs, self.couplings = get_frag_W(self.mf, self, pcoupling=not ("old" 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) diff --git a/vayesta/core/screening/screening_crpa.py b/vayesta/core/screening/screening_crpa.py index 6499ccd9d..5530b8eed 100644 --- a/vayesta/core/screening/screening_crpa.py +++ b/vayesta/core/screening/screening_crpa.py @@ -14,7 +14,7 @@ class cRPAError(RuntimeError): pass -def get_frag_W(mf, fragment, log=None): +def get_frag_W(mf, fragment, pcoupling=True, 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. @@ -37,7 +37,7 @@ def get_frag_W(mf, fragment, log=None): log.info("Generating screened interaction via frequency dependent cRPA.") try: - l_a, l_b, crpa = set_up_W_crpa(mf, fragment, log) + l_a, l_b, crpa = set_up_W_crpa(mf, fragment, pcoupling, log=log) except cRPAError as e: freqs = np.zeros((0,)) nmo = mf.mo_coeff.shape[1] @@ -50,7 +50,7 @@ def get_frag_W(mf, fragment, log=None): return freqs, couplings -def get_frag_deltaW(mf, fragment, log=None): +def get_frag_deltaW(mf, fragment, pcoupling=True, 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. @@ -71,7 +71,7 @@ def get_frag_deltaW(mf, fragment, log=None): log.info("Generating screened interaction via static limit of cRPA.") try: - l_a, l_b, crpa = set_up_W_crpa(mf, fragment, log) + l_a, l_b, crpa = set_up_W_crpa(mf, fragment, pcoupling, log=log) except cRPAError as e: nmo = mf.mo_coeff.shape[1] delta_w = tuple([np.zeros([nmo] * 4)] * 4) @@ -79,20 +79,20 @@ def get_frag_deltaW(mf, fragment, log=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 = - 2.0 * (crpa.freqs_ss ** (-1)) + static_fac = - 1.0 * (crpa.freqs_ss ** (-1)) delta_w = ( - einsum("npq,n,nrs->pqrs", l_a, static_fac, l_a), - einsum("npq,n,nrs->pqrs", l_a, static_fac, l_b), - einsum("npq,n,nrs->pqrs", l_b, static_fac, l_b), + 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, log=None): +def set_up_W_crpa(mf, fragment, pcoupling=True, 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 = get_crpa(mf, fragment, log) + 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) @@ -127,12 +127,38 @@ def set_up_W_crpa(mf, fragment, log=None): 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) + return l_a, l_b, crpa def get_crpa(orig_mf, f, log): - def construct_crpa_rot(f): + 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]") @@ -147,16 +173,17 @@ def construct_crpa_rot(f): 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 scipy.linalg.null_space(rot_ova).T, scipy.linalg.null_space(rot_ovb).T + return rot_ova, rot_ovb - rot_ov = construct_crpa_rot(f) + 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 + return crpa, rot_loc, rot_ov def ao2mo(mf, mo_coeff=None, ijslice=None): diff --git a/vayesta/rpa/ssrpa.py b/vayesta/rpa/ssrpa.py index 829c7146d..0dd042d00 100644 --- a/vayesta/rpa/ssrpa.py +++ b/vayesta/rpa/ssrpa.py @@ -204,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 :] @@ -223,6 +222,13 @@ def _gen_arrays(self, xc_kernel=None, alpha=1.0): 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([epsa, epsb]) diff --git a/vayesta/rpa/ssurpa.py b/vayesta/rpa/ssurpa.py index a8221ef8a..0c605cdc5 100644 --- a/vayesta/rpa/ssurpa.py +++ b/vayesta/rpa/ssurpa.py @@ -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]) diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index fd1a3cac6..7327b4489 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -431,7 +431,8 @@ def spin_integrate_and_report(m, warn_threshold=1e-6): 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.mf, self._fragment, self.log) + delta, crpa = screening_crpa.get_frag_deltaW(self.mf, self._fragment, + pcoupling=not ("old" 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 @@ -755,8 +756,9 @@ 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): From b3d8e275e837ddc36aab92d95e31208dfa5f62dd Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Mon, 17 Oct 2022 14:01:50 +0100 Subject: [PATCH 16/33] Added option to only use cRPA screening in the ov space, and changed default to not include irreducible polarisability. --- vayesta/core/qemb/fragment.py | 3 ++- vayesta/core/screening/screening_crpa.py | 19 ++++++++++++++----- vayesta/solver/hamiltonian.py | 4 +++- 3 files changed, 19 insertions(+), 7 deletions(-) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index c0b0cb807..aa4ebbf65 100644 --- a/vayesta/core/qemb/fragment.py +++ b/vayesta/core/qemb/fragment.py @@ -1040,7 +1040,8 @@ def get_solver(self, solver=None): 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=not ("old" 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, diff --git a/vayesta/core/screening/screening_crpa.py b/vayesta/core/screening/screening_crpa.py index 5530b8eed..69f92c608 100644 --- a/vayesta/core/screening/screening_crpa.py +++ b/vayesta/core/screening/screening_crpa.py @@ -14,7 +14,7 @@ class cRPAError(RuntimeError): pass -def get_frag_W(mf, fragment, pcoupling=True, log=None): +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. @@ -37,7 +37,7 @@ def get_frag_W(mf, fragment, pcoupling=True, log=None): log.info("Generating screened interaction via frequency dependent cRPA.") try: - l_a, l_b, crpa = set_up_W_crpa(mf, fragment, pcoupling, log=log) + 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] @@ -50,7 +50,7 @@ def get_frag_W(mf, fragment, pcoupling=True, log=None): return freqs, couplings -def get_frag_deltaW(mf, fragment, pcoupling=True, log=None): +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. @@ -71,7 +71,7 @@ def get_frag_deltaW(mf, fragment, pcoupling=True, log=None): log.info("Generating screened interaction via static limit of cRPA.") try: - l_a, l_b, crpa = set_up_W_crpa(mf, fragment, pcoupling, log=log) + 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) @@ -88,7 +88,7 @@ def get_frag_deltaW(mf, fragment, pcoupling=True, log=None): ) return delta_w, crpa -def set_up_W_crpa(mf, fragment, pcoupling=True, log=None): +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.") @@ -153,6 +153,15 @@ def set_up_W_crpa(mf, fragment, pcoupling=True, log=None): 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 diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 7327b4489..7a22c6375 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -432,7 +432,9 @@ def spin_integrate_and_report(m, warn_threshold=1e-6): elif self.opts.screening[:4] == "crpa": bare_eris = self.get_eris_bare() delta, crpa = screening_crpa.get_frag_deltaW(self.mf, self._fragment, - pcoupling=not ("old" in self.opts.screening), log=self.log) + 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 From fa4fb1fd6204a38b9976af1e17864477d4728dc5 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Tue, 1 Nov 2022 22:25:38 +0000 Subject: [PATCH 17/33] Added RPA external correction for ewf. --- vayesta/core/qemb/fragment.py | 6 +- vayesta/core/qemb/qemb.py | 20 ++++-- vayesta/core/screening/screening_moment.py | 20 ++---- vayesta/ewf/ewf.py | 1 + vayesta/ewf/fragment.py | 4 +- vayesta/solver/hamiltonian.py | 78 ++++++++++------------ 6 files changed, 66 insertions(+), 63 deletions(-) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index aa4ebbf65..f8fa353ed 100644 --- a/vayesta/core/qemb/fragment.py +++ b/vayesta/core/qemb/fragment.py @@ -74,6 +74,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 @@ -1035,7 +1036,10 @@ def get_solver(self, solver=None): cluster_solver = solver_cls(cl_ham, **solver_opts) if self.opts.screening is not None: cluster_solver.hamil.add_screening(self._seris_ov) - return cluster_solver + e_loc_rpa = None + if self.base.opts.ext_rpa_correction: + e_loc_rpa = cl_ham.calc_loc_erpa(*self._seris_ov[1:]) + return cluster_solver, e_loc_rpa def get_frag_hamil(self): if self.opts.screening is not None: diff --git a/vayesta/core/qemb/qemb.py b/vayesta/core/qemb/qemb.py index 4424c0713..ee79e6b63 100644 --- a/vayesta/core/qemb/qemb.py +++ b/vayesta/core/qemb/qemb.py @@ -114,6 +114,7 @@ class Options(OptionsBase): 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 + ext_rpa_correction: bool = False class Embedding: """Base class for quantum embedding methods. @@ -518,6 +519,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 not self.opts.ext_rpa_correction: + 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,17 +762,19 @@ 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: + nmomscr = len([x.opts.screening for x in self.fragments if x.opts.screening == "mrpa"]) + if self.opts.ext_rpa_correction: + if nmomscr < self.nfrag: + raise NotImplementedError("External dRPA correction currently requires all fragments use mrpa screening.") # 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: + self.e_rpa, energy_error = rpa.kernel_energy(correction='linear') + 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, **kwargs) + build_screened_eris(self, *args, store_m0=self.opts.ext_rpa_correction, **kwargs) # Symmetry between fragments # -------------------------- @@ -1565,6 +1574,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/screening/screening_moment.py b/vayesta/core/screening/screening_moment.py index d5063b311..5ad5844f6 100644 --- a/vayesta/core/screening/screening_moment.py +++ b/vayesta/core/screening/screening_moment.py @@ -7,7 +7,7 @@ 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_delta_e : 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 @@ -51,7 +51,7 @@ def build_screened_eris(emb, fragments=None, cderi_ov=None, calc_e=True, npoints r_virs = [f.get_overlap('mo[vir]|cluster[vir]') for f in fragments] target_rots, ovs_active = _get_target_rot(r_occs, r_virs) - local_moments, erpa = calc_moms_RIRPA(emb.mf, target_rots, ovs_active, log, cderi_ov, calc_e, npoints) + 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) @@ -105,18 +105,12 @@ def get_eps_singlespin(no_, nv_, mo_energy): kcbb = kc[ova:, ova:].reshape((no[1], nv[1], no[1], nv[1])) kc = (kcaa, kcab, kcbb) - f._seris_ov = (kc, mom, amb) + 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, calc_e, npoints): +def calc_moms_RIRPA(mf, target_rots, ovs_active, log, cderi_ov, npoints): rpa = ssRIRPA(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: @@ -138,7 +132,7 @@ def calc_moms_RIRPA(mf, target_rots, ovs_active, log, cderi_ov, calc_e, npoints) local_moments += [mom] n += sum(nov) - return local_moments, erpa + return local_moments def calc_moms_RPA(mf, target_rots, ovs_active, log, cderi_ov, calc_e, npoints): rpa = ssRPA(mf, log=log) diff --git a/vayesta/ewf/ewf.py b/vayesta/ewf/ewf.py index 11fee6151..adabbd690 100644 --- a/vayesta/ewf/ewf.py +++ b/vayesta/ewf/ewf.py @@ -291,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..4126f74ba 100644 --- a/vayesta/ewf/fragment.py +++ b/vayesta/ewf/fragment.py @@ -189,7 +189,7 @@ def kernel(self, solver=None, init_guess=None): init_guess = self.get_init_guess(init_guess, solver, cluster) # Create solver object - cluster_solver = self.get_solver(solver) + cluster_solver, e_corr_rpa = self.get_solver(solver) # --- Chemical potential cpt_frag = self.base.opts.global_frag_chempot @@ -253,7 +253,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/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 7a22c6375..1fb2244a5 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -360,47 +360,38 @@ 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) - - # Use bare coulomb interaction from hamiltonian; this could well be cached in future. - bare_eris = self.get_eris_bare() - - self._seris = screening_moment.get_screened_eris_full(bare_eris, seris_intermed) - - elif self.opts.screening == "crpa": - raise NotImplementedError() - else: - raise ValueError("Unknown cluster screening protocol: %s" % self.opts.screening) - - def calc_loc_erpa(self): - - 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 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 get_product_projector(p_o, p_v, no, nv): - return einsum("ij,ab->iajb", p_o, p_v).reshape((no * nv, no * nv)) - - 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) - - proj = get_product_projector() - eloc = 0.5 * einsum("pq,qr,rp->", proj, m0, ApB) - einsum("pq,qp->", proj, ApB + AmB) - return eloc + def calc_loc_erpa(self, m0, amb): + + 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) + 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): """Add screened interactions into the Hamiltonian.""" @@ -424,7 +415,7 @@ def spin_integrate_and_report(m, warn_threshold=1e-6): 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; this could well be cached in future. + # 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: @@ -794,3 +785,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() From d0f6242a37bc60377d91fcd6d29406842cbae15f Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Thu, 3 Nov 2022 16:18:52 +0000 Subject: [PATCH 18/33] Added RPA external correction via cumulant 2rdm contribution. --- vayesta/core/qemb/fragment.py | 5 +++- vayesta/core/qemb/qemb.py | 24 ++++++++++++++---- vayesta/solver/hamiltonian.py | 47 +++++++++++++++++++++++------------ 3 files changed, 54 insertions(+), 22 deletions(-) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index f8fa353ed..b23a392df 100644 --- a/vayesta/core/qemb/fragment.py +++ b/vayesta/core/qemb/fragment.py @@ -1038,7 +1038,10 @@ def get_solver(self, solver=None): cluster_solver.hamil.add_screening(self._seris_ov) e_loc_rpa = None if self.base.opts.ext_rpa_correction: - e_loc_rpa = cl_ham.calc_loc_erpa(*self._seris_ov[1:]) + 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 = cl_ham.calc_loc_erpa(*self._seris_ov[1:], cumulant) return cluster_solver, e_loc_rpa def get_frag_hamil(self): diff --git a/vayesta/core/qemb/qemb.py b/vayesta/core/qemb/qemb.py index ee79e6b63..0aa314940 100644 --- a/vayesta/core/qemb/qemb.py +++ b/vayesta/core/qemb/qemb.py @@ -114,7 +114,7 @@ class Options(OptionsBase): 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 - ext_rpa_correction: bool = False + ext_rpa_correction: Optional[str] = None class Embedding: """Base class for quantum embedding methods. @@ -521,7 +521,7 @@ def e_nuc(self): @property def e_nonlocal(self): - if not self.opts.ext_rpa_correction: + 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)]) @@ -764,11 +764,25 @@ def get_eris_object(self, postscf, fock=None): def 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.") - # Calculate total dRPA energy in N^4 time; this is cheaper than screening calculations. - rpa = ssRIRPA(self.mf, log=self.log) - self.e_rpa, energy_error = rpa.kernel_energy(correction='linear') + + 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] + # Just need to project the RHS and we're done. + self.e_rpa = - 0.5 * einsum("pq,pq->", m0, 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") diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 1fb2244a5..8907a8009 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -360,7 +360,7 @@ 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 calc_loc_erpa(self, m0, amb): + def calc_loc_erpa(self, m0, amb, only_cumulant=False): no, nv = self.cluster.nocc_active, self.cluster.nvir_active nov = no * nv @@ -373,21 +373,36 @@ 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) - 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 + + + if only_cumulant: + + 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), v) + self.log.info("Computed fragment RPA cumulant energy contribution for cluster %s as %s", self._fragment.id, + energy_string(erpa)) + return erpa + + else: + 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)) From c84d8d9f70d143e1e1074a78e7cc20e25f18ce3a Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Mon, 7 Nov 2022 15:43:49 +0000 Subject: [PATCH 19/33] Added fixed expression for nonlocal cumulant correction. --- vayesta/core/qemb/qemb.py | 4 +-- vayesta/solver/hamiltonian.py | 46 ++++++++++++++++++----------------- 2 files changed, 26 insertions(+), 24 deletions(-) diff --git a/vayesta/core/qemb/qemb.py b/vayesta/core/qemb/qemb.py index 0aa314940..6143833c3 100644 --- a/vayesta/core/qemb/qemb.py +++ b/vayesta/core/qemb/qemb.py @@ -777,8 +777,8 @@ def build_screened_eris(self, *args, **kwargs): l_ = np.concatenate([l_, l_], axis=1) m0 = rpa.kernel_moms(0, target_rot=l_)[0][0] - # Just need to project the RHS and we're done. - self.e_rpa = - 0.5 * einsum("pq,pq->", m0, l_) + # 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') diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 8907a8009..8415d46cc 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -4,13 +4,11 @@ import numpy as np import pyscf.lib import pyscf.scf -import scipy.linalg from vayesta.core.util import dot, einsum, OptionsBase, break_into_lines, log_time -from typing import Optional -from vayesta.core.types import Orbitals from vayesta.core.screening import screening_moment, screening_crpa -from vayesta.rpa import ssRPA +from vayesta.core.types import Orbitals +from typing import Optional def is_ham(ham): @@ -374,29 +372,31 @@ def gen_spin_components(mat): m0_aa, m0_ab, m0_bb = gen_spin_components(m0) - if only_cumulant: def compute_e_rrpa(proj): def pr(m): - m = m.reshape((no, nv, no*nv)) + 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), v) - self.log.info("Computed fragment RPA cumulant energy contribution for cluster %s as %s", self._fragment.id, + 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: d_aa, d_ab, d_bb = gen_spin_components(amb) # This should be zero. - assert(abs(d_ab).max() < 1e-10) + assert (abs(d_ab).max() < 1e-10) def compute_e_rrpa(proj): def pr(m): - m = m.reshape((no, nv, no*nv)) + m = m.reshape((no, nv, no * nv)) m = np.tensordot(proj, m, axes=[0, 0]) - return m.reshape((no*nv, no*nv)) + 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()) @@ -425,11 +425,12 @@ def spin_integrate_and_report(m, warn_threshold=1e-6): self.log.info("Largest change in screened interactions due to spin integration: %e", max(dev)) return spat + seris = None 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) + 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) @@ -439,7 +440,7 @@ def spin_integrate_and_report(m, warn_threshold=1e-6): bare_eris = self.get_eris_bare() delta, crpa = screening_crpa.get_frag_deltaW(self.mf, self._fragment, pcoupling=("pcoupled" in self.opts.screening), - only_ov_screened= ("ov" 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!") @@ -459,15 +460,15 @@ def spin_integrate_and_report(m, warn_threshold=1e-6): raise ValueError("Unknown cluster screening protocol: %s" % self.opts.screening) def report_screening(screened, bare, spins): - maxidx = np.unravel_index(np.argmax(abs(screened-bare)), bare.shape) + 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]) + 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( + 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)) @@ -480,13 +481,13 @@ def report_screening(screened, bare, spins): 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) + 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)]) + 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) @@ -766,7 +767,7 @@ def get_polaritonic_shifted_couplings(self): self.log.critical("Polaritonic shifted bosonic fermion-boson couplings break cluster spin symmetry; please" " use an unrestricted formalism.") raise RuntimeError("Polaritonic shifted bosonic fermion-boson couplings break cluster spin symmetry; please" - " use an unrestricted formalism.") + " use an unrestricted formalism.") return couplings[0] def get_eb_dm_polaritonic_shift(self, dm1): @@ -775,6 +776,7 @@ def get_eb_dm_polaritonic_shift(self, dm1): def _add_screening(self, seris_intermed=None, spin_integrate=True): return self.get_eris_bare() + class EB_UClusterHamiltonian(UClusterHamiltonian, EB_RClusterHamiltonian): @dataclasses.dataclass class Options(EB_RClusterHamiltonian.Options): From 0c15bfbcacb01e12cf94c4ac572405f0a5c34310 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Mon, 28 Nov 2022 18:21:19 +0000 Subject: [PATCH 20/33] Added option to match cluster fock matrix, rather than one-body hamiltonian. --- vayesta/core/qemb/fragment.py | 5 +++-- vayesta/core/qemb/qemb.py | 1 + vayesta/solver/hamiltonian.py | 15 +++++++++++---- 3 files changed, 15 insertions(+), 6 deletions(-) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index b23a392df..995683c17 100644 --- a/vayesta/core/qemb/fragment.py +++ b/vayesta/core/qemb/fragment.py @@ -46,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 @@ -1052,7 +1053,7 @@ def get_frag_hamil(self): 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 + raise AbstractMethodError \ No newline at end of file diff --git a/vayesta/core/qemb/qemb.py b/vayesta/core/qemb/qemb.py index 6143833c3..fb0baf742 100644 --- a/vayesta/core/qemb/qemb.py +++ b/vayesta/core/qemb/qemb.py @@ -115,6 +115,7 @@ class Options(OptionsBase): symmetry_mf_tol: float = 1e-5 # Tolerance for mean-field solution screening: Optional[str] = None ext_rpa_correction: Optional[str] = None + match_cluster_fock: bool = False class Embedding: """Base class for quantum embedding methods. diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 8415d46cc..3967cb53f 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -59,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): @@ -111,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: @@ -132,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 @@ -324,12 +329,14 @@ 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( From b26fdc14d7dcb372608acdd5c12aa4c282c3bee9 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Thu, 25 May 2023 12:03:11 +0100 Subject: [PATCH 21/33] Update crpa and add tests for cRPA calculation. --- vayesta/solver/hamiltonian.py | 4 ++-- vayesta/tests/ewf/test_screening.py | 30 +++++++++++++++++++---------- 2 files changed, 22 insertions(+), 12 deletions(-) diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 3967cb53f..433dee917 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -5,7 +5,7 @@ import pyscf.lib import pyscf.scf -from vayesta.core.util import dot, einsum, OptionsBase, break_into_lines, log_time +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 typing import Optional @@ -445,7 +445,7 @@ def spin_integrate_and_report(m, warn_threshold=1e-6): 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.mf, self._fragment, + 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) 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__) From 315cfae42babbca076964e2e04a5e140bc4efde4 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Wed, 31 May 2023 16:54:19 +0100 Subject: [PATCH 22/33] Added Exception when requesting RPA external correction with unrestricted reference. --- vayesta/core/qemb/uqemb.py | 7 +++++++ 1 file changed, 7 insertions(+) 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 From c99489cf2127460b47346e63f805bf2f5aa40918 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Fri, 2 Jun 2023 14:33:05 +0100 Subject: [PATCH 23/33] Remove duplicate log message. --- vayesta/solver/hamiltonian.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 433dee917..8d632346b 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -428,12 +428,8 @@ def spin_integrate_and_report(m, warn_threshold=1e-6): 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)) - else: - self.log.info("Largest change in screened interactions due to spin integration: %e", max(dev)) - return spat - seris = None if self.opts.screening is None: raise ValueError("Attempted to add screening to fragment with no screening protocol specified.") if self.opts.screening == "mrpa": From fd05d9c8dffe8a4af5097eeb354e06c00f2e95b0 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Fri, 2 Jun 2023 14:46:37 +0100 Subject: [PATCH 24/33] Added documentation on what screened interaction options do in the Hamiltonian class. --- vayesta/core/qemb/qemb.py | 2 +- vayesta/solver/hamiltonian.py | 11 ++++++++++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/vayesta/core/qemb/qemb.py b/vayesta/core/qemb/qemb.py index fb0baf742..d0a1dd8c3 100644 --- a/vayesta/core/qemb/qemb.py +++ b/vayesta/core/qemb/qemb.py @@ -113,7 +113,7 @@ 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 diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 8d632346b..153ddd9e8 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -416,11 +416,20 @@ def pr(m): return compute_e_rrpa(po) def add_screening(self, seris_intermed=None): - """Add screened interactions into the Hamiltonian.""" + """ + 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. + + seris_intermed is only required for mRPA interactions. + """ self._seris = self._add_screening(seris_intermed, spin_integrate=True) def _add_screening(self, seris_intermed=None, spin_integrate=True): + 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 From b6c92903095b27df552389c3bd0d88bd12ce3cbf Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Fri, 9 Jun 2023 14:37:03 +0100 Subject: [PATCH 25/33] Restored missing import in ssURPA. --- vayesta/rpa/ssurpa.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vayesta/rpa/ssurpa.py b/vayesta/rpa/ssurpa.py index 0c605cdc5..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): From 0cd658b5361536f179b26199474bc4a1b98ede40 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Fri, 28 Jul 2023 14:33:48 +0100 Subject: [PATCH 26/33] Moved calculation of external RPA corrections to separate instance method. --- vayesta/core/qemb/fragment.py | 15 +++++++++++++-- vayesta/ewf/fragment.py | 5 +++-- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index 995683c17..353ed2345 100644 --- a/vayesta/core/qemb/fragment.py +++ b/vayesta/core/qemb/fragment.py @@ -1037,16 +1037,27 @@ def get_solver(self, solver=None): cluster_solver = solver_cls(cl_ham, **solver_opts) if self.opts.screening is not None: cluster_solver.hamil.add_screening(self._seris_ov) - e_loc_rpa = None if self.base.opts.ext_rpa_correction: 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 = cl_ham.calc_loc_erpa(*self._seris_ov[1:], cumulant) - return cluster_solver, e_loc_rpa + return cluster_solver + + def get_ext_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 self + 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), diff --git a/vayesta/ewf/fragment.py b/vayesta/ewf/fragment.py index 4126f74ba..96e5e4399 100644 --- a/vayesta/ewf/fragment.py +++ b/vayesta/ewf/fragment.py @@ -189,8 +189,9 @@ def kernel(self, solver=None, init_guess=None): init_guess = self.get_init_guess(init_guess, solver, cluster) # Create solver object - cluster_solver, e_corr_rpa = self.get_solver(solver) - + cluster_solver = self.get_solver(solver) + # Calculate cluster energy at the level of RPA. + e_corr_rpa = self.get_ext_rpa_correction(cluster_solver.hamil) # --- Chemical potential cpt_frag = self.base.opts.global_frag_chempot if self.opts.nelectron_target is not None: From aead9666a3a4e868e346e73e01fb0c83b2df8627 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Fri, 28 Jul 2023 14:37:41 +0100 Subject: [PATCH 27/33] Avoided clash in e_nonlocal namespace. --- vayesta/edmet/edmet.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/vayesta/edmet/edmet.py b/vayesta/edmet/edmet.py index 5e2dda8e9..a399bc16c 100644 --- a/vayesta/edmet/edmet.py +++ b/vayesta/edmet/edmet.py @@ -54,6 +54,17 @@ def eps(self): eps = eps.reshape(-1) return eps, eps + @property + def e_nonlocal(self): + try: + return self._e_nonlocal + except AttributeError: + return 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 From 13366d3b6f142717a49af499b43265fdde698246 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Fri, 28 Jul 2023 15:02:04 +0100 Subject: [PATCH 28/33] Remove defunct code. --- vayesta/core/qemb/fragment.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index 353ed2345..aefc32eae 100644 --- a/vayesta/core/qemb/fragment.py +++ b/vayesta/core/qemb/fragment.py @@ -1037,11 +1037,6 @@ def get_solver(self, solver=None): cluster_solver = solver_cls(cl_ham, **solver_opts) if self.opts.screening is not None: cluster_solver.hamil.add_screening(self._seris_ov) - if self.base.opts.ext_rpa_correction: - 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 = cl_ham.calc_loc_erpa(*self._seris_ov[1:], cumulant) return cluster_solver def get_ext_rpa_correction(self, hamil=None): From fd94a198016b4ee97f541ca419fef4254172c99c Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Fri, 28 Jul 2023 15:18:38 +0100 Subject: [PATCH 29/33] Added warnings for calculations using cRPA screened interactions. --- vayesta/core/screening/screening_crpa.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/vayesta/core/screening/screening_crpa.py b/vayesta/core/screening/screening_crpa.py index 69f92c608..0de9fd65d 100644 --- a/vayesta/core/screening/screening_crpa.py +++ b/vayesta/core/screening/screening_crpa.py @@ -70,6 +70,8 @@ def get_frag_deltaW(mf, fragment, pcoupling=True, only_ov_screened=False, log=No """ 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: From d621499e51962359f6a07df017939b97ae6dddf1 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Fri, 28 Jul 2023 16:48:13 +0100 Subject: [PATCH 30/33] Added tests for RPA external corrections. --- vayesta/tests/ewf/test_rpa_correction.py | 47 ++++++++++++++++++++++++ vayesta/tests/testsystems.py | 1 + 2 files changed, 48 insertions(+) create mode 100644 vayesta/tests/ewf/test_rpa_correction.py diff --git a/vayesta/tests/ewf/test_rpa_correction.py b/vayesta/tests/ewf/test_rpa_correction.py new file mode 100644 index 000000000..d0d79acf8 --- /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) \ No newline at end of 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") From c825f43b7e9c786a3d3b147adcd46214033c0dd1 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Fri, 28 Jul 2023 16:48:36 +0100 Subject: [PATCH 31/33] Rename function to generate local RPA correlation energy. --- vayesta/core/qemb/fragment.py | 2 +- vayesta/ewf/fragment.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index aefc32eae..e9c6779b5 100644 --- a/vayesta/core/qemb/fragment.py +++ b/vayesta/core/qemb/fragment.py @@ -1039,7 +1039,7 @@ def get_solver(self, solver=None): cluster_solver.hamil.add_screening(self._seris_ov) return cluster_solver - def get_ext_rpa_correction(self, hamil=None): + def get_local_rpa_correction(self, hamil=None): e_loc_rpa = None if self.base.opts.ext_rpa_correction: hamil = hamil or self.hamil diff --git a/vayesta/ewf/fragment.py b/vayesta/ewf/fragment.py index 96e5e4399..e0de257f3 100644 --- a/vayesta/ewf/fragment.py +++ b/vayesta/ewf/fragment.py @@ -191,7 +191,7 @@ 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_ext_rpa_correction(cluster_solver.hamil) + 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: From ed6d8381061bfc8fab776cfac24fd15773a980bc Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Fri, 28 Jul 2023 17:16:59 +0100 Subject: [PATCH 32/33] Added example for screening and RPA external corrections. --- .../molecules/08-screening-rpa-corrections.py | 49 +++++++++++++++++++ 1 file changed, 49 insertions(+) create mode 100644 examples/ewf/molecules/08-screening-rpa-corrections.py 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)) From 0270c476e3a191ead914f8bc74904bfa0464e397 Mon Sep 17 00:00:00 2001 From: Charles Scott Date: Thu, 3 Aug 2023 15:44:14 +0100 Subject: [PATCH 33/33] Response to Ollie's review. --- vayesta/core/qemb/fragment.py | 6 ++---- vayesta/edmet/edmet.py | 9 +++++---- vayesta/tests/ewf/test_rpa_correction.py | 2 +- 3 files changed, 8 insertions(+), 9 deletions(-) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index e9c6779b5..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 @@ -1051,8 +1051,6 @@ def get_local_rpa_correction(self, hamil=None): def get_frag_hamil(self): if self.opts.screening is not None: - #if self - 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), @@ -1062,4 +1060,4 @@ def get_frag_hamil(self): cache_eris=self.opts.store_eris, match_fock=self.opts.match_cluster_fock) def get_solver_options(self, *args, **kwargs): - raise AbstractMethodError \ No newline at end of file + raise AbstractMethodError diff --git a/vayesta/edmet/edmet.py b/vayesta/edmet/edmet.py index a399bc16c..6e73338ba 100644 --- a/vayesta/edmet/edmet.py +++ b/vayesta/edmet/edmet.py @@ -56,10 +56,7 @@ def eps(self): @property def e_nonlocal(self): - try: - return self._e_nonlocal - except AttributeError: - return 0.0 + return self._e_nonlocal or 0.0 @e_nonlocal.setter def e_nonlocal(self, value): @@ -386,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/tests/ewf/test_rpa_correction.py b/vayesta/tests/ewf/test_rpa_correction.py index d0d79acf8..63e9265db 100644 --- a/vayesta/tests/ewf/test_rpa_correction.py +++ b/vayesta/tests/ewf/test_rpa_correction.py @@ -44,4 +44,4 @@ def test_rpa_correction(self): def test_cumulant_correction(self): enl, erpa = self.get_nl_energies("cumulant", "full") - self.assertAlmostEqual(enl, 0.0) \ No newline at end of file + self.assertAlmostEqual(enl, 0.0)