Merge branch 'master' into variational_embedding

.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

docs/source/quickstart/edmetfinite.py

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

examples/ewf/molecules/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)

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))

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))

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()
 

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)

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)

vayesta/core/bath/dmet.py

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

vayesta/core/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