diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index a0f2f1e7e..a9513f827 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -34,6 +34,8 @@ 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 .[dmet,ebcc] --user - name: Run unit tests run: | diff --git a/examples/ewf/molecules/40-spectral-moments.py b/examples/ewf/molecules/40-spectral-moments.py new file mode 100644 index 000000000..f7ba0a437 --- /dev/null +++ b/examples/ewf/molecules/40-spectral-moments.py @@ -0,0 +1,78 @@ +import pyscf +import pyscf.gto +import pyscf.scf +import pyscf.cc +import vayesta +import vayesta.ewf +import vayesta.misc.molecules +import numpy as np + + +mol = pyscf.gto.Mole() +# mol.atom = """ +# H 0.0000 0.0000 0.0000 +# H 0.0000 0.0000 0.7444 +# """ + +mol.atom = """ +O 0.0000 0.0000 0.1173 +H 0.0000 0.7572 -0.4692 +H 0.0000 -0.7572 -0.4692 +""" + +mol.atom = vayesta.misc.molecules.alkane(4) +print(mol.atom) +mol.basis = 'sto-6g' +mol.output = 'pyscf.out' +mol.build() + +# Hartree-Fock +mf = pyscf.scf.RHF(mol) +mf.kernel() + +niter = (5,0) + +# Embedded CCSD +emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=-1), solver_options=dict(solve_lambda=True)) +emb.opts.nmoments= 2*niter[0] + 2 + +with emb.iao_fragmentation() as f: + n = len(mol.atom) + f.add_atomic_fragment([0,1,2,3]) + for i in range(4, n-4, 3): + f.add_atomic_fragment([i,i+1, i+2]) + f.add_atomic_fragment([n-1,n-2,n-3,n-4]) + + +emb.kernel() + +# Reference full system CCSD: +cc = pyscf.cc.CCSD(mf) +cc.kernel() +cc.solve_lambda() + +gfcc = pyscf.cc.gfccsd.GFCCSD(cc, niter=niter) +gfcc.kernel() +ip = gfcc.ipgfccsd(nroots=1)[0] + + +th = gfcc.build_hole_moments() +print(len(th)) +#tp = gfcc.build_part_moments() + +moms = vayesta.ewf.moments.make_ccsdgf_moms(emb) + +print(th) +print(moms[0]) + +for i in range(len(th)): + th[i] = (th[i] + th[i].T)/2 + moms[0][i] = (moms[0][i] + moms[0][i].T)/2 + +for i in range(len(th)): + mask = np.abs(th[i]) > 1e-10 + print("%d mom: norm = %e maxerr = %e"%(i, np.linalg.norm(th[i]-moms[0][i]), (np.abs(th[i]-moms[0][i])).max())) + + +mask = np.abs(th) > 1e-10 +print(np.abs(th-moms[0])) diff --git a/examples/ewf/molecules/41-ccsd-incluster-moments.py b/examples/ewf/molecules/41-ccsd-incluster-moments.py new file mode 100644 index 000000000..d2120f14b --- /dev/null +++ b/examples/ewf/molecules/41-ccsd-incluster-moments.py @@ -0,0 +1,34 @@ +import pyscf +import pyscf.gto +import pyscf.scf +import pyscf.cc +import vayesta +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-pVDZ' +mol.output = 'pyscf.txt' +mol.build() + +# Hartree-Fock +mf = pyscf.scf.RHF(mol) +mf.kernel() + +# Embedded CCSD +emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=1e-6), solver_options=dict(n_moments=(2,4), solve_lambda=True)) +emb.kernel() + +# Reference full system CCSD: +cc = pyscf.cc.CCSD(mf) +cc.kernel() + +print("E(HF)= %+16.8f Ha" % mf.e_tot) +print("E(CCSD)= %+16.8f Ha" % cc.e_tot) +print("E(Emb. CCSD)[WF]= %+16.8f Ha" % emb.e_tot) +print("E(Emb. CCSD)[DM]= %+16.8f Ha" % emb.get_dm_energy()) diff --git a/examples/ewf/molecules/42-fci-incluster-moments.py b/examples/ewf/molecules/42-fci-incluster-moments.py new file mode 100644 index 000000000..6d29e0158 --- /dev/null +++ b/examples/ewf/molecules/42-fci-incluster-moments.py @@ -0,0 +1,34 @@ +import pyscf +import pyscf.gto +import pyscf.scf +import pyscf.cc +import vayesta +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-pVDZ' +mol.output = 'pyscf.txt' +mol.build() + +# Hartree-Fock +mf = pyscf.scf.RHF(mol) +mf.kernel() + +# Embedded CCSD +emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=1e-6), solver_options=dict(n_moments=(5,4), solve_lambda=True)) +emb.kernel() + +# Reference full system CCSD: +cc = pyscf.cc.CCSD(mf) +cc.kernel() + +print("E(HF)= %+16.8f Ha" % mf.e_tot) +print("E(CCSD)= %+16.8f Ha" % cc.e_tot) +print("E(Emb. CCSD)[WF]= %+16.8f Ha" % emb.e_tot) +print("E(Emb. CCSD)[DM]= %+16.8f Ha" % emb.get_dm_energy()) diff --git a/vayesta/__init__.py b/vayesta/__init__.py index a47ce82a9..a464c5b10 100644 --- a/vayesta/__init__.py +++ b/vayesta/__init__.py @@ -93,6 +93,7 @@ def import_package(name, required=True): # Optional import_package('mpi4py', False) import_package('cvxpy', False) +dyson = import_package('dyson', False) ebcc = import_package('ebcc', False) # --- Git hashes @@ -113,6 +114,14 @@ def get_git_hash(dir): pdir = os.path.dirname(os.path.dirname(pyscf.__file__)) phash = get_git_hash(pdir) log.debug(" * PySCF: %s", phash) +if dyson is not None: + ddir = os.path.dirname(os.path.dirname(dyson.__file__)) + dhash = get_git_hash(ddir) + log.debug(" * Dyson: %s", dhash) +if ebcc is not None: + edir = os.path.dirname(os.path.dirname(ebcc.__file__)) + ehash = get_git_hash(edir) + log.debug(" * EBCC: %s", ehash) # --- System information log.debug('System: node= %s processor= %s' % (platform.node(), platform.processor())) diff --git a/vayesta/core/qemb/fragment.py b/vayesta/core/qemb/fragment.py index 0c9023a0a..a6a152de5 100644 --- a/vayesta/core/qemb/fragment.py +++ b/vayesta/core/qemb/fragment.py @@ -76,6 +76,7 @@ class Results: # --- Wave-function wf: WaveFunction = None # WaveFunction object (MP2, CCSD,...) pwf: WaveFunction = None # Fragment-projected wave function + moms: tuple = None def __init__(self, base, fid, name, c_frag, c_env, diff --git a/vayesta/core/qemb/qemb.py b/vayesta/core/qemb/qemb.py index be2810dcd..d28b89d13 100644 --- a/vayesta/core/qemb/qemb.py +++ b/vayesta/core/qemb/qemb.py @@ -94,6 +94,7 @@ class Options(OptionsBase): solver_options: dict = OptionsBase.dict_with_defaults( # General conv_tol=None, + n_moments=None, # CCSD solve_lambda=True, conv_tol_normt=None, level_shift=None, diis_space=None, diis_start_cycle=None, diff --git a/vayesta/ewf/fragment.py b/vayesta/ewf/fragment.py index 689c7d01d..e9aacfffc 100644 --- a/vayesta/ewf/fragment.py +++ b/vayesta/ewf/fragment.py @@ -72,6 +72,7 @@ class Results(BaseFragment.Results): n_active: int = None ip_energy: np.ndarray = None ea_energy: np.ndarray = None + moms: tuple = None @property def dm1(self): @@ -246,9 +247,13 @@ def kernel(self, solver=None, init_guess=None): proj = self.get_overlap('proj|cluster-occ') pwf = pwf.project(proj, inplace=False) + # Moments + + moms = cluster_solver.hole_moments, cluster_solver.particle_moments + # --- 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) + converged=cluster_solver.converged, wf=wf, pwf=pwf, moms=moms) self.hamil = cluster_solver.hamil diff --git a/vayesta/solver/ccsd.py b/vayesta/solver/ccsd.py index bb6d4c19a..3555c8284 100644 --- a/vayesta/solver/ccsd.py +++ b/vayesta/solver/ccsd.py @@ -23,9 +23,12 @@ class Options(ClusterSolver.Options): level_shift: float = None # Level shift for virtual orbitals to stabilize ccsd iterations init_guess: str = 'MP2' # ['MP2', 'CISD'] solve_lambda: bool = True # If false, use Lambda=T approximation + n_moments: tuple = None # Self-consistent mode sc_mode: int = None + + def kernel(self, t1=None, t2=None, l1=None, l2=None, coupled_fragments=None, t_diagnostic=True): mf_clus, frozen = self.hamil.to_pyscf_mf(allow_dummy_orbs=True, allow_df=True) solver_cls = self.get_solver_class(mf_clus) @@ -70,6 +73,35 @@ def kernel(self, t1=None, t2=None, l1=None, l2=None, coupled_fragments=None, t_d self.wf = CCSD_WaveFunction(self.hamil.mo, mycc.t1, mycc.t2, l1=l1, l2=l2) + # In-cluster Moments + nmom = self.opts.n_moments + if nmom is not None: + try: + from dyson.expressions import CCSD + except ImportError: + self.log.error("Dyson not found - required for moment calculations") + self.log.info("Skipping in-cluster moment calculations") + return + self.log.info("Calculating in-cluster CCSD moments %s"%str(nmom)) + # expr = CCSD["1h"](mf_clus, t1=mycc.t1, t2=mycc.t2, l1=l1, l2=l2) + # vecs_bra = expr.build_gf_vectors(nmom[0], left=True) + # amps_bra = [expr.eom.vector_to_amplitudes(amps[n,p], ccm.nmo, ccm.nocc) for p in range(ccm.nmo) for n in range(nmom)] + # vecs_ket = expr.build_gf_vectors(nmom[0], left=False) + # amps_ket = [expr.eom.vector_to_amplitudes(amps[n,p], ccm.nmo, ccm.nocc) for p in range(ccm.nmo) for n in range(nmom)] + # self.ip_moment_amplitudes = (amps_bra, amps_ket) + + # expr = CCSD["1p"](mf_clus, t1=mycc.t1, t2=mycc.t2, l1=l1, l2=l2) + # vecs_bra = expr.build_gf_vectors(nmom[0], left=True) + # amps_bra = [expr.eom.vector_to_amplitudes(amps[n,p], ccm.nmo, ccm.nocc) for p in range(ccm.nmo) for n in range(nmom)] + # vecs_ket = expr.build_gf_vectors(nmom[0], left=False) + # amps_ket = [expr.eom.vector_to_amplitudes(amps[n,p], ccm.nmo, ccm.nocc) for p in range(ccm.nmo) for n in range(nmom)] + # self.ea_moment_amplitudes = (amps_bra, amps_ket) + + expr = CCSD["1h"](mf_clus, t1=mycc.t1, t2=mycc.t2, l1=l1, l2=l2) + self.hole_moments = expr.build_gf_moments(nmom[0]) + expr = CCSD["1p"](mf_clus, t1=mycc.t1, t2=mycc.t2, l1=l1, l2=l2) + self.particle_moments = expr.build_gf_moments(nmom[1]) + def get_solver_class(self, mf): if hasattr(mf, "with_df") and mf.with_df is not None: return pyscf.cc.dfccsd.RCCSD diff --git a/vayesta/solver/fci.py b/vayesta/solver/fci.py index 09c81d59e..cd2308b69 100644 --- a/vayesta/solver/fci.py +++ b/vayesta/solver/fci.py @@ -24,7 +24,7 @@ class Options(ClusterSolver.Options): davidson_only: bool = True init_guess: str = 'default' init_guess_noise: float = 1e-5 - + n_moments: tuple = None cisd_solver = RCISD_Solver def __init__(self, *args, **kwargs): @@ -78,6 +78,22 @@ def kernel(self, ci=None): self.converged = self.solver.converged self.wf = FCI_WaveFunction(self.hamil.mo, self.civec) + # In-cluster Moments + nmom = self.opts.n_moments + if nmom is not None: + try: + from dyson.expressions import FCI + except ImportError: + self.log.error("Dyson not found - required for moment calculations") + self.log.info("Skipping in-cluster moment calculations") + return + self.log.info("Calculating in-cluster FCI moments %s"%str(nmom)) + mf_clus, frozen = self.hamil.to_pyscf_mf(allow_dummy_orbs=True, allow_df=True) + expr = FCI["1h"](mf_clus, c_ci=self.civec) + self.hole_moments = expr.build_gf_moments(nmom[0]) + expr = FCI["1p"](mf_clus, c_ci=self.civec) + self.particle_moments = expr.build_gf_moments(nmom[1]) + class UFCI_Solver(UClusterSolver, FCI_Solver): @dataclasses.dataclass diff --git a/vayesta/solver/solver.py b/vayesta/solver/solver.py index 1cadb15df..09c4b732b 100644 --- a/vayesta/solver/solver.py +++ b/vayesta/solver/solver.py @@ -33,6 +33,8 @@ def __init__(self, hamil, log=None, **kwargs): self.wf = None self.dm1 = None self.dm2 = None + self.hole_moments = None + self.particle_moments = None @property def v_ext(self): diff --git a/vayesta/tests/ewf/test_moments.py b/vayesta/tests/ewf/test_moments.py new file mode 100644 index 000000000..68003d934 --- /dev/null +++ b/vayesta/tests/ewf/test_moments.py @@ -0,0 +1,90 @@ +import unittest +import pytest +import numpy as np + +import pyscf +import pyscf.cc + +import vayesta +import vayesta.ewf +from vayesta.tests.common import TestCase +from vayesta.tests import testsystems + +class Test_RFCI(TestCase): + + def test(self): + + #RHF + mf = testsystems.h6_sto6g.rhf() + + try: + from dyson.expressions import FCI + except ImportError: + pytest.skip("Could not import Dyson. Skipping FCI moment tests.") + return + + fci = FCI["1h"](mf) + fci_ip = fci.build_gf_moments(4) + + fci = FCI["1p"](mf) + fci_ea = fci.build_gf_moments(4) + + #Full bath EWF + ewf = vayesta.ewf.EWF(mf, bath_options=dict(bathtype='full'), solver_options=dict(n_moments=(4,4)), solver='FCI') + ewf.kernel() + + for f in ewf.fragments: + ip, ea = f.results.moms + + cx = f.get_overlap('mo|cluster') + ip = np.einsum('pP,qQ,nPQ->npq', cx, cx, ip) + ea = np.einsum('pP,qQ,nPQ->npq', cx, cx, ea) + + self.assertTrue(np.allclose(ip, fci_ip, atol=1e-14)) + self.assertTrue(np.allclose(ea, fci_ea, atol=1e-14)) + + +# class Test_RCCSD(TestCase): + +# def test(self): + +# #RHF +# mf = testsystems.water_sto3g.rhf() + +# try: +# from dyson.expressions import CCSD +# except ImportError: +# pytest.skip("Could not import Dyson. Skipping CCSD moment tests.") +# return + +# cc = CCSD["1h"](mf) +# cc_ip = cc.build_gf_moments(4) + +# cc = CCSD["1p"](mf) +# cc_ea = cc.build_gf_moments(4) + +# #Full bath EWF +# ewf = vayesta.ewf.EWF(mf, bath_type='full', solver_options=dict(n_moments=(4,4)), solver='CCSD') +# ewf.kernel() + +# for f in ewf.fragments: +# ip, ea = f.results.moms + +# cx = f.get_overlap('mo|cluster') +# ip = np.einsum('pP,qQ,nPQ->npq', cx, cx, ip) +# ea = np.einsum('pP,qQ,nPQ->npq', cx, cx, ea) + +# for i in range(len(ip)): +# print(i, np.trace(np.dot(ip[i], ip[i])), np.trace(np.dot(cc_ip[i], cc_ip[i])), abs(np.trace(np.dot(ip[i], ip[i])) - np.trace(np.dot(cc_ip[i], cc_ip[i])))) +# print('Norm %f'%np.linalg.norm(ip- cc_ip)) +# print('IP %f'%np.linalg.norm(ea- cc_ea)) +# self.assertAlmostEqual(np.linalg.norm(ip- cc_ip), 0.0) +# self.assertAlmostEqual(np.linalg.norm(ea- cc_ea), 0.0) + +# # self.assertTrue(np.allclose(ip, cc_ip)) +# # self.assertTrue(np.allclose(ea, cc_ea)) + +if __name__ == '__main__': + print("Running %s" % __file__) + unittest.main() + diff --git a/vayesta/tests/ewf/test_rdm_energy.py b/vayesta/tests/ewf/test_rdm_energy.py index 8cf69c8f5..d42bb5ace 100644 --- a/vayesta/tests/ewf/test_rdm_energy.py +++ b/vayesta/tests/ewf/test_rdm_energy.py @@ -31,6 +31,7 @@ def test(self): self.assertAlmostEqual(gl, cc.e_corr) self.assertAlmostEqual(lg, cc.e_corr) self.assertAlmostEqual(gg, cc.e_corr) + self.assertAlmostEqual(ewf.get_dm_energy(), cc.e_tot) class Test_UHF(TestCase): @@ -56,6 +57,30 @@ def test(self): self.assertAlmostEqual(gl, cc.e_corr) self.assertAlmostEqual(lg, cc.e_corr) self.assertAlmostEqual(gg, cc.e_corr) + self.assertAlmostEqual(ewf.get_dm_energy(), cc.e_tot) + +# def test_h2_solid(self): +# +# #RHF +# mf = testsystems.h2_sto3g_331_2d.rhf() +# +# #CCSD +# cc = pyscf.cc.CCSD(mf) +# cc.kernel() +# +# #Full bath EWF +# ewf = vayesta.ewf.EWF(mf, bno_threshold=-1) +# ewf.kernel() +# +# ll = ewf.get_rdm2_corr_energy(global_dm1=False, global_dm2=False) +# gl = ewf.get_rdm2_corr_energy(global_dm1=True, global_dm2=False) +# lg = ewf.get_rdm2_corr_energy(global_dm1=False, global_dm2=True) +# gg = ewf.get_rdm2_corr_energy(global_dm1=True, global_dm2=True) +# +# self.assertAlmostEqual(ll, cc.e_corr) +# self.assertAlmostEqual(gl, cc.e_corr) +# self.assertAlmostEqual(lg, cc.e_corr) +# self.assertAlmostEqual(gg, cc.e_corr) if __name__ == '__main__':