Merge pull request #168 from BoothGroup/callback_solver

Callback solver - CISD amplitudes

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)

vayesta/ewf/ewf.py

@@ -179,13 +179,13 @@ def kernel(self):
             self.log.error("Some fragments did not converge!")
         self.converged = conv
 
-        if self.solver.lower() == "callback":
-            self.log.info("Total wall time:  %s", time_string(timer() - t_start))
-            return
         # --- Evaluate correlation energy and log information
-        self.e_corr = self.get_e_corr()
-
-
+        try:
+            self.e_corr = self.get_e_corr()
+        except AttributeError as e:
+            self.log.error("Could not calculate correlation energy")
+            self.e_corr = np.nan
+
         self.log.output("E(MF)=   %s", energy_string(self.e_mf))
         self.log.output("E(corr)= %s", energy_string(self.e_corr))
         self.log.output("E(tot)=  %s", energy_string(self.e_tot))

vayesta/tests/solver/test_callback.py

@@ -49,12 +49,14 @@ class TestSolvers(TestCase):
     def _test(self, key):
         mf = getattr(getattr(testsystems, key[0]), key[1])()
         bath_opts = dict(bathtype="dmet")
-        emb = vayesta.ewf.EWF(mf, solver=key[2],  energy_functional='dm', bath_options=bath_opts)
+        emb = vayesta.ewf.EWF(mf, solver=key[2],  energy_functional='wf', bath_options=bath_opts)
         emb.kernel()
 
         for dm in [True, False]:
+            if key[2] == 'CISD' and dm:
+                return
             callback = lambda mf: callbacks[key[2]](mf, dm=dm)
-            emb_callback = vayesta.ewf.EWF(mf, solver="CALLBACK",  energy_functional='dm', bath_options=bath_opts, solver_options=dict(callback=callback))
+            emb_callback = vayesta.ewf.EWF(mf, solver="CALLBACK",  energy_functional='wf', bath_options=bath_opts, solver_options=dict(callback=callback))
             emb_callback.kernel()
 
             for a, b in [(False, False), (True, False), (True, True)]:
@@ -67,8 +69,8 @@ def test_rccsd_water(self):
     def test_fci_h6(self):
         self._test(("h6_sto6g", "rhf", "FCI"))
 
-    # def test_cisd_lih(self):
-    #     self._test(("lih_ccpvdz", "rhf", "CISD"))
+    def test_cisd_lih(self):
+        self._test(("lih_ccpvdz", "rhf", "CISD"))
 
 if __name__ == "__main__":
     print("Running %s" % __file__)