diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 2f401baed..38182ad3a 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -34,12 +34,10 @@ jobs: run: | python -m pip install wheel --user python -m pip install setuptools --upgrade --user - python -m pip install https://github.com/BoothGroup/dyson/archive/master.zip python -m pip install pybind11[global] export pybind11_DIR=$(python -m pybind11 --cmakedir) python -m pip install -e .[dmet,ebcc,pygnme] --user - - name: Run unit tests run: | export LD_LIBRARY_PATH=$(dirname $(python -c 'import pygnme;print(pygnme.__path__[0])') ):$LD_LIBRARY_PATH diff --git a/docs/source/quickstart/edmetfinite.py b/docs/source/quickstart/edmetfinite.py index 091e29e0b..180e7d7de 100644 --- a/docs/source/quickstart/edmetfinite.py +++ b/docs/source/quickstart/edmetfinite.py @@ -26,7 +26,7 @@ # The KS object needs to be converted to a HF object: lda = lda.to_rhf() -emb_lda = vayesta.edmet.EDMET(lda, dmet_threshold=1e-12, solver="EBCCSD", +emb_lda = vayesta.edmet.EDMET(lda, bath_options=dict(dmet_threshold=1e-12), solver="CCSD-S-1-1", oneshot=True, make_dd_moments=False) with emb_lda.iao_fragmentation() as f: @@ -42,7 +42,7 @@ # The KS opbject needs to be converted to a HF object: b3lyp = b3lyp.to_rhf() -emb_b3lyp = vayesta.edmet.EDMET(b3lyp, dmet_threshold=1e-12, solver="EBCCSD", oneshot=True, make_dd_moments=False) +emb_b3lyp = vayesta.edmet.EDMET(b3lyp, bath_options=dict(dmet_threshold=1e-12), solver="CCSD-S-1-1", oneshot=True, make_dd_moments=False) with emb_b3lyp.iao_fragmentation() as f: f.add_all_atomic_fragments() emb_b3lyp.kernel() @@ -51,7 +51,7 @@ hf = pyscf.scf.RHF(mol).density_fit() hf.kernel() -emb_hf = vayesta.edmet.EDMET(hf, dmet_threshold=1e-12, solver="EBCCSD", +emb_hf = vayesta.edmet.EDMET(hf, bath_options=dict(dmet_threshold=1e-12), solver="CCSD-S-1-1", oneshot=True, make_dd_moments=False) with emb_hf.iao_fragmentation() as f: @@ -65,7 +65,7 @@ print("E(LDA)= %+16.8f Ha" % lda.e_tot) print("E(Emb. CCSD @LDA)= %+16.8f Ha" % emb_lda.e_tot) print("E(HF @LDA)= %+16.8f Ha" % emb_lda.e_mf) -print("E(CCSD)= %+16.8f Ha" % cc.) +print("E(CCSD)= %+16.8f Ha" % cc.e_tot) print("E(B3LYP)= %+16.8f Ha" % b3lyp.e_tot) print("E(HF)= %+16.8f Ha" % hf.e_tot) print("E(HF @B3LYP)= %+16.8f Ha" % emb_b3lyp.e_mf) diff --git a/examples/ewf/molecules/26-ebcc-solvers.py b/examples/ewf/molecules/26-ebcc-solvers.py new file mode 100644 index 000000000..90f2d7062 --- /dev/null +++ b/examples/ewf/molecules/26-ebcc-solvers.py @@ -0,0 +1,53 @@ +import pyscf +import pyscf.gto +import pyscf.scf +import pyscf.mcscf + +import vayesta +import vayesta.ewf + +mol = pyscf.gto.Mole() +mol.atom = ['N 0 0 0', 'N 0 0 2'] +mol.basis = 'aug-cc-pvdz' +mol.output = 'pyscf.out' +mol.build() + +# Hartree-Fock +mf = pyscf.scf.RHF(mol) +mf.kernel() + +# Reference CASCI +casci = pyscf.mcscf.CASCI(mf, 8, 10) +casci.kernel() + +# Reference CASSCF +casscf = pyscf.mcscf.CASSCF(mf, 8, 10) +casscf.kernel() + +def get_emb_result(ansatz, bathtype='full'): + # Uses fastest available solver for given ansatz; PySCF if available, otherwise ebcc. + emb = vayesta.ewf.EWF(mf, solver=ansatz, bath_options=dict(bathtype=bathtype), + solver_options=dict(solve_lambda=False)) + # Both these alternative specifications will always use an ebcc solver. + # Note that the capitalization of the solver name other than the ansatz is arbitrary. + #emb = vayesta.ewf.EWF(mf, solver=f'EB{ansatz}', bath_options=dict(bathtype=bathtype), + # solver_options=dict(solve_lambda=False)) + #emb = vayesta.ewf.EWF(mf, solver='ebcc', bath_options=dict(bathtype=bathtype), + # solver_options=dict(solve_lambda=False, ansatz=ansatz)) + + with emb.iao_fragmentation() as f: + with f.rotational_symmetry(2, "y", center=(0, 0, 1)): + f.add_atomic_fragment(0) + emb.kernel() + return emb.e_tot + +e_ccsd = get_emb_result('CCSD', 'full') +e_ccsdt = get_emb_result('CCSDT', 'dmet') +e_ccsdtprime = get_emb_result("CCSDt'", 'full') + +print("E(HF)= %+16.8f Ha" % mf.e_tot) +print("E(CASCI)= %+16.8f Ha" % casci.e_tot) +print("E(CASSCF)= %+16.8f Ha" % casscf.e_tot) +print("E(CCSD, complete)= %+16.8f Ha" % e_ccsd) +print("E(emb. CCSDT, DMET CAS)= %+16.8f Ha" % e_ccsdt) +print("E(emb. CCSDt', complete+DMET active space)= %+16.8f Ha" % e_ccsdtprime) diff --git a/examples/ewf/molecules/35-screening-rpa-corrections.py b/examples/ewf/molecules/35-screening-rpa-corrections.py new file mode 100644 index 000000000..17375d4ef --- /dev/null +++ b/examples/ewf/molecules/35-screening-rpa-corrections.py @@ -0,0 +1,49 @@ +import pyscf.gto +import pyscf.scf +import pyscf.cc +import vayesta.ewf + + +mol = pyscf.gto.Mole() +mol.atom = """ +O 0.0000 0.0000 0.1173 +H 0.0000 0.7572 -0.4692 +H 0.0000 -0.7572 -0.4692 +""" +mol.basis = 'cc-pVTZ' +mol.output = 'pyscf.out' +mol.build() + +# Hartree-Fock +mf = pyscf.scf.RHF(mol).density_fit() +mf.kernel() + +# Reference full system CCSD: +cc = pyscf.cc.CCSD(mf) +cc.kernel() + +eta = 1e-6 + +# Embedded CCSD calculation with bare interactions and no energy correction. +emb_bare = vayesta.ewf.EWF(mf, bath_options=dict(threshold=eta)) +emb_bare.kernel() + +# Embedded CCSD with mRPA screened interactions and RPA cumulant approximation for nonlocal correlations. +emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=eta), screening="mrpa", ext_rpa_correction="cumulant") +emb.kernel() +e_nonlocal_cumulant = emb.e_nonlocal +# Embedded CCSD with mRPA screened interactions and delta RPA approximation for nonlocal correlations. +emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=eta), screening="mrpa", ext_rpa_correction="erpa") +emb.kernel() +e_nonlocal_erpa = emb.e_nonlocal + +# Note that mRPA screening and external corrections often cancel with each other in the case of the energy. +print("E(CCSD)= %+16.8f Ha" % cc.e_tot) +print("E(RPA)= %+16.8f Ha (error= %+.8f Ha)" % (emb.e_mf + emb.e_rpa, + emb.e_mf + emb.e_rpa - cc.e_tot)) +print("E(Emb. CCSD)= %+16.8f Ha (error= %+.8f Ha)" % (emb_bare.e_tot, emb_bare.e_tot-cc.e_tot)) +print("E(Emb. Screened CCSD)= %+16.8f Ha (error= %+.8f Ha)" % (emb.e_tot, emb.e_tot-cc.e_tot)) +print("E(Emb. Screened CCSD + \Delta E_k)= %+16.8f Ha (error= %+.8f Ha)" % (emb.e_tot+e_nonlocal_cumulant, + emb.e_tot+e_nonlocal_cumulant-cc.e_tot)) +print("E(Emb. Screened CCSD + \Delta RPA)= %+16.8f Ha (error= %+.8f Ha)" % (emb.e_tot+e_nonlocal_erpa, + emb.e_tot+e_nonlocal_erpa-cc.e_tot)) diff --git a/examples/ewf/molecules/36-crpa-cas-screening.py b/examples/ewf/molecules/36-crpa-cas-screening.py new file mode 100644 index 000000000..c15022df0 --- /dev/null +++ b/examples/ewf/molecules/36-crpa-cas-screening.py @@ -0,0 +1,45 @@ +import pyscf.gto +import pyscf.scf +import pyscf.cc +import pyscf.mcscf +import vayesta.ewf +from vayesta.misc import molecules + +mol = pyscf.gto.Mole() +mol.atom = molecules.arene(6) +mol.basis = 'cc-pVDZ' +mol.output = 'pyscf.out' +mol.build() + +mf = pyscf.scf.RHF(mol).density_fit() +mf.kernel() + +ncas = (4,4) + +# Reference CASCI +mycasci = pyscf.mcscf.CASCI(mf, *ncas) +mycasci.kernel() + +# Equivalent CAS calculation with bare interactions +casci_bare = vayesta.ewf.EWF(mf, bath_options=dict(bathtype="dmet"), screening=None, solver="FCI") +with casci_bare.cas_fragmentation() as f: + f.add_cas_fragment(*ncas) +casci_bare.kernel() + +# CAS calculation with cRPA screening +casci_crpa = vayesta.ewf.EWF(mf, bath_options=dict(bathtype="dmet"), screening="crpa", solver="FCI") +with casci_crpa.cas_fragmentation() as f: + f.add_cas_fragment(*ncas) +casci_crpa.kernel() + +# CAS calculation with mRPA screening +casci_mrpa = vayesta.ewf.EWF(mf, bath_options=dict(bathtype="dmet"), screening="mrpa", solver="FCI") +with casci_mrpa.cas_fragmentation() as f: + f.add_cas_fragment(*ncas) +casci_mrpa.kernel() + +print("E(HF)= %+16.8f Ha" % mf.e_tot) +print("ΔE(CASCI)= %+16.8f Ha" % (mycasci.e_tot-mf.e_tot)) +print("ΔE(Emb. CASCI, bare interactions)= %+16.8f Ha (diff= %+.8f Ha)" % (casci_bare.e_tot-mf.e_tot, casci_bare.e_tot-mycasci.e_tot)) +print("ΔE(Emb. CASCI, cRPA screening)= %+16.8f Ha (diff= %+.8f Ha)" % (casci_crpa.e_tot-mf.e_tot, casci_crpa.e_tot-mycasci.e_tot)) +print("ΔE(Emb. CASCI, mRPA screening)= %+16.8f Ha (diff= %+.8f Ha)" % (casci_mrpa.e_tot-mf.e_tot, casci_mrpa.e_tot-mycasci.e_tot)) diff --git a/examples/ewf/solids/01-simple-sym.py b/examples/ewf/solids/01-simple-sym.py index 72d974a44..51a137a3a 100644 --- a/examples/ewf/solids/01-simple-sym.py +++ b/examples/ewf/solids/01-simple-sym.py @@ -29,7 +29,8 @@ emb.kernel() # Some versions of pyscf may encounter an error when generating IAOs with space_group_symmetry=True. -# To avoid this, turn off space group symmetry before generating IAOs, as below: +# To avoid this, either update to a version more recent that 2.3, or turn off space group symmetry before generating +# IAOs, as in the lines below. # emb.mf.mol.space_group_symmetry = False # emb.kernel() diff --git a/examples/ewf/solids/20-mRPA-interactions.py b/examples/ewf/solids/20-mRPA-interactions.py new file mode 100644 index 000000000..f4cc9270f --- /dev/null +++ b/examples/ewf/solids/20-mRPA-interactions.py @@ -0,0 +1,44 @@ +import pyscf +import pyscf.pbc +import pyscf.pbc.scf +import pyscf.pbc.dft +import pyscf.pbc.cc + +import vayesta +import vayesta.ewf +from vayesta.misc import solids + +cell = pyscf.pbc.gto.Cell() +cell.a, cell.atom = solids.graphene() +cell.basis = 'sto-3g' +cell.output = 'pyscf.out' +cell.dimension = 2 +cell.space_group_symmetry = True +cell.symmorphic = True +cell.build() + +# HF +kmesh = [2,2,1] +kpts = cell.make_kpts(kmesh, space_group_symmetry=True, time_reversal_symmetry=True, symmorphic=True) +hf = pyscf.pbc.scf.KRHF(cell, cell.make_kpts([2,2,1])) +hf = hf.density_fit() +hf.kernel() +# This may be required to avoid issues discussed in 01-simple-sym.py. +hf.mol.space_group_symmetry = False +# Run embedded +emb_bare = vayesta.ewf.EWF(hf, bath_options=dict(bathtype="rpa", threshold=1e-2), screening=None) +emb_bare.kernel() +# Run calculation using screened interactions and cumulant correction for nonlocal energy. +emb_mrpa = vayesta.ewf.EWF(hf, bath_options=dict(bathtype="rpa", threshold=1e-2), screening="mrpa", + ext_rpa_correction="cumulant") +emb_mrpa.kernel() + +# Reference full system CCSD: +cc = pyscf.pbc.cc.KCCSD(hf) +cc.kernel() + +print("Error(HF)= %+16.8f Ha" % (hf.e_tot - cc.e_tot)) +print("Error(Emb. bare CCSD)= %+16.8f Ha" % (emb_bare.e_tot - cc.e_tot)) +print("Error(Emb. mRPA CCSD)= %+16.8f Ha" % (emb_mrpa.e_tot - cc.e_tot)) +print("Error(Emb. mRPA CCSD + ΔE_k)= %+16.8f Ha" % (emb_mrpa.e_tot+emb_mrpa.e_nonlocal-cc.e_tot)) +print("E(CCSD)= %+16.8f Ha" % cc.e_tot) diff --git a/vayesta/core/bath/__init__.py b/vayesta/core/bath/__init__.py index 34958b76d..fb76ddc26 100644 --- a/vayesta/core/bath/__init__.py +++ b/vayesta/core/bath/__init__.py @@ -10,6 +10,8 @@ from vayesta.core.bath.bno import MP2_BNO_Bath as MP2_Bath_RHF from vayesta.core.bath.bno import UMP2_BNO_Bath as MP2_Bath_UHF +from vayesta.core.bath.rpa import RPA_BNO_Bath + from vayesta.core.bath.r2bath import R2_Bath_RHF def DMET_Bath(fragment, *args, **kwargs): @@ -30,6 +32,12 @@ def MP2_Bath(fragment, *args, **kwargs): if fragment.base.is_uhf: return MP2_Bath_UHF(fragment, *args, **kwargs) +def RPA_Bath(fragment, *args, **kwargs): + if fragment.base.is_rhf: + return RPA_BNO_Bath(fragment, *args, **kwargs) + if fragment.base.is_uhf: + raise NotImplementedError + def R2_Bath(fragment, *args, **kwargs): if fragment.base.is_rhf: return R2_Bath_RHF(fragment, *args, **kwargs) diff --git a/vayesta/core/bath/dmet.py b/vayesta/core/bath/dmet.py index e4a1b0398..fbecb3856 100644 --- a/vayesta/core/bath/dmet.py +++ b/vayesta/core/bath/dmet.py @@ -47,7 +47,8 @@ def kernel(self): # --- Separate occupied and virtual in cluster cluster = [self.c_frag, c_dmet] self.base._check_orthonormal(*cluster, mo_name='cluster MO') - c_cluster_occ, c_cluster_vir = self.fragment.diagonalize_cluster_dm(*cluster, tol=2*self.dmet_threshold) + tol = self.base.opts.bath_options['occupation_tolerance'] + c_cluster_occ, c_cluster_vir = self.fragment.diagonalize_cluster_dm(*cluster, tol=tol) # Canonicalize c_cluster_occ = self.fragment.canonicalize_mo(c_cluster_occ)[0] c_cluster_vir = self.fragment.canonicalize_mo(c_cluster_vir)[0] diff --git a/vayesta/core/bath/rpa.py b/vayesta/core/bath/rpa.py new file mode 100644 index 000000000..be9f7d55f --- /dev/null +++ b/vayesta/core/bath/rpa.py @@ -0,0 +1,96 @@ +from vayesta.core.bath.bno import BNO_Threshold, BNO_Bath + +from vayesta.core.util import AbstractMethodError, brange, dot, einsum, fix_orbital_sign, hstack, time_string, timer +from vayesta.core.types import Cluster + +from vayesta.rpa.rirpa import ssRIdRRPA + +from vayesta.core.eris import get_cderi +from vayesta.core import spinalg +from vayesta.core.bath import helper + + +import numpy as np + +class RPA_BNO_Bath(BNO_Bath): + + def __init__(self, *args, project_dmet_order=0, project_dmet_mode='full', project_dmet=None, **kwargs): + self.project_dmet_order = project_dmet_order + self.project_dmet_mode = project_dmet_mode + super().__init__(*args, **kwargs) + + + def make_bno_coeff(self, cderis=None): + """Construct RPA bath natural orbital coefficients and occupation numbers. + + This routine works for both for spin-restricted and unrestricted. + + Parameters + ---------- + cderis: cderis in the particle-hole space. + + Returns + ------- + c_bno: (n(AO), n(BNO)) array + Bath natural orbital coefficients. + n_bno: (n(BNO)) array + Bath natural orbital occupation numbers. + """ + + t_init = timer() + + if cderis is None: + cderis = get_cderi(self.base, (self.base.mo_coeff_occ, self.base.mo_coeff_vir), compact=False) + + if self.occtype == "occupied": + proj = dot(self.dmet_bath.c_cluster_vir.T, self.base.get_ovlp(), self.fragment.c_frag, + self.fragment.c_frag.T, self.base.get_ovlp(), self.dmet_bath.c_cluster_vir) + + rot_vir = dot(self.dmet_bath.c_cluster_vir.T, self.base.get_ovlp(), self.base.mo_coeff_vir) + rot_occ = np.eye(self.base.nocc) + else: + proj = dot(self.dmet_bath.c_cluster_occ.T, self.base.get_ovlp(), self.fragment.c_frag, self.fragment.c_frag.T, + self.base.get_ovlp(), self.dmet_bath.c_cluster_occ) + rot_occ = dot(self.dmet_bath.c_cluster_occ.T, self.base.get_ovlp(), self.base.mo_coeff_occ) + rot_vir = np.eye(self.base.nvir) + + loc_excit_shape = (rot_occ.shape[0], rot_vir.shape[0]) + + # Get target rotation in particle-hole excitation space. + # This is of size O(N), so this whole procedure scales as O(N^4) + + target_rot = einsum("ij,ab->iajb", rot_occ, rot_vir) + target_rot = target_rot.reshape(np.product(target_rot.shape[:2]), np.product(target_rot.shape[2:])) + + t0 = timer() + myrpa = ssRIdRRPA(self.base.mf, lov=cderis) + # This initially calculates the spin-summed zeroth moment, then deducts the spin-dependent component and + # accounts for factor of two from different spin channels. + m0 = (myrpa.kernel_moms(0, target_rot=target_rot, return_spatial=True)[0][0] - target_rot) / 2.0 + m0 = -dot(m0, target_rot.T).reshape(loc_excit_shape + loc_excit_shape) + if self.occtype == "occupied": + corr_dm = einsum("iajb,ab->ij", m0, proj) + else: + corr_dm = einsum("iajb,ij->ab", m0, proj) + t_eval = timer()-t0 + + corr_dm = (corr_dm + corr_dm.T) / 2 + + # --- Diagonalize environment-environment block + if self.occtype == 'occupied': + corr_dm = self._rotate_dm(corr_dm, dot(self.dmet_bath.c_env_occ.T, self.base.get_ovlp(), self.base.mo_coeff_occ)) + elif self.occtype == 'virtual': + corr_dm = self._rotate_dm(corr_dm, dot(self.dmet_bath.c_env_vir.T, self.base.get_ovlp(), self.base.mo_coeff_vir)) + t0 = timer() + r_bno, n_bno = self._diagonalize_dm(corr_dm) + t_diag = timer()-t0 + c_bno = spinalg.dot(self.c_env, r_bno) + c_bno = fix_orbital_sign(c_bno)[0] + + self.log.timing("Time RPA bath: evaluation= %s diagonal.= %s total= %s", + *map(time_string, (t_eval, t_diag, (timer()-t_init)))) + + if min(n_bno) < 0.0: + self.log.critical("Negative bath occupation number encountered: %s", n_bno) + + return c_bno, n_bno, 0.0 diff --git a/vayesta/core/eris.py b/vayesta/core/eris.py new file mode 100644 index 000000000..39a7015db --- /dev/null +++ b/vayesta/core/eris.py @@ -0,0 +1,136 @@ +from pyscf.mp.mp2 import _mo_without_core +from vayesta.core.ao2mo import kao2gmo_cderi +from vayesta.core.ao2mo import postscf_ao2mo +from vayesta.core.ao2mo import postscf_kao2gmo + +import numpy as np +from vayesta.core.util import * +import pyscf.lib + + +def get_cderi(emb, mo_coeff, compact=False, blksize=None): + """Get density-fitted three-center integrals in MO basis.""" + if compact: + raise NotImplementedError() + if emb.kdf is not None: + return kao2gmo_cderi(emb.kdf, mo_coeff) + + if np.ndim(mo_coeff[0]) == 1: + mo_coeff = (mo_coeff, mo_coeff) + + nao = emb.mol.nao + try: + naux = (emb.df.auxcell.nao if hasattr(emb.df, 'auxcell') else emb.df.auxmol.nao) + except AttributeError: + naux = emb.df.get_naoaux() + + cderi = np.zeros((naux, mo_coeff[0].shape[-1], mo_coeff[1].shape[-1])) + cderi_neg = None + if blksize is None: + blksize = int(1e9 / naux * nao * nao * 8) + # PBC: + if hasattr(emb.df, 'sr_loop'): + blk0 = 0 + for labr, labi, sign in emb.df.sr_loop(compact=False, blksize=blksize): + assert np.allclose(labi, 0) + assert (cderi_neg is None) # There should be only one block with sign -1 + labr = labr.reshape(-1, nao, nao) + if (sign == 1): + blk1 = (blk0 + labr.shape[0]) + blk = np.s_[blk0:blk1] + blk0 = blk1 + cderi[blk] = einsum('Lab,ai,bj->Lij', labr, mo_coeff[0], mo_coeff[1]) + elif (sign == -1): + cderi_neg = einsum('Lab,ai,bj->Lij', labr, mo_coeff[0], mo_coeff[1]) + return cderi, cderi_neg + # No PBC: + blk0 = 0 + for lab in emb.df.loop(blksize=blksize): + blk1 = (blk0 + lab.shape[0]) + blk = np.s_[blk0:blk1] + blk0 = blk1 + lab = pyscf.lib.unpack_tril(lab) + cderi[blk] = einsum('Lab,ai,bj->Lij', lab, mo_coeff[0], mo_coeff[1]) + return cderi, None + + +@log_method() +def get_eris_array(emb, mo_coeff, compact=False): + """Get electron-repulsion integrals in MO basis as a NumPy array. + + Parameters + ---------- + mo_coeff: [list(4) of] (n(AO), n(MO)) array + MO coefficients. + + Returns + ------- + eris: (n(MO), n(MO), n(MO), n(MO)) array + Electron-repulsion integrals in MO basis. + """ + # PBC with k-points: + if emb.kdf is not None: + if compact: + raise NotImplementedError + if np.ndim(mo_coeff[0]) == 1: + mo_coeff = 4 * [mo_coeff] + cderi1, cderi1_neg = kao2gmo_cderi(emb.kdf, mo_coeff[:2]) + if (mo_coeff[0] is mo_coeff[2]) and (mo_coeff[1] is mo_coeff[3]): + cderi2, cderi2_neg = cderi1, cderi1_neg + else: + cderi2, cderi2_neg = kao2gmo_cderi(emb.kdf, mo_coeff[2:]) + eris = einsum('Lij,Lkl->ijkl', cderi1.conj(), cderi2) + if cderi1_neg is not None: + eris -= einsum('Lij,Lkl->ijkl', cderi1_neg.conj(), cderi2_neg) + return eris + # Molecules and Gamma-point PBC: + if hasattr(emb.mf, 'with_df') and emb.mf.with_df is not None: + eris = emb.mf.with_df.ao2mo(mo_coeff, compact=compact) + elif emb.mf._eri is not None: + eris = pyscf.ao2mo.kernel(emb.mf._eri, mo_coeff, compact=compact) + else: + eris = emb.mol.ao2mo(mo_coeff, compact=compact) + if not compact: + if isinstance(mo_coeff, np.ndarray) and mo_coeff.ndim == 2: + shape = 4 * [mo_coeff.shape[-1]] + else: + shape = [mo.shape[-1] for mo in mo_coeff] + eris = eris.reshape(shape) + return eris + + +@log_method() +def get_eris_object(emb, postscf, fock=None): + """Get ERIs for post-SCF methods. + + For folded PBC calculations, this folds the MO back into k-space + and contracts with the k-space three-center integrals.. + + Parameters + ---------- + postscf: one of the following PySCF methods: MP2, CCSD, RCCSD, DFCCSD + Post-SCF method with attribute mo_coeff set. + + Returns + ------- + eris: _ChemistsERIs + ERIs which can be used for the respective post-scf method. + """ + if fock is None: + if isinstance(postscf, pyscf.mp.mp2.MP2): + fock = emb.get_fock() + elif isinstance(postscf, (pyscf.ci.cisd.CISD, pyscf.cc.ccsd.CCSD)): + fock = emb.get_fock(with_exxdiv=False) + else: + raise ValueError("Unknown post-SCF method: %r", type(postscf)) + # For MO energies, always use get_fock(): + mo_act = _mo_without_core(postscf, postscf.mo_coeff) + mo_energy = einsum('ai,ab,bi->i', mo_act, emb.get_fock(), mo_act) + e_hf = emb.mf.e_tot + + # Fold MOs into k-point sampled primitive cell, to perform efficient AO->MO transformation: + if emb.kdf is not None: + return postscf_kao2gmo(postscf, emb.kdf, fock=fock, mo_energy=mo_energy, e_hf=e_hf) + # Regular AO->MO transformation + eris = postscf_ao2mo(postscf, fock=fock, mo_energy=mo_energy, e_hf=e_hf) + return eris diff --git a/vayesta/core/helper.py b/vayesta/core/helper.py index ff261e2eb..ade0b6dab 100644 --- a/vayesta/core/helper.py +++ b/vayesta/core/helper.py @@ -94,14 +94,7 @@ def unpack_arrays(packed, dtype=float, maxdim=8): #np.asarray([True, False, False]) ] pack = pack_arrays(*arrays_in) - - print(np.sum([sys.getsizeof(x) for x in arrays_in])) - print(sys.getsizeof(pack)) - arrays_out = unpack_arrays(pack) - print(len(arrays_out)) - print(arrays_out) - assert len(arrays_in) == len(arrays_out) for i, x in enumerate(arrays_in): if x is None: diff --git a/vayesta/core/qemb/corrfunc.py b/vayesta/core/qemb/corrfunc.py index d3304dea7..1b219e7fa 100644 --- a/vayesta/core/qemb/corrfunc.py +++ b/vayesta/core/qemb/corrfunc.py @@ -7,7 +7,7 @@ from vayesta.misc import corrfunc -def get_corrfunc_mf(emb, kind, dm1=None, atoms=None, projection='sao'): +def get_corrfunc_mf(emb, kind, dm1=None, atoms=None, projection='sao', orbital_filter=None): """dm1 in MO basis""" if emb.spinsym == 'unrestricted' and kind.lower() in ('n,n', 'dn,dn'): raise NotImplementedError @@ -34,7 +34,7 @@ def get_corrfunc_mf(emb, kind, dm1=None, atoms=None, projection='sao'): func = funcs.get(kind.lower()) if func is None: raise ValueError(kind) - atoms1, atoms2, proj = emb._get_atom_projectors(atoms, projection) + atoms1, atoms2, proj = emb._get_atom_projectors(atoms, projection, orbital_filter=orbital_filter) corr = np.zeros((len(atoms1), len(atoms2))) for a, atom1 in enumerate(atoms1): for b, atom2 in enumerate(atoms2): diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index a6a152de5..2e47cd446 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 @@ -24,10 +24,13 @@ from vayesta.core.bath import MP2_Bath from vayesta.core.bath import Full_Bath from vayesta.core.bath import R2_Bath +from vayesta.core.bath import RPA_Bath + # Other from vayesta.misc.cubefile import CubeFile from vayesta.mpi import mpi from vayesta.solver import get_solver_class, check_solver_config, ClusterHamiltonian +from vayesta.core.screening.screening_crpa import get_frag_W # Get MPI rank of fragment get_fragment_mpi_rank = lambda *args : args[0].mpi_rank @@ -45,6 +48,7 @@ class Options(OptionsBase): store_eris: bool = None # If True, ERIs will be cached in Fragment.hamil dm_with_frozen: bool = None # TODO: is still used? screening: typing.Optional[str] = None + match_cluster_fock: bool = None # Fragment specific # ----------------- # Auxiliary fragments are treated before non-auxiliary fragments, but do not contribute to expectation values @@ -73,6 +77,7 @@ class Results: converged: bool = None # True, if solver reached convergence criterion or no convergence required (eg. MP2 solver) # --- Energies e_corr: float = None # Fragment correlation energy contribution + e_corr_rpa: float = None # Fragment correlation energy contribution at the level of RPA. # --- Wave-function wf: WaveFunction = None # WaveFunction object (MP2, CCSD,...) pwf: WaveFunction = None # Fragment-projected wave function @@ -517,8 +522,8 @@ def diagonalize_cluster_dm(self, *mo_coeff, dm1=None, norm=2, tol=1e-4): sc = np.dot(self.base.get_ovlp(), c_cluster) dm = dot(sc.T, dm1, sc) e, r = np.linalg.eigh(dm) - if tol and not np.allclose(np.fmin(abs(e), abs(e-norm)), 0, atol=tol, rtol=0): - raise RuntimeError("Eigenvalues of cluster-DM not all close to 0 or %d:\n%s" % (norm, e)) + if tol and not np.allclose(np.fmin(abs(e), abs(e-norm)), 0, atol=2*tol, rtol=0): + self.log.warn("Eigenvalues of cluster-DM not all close to 0 or %d:\n%s" % (norm, e)) e, r = e[::-1], r[:,::-1] c_cluster = np.dot(c_cluster, r) c_cluster = fix_orbital_sign(c_cluster)[0] @@ -818,6 +823,11 @@ def get_bath(occtype): c_buffer = None return MP2_Bath(self, dmet_bath=dmet, occtype=occtype, c_buffer=c_buffer, project_dmet_order=project_dmet_order, project_dmet_mode=project_dmet_mode) + if btype == 'rpa': + project_dmet_order = self._get_bath_option('project_dmet_order', occtype) + project_dmet_mode = self._get_bath_option('project_dmet_mode', occtype) + return RPA_Bath(self, dmet_bath=dmet, occtype=occtype, c_buffer=None, + project_dmet_order=project_dmet_order, project_dmet_mode=project_dmet_mode) raise NotImplementedError('bathtype= %s' % btype) self._bath_factory_occ = get_bath(occtype='occupied') self._bath_factory_vir = get_bath(occtype='virtual') @@ -839,7 +849,7 @@ def get_orbitals(occtype): elif btype == 'ewdmet': order = self._get_bath_option('order', occtype) c_bath, c_frozen = factory.get_bath(order) - elif btype == 'mp2': + elif btype in ['mp2', 'rpa', 'rpacorr']: threshold = self._get_bath_option('threshold', occtype) truncation = self._get_bath_option('truncation', occtype) bno_threshold = BNO_Threshold(truncation, threshold) @@ -863,7 +873,7 @@ def get_orbitals(occtype): def check_occup(mo_coeff, expected): occup = self.get_mo_occupation(mo_coeff) # RHF - atol = self.opts.bath_options['dmet_threshold'] + atol = self.opts.bath_options['occupation_tolerance'] if np.ndim(occup[0]) == 0: assert np.allclose(occup, 2*expected, rtol=0, atol=2*atol) else: @@ -1036,10 +1046,25 @@ def get_solver(self, solver=None): cluster_solver.hamil.add_screening(self._seris_ov) return cluster_solver + def get_local_rpa_correction(self, hamil=None): + e_loc_rpa = None + if self.base.opts.ext_rpa_correction: + hamil = hamil or self.hamil + cumulant = self.base.opts.ext_rpa_correction == "cumulant" + if self.base.opts.ext_rpa_correction not in ["erpa", "cumulant"]: + raise ValueError("Unknown external rpa correction %s specified.") + e_loc_rpa = hamil.calc_loc_erpa(*self._seris_ov[1:], cumulant) + return e_loc_rpa + def get_frag_hamil(self): + if self.opts.screening is not None: + if "crpa_full" in self.opts.screening: + self.bos_freqs, self.couplings = get_frag_W(self.mf, self, pcoupling=("pcoupled" in self.opts.screening), + only_ov_screened=("ov" in self.opts.screening), + log=self.log) # This detects based on fragment what kind of Hamiltonian is appropriate (restricted and/or EB). return ClusterHamiltonian(self, self.mf, self.log, screening=self.opts.screening, - cache_eris=self.opts.store_eris) + cache_eris=self.opts.store_eris, match_fock=self.opts.match_cluster_fock) def get_solver_options(self, *args, **kwargs): raise AbstractMethodError diff --git a/vayesta/core/qemb/qemb.py b/vayesta/core/qemb/qemb.py index d28b89d13..472c4911f 100644 --- a/vayesta/core/qemb/qemb.py +++ b/vayesta/core/qemb/qemb.py @@ -16,16 +16,14 @@ import pyscf.pbc import pyscf.pbc.tools import pyscf.lib + from pyscf.mp.mp2 import _mo_without_core from vayesta.core.foldscf import FoldedSCF, fold_scf from vayesta.core.util import (OptionsBase, OrthonormalityError, SymmetryError, dot, einsum, energy_string, getattr_recursive, hstack, log_method, log_time, with_doc) -from vayesta.core import spinalg -from vayesta.core.ao2mo import kao2gmo_cderi -from vayesta.core.ao2mo import postscf_ao2mo -from vayesta.core.ao2mo import postscf_kao2gmo +from vayesta.core import spinalg, eris from vayesta.core.scmf import PDMET, Brueckner -from vayesta.core.qemb.scrcoulomb import build_screened_eris +from vayesta.core.screening.screening_moment import build_screened_eris from vayesta.mpi import mpi from vayesta.core.qemb.register import FragmentRegister from vayesta.rpa import ssRIRPA @@ -67,7 +65,7 @@ class Options(OptionsBase): # --- Bath options bath_options: dict = OptionsBase.dict_with_defaults( # General - bathtype='dmet', canonicalize=True, + bathtype='dmet', canonicalize=True, occupation_tolerance=1e-8, # DMET bath dmet_threshold=1e-8, # R2 bath @@ -105,7 +103,7 @@ class Options(OptionsBase): # EBFCI/EBCCSD max_boson_occ=2, # EBCC - ansatz=None, + ansatz=None, store_as_ccsd=None, fermion_wf=False, # Dump dumpfile='clusters.h5', # MP2 @@ -113,7 +111,9 @@ class Options(OptionsBase): # --- Other symmetry_tol: float = 1e-6 # Tolerance (in Bohr) for atomic positions symmetry_mf_tol: float = 1e-5 # Tolerance for mean-field solution - screening: Optional[str] = None + screening: Optional[str] = None # What form of screening to use in clusters. + ext_rpa_correction: Optional[str] = None + match_cluster_fock: bool = False class Embedding: """Base class for quantum embedding methods. @@ -518,6 +518,13 @@ def e_nuc(self): """Nuclear-repulsion energy per unit cell (not folded supercell).""" return self.mol.energy_nuc()/self.ncells + @property + def e_nonlocal(self): + if self.opts.ext_rpa_correction is None: + return 0.0 + e_local = sum([x.results.e_corr_rpa * x.symmetry_factor for x in self.get_fragments(sym_parent=None)]) + return self.e_rpa - (e_local / self.ncells) + # Embedding properties @property @@ -631,137 +638,53 @@ def get_ovlp_power(self, power): spow = pyscf.pbc.tools.k2gamma.to_supercell_ao_integrals(self.kcell, self.kpts, spowk) return spow - def get_cderi(self, mo_coeff, compact=False, blksize=None): - """Get density-fitted three-center integrals in MO basis.""" - if compact: - raise NotImplementedError() - if self.kdf is not None: - return kao2gmo_cderi(self.kdf, mo_coeff) - - if np.ndim(mo_coeff[0]) == 1: - mo_coeff = (mo_coeff, mo_coeff) - - nao = self.mol.nao - naux = (self.df.auxcell.nao if hasattr(self.df, 'auxcell') else self.df.auxmol.nao) - cderi = np.zeros((naux, mo_coeff[0].shape[-1], mo_coeff[1].shape[-1])) - cderi_neg = None - if blksize is None: - blksize = int(1e9 / naux*nao*nao * 8) - # PBC: - if hasattr(self.df, 'sr_loop'): - blk0 = 0 - for labr, labi, sign in self.df.sr_loop(compact=False, blksize=blksize): - assert np.allclose(labi, 0) - assert (cderi_neg is None) # There should be only one block with sign -1 - labr = labr.reshape(-1, nao, nao) - if (sign == 1): - blk1 = (blk0 + labr.shape[0]) - blk = np.s_[blk0:blk1] - blk0 = blk1 - cderi[blk] = einsum('Lab,ai,bj->Lij', labr, mo_coeff[0], mo_coeff[1]) - elif (sign == -1): - cderi_neg = einsum('Lab,ai,bj->Lij', labr, mo_coeff[0], mo_coeff[1]) - return cderi, cderi_neg - # No PBC: - blk0 = 0 - for lab in self.df.loop(blksize=blksize): - blk1 = (blk0 + lab.shape[0]) - blk = np.s_[blk0:blk1] - blk0 = blk1 - lab = pyscf.lib.unpack_tril(lab) - cderi[blk] = einsum('Lab,ai,bj->Lij', lab, mo_coeff[0], mo_coeff[1]) - return cderi, None - - @log_method() - def get_eris_array(self, mo_coeff, compact=False): - """Get electron-repulsion integrals in MO basis as a NumPy array. - - Parameters - ---------- - mo_coeff: [list(4) of] (n(AO), n(MO)) array - MO coefficients. - - Returns - ------- - eris: (n(MO), n(MO), n(MO), n(MO)) array - Electron-repulsion integrals in MO basis. - """ - # PBC with k-points: - if self.kdf is not None: - if compact: - raise NotImplementedError - if np.ndim(mo_coeff[0]) == 1: - mo_coeff = 4*[mo_coeff] - cderi1, cderi1_neg = kao2gmo_cderi(self.kdf, mo_coeff[:2]) - if (mo_coeff[0] is mo_coeff[2]) and (mo_coeff[1] is mo_coeff[3]): - cderi2, cderi2_neg = cderi1, cderi1_neg - else: - cderi2, cderi2_neg = kao2gmo_cderi(self.kdf, mo_coeff[2:]) - eris = einsum('Lij,Lkl->ijkl', cderi1.conj(), cderi2) - if cderi1_neg is not None: - eris -= einsum('Lij,Lkl->ijkl', cderi1_neg.conj(), cderi2_neg) - return eris - # Molecules and Gamma-point PBC: - if hasattr(self.mf, 'with_df') and self.mf.with_df is not None: - eris = self.mf.with_df.ao2mo(mo_coeff, compact=compact) - elif self.mf._eri is not None: - eris = pyscf.ao2mo.kernel(self.mf._eri, mo_coeff, compact=compact) - else: - eris = self.mol.ao2mo(mo_coeff, compact=compact) - if not compact: - if isinstance(mo_coeff, np.ndarray) and mo_coeff.ndim == 2: - shape = 4*[mo_coeff.shape[-1]] - else: - shape = [mo.shape[-1] for mo in mo_coeff] - eris = eris.reshape(shape) - return eris - - @log_method() - def get_eris_object(self, postscf, fock=None): - """Get ERIs for post-SCF methods. + get_cderi = eris.get_cderi - For folded PBC calculations, this folds the MO back into k-space - and contracts with the k-space three-center integrals.. + get_eris_array = eris.get_eris_array - Parameters - ---------- - postscf: one of the following PySCF methods: MP2, CCSD, RCCSD, DFCCSD - Post-SCF method with attribute mo_coeff set. - - Returns - ------- - eris: _ChemistsERIs - ERIs which can be used for the respective post-scf method. - """ - if fock is None: - if isinstance(postscf, pyscf.mp.mp2.MP2): - fock = self.get_fock() - elif isinstance(postscf, (pyscf.ci.cisd.CISD, pyscf.cc.ccsd.CCSD)): - fock = self.get_fock(with_exxdiv=False) - else: - raise ValueError("Unknown post-SCF method: %r", type(postscf)) - # For MO energies, always use get_fock(): - mo_act = _mo_without_core(postscf, postscf.mo_coeff) - mo_energy = einsum('ai,ab,bi->i', mo_act, self.get_fock(), mo_act) - e_hf = self.mf.e_tot - - # Fold MOs into k-point sampled primitive cell, to perform efficient AO->MO transformation: - if self.kdf is not None: - return postscf_kao2gmo(postscf, self.kdf, fock=fock, mo_energy=mo_energy, e_hf=e_hf) - # Regular AO->MO transformation - eris = postscf_ao2mo(postscf, fock=fock, mo_energy=mo_energy, e_hf=e_hf) - return eris + get_eris_object = eris.get_eris_object @log_method() @with_doc(build_screened_eris) def build_screened_eris(self, *args, **kwargs): - scrfrags = [x.opts.screening for x in self.fragments if x.opts.screening is not None] - if len(scrfrags) > 0: - # Calculate total dRPA energy in N^4 time; this is cheaper than screening calculations. - rpa = ssRIRPA(self.mf, log=self.log) - self.e_nonlocal, energy_error = rpa.kernel_energy(correction='linear') - if scrfrags.count("mrpa") > 0: - build_screened_eris(self, *args, **kwargs) + nmomscr = len([x.opts.screening for x in self.fragments if x.opts.screening == "mrpa"]) + lov = None + if self.opts.ext_rpa_correction: + cumulant = self.opts.ext_rpa_correction == "cumulant" + if nmomscr < self.nfrag: + raise NotImplementedError("External dRPA correction currently requires all fragments use mrpa screening.") + + if self.opts.ext_rpa_correction not in ["erpa", "cumulant"]: + raise ValueError("Unknown external rpa correction %s specified.") + lov = self.get_cderi((self.mo_coeff_occ, self.mo_coeff_vir)) + rpa = ssRIRPA(self.mf, log=self.log, lov=lov) + if cumulant: + lp = lov[0] + lp = lp.reshape((lp.shape[0], -1)) + l_ = l_2 = lp + + if lov[1] is not None: + ln = lov[1].reshape((lov[1].shape[0], -1)) + l_ = np.concatenate([lp, ln], axis=0) + l_2 = np.concatenate([lp, -ln], axis=0) + + l_ = np.concatenate([l_, l_], axis=1) + l_2 = np.concatenate([l_2, l_2], axis=1) + + m0 = rpa.kernel_moms(0, target_rot=l_)[0][0] + # Deduct effective mean-field contribution and project the RHS and we're done. + self.e_rpa = 0.5 * einsum("pq,pq->", m0 - l_, l_2) + 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.e_rpa = self.e_rpa / self.ncells + self.log.info("Set total RPA correlation energy contribution as %s", energy_string(self.e_rpa)) + if nmomscr > 0: + self.log.info("") + self.log.info("SCREENED INTERACTION SETUP") + self.log.info("==========================") + with log_time(self.log.timing, "Time for screened interation setup: %s"): + build_screened_eris(self, *args, store_m0=self.opts.ext_rpa_correction, cderi_ov=lov, **kwargs) # Symmetry between fragments # -------------------------- @@ -1561,6 +1484,7 @@ def _reset_fragments(self, *args, **kwargs): def _reset(self): self.e_corr = None self.converged = False + self.e_rpa = None def reset(self, *args, **kwargs): self._reset() diff --git a/vayesta/core/qemb/uqemb.py b/vayesta/core/qemb/uqemb.py index fba6e8244..46cce67f9 100644 --- a/vayesta/core/qemb/uqemb.py +++ b/vayesta/core/qemb/uqemb.py @@ -196,6 +196,13 @@ def get_eris_object(self, postscf, fock=None): eris = postscf_ao2mo(postscf, fock=fock, mo_energy=mo_energy, e_hf=e_hf) return eris + @log_method() + @with_doc(Embedding.build_screened_eris) + def build_screened_eris(self, *args, **kwargs): + if self.opts.ext_rpa_correction is not None: + raise NotImplementedError("External RPA correlation energy only implemented for restricted references!") + super().build_screened_eris(*args, **kwargs) + def update_mf(self, mo_coeff, mo_energy=None, veff=None): """Update underlying mean-field object.""" # Chech orthonormal MOs diff --git a/vayesta/core/screening/__init__.py b/vayesta/core/screening/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/vayesta/core/screening/screening_crpa.py b/vayesta/core/screening/screening_crpa.py new file mode 100644 index 000000000..0de9fd65d --- /dev/null +++ b/vayesta/core/screening/screening_crpa.py @@ -0,0 +1,223 @@ +import scipy.linalg + +from vayesta.rpa import ssRPA +from .screening_moment import _get_target_rot +import copy +from vayesta.core.util import * +import numpy as np +from pyscf import lib +from pyscf.lib import logger +from pyscf.ao2mo import _ao2mo +from pyscf import __config__ + + +class cRPAError(RuntimeError): + pass + +def get_frag_W(mf, fragment, pcoupling=True, only_ov_screened=False, log=None): + """Generates screened coulomb interaction due to screening at the level of cRPA. + Note that this currently scales as O(N_frag N^6), so is not practical without further refinement. + + Parameters + ---------- + mf : pyscf.scf object + Mean-field instance. + fragment : vayesta.qemb.Fragment subclass + Fragments for the calculation, used to define local interaction space. + log : logging.Logger, optional + Logger object. If None, the logger of the `emb` object is used. Default: None. + + Returns + ------- + freqs : np.array + Effective bosonic frequencies for cRPA screening. + couplings : np.array + Effective bosonic couplings for cRPA screening. + """ + + log.info("Generating screened interaction via frequency dependent cRPA.") + try: + l_a, l_b, crpa = set_up_W_crpa(mf, fragment, pcoupling, only_ov_screened, log=log) + except cRPAError as e: + freqs = np.zeros((0,)) + nmo = mf.mo_coeff.shape[1] + couplings = (np.zeros((0, nmo, nmo)), np.zeros((0, nmo, nmo))) + else: + freqs = crpa.freqs_ss + couplings = (l_a, l_b) + log.info("cRPA resulted in %d poles", len(freqs)) + + return freqs, couplings + + +def get_frag_deltaW(mf, fragment, pcoupling=True, only_ov_screened=False, log=None): + """Generates change in coulomb interaction due to screening at the level of static limit of cRPA. + Note that this currently scales as O(N_frag N^6), so is not practical without further refinement. + + Parameters + ---------- + mf : pyscf.scf object + Mean-field instance. + fragment : vayesta.qemb.Fragment subclass + Fragments for the calculation, used to define local interaction space. + log : logging.Logger, optional + Logger object. If None, the logger of the `emb` object is used. Default: None. + + Returns + ------- + deltaW : np.array + Change in cluster local coulomb interaction due to cRPA screening. + """ + + log.info("Generating screened interaction via static limit of cRPA.") + log.warning("This is poorly defined for non-CAS fragmentations.") + log.warning("This implementation is expensive, with O(N^6) computational cost per cluster.") + try: + l_a, l_b, crpa = set_up_W_crpa(mf, fragment, pcoupling, only_ov_screened=only_ov_screened, log=log) + except cRPAError as e: + nmo = mf.mo_coeff.shape[1] + delta_w = tuple([np.zeros([nmo] * 4)] * 4) + crpa = None + else: + # Have a factor of -2 due to negative value of RPA dd response, and summation of + # the excitation and deexcitation branches of the dd response. + static_fac = - 1.0 * (crpa.freqs_ss ** (-1)) + + delta_w = ( + einsum("npq,n,nrs->pqrs", l_a, static_fac, l_a) + einsum("nqp,n,nsr->pqrs", l_a, static_fac, l_a), + einsum("npq,n,nrs->pqrs", l_a, static_fac, l_b) + einsum("nqp,n,nsr->pqrs", l_a, static_fac, l_b), + einsum("npq,n,nrs->pqrs", l_b, static_fac, l_b) + einsum("nqp,n,nsr->pqrs", l_b, static_fac, l_b), + ) + return delta_w, crpa + +def set_up_W_crpa(mf, fragment, pcoupling=True, only_ov_screened=False, log=None): + is_rhf = np.ndim(mf.mo_coeff[1]) == 1 + if not hasattr(mf, "with_df"): + raise NotImplementedError("Screened interactions require density-fitting.") + crpa, rot_loc, rot_crpa = get_crpa(mf, fragment, log) + # Now need to calculate interactions. + nmo = mf.mo_coeff.shape[1] + nocc = sum(mf.mo_occ.T > 0) + c = mf.mo_coeff + if is_rhf: + nocc = (nocc, nocc) + c = (c, c) + + # First calculate alpha contribution. + lov_a = ao2mo(mf, mo_coeff=c[0], ijslice=(0, nocc[0], nocc[0], nmo)).reshape((-1, crpa.ova)) + lov_a = dot(lov_a, crpa.ov_rot[0].T) + l_aux = dot(lov_a, crpa.XpY_ss[0]) + del lov_a + lov_b = ao2mo(mf, mo_coeff=c[1], ijslice=(0, nocc[1], nocc[1], nmo)).reshape((-1, crpa.ovb)) + lov_b = dot(lov_b, crpa.ov_rot[1].T) + + # This is a decomposition of the cRPA reducible dd response in the auxilliary basis + l_aux += dot(lov_b, crpa.XpY_ss[1]) + del lov_b + + # This is expensive, and we'd usually want to avoid doing it twice unnecessarily, but other things will be worse. + # Now calculate the coupling back to the fragment itself. + c_act = fragment.cluster.c_active + if is_rhf: + c_act = (c_act, c_act) + + # Calculating this quantity scales as O(N^3), rather than O(N^4) if we explicitly calculated the cRPA dd response + # in the auxiliary basis. + lpqa_loc = ao2mo(mf, mo_coeff=c_act[0]) + l_a = einsum("npq,nm->mpq", lpqa_loc, l_aux) + del lpqa_loc + lpqb_loc = ao2mo(mf, mo_coeff=c_act[1]) + l_b = einsum("npq,nm->mpq", lpqb_loc, l_aux) + del lpqb_loc + + if pcoupling: + # Need to calculate additional contribution to the couplings resulting from rotation of the irreducible + # polarisability. + # Generate the full-system matrix of orbital energy differences. + eps = ssRPA(mf)._gen_eps() + # First, generate epsilon couplings between cluster and crpa spaces. + eps_fb = [einsum("p,qp,rp->qr", e, l, nl) for e, l, nl in zip(eps, rot_loc, crpa.ov_rot)] + # Then generate X and Y values for this correction. + x_crpa = [(p + m)/2 for p, m in zip(crpa.XpY_ss, crpa.XmY_ss)] + y_crpa = [(p - m)/2 for p, m in zip(crpa.XpY_ss, crpa.XmY_ss)] + # Contract with epsilon values + a_fb = [dot(e, x) for x, e in zip(x_crpa, eps_fb)] + b_fb = [dot(e, y) for y, e in zip(y_crpa, eps_fb)] + no = fragment.cluster.nocc_active + if isinstance(no, int): + no = (no, no) + nv = fragment.cluster.nvir_active + if isinstance(nv, int): + nv = (nv, nv) + l_a[:, :no[0], no[0]:] += a_fb[0].T.reshape((a_fb[0].shape[-1], no[0], nv[0])) + l_b[:, :no[1], no[1]:] += a_fb[1].T.reshape((a_fb[1].shape[-1], no[1], nv[1])) + + l_a[:, no[0]:, :no[0]] += b_fb[0].T.reshape((b_fb[0].shape[-1], no[0], nv[0])).transpose(0,2,1) + l_b[:, no[1]:, :no[1]] += b_fb[1].T.reshape((b_fb[1].shape[-1], no[1], nv[1])).transpose(0,2,1) + + if only_ov_screened: + # Zero out all contributions screening oo or vv contributions. + no = fragment.cluster.nocc_active + if isinstance(no, int): + no = (no, no) + l_a[:, no[0]:, no[0]:] = 0.0 + l_a[:, :no[0], :no[0]] = 0.0 + l_b[:, no[1]:, no[1]:] = 0.0 + l_b[:, :no[1], :no[1]] = 0.0 + return l_a, l_b, crpa + + +def get_crpa(orig_mf, f, log): + + def construct_loc_rot(f): + """Constructs the rotation of the overall mean-field space into which """ + ro = f.get_overlap("cluster[occ]|mo[occ]") + rv = f.get_overlap("cluster[vir]|mo[vir]") + + if isinstance(ro, np.ndarray): + ro = (ro, ro) + if isinstance(rv, np.ndarray): + rv = (rv, rv) + + rot_ova = einsum("Ij,Ab->IAjb", ro[0], rv[0]) + rot_ova = rot_ova.reshape((rot_ova.shape[0] * rot_ova.shape[1], -1)) + + rot_ovb = einsum("Ij,Ab->IAjb", ro[1], rv[1]) + rot_ovb = rot_ovb.reshape((rot_ovb.shape[0] * rot_ovb.shape[1], -1)) + return rot_ova, rot_ovb + + rot_loc = construct_loc_rot(f) + rot_ov = scipy.linalg.null_space(rot_loc[0]).T, scipy.linalg.null_space(rot_loc[1]).T + if rot_ov[0].shape[0] == 0 and rot_ov[1].shape[0] == 0: + log.warning("cRPA space contains no excitations! Interactions will be unscreened.") + raise cRPAError("cRPA space contains no excitations!") + # RPA calculation and new_mf will contain all required information for the response. + crpa = ssRPA(orig_mf, ov_rot=rot_ov) + crpa.kernel() + return crpa, rot_loc, rot_ov + + +def ao2mo(mf, mo_coeff=None, ijslice=None): + """Get MO basis density-fitted integrals. + """ + + if mo_coeff is None: + mo_coeff = mf.mo_coeff + nmo = mo_coeff.shape[1] + naux = mf.with_df.get_naoaux() + mem_incore = (2 * nmo ** 2 * naux) * 8 / 1e6 + mem_now = lib.current_memory()[0] + + mo = np.asarray(mo_coeff, order='F') + if ijslice is None: + ijslice = (0, nmo, 0, nmo) + + finshape = (naux, ijslice[1] - ijslice[0], ijslice[3] - ijslice[2]) + + Lpq = None + if (mem_incore + mem_now < 0.99 * mf.max_memory) or mf.mol.incore_anyway: + Lpq = _ao2mo.nr_e2(mf.with_df._cderi, mo, ijslice, aosym='s2', out=Lpq) + return Lpq.reshape(*finshape) + else: + logger.warn(mf, 'Memory may not be enough!') + raise NotImplementedError diff --git a/vayesta/core/qemb/scrcoulomb.py b/vayesta/core/screening/screening_moment.py similarity index 78% rename from vayesta/core/qemb/scrcoulomb.py rename to vayesta/core/screening/screening_moment.py index dacbefbc3..7d6ce7a61 100644 --- a/vayesta/core/qemb/scrcoulomb.py +++ b/vayesta/core/screening/screening_moment.py @@ -2,12 +2,12 @@ import numpy as np import scipy import scipy.linalg -from vayesta.rpa import ssRIRPA +from vayesta.rpa import ssRIRPA, ssRPA from vayesta.core.util import dot, einsum from vayesta.mpi import mpi -def build_screened_eris(emb, fragments=None, cderi_ov=None, calc_e=True, npoints=48, log=None): +def build_screened_eris(emb, fragments=None, cderi_ov=None, store_m0=True, npoints=48, log=None): """Generates renormalised coulomb interactions for use in local cluster calculations. Currently requires unrestricted system. @@ -21,8 +21,8 @@ def build_screened_eris(emb, fragments=None, cderi_ov=None, calc_e=True, npoints cderi_ov : np.array or tuple of np.array, optional. Cholesky-decomposed ERIs in the particle-hole basis of mf. If mf is unrestricted this should be a list of arrays corresponding to the different spin channels. - calc_ecorrection : bool, optional. - Whether to calculate a nonlocal energy correction at the level of RPA + store_m0 : bool, optional. + Whether to store the local zeroth moment in the fragment class for use later. npoints : int, optional Number of points for numerical integration. Default: 48. log : logging.Logger, optional @@ -45,51 +45,35 @@ def build_screened_eris(emb, fragments=None, cderi_ov=None, calc_e=True, npoints if fragments is None: fragments = emb.get_fragments(active=True, sym_parent=None, mpi_rank=mpi.rank) fragments = [f for f in fragments if f.opts.screening == 'mrpa'] - if emb.spinsym != 'unrestricted': - raise NotImplementedError("Screened interactions require a spin-unrestricted formalism.") if emb.df is None: raise NotImplementedError("Screened interactions require density-fitting.") r_occs = [f.get_overlap('mo[occ]|cluster[occ]') for f in fragments] r_virs = [f.get_overlap('mo[vir]|cluster[vir]') for f in fragments] target_rots, ovs_active = _get_target_rot(r_occs, r_virs) - rpa = ssRIRPA(emb.mf, log=log, Lpq=cderi_ov) - if calc_e: - # This scales as O(N^4) - erpa, energy_error = rpa.kernel_energy(correction='linear') - else: - erpa = None - - tr = np.concatenate(target_rots, axis=0) - if sum(sum(ovs_active)) > 0: - # Computation scales as O(N^4) - moms_interact, est_errors = rpa.kernel_moms(0, tr, npoints=npoints) - momzero_interact = moms_interact[0] - else: - momzero_interact = np.zeros_like(np.concatenate(tr, axis=0)) - - # Now need to separate into local contributions - n = 0 - local_moments = [] - for nov, rot in zip(ovs_active, target_rots): - # Computation costs O(N^2 N_clus^2) - # Get corresponding section of overall moment, then project to just local contribution. - local_moments += [dot(momzero_interact[n:n+sum(nov)], rot.T)] - n += sum(nov) + local_moments = calc_moms_RIRPA(emb.mf, target_rots, ovs_active, log, cderi_ov, npoints) + # Could generate moments using N^6 moments instead, but just for debugging. + #local_moments, erpa = calc_moms_RPA(emb.mf, target_rots, ovs_active, log, cderi_ov, calc_e, npoints) # Then construct the RPA coupling matrix A-B, given by the diagonal matrix of energy differences. no = np.array(sum(emb.mf.mo_occ.T > 0)) - norb = np.array(emb.mo_coeff).shape[1] + if no.size == 1: no = np.array([int(no), int(no)]) + norb = emb.mo_coeff[0].shape[0] nv = norb - no + mo_e = emb.mf.mo_energy + + if isinstance(mo_e[0], float): + mo_e = np.array([mo_e, mo_e]) + def get_eps_singlespin(no_, nv_, mo_energy): eps = np.zeros((no_, nv_)) eps = eps + mo_energy[no_:] eps = (eps.T - mo_energy[:no_]).T eps = eps.reshape(-1) return eps - eps = np.concatenate([get_eps_singlespin(no[0], nv[0], emb.mf.mo_energy[0]), - get_eps_singlespin(no[1], nv[1], emb.mf.mo_energy[1])]) + eps = np.concatenate([get_eps_singlespin(no[0], nv[0], mo_e[0]), + get_eps_singlespin(no[1], nv[1], mo_e[1])]) # And use this to perform inversion to calculate interaction in cluster. seris_ov = [] @@ -97,7 +81,11 @@ def get_eps_singlespin(no_, nv_, mo_energy): amb = einsum("pn,qn,n->pq", rot, rot, eps) # O(N^2 N_clus^4) # Everything from here on is independent of system size, scaling at most as O(N_clus^6) # (arrays have side length equal to number of cluster single-particle excitations). - mominv = np.linalg.inv(mom) + e, c = np.linalg.eigh(mom) + if min(e) < 1e-4: + log.warning("Small eigenvalue of local rpa moment in %s: %e", f.name, min(e)) + + mominv = einsum("pn,n,qn->pq", c, e**(-1), c) apb = dot(mominv, amb, mominv) # This is the renormalised coulomb kernel in the cluster. @@ -108,23 +96,57 @@ def get_eps_singlespin(no_, nv_, mo_energy): # indices. no = f.cluster.nocc_active nv = f.cluster.nvir_active + if isinstance(no, int): + no = (no, no) + nv = (nv, nv) + kcaa = kc[:ova, :ova].reshape((no[0], nv[0], no[0], nv[0])) kcab = kc[:ova, ova:].reshape((no[0], nv[0], no[1], nv[1])) kcbb = kc[ova:, ova:].reshape((no[1], nv[1], no[1], nv[1])) - if kcaa.shape == kcbb.shape: - # This is not that meaningful, since the alpha and beta basis themselves may be different: - log.info("Screened interations in %s: spin-symmetry= %.3e spin-dependence= %.3e", - f.name, abs(kcaa-kcbb).max(), abs((kcaa+kcbb)/2-kcab).max()) kc = (kcaa, kcab, kcbb) - f._seris_ov = kc + f._seris_ov = (kc, mom, amb) if store_m0 else (kc,) seris_ov.append(kc) - return seris_ov, erpa + +def calc_moms_RIRPA(mf, target_rots, ovs_active, log, cderi_ov, npoints): + rpa = ssRIRPA(mf, log=log, lov=cderi_ov) + tr = np.concatenate(target_rots, axis=0) + if sum(sum(ovs_active)) > 0: + # Computation scales as O(N^4) + moms_interact, est_errors = rpa.kernel_moms(0, tr, npoints=npoints) + momzero_interact = moms_interact[0] + else: + momzero_interact = np.zeros_like(np.concatenate(tr, axis=0)) + + # Now need to separate into local contributions + n = 0 + local_moments = [] + for nov, rot in zip(ovs_active, target_rots): + # Computation costs O(N^2 N_clus^2) + # Get corresponding section of overall moment, then project to just local contribution. + mom = dot(momzero_interact[n:n+sum(nov)], rot.T) + # This isn't exactly symmetric due to numerical integration, so enforce here. + mom = (mom + mom.T) / 2 + local_moments += [mom] + n += sum(nov) + + return local_moments + +def calc_moms_RPA(mf, target_rots, ovs_active, log, cderi_ov, calc_e, npoints): + rpa = ssRPA(mf, log=log) + erpa = rpa.kernel() + mom0 = rpa.gen_moms(0)[0] + local_moments = [] + for rot in target_rots: + local_moments += [dot(rot, mom0, rot.T)] + return local_moments, erpa def get_screened_eris_full(eris, seris_ov, copy=True, log=None): """Build full array of screened ERIs, given the bare ERIs and screening.""" + log.info("Generating screened interaction to conserve zeroth moment of the dd response.") + def replace_ov(full, ov, spins): out = full.copy() if copy else full no1, no2 = ov.shape[0], ov.shape[2] @@ -134,12 +156,11 @@ def replace_ov(full, ov, spins): out[v1,o1,o2,v2] = ov.transpose([1, 0, 2, 3]) out[o1,v1,v2,o2] = ov.transpose([0, 1, 3, 2]) out[v1,o1,v2,o2] = ov.transpose([1, 0, 3, 2]) - if log: - maxidx = np.unravel_index(np.argmax(abs(full-out)), full.shape) - log.info("Maximally screened element of W(%2s|%2s): V= %.3e -> W= %.3e (delta= %.3e)", - 2*spins[0], 2*spins[1], full[maxidx], out[maxidx], out[maxidx]-full[maxidx]) return out + if isinstance(eris, np.ndarray): + eris = (eris, eris, eris) + seris = (replace_ov(eris[0], seris_ov[0], 'aa'), replace_ov(eris[1], seris_ov[1], 'ab'), replace_ov(eris[2], seris_ov[2], 'bb')) @@ -226,6 +247,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/core/types/wf/ccsd.py b/vayesta/core/types/wf/ccsd.py index 45516f5bd..a4ac69082 100644 --- a/vayesta/core/types/wf/ccsd.py +++ b/vayesta/core/types/wf/ccsd.py @@ -5,12 +5,13 @@ # import pyscf # import pyscf.cc from vayesta.core import spinalg -from vayesta.core.util import NotCalculatedError, Object, callif, einsum +from vayesta.core.util import NotCalculatedError, Object, callif, einsum, dot from vayesta.core.types import wf as wf_types from vayesta.core.types.orbitals import SpatialOrbitals from vayesta.core.types.wf.project import (project_c1, project_c2, project_uc1, project_uc2, symmetrize_c2, - symmetrize_uc2) + symmetrize_uc2, transform_c1, transform_c2, transform_uc1, transform_uc2) from vayesta.core.helper import pack_arrays, unpack_arrays +from scipy.linalg import block_diag def CCSD_WaveFunction(mo, t1, t2, **kwargs): @@ -170,6 +171,42 @@ def as_ccsd(self): def as_fci(self): raise NotImplementedError + def rotate(self, t, inplace=False): + """Rotate wavefunction representation to another basis. + Only rotations which don't mix occupied and virtual orbitals are supported. + Assumes rotated orbitals have same occupancy ordering as originals. + """ + o = self.mo.occ > 0 + v = self.mo.occ == 0 + + to = t[np.ix_(o, o)] + tv = t[np.ix_(v, v)] + tov = t[np.ix_(o, v)] + tvo = t[np.ix_(v, o)] + if abs(tov).max() > 1e-12 or abs(tvo).max() > 1e-12: + raise ValueError("Provided rotation mixes occupied and virtual orbitals.") + return self.rotate_ov(to, tv, inplace=inplace) + + def rotate_ov(self, to, tv, inplace=False): + """Rotate wavefunction representation to another basis. + Only rotations which don't mix occupied and virtual orbitals are supported. + + Parameters + ---------- + to : new occupied orbital coefficients in terms of current ones. + tv : new virtual orbital coefficients in terms of current ones. + inplace : Whether to transform in-place or return a new object. + """ + wf = self if inplace else self.copy() + wf.mo.basis_transform(lambda c: dot(c, block_diag(to, tv)), inplace=True) + wf.t1 = transform_c1(wf.t1, to, tv) + wf.t2 = transform_c2(wf.t2, to, tv) + if wf.l1 is not None: + wf.l1 = transform_c1(wf.l1, to, tv) + if wf.l2 is not None: + wf.l2 = transform_c2(wf.l2, to, tv) + return wf + class UCCSD_WaveFunction(RCCSD_WaveFunction): @@ -296,6 +333,56 @@ def multiply(self, factor): if self.l2 is not None: self.l2 = spinalg.multiply(self.l2, len(self.l2)*[factor]) + def rotate(self, t, inplace=False): + """Rotate wavefunction representation to another basis. + Only rotations which don't mix occupied and virtual orbitals are supported. + Assumes rotated orbitals have same occupancy ordering as originals. + """ + # Allow support for same rotation for alpha and beta. + if isinstance(t, np.ndarray) and t.ndim == 2: + t = (t, t) + def get_components(tsp, occ): + o = occ > 0 + v = occ == 0 + tspo = tsp[np.ix_(o, o)] + tspv = tsp[np.ix_(v, v)] + tspov = tsp[np.ix_(o, v)] + tspvo = tsp[np.ix_(v, o)] + if abs(tspov).max() > 1e-12 or abs(tspvo).max() > 1e-12: + raise ValueError("Provided rotation mixes occupied and virtual orbitals.") + return tspo, tspv + + toa, tva = get_components(t[0], self.mo.alpha.occ) + tob, tvb = get_components(t[1], self.mo.beta.occ) + + return self.rotate_ov((toa, tob), (tva, tvb), inplace=inplace) + + def rotate_ov(self, to, tv, inplace=False): + """Rotate wavefunction representation to another basis. + Only rotations which don't mix occupied and virtual orbitals are supported. + + Parameters + ---------- + to : new occupied orbital coefficients in terms of current ones. + tv : new virtual orbital coefficients in terms of current ones. + inplace : Whether to transform in-place or return a new object. + """ + wf = self if inplace else self.copy() + if isinstance(to, np.ndarray) and len(to) == 2: + assert(isinstance(tv, np.ndarray) and len(tv) == 2) + trafo = lambda c: dot(c, block_diag(to, tv)) + else: + trafo = [lambda c: dot(c, x) for x in (block_diag(to[0], tv[0]), block_diag(to[1], tv[1]))] + + wf.mo.basis_transform(trafo, inplace=True) + wf.t1 = transform_uc1(wf.t1, to, tv) + wf.t2 = transform_uc2(wf.t2, to, tv) + if wf.l1 is not None: + wf.l1 = transform_uc1(wf.l1, to, tv) + if wf.l2 is not None: + wf.l2 = transform_uc2(wf.l2, to, tv) + return wf + #def pack(self, dtype=float): # """Pack into a single array of data type `dtype`. diff --git a/vayesta/core/types/wf/project.py b/vayesta/core/types/wf/project.py index c1137f948..d697e25f9 100644 --- a/vayesta/core/types/wf/project.py +++ b/vayesta/core/types/wf/project.py @@ -1,6 +1,7 @@ """Utility functions for projection of wave functions.""" import numpy as np +from vayesta.core.util import einsum, dot def project_c1(c1, p): @@ -45,3 +46,27 @@ def symmetrize_uc2(c2, inplace=True): if len(c2) == 4: c2ab = (c2ab + c2[2].transpose(1,0,3,2))/2 return (c2aa, c2ab, c2bb) + +def transform_c1(c1, to, tv): + if c1 is None: return None + return dot(to.T, c1, tv) + +def transform_c2(c2, to, tv, to2=None, tv2=None): + if c2 is None: return None + if to2 is None: to2 = to + if tv2 is None: tv2 = tv + # Use einsum for now- tensordot would be faster but less readable. + return einsum("ijab,iI,jJ,aA,bB->IJAB", c2, to, to2, tv, tv2) + +def transform_uc1(c1, to, tv): + if c1 is None: return None + return (transform_c1(c1[0], to[0], tv[0]), + transform_c1(c1[1], to[1], tv[1])) + +def transform_uc2(c2, to, tv): + if c2 is None: return None + c2ba = (c2[2] if len(c2) == 4 else c2[1].transpose(1,0,3,2)) + return (transform_c2(c2[0], to[0], tv[0]), + transform_c2(c2[1], to[0], tv[0], to[1], tv[1]), + transform_c2(c2ba, to[1], tv[1], to[0], tv[0]), + transform_c2(c2[-1], to[1], tv[1])) diff --git a/vayesta/core/types/wf/wf.py b/vayesta/core/types/wf/wf.py index c9776b1b5..81b1bab79 100644 --- a/vayesta/core/types/wf/wf.py +++ b/vayesta/core/types/wf/wf.py @@ -65,6 +65,12 @@ def make_rdm1(self, *args, **kwargs): def make_rdm2(self, *args, **kwargs): raise AbstractMethodError + def rotate_ov(self, *args, **kwargs): + raise AbstractMethodError + + def rotate(self, *args, **kwargs): + raise AbstractMethodError + @staticmethod def from_pyscf(obj, **kwargs): # HF diff --git a/vayesta/dmet/dmet.py b/vayesta/dmet/dmet.py index 911268450..131b6abd0 100644 --- a/vayesta/dmet/dmet.py +++ b/vayesta/dmet/dmet.py @@ -128,11 +128,10 @@ def kernel(self): if type(nelec_mf) == tuple: nelec_mf = sum(nelec_mf) - if self.opts.screening == 'mrpa': - for f in self.get_fragments(sym_parent=None): - f.make_bath() - f.make_cluster() - self.build_screened_eris() + for f in self.get_fragments(sym_parent=None): + f.make_bath() + f.make_cluster() + self.build_screened_eris() def electron_err(cpt, construct_bath=False): err = self.calc_electron_number_defect(cpt, nelec_mf, sym_parents, nsym, construct_bath) diff --git a/vayesta/dmet/fragment.py b/vayesta/dmet/fragment.py index ea57d83b9..cff53c87b 100644 --- a/vayesta/dmet/fragment.py +++ b/vayesta/dmet/fragment.py @@ -67,6 +67,7 @@ def kernel(self, solver=None, init_guess=None, seris_ov=None, construct_bath=Tru return None cluster_solver = self.get_solver(solver) + # Chemical potential if chempot is not None: px = self.get_fragment_projector(self.cluster.c_active) diff --git a/vayesta/edmet/edmet.py b/vayesta/edmet/edmet.py index 5e2dda8e9..6e73338ba 100644 --- a/vayesta/edmet/edmet.py +++ b/vayesta/edmet/edmet.py @@ -54,6 +54,14 @@ def eps(self): eps = eps.reshape(-1) return eps, eps + @property + def e_nonlocal(self): + return self._e_nonlocal or 0.0 + + @e_nonlocal.setter + def e_nonlocal(self, value): + self._e_nonlocal = value + def check_solver(self, solver): is_uhf = np.ndim(self.mo_coeff[1]) == 2 is_eb = True @@ -375,4 +383,8 @@ def run_exact_full_ac(self, xc_kernel=None, deg=5, calc_local=False, cluster_con else: raise NotImplementedError + def _reset(self): + super()._reset() + self._e_nonlocal = None + REDMET = EDMET diff --git a/vayesta/ewf/ewf.py b/vayesta/ewf/ewf.py index cf81c3c05..e2c0f687f 100644 --- a/vayesta/ewf/ewf.py +++ b/vayesta/ewf/ewf.py @@ -27,6 +27,8 @@ class Options(Embedding.Options): # --- Bath settings bath_options: dict = Embedding.Options.change_dict_defaults('bath_options', bathtype='mp2', threshold=1e-8) + solver_options: dict = Embedding.Options.change_dict_defaults('solver_options', + fermion_wf=True) #ewdmet_max_order: int = 1 # If multiple bno thresholds are to be calculated, we can project integrals and amplitudes from a previous larger cluster: project_eris: bool = False # Project ERIs from a pervious larger cluster (corresponding to larger eta), can result in a loss of accuracy especially for large basis sets! @@ -139,12 +141,7 @@ def kernel(self): self.communicate_clusters() # --- Screened Coulomb interaction - if any(x.opts.screening is not None for x in fragments): - self.log.info("") - self.log.info("SCREENING INTERACTIONS") - self.log.info("======================") - with log_time(self.log.timing, "Time for screened interations: %s"): - self.build_screened_eris() + self.build_screened_eris() # --- Loop over fragments with no symmetry parent and with own MPI rank self.log.info("") @@ -232,7 +229,7 @@ def t1_diagnostic(self, warntol=0.02): # Defaults def make_rdm1(self, *args, **kwargs): - if self.solver.lower() == 'ccsd': + if "cc" in self.solver.lower(): return self._make_rdm1_ccsd_global_wf(*args, **kwargs) if self.solver.lower() == 'mp2': return self._make_rdm1_mp2_global_wf(*args, **kwargs) @@ -294,6 +291,7 @@ def _make_rdm2_ccsd_proj_lambda(self, *args, **kwargs): def get_e_corr(self, functional=None, **kwargs): functional = (functional or self.opts.energy_functional) + if functional == 'projected': self.log.warning("functional='projected' is deprecated; use functional='wf' instead.") functional = 'wf' diff --git a/vayesta/ewf/fragment.py b/vayesta/ewf/fragment.py index e9aacfffc..a0cf8d123 100644 --- a/vayesta/ewf/fragment.py +++ b/vayesta/ewf/fragment.py @@ -100,8 +100,13 @@ def __init__(self, *args, **kwargs): # For self-consistent mode self.solver_results = None + def _reset(self, *args, **kwargs): + super()._reset(*args, **kwargs) + # Need to unset these so can be regenerated each iteration. + self.opts.c_cas_occ = self.opts.c_cas_vir = None + def set_cas(self, iaos=None, c_occ=None, c_vir=None, minao='auto', dmet_threshold=None): - """Set complete active space for tailored CCSD""" + """Set complete active space for tailored CCSD and active-space CC methods.""" if dmet_threshold is None: dmet_threshold = 2*self.opts.bath_options['dmet_threshold'] if iaos is not None: @@ -114,7 +119,8 @@ def set_cas(self, iaos=None, c_occ=None, c_vir=None, minao='auto', dmet_threshol c_env = frag.get_env_coeff(indices) bath = DMET_Bath(self, dmet_threshold=dmet_threshold) c_dmet = bath.make_dmet_bath(c_env)[0] - c_iao_occ, c_iao_vir = self.diagonalize_cluster_dm(c_iao, c_dmet, tol=2*dmet_threshold) + tol = self.opts.bath_options['occupation_tolerance'] + c_iao_occ, c_iao_vir = self.diagonalize_cluster_dm(c_iao, c_dmet, tol=2*tol) else: c_iao_occ = c_iao_vir = None @@ -190,7 +196,8 @@ def kernel(self, solver=None, init_guess=None): # Create solver object cluster_solver = self.get_solver(solver) - + # Calculate cluster energy at the level of RPA. + e_corr_rpa = self.get_local_rpa_correction(cluster_solver.hamil) # --- Chemical potential cpt_frag = self.base.opts.global_frag_chempot if self.opts.nelectron_target is not None: @@ -253,7 +260,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 @@ -280,8 +287,9 @@ def get_solver_options(self, solver): self.log.debugv("Passing fragment option %s to solver.", attr) solver_opts[attr] = getattr(self.opts, attr) - if solver.upper() == 'TCCSD': - solver_opts['tcc'] = True + has_actspace = ((solver == "TCCSD") or ("CCSDt'" in solver) or + (solver.upper() == "EBCC" and self.opts.solver_options['ansatz'] == "CCSDt'")) + if has_actspace: # Set CAS orbitals if self.opts.c_cas_occ is None: self.log.warning("Occupied CAS orbitals not set. Setting to occupied DMET cluster orbitals.") @@ -291,7 +299,8 @@ def get_solver_options(self, solver): self.opts.c_cas_vir = self._dmet_bath.c_cluster_vir solver_opts['c_cas_occ'] = self.opts.c_cas_occ solver_opts['c_cas_vir'] = self.opts.c_cas_vir - solver_opts['tcc_fci_opts'] = self.opts.tcc_fci_opts + if solver == "TCCSD": + solver_opts['tcc_fci_opts'] = self.opts.tcc_fci_opts elif solver.upper() == 'DUMP': solver_opts['filename'] = self.opts.solver_options['dumpfile'] solver_opts['external_corrections'] = self.flags.external_corrections diff --git a/vayesta/ewf/ufragment.py b/vayesta/ewf/ufragment.py index 4006d8153..6d4fb381d 100644 --- a/vayesta/ewf/ufragment.py +++ b/vayesta/ewf/ufragment.py @@ -11,8 +11,14 @@ class Fragment(RFragment, BaseFragment): - def set_cas(self, *args, **kwargs): - raise NotImplementedError() + def set_cas(self, iaos=None, c_occ=None, c_vir=None, minao='auto', dmet_threshold=None): + """Set complete active space for tailored CCSD and active-space CC methods.""" + if iaos is not None: + raise NotImplementedError("Unrestricted IAO-based CAS not implemented yet.") + + self.opts.c_cas_occ = c_occ + self.opts.c_cas_vir = c_vir + return c_occ, c_vir def get_fragment_energy(self, c1, c2, hamil=None, fock=None, axis1='fragment', c2ba_order='ba'): """Calculate fragment correlation energy contribution from projected C1, C2. diff --git a/vayesta/rpa/rirpa/NI_eval.py b/vayesta/rpa/rirpa/NI_eval.py index a88b545b1..96ece7d26 100644 --- a/vayesta/rpa/rirpa/NI_eval.py +++ b/vayesta/rpa/rirpa/NI_eval.py @@ -181,9 +181,8 @@ def find_good_start(ainit=1e-6, scale_fac=10.0, maxval=1e8, relevance_factor=5): optval = fvals[optarg] # Now find the values which are within reach of lowest value relevant = np.where(fvals < relevance_factor * optval)[0] - minarg = min(relevant[0], optarg - 1) - maxarg = min(relevant[-1], optarg + 1) + maxarg = max(relevant[-1], optarg + 1) return [ainit * scale_fac**x for x in (optarg, minarg, maxarg)] solve = 1 @@ -216,6 +215,11 @@ def find_good_start(ainit=1e-6, scale_fac=10.0, maxval=1e8, relevance_factor=5): solve, res.fun, ) + else: + self.log.info( + "Used newton optimisation to find optimal quadrature grid: a= %.2e", + solve, + ) return solve def fix_params(self): diff --git a/vayesta/rpa/rirpa/RIRPA.py b/vayesta/rpa/rirpa/RIRPA.py index 167c49f1a..622a51866 100644 --- a/vayesta/rpa/rirpa/RIRPA.py +++ b/vayesta/rpa/rirpa/RIRPA.py @@ -5,9 +5,12 @@ import pyscf.lib from vayesta.core.util import dot, einsum, time_string, timer from vayesta.rpa.rirpa import momzero_NI, energy_NI +from vayesta.core.eris import get_cderi +memory_string = lambda: "Memory usage: %.2f GB" % (pyscf.lib.current_memory()[0] / 1e3) -class ssRIRPA: + +class ssRIRRPA: """Approach based on equations expressed succinctly in the appendix of Furche, F. (2001). PRB, 64(19), 195120. https://doi.org/10.1103/PhysRevB.64.195120 WARNING: Should only be used with canonical mean-field orbital coefficients in mf.mo_coeff and RHF. @@ -27,8 +30,8 @@ class ssRIRPA: svd_tol : float, optional Threshold defining negligible singular values when compressing various decompositions. Default value is 1e-12. - Lpq : np.ndarray, optional - CDERIs in mo basis of provided mean field. + lov : np.ndarray or tuple of np.ndarray, optional + occupied-virtual CDERIs in mo basis of provided mean field. If None recalculated from AOs. Default value is None. compress : int, optional How thoroughly to attempt compression of the low-rank representations of various matrices. @@ -48,7 +51,7 @@ def __init__( log=None, err_tol=1e-6, svd_tol=1e-12, - Lpq=None, + lov=None, compress=0, ): self.mf = dfmf @@ -57,10 +60,27 @@ def __init__( self.err_tol = err_tol self.svd_tol = svd_tol self.e_corr_ss = None - self.Lpq = Lpq + self.lov = lov # Determine how many times to attempt compression of low-rank expressions for various matrices. self.compress = compress + @property + def mol(self): + return self.mf.mol + + @property + def df(self): + return self.mf.with_df + + @property + def kdf(self): + if hasattr(self.mf, "kmf"): + return self.mf.kmf.with_df + if hasattr(self.mf, "kpts"): + if self.mf.kpts is not None: + return self.mf.with_df + return None + @property def nocc(self): return sum(self.mf.mo_occ > 0) @@ -120,11 +140,17 @@ def e_tot(self): return self.mf.e_tot + self.e_corr @property - def D(self): + def eps(self): eps = np.zeros((self.nocc, self.nvir)) eps = eps + self.mo_energy_vir eps = (eps.T - self.mo_energy_occ).T eps = eps.reshape((self.ov,)) + return eps + + + @property + def D(self): + eps = self.eps D = np.concatenate([eps, eps]) return D @@ -203,6 +229,7 @@ def kernel_moms(self, max_moment, target_rot=None, **kwargs): moments[i] = einsum("pq,q->pq", moments[i - 2], D**2) + dot( moments[i - 2], ri_mp[1].T, ri_mp[0] ) + self.record_memory() if max_moment > 0: self.log.info( "RIRPA Higher Moments wall time: %s", @@ -267,12 +294,10 @@ def _kernel_mom0( # Evaluate (MP)^{1/2} - D, niworker = momzero_NI.MomzeroDeductD(*inputs) integral_offset = einsum("lp,p->lp", target_rot, self.D) - moment_offset = np.zeros_like(target_rot) elif integral_deduct is None: # Explicitly evaluate (MP)^{1/2}, with no offsets. niworker = momzero_NI.MomzeroDeductNone(*inputs) integral_offset = np.zeros_like(target_rot) - moment_offset = np.zeros_like(target_rot) elif integral_deduct == "HO": niworker = momzero_NI.MomzeroDeductHigherOrder(*inputs) offset_niworker = momzero_NI.MomzeroOffsetCalcGaussLag(*inputs) @@ -284,26 +309,22 @@ def _kernel_mom0( # estval2 = einsum("rp,pq,np,nq->rq", target_rot, mat ** (-1), ri_mp[0], ri_mp[1]) # self.log.info("Error in numerical Offset Approximation=%6.4e",abs(estval - estval2).max()) integral_offset = einsum("lp,p->lp", target_rot, self.D) + estval - moment_offset = np.zeros_like(target_rot) else: raise ValueError("Unknown integral offset specification.`") if return_niworker: return niworker, offset_niworker - # niworker.test_diag_derivs(4.0) if adaptive_quad: # Can also make use of scipy adaptive quadrature routines; this is more expensive but a good sense-check. integral, upper_bound = 2 * niworker.kernel_adaptive() else: integral, upper_bound = niworker.kernel(a=ainit, opt_quad=opt_quad) - # Need to construct RI representation of P^{-1} + ri_apb_inv = construct_inverse_RI(self.D, ri_apb) - if self.compress > 5: - ri_apb_inv = self.compress_low_rank(*ri_apb_inv, name="(A+B)^-1") - mom0 = einsum("pq,q->pq", integral + integral_offset, self.D ** (-1)) - np.dot( - np.dot(integral + integral_offset, ri_apb_inv[0].T), ri_apb_inv[1] - ) + + mom0 = self.mult_apbinv(integral + integral_offset, ri_apb_inv) + # Also need to convert error estimate of the integral into one for the actual evaluated quantity. # Use Cauchy-Schwartz to both obtain an upper bound on resulting mom0 error, and efficiently obtain upper bound # on norm of low-rank portion of P^{-1}. @@ -317,7 +338,7 @@ def _kernel_mom0( mom0_ub = None mom_lb = ( - self.test_eta0_error(mom0 + moment_offset, target_rot, ri_apb, ri_amb) + self.test_eta0_error(mom0, target_rot, ri_apb, ri_amb) if analytic_lower_bound else None ) @@ -326,7 +347,14 @@ def _kernel_mom0( "RIRPA Zeroth Moment wall time: %s", time_string(timer() - t_start) ) - return mom0 + moment_offset, (mom0_ub, mom_lb) + return mom0, (mom0_ub, mom_lb) + + def mult_apbinv(self, integral, ri_apb_inv): + if self.compress > 5: + ri_apb_inv = self.compress_low_rank(*ri_apb_inv, name="(A+B)^-1") + mom0 = integral * (self.D ** (-1))[None] + mom0 -= dot(dot(integral, ri_apb_inv[0].T), ri_apb_inv[1]) + return mom0 def test_eta0_error(self, mom0, target_rot, ri_apb, ri_amb): """Test how well our obtained zeroth moment obeys relation used to derive it, namely @@ -378,9 +406,7 @@ def kernel_trMPrt(self, npoints=48, ainit=10): niworker = energy_NI.NITrRootMP(*inputs) integral, err = niworker.kernel(a=ainit, opt_quad=True) # Compute offset; possible analytically in N^3 for diagonal. - offset = sum(self.D) + 0.5 * einsum( - "p,np,np->", self.D ** (-1), ri_mp[0], ri_mp[1] - ) + offset = sum(self.D) + 0.5 * np.tensordot(ri_mp[0] * (self.D ** (-1)), ri_mp[1], ((0, 1), (0, 1))) return integral[0] + offset, err def kernel_energy(self, npoints=48, ainit=10, correction="linear"): @@ -388,9 +414,11 @@ def kernel_energy(self, npoints=48, ainit=10, correction="linear"): t_start = timer() e1, err = self.kernel_trMPrt(npoints, ainit) e2 = 0.0 - ri_apb_eri = self.get_apb_eri_ri() + ri_apb_eri, ri_apb_eri_neg = self.get_apb_eri_ri() # Note that eri contribution to A and B is equal, so can get trace over one by dividing by two e3 = sum(self.D) + einsum("np,np->", ri_apb_eri, ri_apb_eri) / 2 + if ri_apb_eri_neg is not None: + e3 -= einsum("np,np->", ri_apb_eri_neg, ri_apb_eri_neg) / 2 err /= 2 if self.rixc is not None and correction is not None: if correction.lower() == "linear": @@ -435,8 +463,8 @@ def direct_AC_integration( local_rot describes the rotation of the ov-excitations to the space of local excitations within a cluster, while fragment_projectors gives the projector within this local excitation space to the actual fragment.""" # Get the coulomb integrals. - ri_eri = self.get_apb_eri_ri() / np.sqrt(2) - + ri_eri, ri_eri_neg = self.get_apb_eri_ri() / np.sqrt(2) + # TODO use ri_eri_neg here. def get_eta_alpha(alpha, target_rot): newrirpa = self.__class__(self.mf, rixc=self.rixc, log=self.log) moms, errs = newrirpa.kernel_moms( @@ -452,6 +480,10 @@ def run_ac_inter(func, deg=5): weights /= 2 return sum([w * func(p) for w, p in zip(weights, points)]) + if ri_eri_neg is not None: + raise NotImplementedError("Use of negative CDERI contributions with direct integration for the RIRPA " + "correlation energy not yet supported.") + naux_eri = ri_eri.shape[0] if local_rot is None or fragment_projectors is None: lrot = ri_eri @@ -634,52 +666,52 @@ def check_errors(self, error, nelements): def construct_RI_AB(self): """Construct the RI expressions for the deviation of A+B and A-B from D.""" - ri_apb_eri = self.get_apb_eri_ri() + ri_apb_eri, ri_neg_apb_eri = self.get_apb_eri_ri() # Use empty AmB contrib initially; this is the dRPA contrib. - ri_amb_eri = np.zeros((0, self.ov * 2)) + ri_amb_eri = np.zeros((0, ri_apb_eri.shape[1])) + + apb_lhs = [ri_apb_eri] + apb_rhs = [ri_apb_eri] + + if ri_neg_apb_eri is not None: + # Negative is factored in on only one side. + apb_lhs += [ri_neg_apb_eri] + apb_rhs += [-ri_neg_apb_eri] + + amb_lhs = [ri_amb_eri] + amb_rhs = [ri_amb_eri] + if self.rixc is not None: ri_a_xc, ri_b_xc = self.get_ab_xc_ri() - ri_apb_xc = [ - np.concatenate([ri_a_xc[0], ri_b_xc[0]], axis=0), - np.concatenate([ri_a_xc[1], ri_b_xc[1]], axis=0), - ] - ri_amb_xc = [ - np.concatenate([ri_a_xc[0], ri_b_xc[0]], axis=0), - np.concatenate([ri_a_xc[1], -ri_b_xc[1]], axis=0), - ] - else: - ri_apb_xc = [np.zeros((0, self.ov * 2))] * 2 - ri_amb_xc = [np.zeros((0, self.ov * 2))] * 2 + apb_lhs += [np.concatenate([ri_a_xc[0], ri_b_xc[0]], axis=0)] + apb_rhs += [np.concatenate([ri_a_xc[1], ri_b_xc[1]], axis=0)] - ri_apb = [np.concatenate([ri_apb_eri, x], axis=0) for x in ri_apb_xc] - ri_amb = [np.concatenate([ri_amb_eri, x], axis=0) for x in ri_amb_xc] + amb_lhs += [np.concatenate([ri_a_xc[0], ri_b_xc[0]], axis=0)] + amb_rhs += [np.concatenate([ri_a_xc[1], -ri_b_xc[1]], axis=0)] - return ri_apb, ri_amb + return (tuple([np.concatenate(x, axis=0) for x in [apb_lhs, apb_rhs]]), + tuple([np.concatenate(x, axis=0) for x in [amb_lhs, amb_rhs]])) def compress_low_rank(self, ri_l, ri_r, name=None): return compress_low_rank(ri_l, ri_r, tol=self.svd_tol, log=self.log, name=name) def get_apb_eri_ri(self): # Coulomb integrals only contribute to A+B. - # This needs to be optimised, but will do for now. - if self.Lpq is None: - v = self.get_3c_integrals() # pyscf.lib.unpack_tril(self.mf._cderi) - Lov = einsum( - "npq,pi,qa->nia", v, self.mo_coeff_occ, self.mo_coeff_vir - ).reshape((self.naux_eri, self.ov)) - else: - Lov = self.Lpq[:, : self.nocc, self.nocc :].reshape( - (self.naux_eri, self.ov) - ) + lov, lov_neg = self.get_cderi() - ri_apb_eri = np.zeros((self.naux_eri, self.ov_tot)) + lov = lov.reshape((lov.shape[0], -1)) + if lov_neg is not None: + lov_neg = lov_neg.reshape((lov_neg.shape[0], -1)) # Need to include factor of two since eris appear in both A and B. - ri_apb_eri[:, : self.ov] = ri_apb_eri[:, self.ov : 2 * self.ov] = ( - np.sqrt(2) * Lov - ) - return ri_apb_eri + ri_apb_eri = np.sqrt(2) * np.concatenate([lov, lov], axis = 1) + + ri_neg_apb_eri = None + if lov_neg is not None: + ri_neg_apb_eri = np.sqrt(2) * np.concatenate([lov_neg, lov_neg], axis=1) + + return ri_apb_eri, ri_neg_apb_eri def get_ab_xc_ri(self): # Have low-rank representation for interactions over and above coulomb interaction. @@ -714,8 +746,17 @@ def get_ab_xc_ri(self): ri_b_xc = [np.concatenate([x, y], axis=1) for x, y in zip(ri_b_aa, ri_b_bb)] return ri_a_xc, ri_b_xc - def get_3c_integrals(self): - return pyscf.lib.unpack_tril(next(self.mf.with_df.loop(blksize=self.naux_eri))) + def get_cderi(self, blksize=None): + if self.lov is None: + return get_cderi(self, (self.mo_coeff_occ, self.mo_coeff_vir), compact=False, blksize=blksize) + else: + if isinstance(self.lov, tuple): + lov, lov_neg = self.lov + else: + assert self.lov.shape == (self.naux_eri, self.nocc, self.nvir) + lov = self.lov + lov_neg = None + return lov, lov_neg def test_spectral_rep(self, freqs): from vayesta.rpa import ssRPA @@ -757,6 +798,8 @@ def get_log_specvals(freq): return log_qvals, log_specvals, get_log_qval, get_log_specvals + def record_memory(self): + self.log.info(" %s", memory_string()) def construct_product_RI(D, ri_1, ri_2): """Given two matrices expressed as low-rank modifications, cderi_1 and cderi_2, of some full-rank matrix D, @@ -808,14 +851,15 @@ def construct_inverse_RI(D, ri): def compress_low_rank(ri_l, ri_r, tol=1e-12, log=None, name=None): naux_init = ri_l.shape[0] - u, s, v = np.linalg.svd(ri_l, full_matrices=False) - nwant = sum(s > tol) - rot = u[:, :nwant] - ri_l = dot(rot.T, ri_l) - ri_r = dot(rot.T, ri_r) - u, s, v = np.linalg.svd(ri_r, full_matrices=False) - nwant = sum(s > tol) - rot = u[:, :nwant] + + inner_prod = dot(ri_l, ri_r.T) + + e, c = np.linalg.eig(inner_prod) + + want = e > tol + nwant = sum(want) + rot = c[:, want] + ri_l = dot(rot.T, ri_l) ri_r = dot(rot.T, ri_r) if nwant < naux_init and log is not None: @@ -833,6 +877,3 @@ def compress_low_rank(ri_l, ri_r, tol=1e-12, log=None, name=None): nwant, ) return ri_l, ri_r - - -ssRIRRPA = ssRIRPA diff --git a/vayesta/rpa/rirpa/RIURPA.py b/vayesta/rpa/rirpa/RIURPA.py index 5206cb154..2c1f3b4eb 100644 --- a/vayesta/rpa/rirpa/RIURPA.py +++ b/vayesta/rpa/rirpa/RIURPA.py @@ -2,6 +2,11 @@ from vayesta.core.util import einsum from vayesta.rpa.rirpa.RIRPA import ssRIRRPA +import pyscf.lib +from vayesta.core.util import einsum +from vayesta.core.eris import get_cderi + +from .RIRPA import ssRIRRPA class ssRIURPA(ssRIRRPA): @@ -90,56 +95,28 @@ def D(self): D = np.concatenate([epsa, epsb]) return D - def construct_RI_AB(self): - """Construct the RI expressions for the deviation of A+B and A-B from D.""" - ri_apb_eri = self.get_apb_eri_ri() - # Use empty AmB contrib initially; this is the dRPA contrib. - ri_amb_eri = np.zeros((0, self.ov_tot)) - if self.rixc is not None: - ri_a_xc, ri_b_xc = self.get_ab_xc_ri() - - ri_apb_xc = [ - np.concatenate([ri_a_xc[0], ri_b_xc[0]], axis=0), - np.concatenate([ri_a_xc[1], ri_b_xc[1]], axis=0), - ] - ri_amb_xc = [ - np.concatenate([ri_a_xc[0], ri_b_xc[0]], axis=0), - np.concatenate([ri_a_xc[1], -ri_b_xc[1]], axis=0), - ] - else: - ri_apb_xc = [np.zeros((0, self.ov_tot))] * 2 - ri_amb_xc = [np.zeros((0, self.ov_tot))] * 2 - - ri_apb = [np.concatenate([ri_apb_eri, x], axis=0) for x in ri_apb_xc] - ri_amb = [np.concatenate([ri_amb_eri, x], axis=0) for x in ri_amb_xc] - - return ri_apb, ri_amb - def get_apb_eri_ri(self): # Coulomb integrals only contribute to A+B. # This needs to be optimised, but will do for now. - if self.Lpq is None: - v = self.get_3c_integrals() - Lov_a = einsum( - "npq,pi,qa->nia", v, self.mo_coeff_occ[0], self.mo_coeff_vir[0] - ).reshape((self.naux_eri, self.ov[0])) - Lov_b = einsum( - "npq,pi,qa->nia", v, self.mo_coeff_occ[1], self.mo_coeff_vir[1] - ).reshape((self.naux_eri, self.ov[1])) - else: - Lov_a = self.Lpq[0][:, : self.nocc[0], self.nocc[0] :].reshape( - (self.naux_eri, self.ov[0]) - ) - Lov_b = self.Lpq[1][:, : self.nocc[1], self.nocc[1] :].reshape( - (self.naux_eri, self.ov[1]) - ) + (lova, lovb), (lova_neg, lovb_neg) = self.get_cderi() - ri_apb_eri = np.zeros((self.naux_eri, sum(self.ov))) + lova = lova.reshape((lova.shape[0], -1)) + lovb = lovb.reshape((lovb.shape[0], -1)) + if lova_neg is not None: + if lovb_neg is None: + raise RuntimeError("Encountered negative cderi contribution in only one spin channel." + "Isn't this impossible?") + lova_neg = lova_neg.reshape((lova_neg.shape[0], -1)) + lovb_neg = lovb_neg.reshape((lovb_neg.shape[0], -1)) # Need to include factor of two since eris appear in both A and B. - ri_apb_eri[:, : self.ov[0]] = np.sqrt(2) * Lov_a - ri_apb_eri[:, self.ov[0] : self.ov_tot] = np.sqrt(2) * Lov_b - return ri_apb_eri + ri_apb_eri = np.sqrt(2) * np.concatenate([lova, lovb], axis=1) + + ri_neg_apb_eri = None + if lova_neg is not None: + ri_neg_apb_eri = np.sqrt(2) * np.concatenate([lova_neg, lovb_neg], axis=1) + + return ri_apb_eri, ri_neg_apb_eri def get_ab_xc_ri(self): # Have low-rank representation for interactions over and above coulomb interaction. @@ -179,3 +156,17 @@ def get_ab_xc_ri(self): ri_a_xc = [np.concatenate([x, y], axis=1) for x, y in zip(ri_a_aa, ri_a_bb)] ri_b_xc = [np.concatenate([x, y], axis=1) for x, y in zip(ri_b_aa, ri_b_bb)] return ri_a_xc, ri_b_xc + + def get_cderi(self, blksize=None): + if self.lov is None: + la, la_neg = get_cderi(self, (self.mo_coeff_occ[0], self.mo_coeff_vir[0]), compact=False, blksize=blksize) + lb, lb_neg = get_cderi(self, (self.mo_coeff_occ[1], self.mo_coeff_vir[1]), compact=False, blksize=blksize) + else: + if isinstance(self.lov, tuple): + (la, lb), (la_neg, lb_neg) = self.lov + else: + assert self.lov[0][0].shape == (self.naux_eri, self.nocc[0], self.nvir[0]) + assert self.lov[0][1].shape == (self.naux_eri, self.nocc[1], self.nvir[1]) + la, lb = self.lov + la_neg = lb_neg = None + return (la, lb), (la_neg, lb_neg) diff --git a/vayesta/rpa/rirpa/RIdRRPA.py b/vayesta/rpa/rirpa/RIdRRPA.py new file mode 100644 index 000000000..0aff47c7d --- /dev/null +++ b/vayesta/rpa/rirpa/RIdRRPA.py @@ -0,0 +1,265 @@ +import numpy as np + +from vayesta.rpa.rirpa.RIRPA import ssRIRRPA +from vayesta.rpa.rirpa.momzero_NI import MomzeroOffsetCalcGaussLag, MomzeroDeductHigherOrder_dRHF +from vayesta.rpa.rirpa.energy_NI import NITrRootMP_dRHF +from vayesta.core.util import dot, time_string, timer, with_doc + +from pyscf import lib + + +class ssRIdRRPA(ssRIRRPA): + """Class for computing direct RPA correlated quantites with a restricted reference state. + """ + + @with_doc(ssRIRRPA.kernel_moms) + def kernel_moms(self, max_moment, target_rot=None, return_spatial=False, **kwargs): + t_start = timer() + self.log.debug("Running specialised dRPA RHF code.") + + if target_rot is None: + self.log.warning( + "Warning; generating full moment rather than local component. Will scale as O(N^5)." + ) + if return_spatial: + target_rot = np.eye(self.ov) + else: + target_rot = np.eye(self.ov_tot) + + ri_decomps = self.get_compressed_MP() + ri_mp, ri_apb, ri_amb = ri_decomps + # First need to calculate zeroth moment. This only generates the spin-independent contribution in a single + # spin channel; the spin-dependent contribution is just the identity rotated to our target rotation. + mom0, err0 = self._kernel_mom0( + target_rot, ri_decomps=ri_decomps, return_spatial=return_spatial, **kwargs + ) + + moments = np.zeros((max_moment + 1,) + mom0.shape, dtype=mom0.dtype) + + moments[0] = mom0 + + t_start_higher = timer() + + if max_moment == 0: + return moments, err0 + + eps = self.eps + + if return_spatial: + # Must have spatial target rotation. + def gen_new_moment(prev_mom): + return prev_mom * (eps**2)[None] + 2 * dot(prev_mom, ri_mp[1].T, ri_mp[0]) + + moments[1] = target_rot * eps[None] + + else: + def gen_new_moment(prev_mom): + prev_aa, prev_bb = prev_mom[:, :self.ov], prev_mom[:, self.ov:] + spat_vv = dot(prev_aa + prev_bb, ri_mp[1].T, ri_mp[0]) + new_aa = spat_vv + prev_aa * (eps**2)[None] + new_bb = spat_vv + prev_bb * (eps**2)[None] + return np.concatenate((new_aa, new_bb), axis=1) + + if target_rot.shape[1] == self.ov: + moments[1, :, :self.ov] = target_rot * eps[None] + else: + moments[1] = target_rot * self.D[None] + + if max_moment > 1: + for i in range(2, max_moment + 1): + moments[i] = gen_new_moment(moments[i - 2]) + self.record_memory() + if max_moment > 0: + self.log.info( + "RIRPA Higher Moments wall time: %s", + time_string(timer() - t_start_higher), + ) + self.log.info( + "Overall RIRPA Moments wall time: %s", time_string(timer() - t_start) + ) + return moments, err0 + + @with_doc(ssRIRRPA._kernel_mom0) + def _kernel_mom0( + self, + target_rot=None, + npoints=48, + ainit=10, + integral_deduct="HO", + opt_quad=True, + adaptive_quad=False, + alpha=1.0, + ri_decomps=None, + return_niworker=False, + analytic_lower_bound=False, + return_spatial=False + ): + t_start = timer() + if analytic_lower_bound or adaptive_quad or integral_deduct != "HO": + raise NotImplementedError("Only core functionality is implemented in dRPA specific code.") + # If we have a rotation in the spinorbital basis we need to stack the different spin channels as additional + # spatial rotations and then sum at the end. + target_rot, stack_spin = self.check_target_rot(target_rot) + trrot = None + # We can then compress the spatial rotation. + if stack_spin: + if return_spatial: + raise ValueError("Requested spatially integrated calculation, but target_rot is spin-dependent.") + target_rot, trrot = self.compress_target_rot(target_rot) + + if ri_decomps is None: + ri_mp, ri_apb, ri_amb = self.get_compressed_MP(alpha) + else: + ri_mp, ri_apb, ri_amb = ri_decomps + self.record_memory() + + offset_niworker = None + inputs = (self.eps, ri_mp[0], ri_mp[1], target_rot, npoints, self.log) + + # We are computing our values using spatial quantities, but relating our calculations to the spinorbital + # resolved basis. As such, we need to be careful to ensure we know which terms are spin-diagonal and which are + # spin-invariant. + + niworker = MomzeroDeductHigherOrder_dRHF(*inputs) + offset_niworker = MomzeroOffsetCalcGaussLag(*inputs) + + if return_niworker: + return niworker, offset_niworker + self.record_memory() + + integral, upper_bound = niworker.kernel(a=ainit, opt_quad=opt_quad) # This contribution is spin-invariant. + integral += offset_niworker.kernel()[0] # As is this one. + + self.record_memory() + + # Free memory. + del ri_mp, inputs, niworker + self.record_memory() + + # Construct A+B inverse. + eps = self.eps + epsinv = eps ** (-1) + # Factor of two from different spin channels. + u = 2 * dot(ri_apb[0] * epsinv[None], ri_apb[1].T) + u = np.linalg.inv(np.eye(u.shape[0]) + u) + self.record_memory() + + # First, compute contribution which is spin-dependent (diagonal in both the integrated value, and (A+B)^-1). + mom0_spinindependent = integral * epsinv[None] + + # Factor of two from sum over spin channels. + mom0_spinindependent -= dot( + dot( + dot( + target_rot * self.eps[None] + 2 * integral, (ri_apb[1] * epsinv[None]).T + ), + u + ), + ri_apb[0] * epsinv[None] + ) + + + self.log.info( + "RIRPA Zeroth Moment wall time: %s", time_string(timer() - t_start) + ) + + if trrot is not None: + target_rot = dot(trrot, target_rot) + mom0_spinindependent = dot(trrot, mom0_spinindependent) + + if return_spatial: + mom0 = target_rot + 2 * mom0_spinindependent + else: + # Want to return quantities in spin-orbital basis. + n = target_rot.shape[0] + if stack_spin: + # Half of target rot are actually other spin. + n = n// 2 + mom0 = np.zeros((n, self.ov_tot)) + + if stack_spin: + mom0[:, :self.ov] = target_rot[:n] + mom0_spinindependent[:n] + mom0_spinindependent[n:] # Has aa and ab interactions. + mom0[:, self.ov:] = target_rot[n:] + mom0_spinindependent[n:] + mom0_spinindependent[:n] # Has bb and ba interactions. + else: + mom0[:, :self.ov] = target_rot + mom0_spinindependent + mom0[:, self.ov:] = mom0_spinindependent + + return mom0, (None, None) + + def kernel_trMPrt(self, npoints=48, ainit=10): + """Evaluate Tr((MP)^(1/2)).""" + ri_mp, ri_apb, ri_amb = self.get_compressed_MP() + inputs = (self.eps, ri_mp[0], ri_mp[1], npoints, self.log) + + niworker = NITrRootMP_dRHF(*inputs) + integral, err = niworker.kernel(a=ainit, opt_quad=True) + # Compute offset; possible analytically in N^3 for diagonal. + offset = 2 * sum(self.eps) + np.tensordot(ri_mp[0] * self.eps ** (-1), ri_mp[1], ((0, 1), (0, 1))) + return integral[0] + offset, err + + def kernel_energy(self, npoints=48, ainit=10, correction="linear"): + + t_start = timer() + e1, err = self.kernel_trMPrt(npoints, ainit) + e2 = 0.0 + cderi, cderi_neg = self.get_cderi() + cderi = cderi.reshape((-1, self.ov)) + if cderi_neg is not None: + cderi_neg = cderi_neg.reshape((-1, self.ov)) + # Note that eri contribution to A and B is equal, so can get trace over one by dividing by two + e3 = 2 * (sum(self.eps) + np.tensordot(cderi, cderi, ((0, 1), (0, 1)))) + if cderi_neg is not None: + e3 -= 2 * np.tensordot(cderi_neg, cderi_neg, ((0, 1), (0, 1))) + err /= 2 + self.e_corr_ss = 0.5 * (e1 + e2 - e3) + self.log.info( + "Total RIRPA Energy Calculation wall time: %s", + time_string(timer() - t_start), + ) + return self.e_corr_ss, err + + def get_compressed_MP(self, alpha=1.0): + lov, lov_neg = self.get_cderi() + + lov = lov.reshape(lov.shape[0], self.ov) + ri_apb = [lov, lov] + + if lov_neg is not None: + lov_neg = lov_neg.reshape(lov_neg.shape[0], self.ov) + + ri_apb = [np.concatenate([lov, lov_neg], axis=0), + np.concatenate([lov, -lov_neg], axis=0)] + + if self.compress > 3: + ri_apb = self.compress_low_rank(*ri_apb, name="A+B") + + ri_apb = [x * np.sqrt(2) * (alpha ** (0.5)) for x in ri_apb] + ri_amb = [np.zeros((0, lov.shape[1]))] * 2 + + ri_mp = [ri_apb[0] * self.eps[None], ri_apb[1]] + if self.compress > 0: + ri_mp = self.compress_low_rank(*ri_mp, name="(A-B)(A+B)") + return ri_mp, ri_apb, ri_amb + + def check_target_rot(self, target_rot): + stack_spins = False + if target_rot is None: + self.log.warning( + "Warning; generating full moment rather than local component. Will scale as O(N^5)." + ) + target_rot = np.eye(self.ov) + elif target_rot.shape[1] == self.ov_tot: + # Provided rotation is in spinorbital space. We want to convert to spatial, but record that we've done this + # so we can convert back later. + target_rot = np.concatenate([target_rot[:, :self.ov], target_rot[:, self.ov:]], axis=0) + stack_spins = True + return target_rot, stack_spins + + def compress_target_rot(self, target_rot, tol=1e-10): + inner_prod = dot(target_rot, target_rot.T) + e, c = np.linalg.eigh(inner_prod) + want = e > tol + rot = c[:, want] + if rot.shape[1] == target_rot.shape[1]: + return target_rot, None + return dot(rot.T, target_rot), rot diff --git a/vayesta/rpa/rirpa/__init__.py b/vayesta/rpa/rirpa/__init__.py index fdcc7c94d..73a2c762e 100644 --- a/vayesta/rpa/rirpa/__init__.py +++ b/vayesta/rpa/rirpa/__init__.py @@ -1,9 +1,13 @@ from vayesta.rpa.rirpa.RIRPA import ssRIRRPA from vayesta.rpa.rirpa.RIURPA import ssRIURPA +from vayesta.rpa.rirpa.RIdRRPA import ssRIdRRPA import pyscf.scf def ssRIRPA(mf, *args, **kwargs): if isinstance(mf, pyscf.scf.uhf.UHF): return ssRIURPA(mf, *args, **kwargs) - return ssRIRRPA(mf, *args, **kwargs) + if "rixc" in kwargs: + if kwargs["rixc"] is not None: + return ssRIRRPA(mf, *args, **kwargs) + return ssRIdRRPA(mf, *args, **kwargs) diff --git a/vayesta/rpa/rirpa/energy_NI.py b/vayesta/rpa/rirpa/energy_NI.py index 136c80e81..bc9740743 100644 --- a/vayesta/rpa/rirpa/energy_NI.py +++ b/vayesta/rpa/rirpa/energy_NI.py @@ -23,9 +23,9 @@ def __init__(self, D, S_L, S_R, npoints, log): out_shape = (1,) diag_shape = (1,) super().__init__(out_shape, diag_shape, npoints, log, True) - self.diagmat1 = self.D**2 + einsum("np,np->p", self.S_L, self.S_R) - self.diagmat2 = self.D**2 self.diagRI = einsum("np,np->p", self.S_L, self.S_R) + self.diagmat1 = self.D**2 + self.diagRI + self.diagmat2 = self.D**2 @property def n_aux(self): @@ -39,7 +39,7 @@ def get_Q(self, freq): """Efficiently construct Q = S_R (D^{-1} G) S_L^T This is generally the limiting """ - S_L = einsum("np,p->np", self.S_L, self.get_F(freq)) + S_L = self.S_L * self.get_F(freq)[None] return dot(self.S_R, S_L.T) @property @@ -104,6 +104,18 @@ def eval_contrib(self, freq): Q = self.get_Q(freq) F = self.get_F(freq) val_aux = np.linalg.inv(np.eye(self.n_aux) + Q) - np.eye(self.n_aux) - res = einsum("p,np,nm,mp,p->", F, self.S_L, val_aux, self.S_R, F) + lhs = dot(val_aux, self.S_L * F[None]) + res = np.tensordot(lhs, self.S_R * F[None], ((0, 1), (0, 1))) res = (freq**2) * res / np.pi return np.array([res]) + + +class NITrRootMP_dRHF(NITrRootMP): + """All provided quantities are now in spatial orbitals. This actually only requires an additional factor in the + get_Q method.""" + def get_Q(self, freq): + # Have equal contributions from both spin channels. + return 2 * super().get_Q(freq) + + def eval_contrib(self, freq): + return 2 * super().eval_contrib(freq) diff --git a/vayesta/rpa/rirpa/momzero_NI.py b/vayesta/rpa/rirpa/momzero_NI.py index cde9c7b21..28bb00ca6 100644 --- a/vayesta/rpa/rirpa/momzero_NI.py +++ b/vayesta/rpa/rirpa/momzero_NI.py @@ -32,10 +32,10 @@ def get_F(self, freq): return (self.D**2 + freq**2) ** (-1) def get_Q(self, freq): - """Efficiently construct Q = S_R (D^{-1} G) S_L^T - This is generally the limiting + """Efficiently construct Q = S_R F S_L^T + This is generally the limiting step. """ - S_L = einsum("np,p->np", self.S_L, self.get_F(freq)) + S_L = self.S_L * self.get_F(freq)[None] return dot(self.S_R, S_L.T) @property @@ -102,10 +102,10 @@ def eval_contrib(self, freq): Q = self.get_Q(freq) rrot = F - lrot = einsum("lq,q->lq", self.target_rot, rrot) + lrot = self.target_rot * rrot[None] val_aux = np.linalg.inv(np.eye(self.n_aux) + Q) - lres = np.dot(lrot, self.S_L.T) - res = dot(dot(lres, val_aux), einsum("np,p->np", self.S_R, rrot)) + lres = dot(lrot, self.S_L.T) + res = dot(dot(lres, val_aux), self.S_R * rrot[None]) return (self.target_rot + (freq**2) * (res - lrot)) / np.pi @@ -169,7 +169,7 @@ def eval_contrib(self, freq): lrot = einsum("lq,q->lq", self.target_rot, F) val_aux = np.linalg.inv(np.eye(self.n_aux) + Q) - np.eye(self.n_aux) res = dot( - dot(dot(lrot, self.S_L.T), val_aux), einsum("np,p->np", self.S_R, rrot) + dot(dot(lrot, self.S_L.T), val_aux), self.S_R * rrot[None] ) res = (freq**2) * res / np.pi return res @@ -194,9 +194,9 @@ def get_offset(self): def eval_contrib(self, freq): # This should be real currently, so can safely do this. expval = np.exp(-freq * self.D) - lrot = einsum("lp,p->lp", self.target_rot, expval) + lrot = self.target_rot * expval[None] rrot = expval - res = dot(dot(lrot, self.S_L.T), einsum("np,p->np", self.S_R, rrot)) + res = dot(dot(lrot, self.S_L.T), self.S_R * rrot[None]) return res def eval_diag_contrib(self, freq): @@ -238,3 +238,13 @@ def diag_sqrt_grad(D, freq): def diag_sqrt_deriv2(D, freq): M = (D + freq**2) ** (-1) return (-2 * M + 10 * (freq**2) * (M**2) - 8 * (freq**4) * (M**3)) / np.pi + + +# Subclass for performing calculations with RHF quantities. + +class MomzeroDeductHigherOrder_dRHF(MomzeroDeductHigherOrder): + """All provided quantities are now in spatial orbitals. This actually only requires an additional factor in the + get_Q method.""" + def get_Q(self, freq): + # Have equal contributions from both spin channels. + return 2 * super().get_Q(freq) diff --git a/vayesta/rpa/ssrpa.py b/vayesta/rpa/ssrpa.py index 4dbb0d152..0dd042d00 100644 --- a/vayesta/rpa/ssrpa.py +++ b/vayesta/rpa/ssrpa.py @@ -94,8 +94,13 @@ def kernel(self, xc_kernel=None, alpha=1.0): AmBrtinv = einsum("pn,n,qn->pq", c2, e2 ** (-0.5), c2) XpY = np.einsum("n,pq,qn->pn", self.freqs_ss ** (-0.5), AmBrt, c) XmY = np.einsum("n,pq,qn->pn", self.freqs_ss ** (0.5), AmBrtinv, c) - self.XpY_ss = (XpY[: self.ova], XpY[self.ova :]) - self.XmY_ss = (XmY[: self.ova], XmY[self.ova :]) + + nov = self.ova + if self.ov_rot is not None: + nov = self.ov_rot[0].shape[0] + + self.XpY_ss = (XpY[:nov], XpY[nov:]) + self.XmY_ss = (XmY[:nov], XmY[nov:]) self.freqs_sf = None self.log.timing("Time to solve RPA problem: %s", time_string(timer() - t0)) @@ -199,8 +204,7 @@ def get_val_alpha(alpha): return e - def _gen_arrays(self, xc_kernel=None, alpha=1.0): - t0 = timer() + def _gen_eps(self): # Only have diagonal components in canonical basis. eps = np.zeros((self.nocc, self.nvir)) eps = eps + self.mf.mo_energy[self.nocc :] @@ -208,18 +212,30 @@ def _gen_arrays(self, xc_kernel=None, alpha=1.0): eps = eps.reshape((self.ova,)) if self.ov_rot is not None: - eps = einsum("pn,n,qn->pq", self.ov_rot, eps, self.ov_rot) - # want to ensure diagonalise eps in our new basis; nb this also rotates the existing rotation for - # consistency. - eps, c = scipy.linalg.eigh(eps) - self.ov_rot = dot(c.T, self.ov_rot) + if self.ov_rot[0].shape[0] > 0: + epsa = einsum("pn,n,qn->pq", self.ov_rot[0], eps, self.ov_rot[0]) + epsa, ca = scipy.linalg.eigh(epsa) + self.ov_rot = (dot(ca.T, self.ov_rot[0]), self.ov_rot[1]) + if self.ov_rot[1].shape[0] > 0: + epsb = einsum("pn,n,qn->pq", self.ov_rot[1], eps, self.ov_rot[1]) + epsb, cb = scipy.linalg.eigh(epsb) + self.ov_rot = (self.ov_rot[0], dot(cb.T, self.ov_rot[1])) + else: + epsa = epsb = eps + return epsa, epsb + + + def _gen_arrays(self, xc_kernel=None, alpha=1.0): + t0 = timer() + + epsa, epsb = self._gen_eps() - AmB = np.concatenate([eps, eps]) + AmB = np.concatenate([epsa, epsb]) fullv = self.get_k() ApB = 2 * fullv * alpha if self.ov_rot is not None: - fullrot = scipy.linalg.block_diagonal(self.ov_rot, self.ov_rot) + fullrot = scipy.linalg.block_diag(self.ov_rot[0], self.ov_rot[1]) ApB = dot(fullrot, ApB, fullrot.T) # At this point AmB is just epsilon so add in. @@ -239,7 +255,7 @@ def _gen_arrays(self, xc_kernel=None, alpha=1.0): M = dot(AmBrt, ApB, AmBrt) self.log.timing("Time to build RPA arrays: %s", time_string(timer() - t0)) - return M, AmB, ApB, (eps, eps), fullv + return M, AmB, ApB, (epsa, epsb), fullv def get_k(self): eris = self.ao2mo() @@ -319,30 +335,21 @@ def get_xc_contribs(self, xc_kernel, c_o, c_v, alpha=1.0): return ApB * alpha, AmB * alpha def gen_moms(self, max_mom, xc_kernel=None): - res = np.zeros((max_mom + 1, self.ov, self.ov)) + nova = self.ova + nov = self.ov + if self.ov_rot is not None: + nova = self.ov_rot[0].shape[0] + nov = nova + self.ov_rot[1].shape[0] + res = np.zeros((max_mom+1, nov, nov)) for x in range(max_mom + 1): # Have different spin components in general; these are alpha-alpha, alpha-beta and beta-beta. - res[x, : self.ova, : self.ova] = np.einsum( - "pn,n,qn->pq", self.XpY_ss[0], self.freqs_ss**x, self.XpY_ss[0] - ) - res[x, self.ova :, self.ova :] = np.einsum( - "pn,n,qn->pq", self.XpY_ss[1], self.freqs_ss**x, self.XpY_ss[1] - ) - res[x, : self.ova, self.ova :] = np.einsum( - "pn,n,qn->pq", self.XpY_ss[0], self.freqs_ss**x, self.XpY_ss[1] - ) - res[x, self.ova :, : self.ova] = res[x][: self.ova, self.ova :].T + res[x, :nova, :nova] = np.einsum("pn,n,qn->pq", self.XpY_ss[0], self.freqs_ss ** x, self.XpY_ss[0]) + res[x, nova:, nova:] = np.einsum("pn,n,qn->pq", self.XpY_ss[1], self.freqs_ss ** x, self.XpY_ss[1]) + res[x, :nova, nova:] = np.einsum("pn,n,qn->pq", self.XpY_ss[0], self.freqs_ss ** x, self.XpY_ss[1]) + res[x, nova:, :nova] = res[x][:nova, nova:].T return res - @property - def mo_coeff(self): - return self.mf.mo_coeff - - @property - def nao(self): - return self.mf.mol.nao - def ao2mo(self, mo_coeff=None, compact=False): """Get the ERIs in MO basis""" mo_coeff = self.mo_coeff if mo_coeff is None else mo_coeff diff --git a/vayesta/rpa/ssurpa.py b/vayesta/rpa/ssurpa.py index ec6a62383..48a353e20 100644 --- a/vayesta/rpa/ssurpa.py +++ b/vayesta/rpa/ssurpa.py @@ -3,7 +3,7 @@ import scipy.linalg from timeit import default_timer as timer -from vayesta.core.util import dot, time_string +from vayesta.core.util import dot, time_string, einsum class ssURPA(ssRPA): @@ -43,9 +43,7 @@ def mo_coeff_vir(self): na, nb = self.nocc return self.mo_coeff[0][:, na:], self.mo_coeff[1][:, nb:] - def _gen_arrays(self, xc_kernel=None, alpha=1.0): - t0 = timer() - + def _gen_eps(self): nocc_a, nocc_b = self.nocc nvir_a, nvir_b = self.nvir # Only have diagonal components in canonical basis. @@ -58,6 +56,11 @@ def _gen_arrays(self, xc_kernel=None, alpha=1.0): epsb = epsb + self.mf.mo_energy[1][nocc_b:] epsb = (epsb.T - self.mf.mo_energy[1][:nocc_b]).T epsb = epsb.reshape((self.ovb,)) + return epsa, epsb + + def _gen_arrays(self, xc_kernel=None, alpha=1.0): + t0 = timer() + epsa, epsb = self._gen_eps() if self.ov_rot is not None: epsa = einsum("pn,n,qn->pq", self.ov_rot[0], epsa, self.ov_rot[0]) @@ -78,7 +81,7 @@ def _gen_arrays(self, xc_kernel=None, alpha=1.0): fullv = self.get_k() ApB = 2 * fullv * alpha if self.ov_rot is not None: - fullrot = scipy.linalg.block_diagonal(self.ov_rot[0], self.ov_rot[1]) + fullrot = scipy.linalg.block_diag(self.ov_rot[0], self.ov_rot[1]) ApB = dot(fullrot, ApB, fullrot.T) # At this point AmB is just epsilon so add in. diff --git a/vayesta/solver/__init__.py b/vayesta/solver/__init__.py index 703a6ba2d..3c0103ca0 100644 --- a/vayesta/solver/__init__.py +++ b/vayesta/solver/__init__.py @@ -2,7 +2,6 @@ from vayesta.solver.cisd import RCISD_Solver, UCISD_Solver from vayesta.solver.coupled_ccsd import coupledRCCSD_Solver from vayesta.solver.dump import DumpSolver -from vayesta.solver.ebcc import REBCC_Solver, UEBCC_Solver, EB_REBCC_Solver, EB_UEBCC_Solver from vayesta.solver.ebfci import EB_EBFCI_Solver, EB_UEBFCI_Solver from vayesta.solver.ext_ccsd import extRCCSD_Solver, extUCCSD_Solver from vayesta.solver.fci import FCI_Solver, UFCI_Solver @@ -10,6 +9,12 @@ from vayesta.solver.mp2 import RMP2_Solver, UMP2_Solver from vayesta.solver.tccsd import TRCCSD_Solver +try: + from vayesta.solver.ebcc import REBCC_Solver, UEBCC_Solver, EB_REBCC_Solver, EB_UEBCC_Solver +except ImportError: + _has_ebcc = False +else: + _has_ebcc = True def get_solver_class(ham, solver): assert (is_ham(ham)) @@ -36,7 +41,6 @@ def _get_solver_class(is_uhf, is_eb, solver, log): def _get_solver_class_internal(is_uhf, is_eb, solver): - solver = solver.upper() # First check if we have a CC approach as implemented in pyscf. if solver == "CCSD" and not is_eb: # Use pyscf solvers. @@ -48,13 +52,13 @@ def _get_solver_class_internal(is_uhf, is_eb, solver): if is_uhf or is_eb: raise ValueError("TCCSD is not implemented for unrestricted or electron-boson calculations!") return TRCCSD_Solver - if solver == "EXTCCSD": + if solver == "extCCSD": if is_eb: raise ValueError("extCCSD is not implemented for electron-boson calculations!") if is_uhf: return extUCCSD_Solver return extRCCSD_Solver - if solver == "COUPLEDCCSD": + if solver == "coupledCCSD": if is_eb: raise ValueError("coupledCCSD is not implemented for electron-boson calculations!") if is_uhf: @@ -62,7 +66,11 @@ def _get_solver_class_internal(is_uhf, is_eb, solver): return coupledRCCSD_Solver # Now consider general CC ansatzes; these are solved via EBCC. - if "CC" in solver: + # Note that we support all capitalisations of `ebcc`, but need `CC` to be capitalised when also using this to + # specify an ansatz. + if "CC" in solver.upper(): + if not _has_ebcc: + raise ImportError(f"{solver} solver is only accessible via ebcc. Please install ebcc.") if is_uhf: if is_eb: solverclass = EB_UEBCC_Solver @@ -73,23 +81,24 @@ def _get_solver_class_internal(is_uhf, is_eb, solver): solverclass = EB_REBCC_Solver else: solverclass = REBCC_Solver - if solver == "EBCC": + if solver.upper() == "EBCC": # Default to `opts.ansatz`. return solverclass - if solver[:2] == "EB": + if solver[:2].upper() == "EB": solver = solver[2:] if solver == "CCSD" and is_eb: # Need to specify CC level for coupled electron-boson model; throw an error rather than assume. raise ValueError( "Please specify a coupled electron-boson CC ansatz as a solver, for example CCSD-S-1-1," "rather than CCSD") - + # This is just a wrapper to allow us to use the solver option as the ansatz kwarg in this case. def get_right_CC(*args, **kwargs): - - if kwargs.get("ansatz", None) is not None: - raise ValueError( - "Desired CC ansatz specified differently in solver and solver_options.ansatz." - "Please use only specify via one approach, or ensure they agree.") + setansatz = kwargs.get("ansatz", None) + if setansatz is not None: + if setansatz != solver: + raise ValueError( + "Desired CC ansatz specified differently in solver and solver_options.ansatz." + "Please use only specify via one approach, or ensure they agree.") kwargs["ansatz"] = solver return solverclass(*args, **kwargs) diff --git a/vayesta/solver/ebcc.py b/vayesta/solver/ebcc.py index 93a72938e..eadf3f08d 100644 --- a/vayesta/solver/ebcc.py +++ b/vayesta/solver/ebcc.py @@ -5,6 +5,7 @@ from vayesta.core.types import WaveFunction, CCSD_WaveFunction from vayesta.core.util import dot, einsum from vayesta.solver.solver import ClusterSolver, UClusterSolver +import ebcc class REBCC_Solver(ClusterSolver): @@ -16,31 +17,43 @@ class Options(ClusterSolver.Options): max_cycle: int = 200 # Max number of iterations conv_tol: float = None # Convergence energy tolerance conv_tol_normt: float = None # Convergence amplitude tolerance + store_as_ccsd: bool = True # Store results as CCSD_WaveFunction + c_cas_occ: np.array = None # Hacky place to put active space orbitals. + c_cas_vir: np.array = None - def __init__(self, *args, **kwargs): - super().__init__(*args, **kwargs) - try: - import ebcc - except ImportError as e: - raise ImportError("Cannot import ebcc; required to use as a solver.") @property def is_fCCSD(self): - return self.opts.ansatz == "CCSD" + return self.opts.store_as_ccsd or self.opts.ansatz == "CCSD" def kernel(self): - import ebcc # Use pyscf mean-field representation. mf_clus, frozen = self.hamil.to_pyscf_mf(allow_dummy_orbs=False) - mycc = ebcc.EBCC(mf_clus, log=self.log, ansatz=self.opts.ansatz, - **self.get_nonnull_solver_opts()) + + mf_clus.mo_coeff, space = self.get_space(self.hamil.cluster.c_active, mf_clus.mo_occ, frozen=frozen) + + mycc = ebcc.EBCC(mf_clus, log=self.log, ansatz=self.opts.ansatz, space=space, **self.get_nonnull_solver_opts()) mycc.kernel() self.converged = mycc.converged if self.opts.solve_lambda: mycc.solve_lambda() self.converged = self.converged and mycc.converged_lambda # Now just need to wrangle EBCC results into wavefunction format. - self.construct_wavefunction(mycc, self.hamil.mo) + self.construct_wavefunction(mycc, self.hamil.get_mo(mf_clus.mo_coeff)) + + def get_space(self, mo_coeff, mo_occ, frozen=None): + s = self.hamil.orig_mf.get_ovlp() + c_active_occ = self.opts.c_cas_occ + c_active_vir = self.opts.c_cas_vir + if (c_active_occ is not None and c_active_vir is not None): + c_active_occ = dot(mo_coeff.T, s, c_active_occ) + c_active_vir = dot(mo_coeff.T, s, c_active_vir) + elif (c_active_occ is not None or c_active_vir is not None): + raise ValueError( + "Both or neither of the occupied and virtual components of an active space must be specified.") + mo_coeff = dot(mo_coeff.T, s, mo_coeff) + return gen_space(mo_coeff[:, mo_occ > 0], mo_coeff[:, mo_occ == 0], c_active_occ, c_active_vir, + frozen_orbs=frozen) def get_nonnull_solver_opts(self): def add_nonull_opt(d, key, newkey): @@ -60,12 +73,20 @@ def construct_wavefunction(self, mycc, mo): self.wf = CCSD_WaveFunction(mo, mycc.t1, mycc.t2, l1=mycc.l1.T, l2=mycc.l2.transpose(2, 3, 0, 1)) except TypeError: self.wf = CCSD_WaveFunction(mo, mycc.t1, mycc.t2) + self.wf.rotate(t=mycc.mo_coeff.T, inplace=True) else: - # Simply alias required quantities for now; this ensures functionality for arbitrary orders of CC via ebcc. + if not (np.allclose(mycc.mo_coeff, np.eye(mycc.nmo), atol=1e-8)): + raise ValueError("Generic wavefunction construction only works for canonical orbitals, please use " + "`store_as_ccsd=True' in solver_options") + self.log.warning( + "Storing EBCC wavefunction as generic wavefunction object; wavefunction-based estimators will not work.") + # Simply alias required quantities for now; this could allow functionality for arbitrary orders of CC via ebcc. self.wf = WaveFunction(mo) self.wf.make_rdm1 = mycc.make_rdm1_f self.wf.make_rdm2 = mycc.make_rdm2_f + # Need to rotate wavefunction back into original cluster active space. + def _debug_exact_wf(self, wf): assert (self.is_fCCSD) mo = self.hamil.mo @@ -96,7 +117,26 @@ def _debug_random_wf(self): class UEBCC_Solver(UClusterSolver, REBCC_Solver): @dataclasses.dataclass class Options(REBCC_Solver.Options, UClusterSolver.Options): - pass + c_cas_occ: np.array = (None, None) # Hacky place to put active space orbitals. + c_cas_vir: np.array = (None, None) + + def get_space(self, mo_coeff, mo_occ, frozen=None): + + s = self.hamil.orig_mf.get_ovlp() + def _get_space(c, occ, co_cas, cv_cas, fr): + # Express active orbitals in terms of cluster orbitals. + if (co_cas is not None and cv_cas is not None): + co_cas = dot(c.T, s, co_cas) + cv_cas = dot(c.T, s, cv_cas) + elif (co_cas is not None or cv_cas is not None): + raise ValueError("Both or neither of the occupied and virtual components of an active space.") + c = dot(c.T, s, c) + return gen_space(c[:, occ > 0], c[:, occ == 0], co_cas, cv_cas, frozen_orbs=fr) + if frozen is None: + frozen = [None, None] + ca, spacea = _get_space(mo_coeff[0], mo_occ[0], self.opts.c_cas_occ[0], self.opts.c_cas_vir[0], frozen[0]) + cb, spaceb = _get_space(mo_coeff[1], mo_occ[1], self.opts.c_cas_occ[1], self.opts.c_cas_vir[1], frozen[1]) + return (ca, cb), (spacea, spaceb) # This should automatically work other than ensuring spin components are in a tuple. def construct_wavefunction(self, mycc, mo): @@ -105,8 +145,15 @@ def construct_wavefunction(self, mycc, mo): def to_spin_tuple1(x): return x.aa, x.bb + # NB EBCC doesn't antisymmetrise the T2s by default, and has an odd factor of two. + + def antisymmetrize_t2(t2): + t2 = t2 - t2.transpose(1, 0, 2, 3) + t2 = t2 - t2.transpose(0, 1, 3, 2) + return t2 / 2.0 + def to_spin_tuple2(x): - return x.aaaa, x.abab, x.bbbb + return antisymmetrize_t2(x.aaaa), x.abab, antisymmetrize_t2(x.bbbb) try: self.wf = CCSD_WaveFunction(mo, @@ -120,10 +167,16 @@ def to_spin_tuple2(x): to_spin_tuple1(mycc.t1), to_spin_tuple2(mycc.t2) ) + self.wf.rotate(t=[x.T for x in mycc.mo_coeff], inplace=True) else: # Simply alias required quantities for now; this ensures functionality for arbitrary orders of CC. self.wf = WaveFunction(mo) - + if not (np.allclose(mycc.mo_coeff[0], np.eye(mycc.nmo), atol=1e-8) and + np.allclose(mycc.mo_coeff[1], np.eye(mycc.nmo), atol=1e-8)): + raise ValueError("Generic wavefunction construction only works for canonical orbitals, please use " + "`store_as_ccsd=True' in solver_options") + self.log.warning( + "Storing EBCC wavefunction as generic wavefunction object; wavefunction-based estimators will not work.") def make_rdm1(*args, **kwargs): dm = mycc.make_rdm1_f(*args, **kwargs) return (dm.aa, dm.bb) @@ -137,7 +190,7 @@ def make_rdm2(*args, **kwargs): self.wf.make_rdm2 = make_rdm2 def _debug_exact_wf(self, wf): - mo = self.cluster.mo + mo = self.hamil.mo # Project onto cluster: ovlp = self.hamil._fragment.base.get_ovlp() roa = dot(wf.mo.coeff_occ[0].T, ovlp, mo.coeff_occ[0]) @@ -172,9 +225,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): @@ -218,3 +274,70 @@ def make_rdmeb(*args, **kwargs): return dm self.wf.make_rdmeb = make_rdmeb + + +def gen_space(c_occ, c_vir, co_active=None, cv_active=None, frozen_orbs=None): + """Given the occupied and virtual orbital coefficients in the local cluster basis, which orbitals are frozen, and + any active space orbitals in this space generate appropriate coefficients and ebcc.Space inputs for a calculation. + Inputs: + c_occ: occupied orbital coefficients in local cluster basis. + c_vir: virtual orbital coefficients in local cluster basis. + co_active: occupied active space orbitals in local cluster basis. + cv_active: virtual active space orbitals in local cluster basis. + frozen_orbs: indices of orbitals to freeze in local cluster basis. + Outputs: + c: coefficients for the active space orbitals in the local cluster basis. + space: ebcc.Space object defining the resulting active space behaviours. + """ + no, nv = c_occ.shape[1], c_vir.shape[1] + + have_actspace = not (co_active is None and cv_active is None) + + if co_active is None or cv_active is None: + if have_actspace: + raise ValueError("Active space must specify both occupied and virtual orbitals.") + # Default to just using original cluster orbitals. + c = np.hstack([c_occ, c_vir]) + + occupied = np.hstack([np.ones(no, dtype=bool), np.zeros(nv, dtype=bool)]) + frozen = np.zeros(no + nv, dtype=bool) + if frozen_orbs is not None: + frozen[frozen_orbs] = True + + def gen_actspace(c_full, c_act, c_frozen, tol=1e-8): + """Generate orbital rotation to define our space, along with frozen and active identifiers.""" + # Check we don't have any overlap between frozen and active orbitals. + if c_frozen.shape[1] > 0: + af_ovlp = dot(c_act.T, c_frozen) + if not np.allclose(np.linalg.svd(af_ovlp)[1], 0, atol=tol): + raise ValueError("Specified frozen and active orbitals overlap!") + if c_act.shape[1] > 0: + # We could use the portion of the active space orbitals inside the cluster instead in future. + full_active_ovlp = dot(c_full.T, c_act) + if not np.allclose(np.linalg.svd(full_active_ovlp)[1], 1, atol=tol): + raise ValueError("Specified active orbitals are not spanned by the full cluster space!") + # Can now safely generate coefficients; get the space spanned by undetermined orbitals. + d_rest = dot(c_full, c_full.T) - dot(c_act, c_act.T) - dot(c_frozen, c_frozen.T) + e, c = np.linalg.eigh(d_rest) + c_rest = c[:, e > tol] + c = np.hstack([c_frozen, c_rest, c_act]) + # Check that we have the right number of orbitals. + assert(c.shape[1] == c_full.shape[1]) + # Generate information. + frozen_orbs = np.zeros(c.shape[1], dtype=bool) + active = np.zeros(c.shape[1], dtype=bool) + frozen_orbs[:c_frozen.shape[1]] = True + active[-c_act.shape[1]:] = True + return c, frozen_orbs, active + + if have_actspace: + c_occ, frozen_occ, active_occ = gen_actspace(c_occ, co_active, c_occ[:, frozen[:no]]) + c_vir, frozen_vir, active_vir = gen_actspace(c_vir, cv_active, c_vir[:, frozen[no:]]) + + c = np.hstack([c_occ, c_vir[:, ::-1]]) + frozen = np.hstack([frozen_occ, frozen_vir[::-1]]) + active = np.hstack([active_occ, active_vir[::-1]]) + else: + active = np.zeros_like(frozen, dtype=bool) + + return c, ebcc.Space(occupied=occupied, frozen=frozen, active=active) diff --git a/vayesta/solver/hamiltonian.py b/vayesta/solver/hamiltonian.py index 2b7832230..4d3a17a1b 100644 --- a/vayesta/solver/hamiltonian.py +++ b/vayesta/solver/hamiltonian.py @@ -4,12 +4,11 @@ import numpy as np import pyscf.lib import pyscf.scf -import scipy.linalg -from vayesta.core.qemb import scrcoulomb +from vayesta.core.util import dot, einsum, OptionsBase, break_into_lines, log_time, energy_string +from vayesta.core.screening import screening_moment, screening_crpa from vayesta.core.types import Orbitals -from vayesta.core.util import dot, einsum, OptionsBase, break_into_lines, log_time -from vayesta.rpa import ssRPA +from typing import Optional def is_ham(ham): @@ -60,6 +59,7 @@ class RClusterHamiltonian: class Options(OptionsBase): screening: Optional[str] = None cache_eris: bool = True + match_fock: bool = True @property def _scf_class(self): @@ -86,7 +86,7 @@ def cluster(self): @property def mo(self): - return Orbitals(self.cluster.c_active, occ=self.cluster.nocc_active) + return self.get_mo() @property def nelec(self): @@ -104,6 +104,15 @@ def target_space_projector(self, c=None): """Projector to the target fragment space within our cluster.""" return self._fragment.get_fragment_projector(self.cluster.c_active, c) + def get_mo(self, mo_coeff=None): + c = self.cluster.c_active + if mo_coeff is not None: + # This specifies the coefficients in the cluster basis used within the calculation. + # Note that we should then apply the opposite rotation to the resulting wavefunction, since Vayesta + # currently assumes the wavefunction is in the basis of c_active. + c = dot(c, mo_coeff) + return Orbitals(c, occ=self.cluster.nocc_active) + # Integrals for the cluster calculations. def get_integrals(self, bare_eris=None, with_vext=True): @@ -112,7 +121,11 @@ def get_integrals(self, bare_eris=None, with_vext=True): seris = self.get_eris_screened() return heff, seris - def get_fock(self, with_vext=True, use_seris=True, with_exxdiv=False): + def get_fock(self, with_vext=True, use_seris=None, with_exxdiv=False): + if use_seris is None: + # Actual default depends on self.opts.match_fock; to reproduce fock we will effectively modify heff + use_seris = not self.opts.match_fock + c = self.cluster.c_active fock = dot(c.T, self._fragment.base.get_fock(with_exxdiv=with_exxdiv), c) if with_vext and self.v_ext is not None: @@ -133,9 +146,9 @@ def get_fock(self, with_vext=True, use_seris=True, with_exxdiv=False): def get_heff(self, eris=None, fock=None, with_vext=True, with_exxdiv=False): if eris is None: - eris = self.get_eris_bare() + eris = self.get_eris_screened() if fock is None: - fock = self.get_fock(with_vext=False, use_seris=False, with_exxdiv=with_exxdiv) + fock = self.get_fock(with_vext=False, with_exxdiv=with_exxdiv) occ = np.s_[:self.cluster.nocc_active] v_act = 2 * einsum('iipq->pq', eris[occ, occ]) - einsum('iqpi->pq', eris[occ, :, :, occ]) h_eff = fock - v_act @@ -325,17 +338,20 @@ def pad(a, diag_val): else: # Just set up heff using standard bare eris. bare_eris = self.get_eris_bare() - heff = pad_to_match(self.get_heff(eris=bare_eris, with_vext=True), dummy_energy) if force_bare_eris: clusmf._eri = pad_to_match(bare_eris, 0.0) else: clusmf._eri = pad_to_match(self.get_eris_screened()) + heff = pad_to_match(self.get_heff(with_vext=True), dummy_energy) + + clusmf.get_hcore = lambda *args, **kwargs: heff if overwrite_fock: - clusmf.get_fock = lambda *args, **kwargs: pad_to_match(self.get_fock(with_vext=True), dummy_energy) - clusmf.get_veff = lambda *args, **kwargs: np.array(clusmf.get_fock(*args, **kwargs)) - np.array( - clusmf.get_hcore()) + clusmf.get_fock = lambda *args, **kwargs: pad_to_match( + self.get_fock(with_vext=True, use_seris=not force_bare_eris), dummy_energy) + clusmf.get_veff = lambda *args, **kwargs: np.array(clusmf.get_fock(*args, **kwargs)) - \ + np.array(clusmf.get_hcore()) return clusmf, orbs_to_freeze @@ -358,48 +374,145 @@ def get_clus_mf_info(self, ao_basis=False, with_vext=True, with_exxdiv=False): # Functionality for use with screened interactions and external corrections. - def add_screening(self, seris_intermed=None): - """Add screened interactions into the Hamiltonian.""" - if self.opts.screening == "mrpa": - assert (seris_intermed is not None) + def calc_loc_erpa(self, m0, amb, only_cumulant=False): - # Use bare coulomb interaction from hamiltonian; this could well be cached in future. - bare_eris = self.get_eris_bare() + no, nv = self.cluster.nocc_active, self.cluster.nvir_active + nov = no * nv + # Bare coulomb interaction in cluster ov-ov space. + v = self.get_eris_bare()[:no, no:, :no, no:].reshape((nov, nov)) + ro = self._fragment.get_overlap("fragment|cluster-occ") + po = dot(ro.T, ro) + + def gen_spin_components(mat): + return mat[:nov, :nov], mat[:nov, nov:], mat[nov:, nov:] - self._seris = scrcoulomb.get_screened_eris_full(bare_eris, seris_intermed) + m0_aa, m0_ab, m0_bb = gen_spin_components(m0) - elif self.opts.screening == "crpa": - raise NotImplementedError() + 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 - 2 * np.eye(nov)), v) + self.log.info("Computed fragment RPA cumulant energy contribution for cluster %s as %s", + self._fragment.id, + energy_string(erpa)) + return erpa else: - raise ValueError("Unknown cluster screening protocol: %s" % self.opts.screening) + d_aa, d_ab, d_bb = gen_spin_components(amb) + # This should be zero. + assert (abs(d_ab).max() < 1e-10) + + def compute_e_rrpa(proj): + def pr(m): + m = m.reshape((no, nv, no * nv)) + m = np.tensordot(proj, m, axes=[0, 0]) + return m.reshape((no * nv, no * nv)) + + erpa = 0.5 * (einsum("pq,qp->", pr(m0_aa), d_aa) + einsum("pq,qp->", pr(m0_bb), d_bb)) + erpa += einsum("pq,qp->", pr(m0_aa + m0_ab + m0_ab.T + m0_bb), v) + erpa -= 0.5 * (pr(d_aa + v + d_bb + v).trace()) + self.log.info("Computed fragment RPA energy contribution for cluster %s as %s", self._fragment.id, + energy_string(erpa)) + return erpa + + compute_e_rrpa(np.eye(no)) + + return compute_e_rrpa(po) + + def add_screening(self, seris_intermed=None): + """ + Adds appropriate screening according to the value of self.opts.screening. + -`None`: gives bare interactions, but this function shouldn't be called in that case. + -'mrpa': moment-conserving interactions. + -'crpa': gives cRPA interactions. Including 'ov' after 'crpa' will only apply cRPA screening in the o-v channel. + Including 'pcoupling' similarly will account for the polarisability in non-canonical cluster spaces. - def calc_loc_erpa(self): + seris_intermed is only required for mRPA interactions. + """ + self._seris = self._add_screening(seris_intermed, spin_integrate=True) - clusmf = self.to_pyscf_mf(force_bare_eris=True) - clusrpa = ssRPA(clusmf) - M, AmB, ApB, eps, v = clusrpa._gen_arrays() - erpa = clusrpa.kernel() - m0 = clusrpa.gen_moms(0) + def _add_screening(self, seris_intermed=None, spin_integrate=True): - def get_product_projector(): - nocc = self.nelec - nvir = tuple([x - y for x, y in zip(self.ncas, self.nelec)]) - p_occ_frag = self.target_space_projector(self.cluster.c_active_occ) - if (not isinstance(p_occ_frag, tuple)) and np.ndim(p_occ_frag) == 2: - p_occ_frag = (p_occ_frag, p_occ_frag) + def spin_integrate_and_report(m, warn_threshold=1e-6): + spat = (m[0] + m[1] + m[2] + m[1].transpose((2, 3, 0, 1))) / 4.0 - def get_product_projector(p_o, p_v, no, nv): - return einsum("ij,ab->iajb", p_o, p_v).reshape((no * nv, no * nv)) + dev = [abs(x - spat).max() for x in m] + [abs(m[2].transpose(2, 3, 0, 1) - spat).max()] + self.log.info("Largest change in screened interactions due to spin integration: %e", max(dev)) + if max(dev) > warn_threshold: + self.log.warning("Significant change in screened interactions due to spin integration: %e", max(dev)) + return spat - pa = get_product_projector(p_occ_frag[0], np.eye(nvir[0]), nocc[0], nvir[0]) - pb = get_product_projector(p_occ_frag[1], np.eye(nvir[1]), nocc[1], nvir[1]) - return scipy.linalg.block_diag(pa, pb) + if self.opts.screening is None: + raise ValueError("Attempted to add screening to fragment with no screening protocol specified.") + if self.opts.screening == "mrpa": + assert (seris_intermed is not None) + # Use bare coulomb interaction from hamiltonian. + bare_eris = self.get_eris_bare() + seris = screening_moment.get_screened_eris_full(bare_eris, seris_intermed[0], log=self.log) + if spin_integrate: + seris = spin_integrate_and_report(seris) + elif self.opts.screening[:4] == "crpa": + bare_eris = self.get_eris_bare() + delta, crpa = screening_crpa.get_frag_deltaW(self.orig_mf, self._fragment, + pcoupling=("pcoupled" in self.opts.screening), + only_ov_screened=("ov" in self.opts.screening), + log=self.log) + if "store" in self.opts.screening: + self.log.warning("Storing cRPA object in Hamiltonian- O(N^4) memory cost!") + self.crpa = crpa + if "full" in self.opts.screening: + # Add a check just in case. + self.log.critical("Static screening of frequency-dependent interactions not supported") + self.log.critical("This statement should be impossible to reach!") + raise ValueError("Static screening of frequency-dependent interactions not supported") + else: + if spin_integrate: + delta = spin_integrate_and_report(delta) + seris = bare_eris + delta + else: + seris = tuple([x + y for x, y in zip(bare_eris, delta)]) + else: + raise ValueError("Unknown cluster screening protocol: %s" % self.opts.screening) - proj = get_product_projector() - eloc = 0.5 * einsum("pq,qr,rp->", proj, m0, ApB) - einsum("pq,qp->", proj, ApB + AmB) - return eloc + def report_screening(screened, bare, spins): + maxidx = np.unravel_index(np.argmax(abs(screened - bare)), bare.shape) + if spins is None: + wstring = "W" + else: + wstring = "W(%2s|%2s)" % (2 * spins[0], 2 * spins[1]) + self.log.info( + "Maximally screened element of %s: V= %.3e -> W= %.3e (delta= %.3e)", + wstring, bare[maxidx], screened[maxidx], screened[maxidx] - bare[maxidx]) + # self.log.info( + # " Corresponding norms%s: ||V||= %.3e, ||W||= %.3e, ||delta||= %.3e", + # " " * len(wstring), np.linalg.norm(bare), np.linalg.norm(screened), + # np.linalg.norm(screened-bare)) + + if spin_integrate: + report_screening(seris, bare_eris, None) + else: + report_screening(seris[0], bare_eris[0], "aa") + report_screening(seris[1], bare_eris[1], "ab") + report_screening(seris[2], bare_eris[2], "bb") + + def get_sym_breaking(norm_aa, norm_ab, norm_bb): + spinsym = abs(norm_aa - norm_bb) / ((norm_aa + norm_bb) / 2) + spindep = abs((norm_aa + norm_bb) / 2 - norm_ab) / ((norm_aa + norm_bb + norm_ab) / 3) + return spinsym, spindep + + bss, bsd = get_sym_breaking(*[np.linalg.norm(x) for x in bare_eris]) + sss, ssd = get_sym_breaking(*[np.linalg.norm(x) for x in seris]) + dss, dsd = get_sym_breaking(*[np.linalg.norm(x - y) for x, y in zip(bare_eris, seris)]) + + self.log.info("Proportional spin symmetry breaking in norms: V= %.3e, W= %.3e, (W-V= %.3e)", bss, sss, dss) + self.log.info("Proportional spin dependence in norms: V= %.3e, W= %.3e, (W-V= %.3e)", bsd, ssd, dsd) + return seris def assert_equal_spin_channels(self, message=""): na, nb = self.ncas @@ -460,6 +573,12 @@ def ncas(self): def nelec(self): return self.cluster.nocc_active + def get_mo(self, mo_coeff=None): + c = self.cluster.c_active + if mo_coeff is not None: + c = tuple([dot(x, y) for x, y in zip(c, mo_coeff)]) + return Orbitals(c, occ=self.cluster.nocc_active) + # Integrals for the cluster calculations. def get_fock(self, with_vext=True, use_seris=True, with_exxdiv=False): @@ -613,6 +732,10 @@ def __contains__(self, item): fock = self.get_fock(with_vext=with_vext, use_seris=not force_bare, with_exxdiv=with_exxdiv) return DummyERIs(getter, valid_blocks=ValidUHFKeys(), fock=fock, nocc=self.cluster.nocc_active) + def add_screening(self, seris_intermed=None): + """Add screened interactions into the Hamiltonian.""" + self._seris = self._add_screening(seris_intermed, spin_integrate=False) + class EB_RClusterHamiltonian(RClusterHamiltonian): @dataclasses.dataclass @@ -669,13 +792,17 @@ def get_polaritonic_shifted_couplings(self): couplings = tuple([x - einsum("pq,n->npq", np.eye(x.shape[1]), temp) for x in self.unshifted_couplings]) if not np.allclose(couplings[0], couplings[1]): self.log.critical("Polaritonic shifted bosonic fermion-boson couplings break cluster spin symmetry; please" - "use an unrestricted formalism.") - raise RuntimeError() + " use an unrestricted formalism.") + raise RuntimeError("Polaritonic shifted bosonic fermion-boson couplings break cluster spin symmetry; please" + " use an unrestricted formalism.") return couplings[0] def get_eb_dm_polaritonic_shift(self, dm1): return (-einsum("n,pq->pqn", self.polaritonic_shift, dm1 / 2),) * 2 + def _add_screening(self, seris_intermed=None, spin_integrate=True): + return self.get_eris_bare() + class EB_UClusterHamiltonian(UClusterHamiltonian, EB_RClusterHamiltonian): @dataclasses.dataclass @@ -702,3 +829,6 @@ def get_polaritonic_shifted_couplings(self): def get_eb_dm_polaritonic_shift(self, dm1): return tuple([-einsum("n,pq->pqn", self.polaritonic_shift, x) for x in dm1]) + + def calc_loc_erpa(self, m0, amb): + raise NotImplementedError() diff --git a/vayesta/tests/core/bath/test_bath.py b/vayesta/tests/core/bath/test_bath.py index 1328941d9..1577d572b 100644 --- a/vayesta/tests/core/bath/test_bath.py +++ b/vayesta/tests/core/bath/test_bath.py @@ -3,7 +3,7 @@ import numpy as np from vayesta.core.bath import DMET_Bath from vayesta.core.bath import EwDMET_Bath -from vayesta.core.bath import MP2_Bath +from vayesta.core.bath import MP2_Bath, RPA_Bath from vayesta.core.qemb import Embedding from vayesta.core.qemb import UEmbedding from vayesta.tests.common import TestCase @@ -51,7 +51,7 @@ def test_ewdmet_bath(self): print("Testing EwDMET bath: kmax= %d moment= %d" % (kmax, order)) self.assertIsNone(np.testing.assert_allclose(mom_cluster[order], mom_full[order], atol=1e-7, rtol=1e-7)) -class BNO_Test(TestCase): +class MP2_BNO_Test(TestCase): def test_bno_Bath(self): rhf = testsystems.ethanol_ccpvdz.rhf() @@ -145,6 +145,35 @@ def test_project_dmet(self): self.assertAllclose(rbno_bath_vir.occup, ubno_bath_vir.occup[0]) self.assertAllclose(rbno_bath_vir.occup, ubno_bath_vir.occup[1]) +class RPA_Test(TestCase): + + def test_bno_Bath(self): + rhf = testsystems.ethanol_631g_df.rhf() + + remb = Embedding(rhf) + with remb.iao_fragmentation() as f: + rfrag = f.add_atomic_fragment('O') + rdmet_bath = DMET_Bath(rfrag) + rdmet_bath.kernel() + rbno_bath_occ = RPA_Bath(rfrag, rdmet_bath, occtype='occupied') + rbno_bath_vir = RPA_Bath(rfrag, rdmet_bath, occtype='virtual') + + # Check maximum, minimum, and mean occupations + n_occ_max = 0.008502908095347851 + n_occ_min = 0.000007766938370992 + n_occ_mean = 0.002300040902721369 + n_vir_max = 0.056091218249012025 + n_vir_min = 0.000122187516348096 + n_vir_mean = 0.009570044321105821 + + # RHF + self.assertAlmostEqual(np.amax(rbno_bath_occ.occup), n_occ_max) + self.assertAlmostEqual(np.amin(rbno_bath_occ.occup), n_occ_min) + self.assertAlmostEqual(np.mean(rbno_bath_occ.occup), n_occ_mean) + self.assertAlmostEqual(np.amax(rbno_bath_vir.occup), n_vir_max) + self.assertAlmostEqual(np.amin(rbno_bath_vir.occup), n_vir_min) + self.assertAlmostEqual(np.mean(rbno_bath_vir.occup), n_vir_mean) + if __name__ == '__main__': print('Running %s' % __file__) diff --git a/vayesta/tests/core/scmf/test_scmf.py b/vayesta/tests/core/scmf/test_scmf.py index 9b67e1974..e349113b3 100644 --- a/vayesta/tests/core/scmf/test_scmf.py +++ b/vayesta/tests/core/scmf/test_scmf.py @@ -7,7 +7,7 @@ class SCMF_Test(TestCase): - + solver = "CCSD" @classmethod def setUpClass(cls): cls.mf = testsystems.h2_dz.rhf() @@ -19,7 +19,7 @@ def tearDownClass(cls): @classmethod @cache def emb(cls, scmf=None): - emb = ewf.EWF(cls.mf, solver_options=dict(solve_lambda=True), bath_options=dict(bathtype='dmet')) + emb = ewf.EWF(cls.mf, solver=cls.solver, solver_options=dict(solve_lambda=True), bath_options=dict(bathtype='dmet')) with emb.sao_fragmentation() as f: f.add_all_atomic_fragments() if scmf == 'pdmet': @@ -47,12 +47,34 @@ def test_brueckner(self): self.assertTrue(emb.with_scmf.converged) self.assertAllclose(emb.with_scmf.e_tot, -1.1417339799464736) + class SCMF_UHF_Test(SCMF_Test): @classmethod def setUpClass(cls): cls.mf = testsystems.h2_dz.rhf().to_uhf() +class SCMF_TCCSD_Test(SCMF_Test): + solver = "TCCSD" + + def test_pdmet(self): + """Test p-DMET.""" + emb0 = self.emb() + emb = self.emb('pdmet') + self.assertAllclose(emb.with_scmf.e_tot_oneshot, emb0.e_tot) + self.assertAllclose(emb.with_scmf.e_tot_oneshot, -1.1419060814972688) + self.assertTrue(emb.with_scmf.converged) + self.assertAllclose(emb.with_scmf.e_tot, -1.1406317042658112) + + def test_brueckner(self): + """Test Brueckner DMET.""" + emb0 = self.emb() + emb = self.emb('brueckner') + self.assertAllclose(emb.with_scmf.e_tot_oneshot, emb0.e_tot) + self.assertAllclose(emb.with_scmf.e_tot_oneshot, -1.1419060814979487) + self.assertTrue(emb.with_scmf.converged) + self.assertAllclose(emb.with_scmf.e_tot, -1.1348718457034288) + if __name__ == '__main__': print('Running %s' % __file__) diff --git a/vayesta/tests/ewf/test_rpa_bath.py b/vayesta/tests/ewf/test_rpa_bath.py new file mode 100644 index 000000000..64a709712 --- /dev/null +++ b/vayesta/tests/ewf/test_rpa_bath.py @@ -0,0 +1,41 @@ +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 +import pytest + +class TestWaterRHF(TestCase): + system = testsystems.water_631g_df + + @classmethod + def setUpClass(cls): + cls.mf = cls.system.rhf() + + @classmethod + def tearDownClass(cls): + del cls.mf + cls.emb.cache_clear() + + @classmethod + @cache + def emb(cls, bno_threshold, solver): + emb = vayesta.ewf.EWF(cls.mf, bath_options=dict(bathtype="rpa", threshold=bno_threshold), solver=solver, + solver_options=dict(conv_tol=1e-12)) + emb.kernel() + return emb + + def test_ccsd_rpa_1(self): + eta=10**-(1.5) + emb = self.emb(eta, 'CCSD') + emb.kernel() + self.assertAllclose(emb.e_tot, -76.10582744548097) + + @pytest.mark.slow + def test_ccsd_rpa_2(self): + eta=1e-2 + emb = self.emb(eta, 'CCSD') + emb.kernel() + self.assertAllclose(emb.e_tot, -76.12444492796294) diff --git a/vayesta/tests/ewf/test_rpa_correction.py b/vayesta/tests/ewf/test_rpa_correction.py new file mode 100644 index 000000000..63e9265db --- /dev/null +++ b/vayesta/tests/ewf/test_rpa_correction.py @@ -0,0 +1,47 @@ +import unittest +import numpy as np +import vayesta +import vayesta.ewf +from vayesta.core.util import cache +from vayesta.tests import testsystems +from vayesta.tests.common import TestCase + + +class Test_RPA_Corrections_Ethanol_RHF(TestCase): + + system = testsystems.ethanol_631g_df + + @property + def mf(self): + return self.system.rhf() + + def get_nl_energies(self, correction, bathtype="dmet"): + + emb = vayesta.ewf.EWF(self.mf, bath_options={"bathtype":bathtype}, solver="CCSD", + screening="mrpa", ext_rpa_correction=correction) + with emb.iao_fragmentation() as f: + f.add_all_atomic_fragments() + emb.kernel() + return emb.e_nonlocal, emb.e_rpa + + def test_rpa_correction(self): + enl, erpa = self.get_nl_energies("erpa") + self.assertAlmostEqual(erpa, -0.29138256715100397) + self.assertAlmostEqual(enl, -0.3087486792144427) + + def test_cumulant_correction(self): + enl, erpa = self.get_nl_energies("cumulant") + self.assertAlmostEqual(erpa, -0.5145262339186916) + self.assertAlmostEqual(enl, -0.31553721476368396) + +class Test_RPA_Corrections_complete(Test_RPA_Corrections_Ethanol_RHF): + """Tests with a complete bath in all clusters. This should give no nonlocal correction in any case.""" + system = testsystems.water_631g_df + + def test_rpa_correction(self): + enl, erpa = self.get_nl_energies("erpa", "full") + self.assertAlmostEqual(enl, 0.0) + + def test_cumulant_correction(self): + enl, erpa = self.get_nl_energies("cumulant", "full") + self.assertAlmostEqual(enl, 0.0) diff --git a/vayesta/tests/ewf/test_screening.py b/vayesta/tests/ewf/test_screening.py index 7c7a0401b..0eacdf1e4 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,45 @@ 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} + + +class TestTwoHoleRHF(TestTwoElectron): + + @classmethod + def setUpClass(cls): + cls.mf = cls.system.rhf() + if __name__ == '__main__': print('Running %s' % __file__) diff --git a/vayesta/tests/rpa/test_rirpa_molecule.py b/vayesta/tests/rpa/test_rirpa_molecule.py index 74247cab0..8c7dbcf6f 100644 --- a/vayesta/tests/rpa/test_rirpa_molecule.py +++ b/vayesta/tests/rpa/test_rirpa_molecule.py @@ -28,8 +28,12 @@ def test_n2_ccpvdz_dRIRPA(self): known_values = {"e_tot": -109.27376877774732} self._test_energy(orig_rpa, known_values) - - rirpa = rpa.ssRIRPA(testsystems.n2_ccpvdz_df.rhf()) + # Check dRPA specialised RHF code. + rirpa = rpa.rirpa.ssRIdRRPA(testsystems.n2_ccpvdz_df.rhf()) + rirpa.kernel_energy() + self._test_energy(rirpa, known_values) + # Check spin-generic code. + rirpa = rpa.rirpa.ssRIRRPA(testsystems.n2_ccpvdz_df.rhf()) rirpa.kernel_energy() self._test_energy(rirpa, known_values) @@ -37,11 +41,14 @@ def test_n2_ccpvdz_dRIRPA(self): self._test_mom0(orig_rpa, rirpa_moms) def test_n2_ccpvdz_dRIRPA_error_estimates(self): - rirpa = rpa.ssRIRPA(testsystems.n2_ccpvdz_df.rhf()) + rirpa = rpa.rirpa.ssRIRRPA(testsystems.n2_ccpvdz_df.rhf()) # Use number of points where errors will be meaningful. - rirpa_moms, error_est = rirpa.kernel_moms(0, analytic_lower_bound=True, npoints=4) - self.assertAlmostEqual(error_est[0], 0.0600120085568549,) - self.assertAlmostEqual(error_est[1], 0.0003040334773499559) + # Note that with this few points the fact that the quadrature optimisation is not invariant to orbital rotations + # can cause issues, so we can just use a specific grid spacing. + rirpa_moms, error_est = rirpa.kernel_moms(0, analytic_lower_bound=True, npoints=4, ainit=1.5214230673470202, + opt_quad=False) + self.assertAlmostEqual(error_est[0], 0.05074756294730469) + self.assertAlmostEqual(error_est[1], 0.00024838720802440015) if __name__ == '__main__': print('Running %s' % __file__) diff --git a/vayesta/tests/rpa/test_rirpa_solids.py b/vayesta/tests/rpa/test_rirpa_solids.py new file mode 100644 index 000000000..d82672a88 --- /dev/null +++ b/vayesta/tests/rpa/test_rirpa_solids.py @@ -0,0 +1,62 @@ +import pytest + +from vayesta import rpa +from vayesta.tests.common import TestCase +from vayesta.tests import testsystems + + +class DiamondRIRPATest(TestCase): + PLACES = 8 + + @classmethod + def setUpClass(cls): + cls.sys = testsystems.diamond_sto3g_s211 + cls.known_results = dict(e_tot=-149.51936410641733, e_corr=-0.19193623440986585) + + def _test_energy(self, myrpa): + """Test the RPA energy. + """ + self.assertAlmostEqual(myrpa.e_corr, self.known_results["e_corr"], self.PLACES) + self.assertAlmostEqual(myrpa.e_tot, self.known_results["e_tot"], self.PLACES) + + @pytest.mark.slow + def test_energy_rhf_opt(self): + """Tests for diamond with optimised RHF dRPA code. + """ + + rirpa = rpa.rirpa.ssRIdRRPA(self.sys.rhf()) + rirpa.kernel_energy() + self._test_energy(rirpa) + + @pytest.mark.fast + def test_energy_rhf_generic(self): + """Tests for diamond with generic RHF RIRPA code. + """ + + rirpa = rpa.rirpa.ssRIRRPA(self.sys.rhf()) + rirpa.kernel_energy() + self._test_energy(rirpa) + + @pytest.mark.slow + def test_energy_uhf(self): + """Tests for diamond with generic UHF RIRPA code. + """ + + rirpa = rpa.rirpa.ssRIURPA(self.sys.uhf()) + rirpa.kernel_energy() + self._test_energy(rirpa) + + @pytest.mark.fast + def test_rhf_moments(self): + gen_rirpa = rpa.rirpa.ssRIRRPA(self.sys.rhf()) + opt_rirpa = rpa.rirpa.ssRIdRRPA(self.sys.rhf()) + mom0_gen = gen_rirpa.kernel_moms(0)[0] + mom0_opt = opt_rirpa.kernel_moms(0)[0] + self.assertAllclose(mom0_gen, mom0_opt, self.PLACES) + +@pytest.mark.slow +class GrapheneRIRPATest(DiamondRIRPATest): + @classmethod + def setUpClass(cls): + cls.sys = testsystems.graphene_sto3g_s211 + cls.known_results = dict(e_tot=-150.15057360171875, e_corr=-0.17724246753903117) \ No newline at end of file diff --git a/vayesta/tests/solver/test_ebcc.py b/vayesta/tests/solver/test_ebcc.py new file mode 100644 index 000000000..3fe69ad13 --- /dev/null +++ b/vayesta/tests/solver/test_ebcc.py @@ -0,0 +1,105 @@ +import pytest + +import pyscf +import pyscf.cc + +import vayesta +import vayesta.ewf + +from vayesta.tests.common import TestCase +from vayesta.tests import testsystems + +# Note that ebCC currently doesn't support density fitting, so we're just testing non-DF results here. + + +class TestEBCC(TestCase): + @classmethod + def setUpClass(cls): + try: + import ebcc + except ImportError: + pytest.skip("Requires ebcc") + + def _test(self, system, mf, ansatz): + # Test a complete bath calculation with given ansatz reproduces full calculation. + mf = getattr(getattr(testsystems, system), mf)() + + emb = vayesta.ewf.EWF(mf, solver=f'EB{ansatz}', bath_options=dict(bathtype='full'), + solver_options=dict(solve_lambda=False)) + emb.kernel() + import ebcc + + cc = ebcc.EBCC(mf, ansatz=ansatz) + cc.kernel() + + self.assertAlmostEqual(emb.e_corr, cc.e_corr) + self.assertAlmostEqual(emb.e_tot, cc.e_tot) + + @pytest.mark.fast + def test_rccsd_h2(self): + return self._test('h2_ccpvdz', 'rhf', 'CCSD') + + @pytest.mark.fast + def test_rccsd_water_sto3g(self): + return self._test('water_sto3g', 'rhf', 'CCSD') + + @pytest.mark.fast + def test_uccsd_water_sto3g(self): + return self._test('water_sto3g', 'uhf', 'CCSD') + + def test_uccsd_water_cation_sto3g(self): + return self._test('water_cation_sto3g', 'uhf', 'CCSD') + + @pytest.mark.fast + def test_rccsdt_water_sto3g(self): + return self._test('water_sto3g', 'rhf', 'CCSDT') + + @pytest.mark.slow + def test_uccsdt_water_cation_sto3g(self): + return self._test('water_cation_sto3g', 'uhf', 'CCSDT') + + +class TestEBCCActSpace(TestCase): + @classmethod + def setUpClass(cls): + try: + import ebcc + except ImportError: + pytest.skip("Requires ebcc") + + def _test(self, system, mf, actansatz, fullansatz, bathtype='dmet', setcas=False): + # Test that active space calculation with complete active space reproduces equivalent calculation using higher- + # level approach of active space. This defaults to a DMET bath space. + mf = getattr(getattr(testsystems, system), mf)() + + embfull = vayesta.ewf.EWF(mf, solver=f'EB{fullansatz}', bath_options=dict(bathtype=bathtype), + solver_options=dict(solve_lambda=False)) + embfull.kernel() + + embact = vayesta.ewf.EWF(mf, solver=f'EB{actansatz}', bath_options=dict(bathtype=bathtype), + solver_options=dict(solve_lambda=False)) + if setcas: + # Set up fragmentation, then set CAS to complete cluster space in previous calculation. + with embact.iao_fragmentation() as f: + f.add_all_atomic_fragments() + + for ffull, fact in zip(embfull.loop(), embact.loop()): + fact.set_cas(c_occ=ffull.cluster.c_active_occ, c_vir=ffull.cluster.c_active_vir) + embact.kernel() + + self.assertAlmostEqual(embact.e_corr, embfull.e_corr) + self.assertAlmostEqual(embact.e_tot, embfull.e_tot) + + @pytest.mark.fast + def test_rccsdtprime_water_sto3g_dmet(self): + return self._test('water_sto3g', 'rhf', "CCSDt'", 'CCSDT', bathtype='dmet', setcas=False) + + def test_uccsdtprime_water_sto3g_dmet(self): + return self._test('water_sto3g', 'uhf', "CCSDt'", 'CCSDT', bathtype='dmet', setcas=False) + + @pytest.mark.slow + def test_rccsdtprime_h4_sto3g_setcas_full(self): + return self._test('h4_sto3g', 'rhf', "CCSDt'", 'CCSDT', bathtype='full', setcas=True) + + def test_uccsdtprime_h3_sto3g_setcas_full(self): + return self._test('h3_sto3g', 'uhf', "CCSDt'", 'CCSDT', bathtype='full', setcas=True) diff --git a/vayesta/tests/testsystems.py b/vayesta/tests/testsystems.py index 1c2c07617..d2a2f79ee 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") @@ -370,9 +371,12 @@ def uhf(self): h2_ccpvdz = TestMolecule(atom="H1 0 0 0; H2 0 0 1.0", basis="cc-pvdz") h2_ccpvdz_df = TestMolecule(atom="H1 0 0 0; H2 0 0 1.0", basis="cc-pvdz", auxbasis="cc-pvdz-jkfit") +h3_sto3g = TestMolecule(atom="H1 0 0 0; H2 0 0 1.0; H3 0 1.0 0", basis="sto3g", spin=1) h3_ccpvdz = TestMolecule(atom="H1 0 0 0; H2 0 0 1.0; H3 0 1.0 0", basis="cc-pvdz", spin=1) h3_ccpvdz_df = TestMolecule(atom="H1 0 0 0; H2 0 0 1.0; H3 0 1.0 0", basis="cc-pvdz", auxbasis="cc-pvdz-jkfit", spin=1) +h4_sto3g = TestMolecule(atom="H1 0 0 0; H2 0 0 1.0; H3 0 1.0 0; H4 0 1.0 1.0", basis="sto3g", spin=0) + heli_631g = TestMolecule(atom="He 0 0 0; Li 0 0 2.0", basis='6-31G', spin=1) h6_dz = TestMolecule(atom=molecules.ring('H', 6, 1.0), basis='cc-pVDZ') @@ -445,6 +449,12 @@ def uhf(self): he_k321 = TestSolid(a, atom="He 0 0 0", basis="def2-svp", auxbasis="def2-svp-ri", kmesh=(3, 2, 1)) he_s321 = TestSolid(a, atom="He 0 0 0", basis="def2-svp", auxbasis="def2-svp-ri", supercell=(3, 2, 1)) +a, atom = solids.graphene() +opts = dict(basis='sto3g', auxbasis='sto3g', exp_to_discard=0.1, dimension=2) +mesh = (2,1,1) +graphene_sto3g_k211 = TestSolid(a=a, atom=atom, kmesh=mesh, **opts) +graphene_sto3g_s211 = TestSolid(a=a, atom=atom, supercell=mesh, **opts) + # Lattices