Refactors solver/__init__.py

examples/ewdmet/65-callback-solver-dyson.py

@@ -0,0 +1,68 @@
+import numpy as np
+import pyscf
+import pyscf.gto
+import pyscf.scf
+import pyscf.fci
+import vayesta
+import vayesta.ewf
+from vayesta.misc.molecules import ring
+from dyson import FCI
+
+
+# User defined FCI solver - takes pyscf mf as input and returns RDMs
+# The mf argment contains the hamiltonain in the orthonormal cluster basis
+# Pyscf or other solvers may be used to solve the cluster problem and may return RDMs, CISD amplitudes or CCSD amplitudes
+# Returning the cluster Green's function moments is also supported. They are calculated with Dyson in this example.
+def solver(mf):
+    fci_1h = FCI["1h"](mf)
+    fci_1p = FCI["1p"](mf)
+
+    # Use MBLGF
+    nmom_max = 4
+    th = fci_1h.build_gf_moments(nmom_max)
+    tp = fci_1p.build_gf_moments(nmom_max)
+
+    norb = mf.mo_coeff.shape[-1]
+    nelec = mf.mol.nelec
+    civec= fci_1h.c_ci
+    dm1, dm2 = pyscf.fci.direct_spin0.make_rdm12(civec, norb, nelec)
+    results = dict(dm1=dm1, dm2=dm2, hole_moments=th, particle_moments=tp, converged=True)
+    return results
+
+natom = 10
+mol = pyscf.gto.Mole()
+mol.atom = ring("H", natom, 1.5)
+mol.basis = "sto-3g"
+mol.output = "pyscf.out"
+mol.verbose = 5
+mol.symmetry = True
+mol.build()
+
+# Hartree-Fock
+mf = pyscf.scf.RHF(mol)
+mf.kernel()
+
+# Vayesta options
+use_sym = True
+nfrag = 1 
+bath_opts = dict(bathtype="ewdmet", order=1, max_order=1)   
+
+# Run vayesta with user defined solver
+emb = vayesta.ewf.EWF(mf, solver="CALLBACK",  energy_functional='dmet', bath_options=bath_opts, solver_options=dict(callback=solver))
+emb.qpewdmet_scmf(proj=2, maxiter=10)
+# Set up fragments
+with emb.iao_fragmentation() as f:
+    if use_sym:
+        # Add rotational symmetry
+        with f.rotational_symmetry(order=natom//nfrag, axis=[0, 0, 1]):
+            f.add_atomic_fragment(range(nfrag))
+    else:
+        # Add all atoms as separate fragments 
+        f.add_all_atomic_fragments()
+emb.kernel()
+
+print("Hartree-Fock energy          : %s"%mf.e_tot)
+print("DMET energy                  : %s"%emb.get_dmet_energy(part_cumulant=False, approx_cumulant=False))
+print("DMET energy   (part-cumulant): %s"%emb.get_dmet_energy(part_cumulant=True, approx_cumulant=False))
+print("DMET energy (approx-cumulant): %s"%emb.get_dmet_energy(part_cumulant=True, approx_cumulant=True))
+

examples/ewf/molecules/65-callback-solver.py → examples/ewf/molecules/65-callback-solver-dms.py

@@ -8,6 +8,8 @@
 from vayesta.misc.molecules import ring
 
 # User defined FCI solver - takes pyscf mf as input and returns RDMs
+# The mf argment contains the hamiltonain in the orthonormal cluster basis
+# Pyscf or other solvers may be used to solve the cluster problem and may return RDMs, CISD amplitudes or CCSD amplitudes
 def solver(mf):
     h1e = mf.get_hcore()
     h2e = mf._eri
@@ -37,7 +39,7 @@ def solver(mf):
 bath_opts = dict(bathtype="dmet")   
 
 # Run vayesta with user defined solver
-emb = vayesta.ewf.EWF(mf, solver="CALLBACK",  energy_functional='dm', bath_options=bath_opts, solver_options=dict(callback=solver))
+emb = vayesta.ewf.EWF(mf, solver="CALLBACK",  energy_functional='dmet', bath_options=bath_opts, solver_options=dict(callback=solver))
 # Set up fragments
 with emb.iao_fragmentation() as f:
     if use_sym:

examples/ewf/molecules/66-callback-solver-amplitudes.py

@@ -0,0 +1,74 @@
+import numpy as np
+import pyscf
+import pyscf.gto
+import pyscf.scf
+import pyscf.fci
+import vayesta
+import vayesta.ewf
+from vayesta.misc.molecules import ring
+from vayesta.core.types.wf.t_to_c import t1_rhf, t2_rhf
+
+# User defined FCI solver - takes pyscf mf as input and returns RDMs
+# The mf argment contains the hamiltonain in the orthonormal cluster basis
+# Pyscf or other solvers may be used to solve the cluster problem and may return RDMs, CISD amplitudes or CCSD amplitudes
+def solver(mf):
+    ci = pyscf.ci.CISD(mf)
+    energy, civec = ci.kernel()
+    c0, c1, c2 = ci.cisdvec_to_amplitudes(civec)
+
+    # To use CI amplitudues use return the following line and set energy_functional='wf' to use the projected energy in the EWF arguments below
+    # return dict(c0=c0, c1=c1, c2=c2, converged=True, energy=ci.e_corr)
+
+    # Convert CISD amplitudes to CCSD amplitudes to be able to make use of the patitioned cumulant energy functional
+    t1 = t1_rhf(c1/c0) 
+    t2 = t2_rhf(t1, c2/c0)
+    return dict(t1=t1, t2=t2, l1=t1, l2=t2, converged=True, energy=ci.e_corr)
+
+natom = 10
+mol = pyscf.gto.Mole()
+mol.atom = ring("H", natom, 1.5)
+mol.basis = "sto-3g"
+mol.output = "pyscf.out"
+mol.verbose = 5
+mol.symmetry = True
+mol.build()
+
+# Hartree-Fock
+mf = pyscf.scf.RHF(mol)
+mf.kernel()
+
+# CISD
+cisd = pyscf.ci.CISD(mf)
+cisd.kernel()
+
+# CCSD
+ccsd = pyscf.cc.CCSD(mf)
+ccsd.kernel()
+
+# FCI
+fci = pyscf.fci.FCI(mf)
+fci.kernel()
+
+# Vayesta options
+use_sym = True
+nfrag = 1 
+bath_opts = dict(bathtype="dmet")   
+
+# Run vayesta with user defined solver
+emb = vayesta.ewf.EWF(mf, solver="CALLBACK",  energy_functional='dm-t2only', bath_options=bath_opts, solver_options=dict(callback=solver))
+# Set up fragments
+with emb.iao_fragmentation() as f:
+    if use_sym:
+        # Add rotational symmetry
+        with f.rotational_symmetry(order=natom//nfrag, axis=[0, 0, 1]):
+            f.add_atomic_fragment(range(nfrag))
+    else:
+        # Add all atoms as separate fragments 
+        f.add_all_atomic_fragments()
+emb.kernel()
+
+print("Hartree-Fock energy          : %s"%mf.e_tot)
+print("CISD Energy                  : %s"%cisd.e_tot)
+print("CCSD Energy                  : %s"%ccsd.e_tot)
+print("FCI  Energy                  : %s"%fci.e_tot)
+print("Emb. Partitioned Cumulant    : %s"%emb.e_tot)

examples/ewf/solids/67-callback-solver-density-fitting.py

@@ -0,0 +1,58 @@
+import numpy as np
+
+import pyscf
+import pyscf.pbc
+import pyscf.pbc.scf
+import pyscf.pbc.cc
+import pyscf.fci
+
+import vayesta
+import vayesta.ewf
+
+# User defined FCI solver - takes pyscf mf as input and returns RDMs
+# The mf argment contains the hamiltonain in the orthonormal cluster basis
+# Pyscf or other solvers may be used to solve the cluster problem and may return RDMs, CISD amplitudes or CCSD amplitudes
+def solver(mf):
+    h1e = mf.get_hcore()
+    # Need to convert cderis into standard 4-index tensor when using denisty fitting for the mean-field
+    cderi = mf.with_df._cderi
+    cderi = pyscf.lib.unpack_tril(cderi)
+    h2e = np.einsum('Lpq,Lrs->pqrs', cderi, cderi)
+    norb = mf.mo_coeff.shape[-1]
+    nelec = mf.mol.nelec
+    energy, civec = pyscf.fci.direct_spin0.kernel(h1e, h2e, norb, nelec, conv_tol=1.e-14)
+    dm1, dm2 = pyscf.fci.direct_spin0.make_rdm12(civec, norb, nelec)
+    results = dict(dm1=dm1, dm2=dm2, converged=True)
+    return results
+
+cell = pyscf.pbc.gto.Cell()
+cell.a = 3.0 * np.eye(3)
+cell.atom = "He 0 0 0"
+cell.basis = "cc-pvdz"
+#cell.exp_to_discard = 0.1
+cell.build()
+
+kmesh = [3, 3, 3]
+kpts = cell.make_kpts(kmesh)
+
+# --- Hartree-Fock
+kmf = pyscf.pbc.scf.KRHF(cell, kpts)
+kmf = kmf.rs_density_fit()
+kmf.kernel()
+
+# Vayesta options
+nfrag = 1 
+bath_opts = dict(bathtype="mp2", dmet_threshold=1e-15)   
+
+# Run vayesta with user defined solver
+emb = vayesta.ewf.EWF(kmf, solver="CALLBACK",  energy_functional='dmet', bath_options=bath_opts, solver_options=dict(callback=solver))
+# Set up fragments
+with emb.iao_fragmentation() as f:
+        f.add_all_atomic_fragments()
+emb.kernel()
+
+print("Hartree-Fock energy          : %s"%kmf.e_tot)
+print("DMET energy                  : %s"%emb.get_dmet_energy(part_cumulant=False, approx_cumulant=False))
+print("DMET energy   (part-cumulant): %s"%emb.get_dmet_energy(part_cumulant=True, approx_cumulant=False))
+print("DMET energy (approx-cumulant): %s"%emb.get_dmet_energy(part_cumulant=True, approx_cumulant=True))
+

pyproject.toml

@@ -41,7 +41,7 @@ dependencies = [
     "scipy>=1.1.0,<=1.10.0",
     "h5py>=2.7",
     "cvxpy>=1.1",
-    "pyscf @ git+https://github.com/pyscf/pyscf@e90d8342e29446365f8570682af4a65fe16be27c",
+    "pyscf @ git+https://github.com/pyscf/pyscf@master",
 ]
 dynamic = ["version"]
 

vayesta/core/qemb/fragment.py

@@ -1160,7 +1160,7 @@ def check_solver(self, solver):
             is_eb = "crpa_full" in self.opts.screening
         else:
             is_eb = False
-        check_solver_config(is_uhf, is_eb, solver, self.log)
+        check_solver_config(solver, is_uhf, is_eb, self.log)
 
     def get_solver(self, solver=None):
         if solver is None:

vayesta/core/qemb/qemb.py

@@ -342,7 +342,7 @@ def init_mf(self, mf):
         mf = copy.copy(mf)
         self.log.debugv("type(mf)= %r", type(mf))
         # If the mean-field has k-points, automatically fold to the supercell:
-        if getattr(mf, "kpts", None) is not None:
+        if isinstance(mf, pyscf.pbc.scf.khf.KSCF):
             with log_time(self.log.timing, "Time for k->G folding of MOs: %s"):
                 mf = fold_scf(mf)
         if isinstance(mf, FoldedSCF):
@@ -1263,7 +1263,7 @@ def make_rdm1_demo(self, *args, **kwargs):
     def make_rdm2_demo(self, *args, **kwargs):
         return make_rdm2_demo_rhf(self, *args, **kwargs)
 
-    def get_dmet_elec_energy(self, part_cumulant=True, approx_cumulant=True):
+    def get_dmet_elec_energy(self, part_cumulant=True, approx_cumulant=True, mpi_target=None):
         """Calculate electronic DMET energy via democratically partitioned density-matrices.
 
         Parameters
@@ -1275,6 +1275,10 @@ def get_dmet_elec_energy(self, part_cumulant=True, approx_cumulant=True):
         approx_cumulant: bool, optional
             If True, the approximate cumulant, containing (delta 1-DM)-squared terms, is partitioned,
             instead of the true cumulant, if `part_cumulant=True`. Default: True.
+        mpi_target: int or None, optional
+            If set to an integer, the result will only be available at the specified MPI rank.
+            If set to None, an MPI allreduce will be performed and the result will be available
+            at all MPI ranks. Default: None. 
 
         Returns
         -------
@@ -1286,7 +1290,8 @@ def get_dmet_elec_energy(self, part_cumulant=True, approx_cumulant=True):
             wx = x.symmetry_factor
             e_dmet += wx * x.get_fragment_dmet_energy(part_cumulant=part_cumulant, approx_cumulant=approx_cumulant)
         if mpi:
-            mpi.world.allreduce(e_dmet)
+            e_dmet = mpi.nreduce(e_dmet, target=mpi_target, logfunc=self.log.timingv)
+
         if part_cumulant:
             dm1 = self.make_rdm1_demo(ao_basis=True)
             if not approx_cumulant:
@@ -1750,4 +1755,4 @@ def check_solver(self, solver):
             is_eb = "crpa_full" in self.opts.screening
         else:
             is_eb = False
-        check_solver_config(is_uhf, is_eb, solver, self.log)
+        check_solver_config(solver, is_uhf, is_eb, self.log)