Merge pull request #103 from BoothGroup/moments_PR

Adds option to calculate in-cluster CCSD and FCI moments via Dyson

.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: |

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

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

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

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

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,

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,

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
 

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

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

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