Adding incomplete draft of EOM-CC code.

pycc/__init__.py

@@ -13,8 +13,9 @@
 from .ccresponse import ccresponse
 from .ccresponse import pertbar
 from pycc.rt.rtcc import rtcc
+from .cceom import cceom
 
-__all__ = ['ccwfn', 'cchbar', 'cclambda', 'ccdensity', 'ccresponse', 'pertbar', 'rtcc']
+__all__ = ['ccwfn', 'cchbar', 'cclambda', 'ccdensity', 'ccresponse', 'pertbar', 'rtcc', 'cceom']
 
 # Handle versioneer
 from ._version import get_versions

pycc/cceom.py

@@ -0,0 +1,137 @@
+"""
+cceom.py: Computes excited-state CC wave functions and energies
+"""
+
+if __name__ == "__main__":
+    raise Exception("This file cannot be invoked on its own.")
+
+import time
+import numpy as np
+
+
+class cceom(object):
+    """
+    An Equation-of-Motion Coupled Cluster Object.
+
+    """
+
+    def __init__(self, cchbar):
+        """
+        Parameters
+        ----------
+        cchbar : PyCC cchbar object
+
+        Returns
+        -------
+        """
+
+        self.cchbar = cchbar
+
+    def solve_eom(self, nstates=1, e_conv=1e-5, r_conv=1e-5, maxiter=100):
+        """
+
+        """
+
+        time_init = time.time()
+
+        hbar = self.cchbar
+        o = hbar.o
+        v = hbar.v
+        no = hbar.no
+        nv = hbar.nv
+        H = hbar.ccwfn.H
+        contract = hbar.contract
+
+        L = nstates * 3 # size of guess space (varies)
+        maxL = L * 10 # max size of subspace (fixed)
+
+        # Initialize guess vectors
+        _, C1 = self.cis(L, o, v, H, contract)
+        C2 = np.zeros((L, no, no, nv, nv))
+
+        # Compute sigma vectors
+        s1 = np.zeros_like(C1)
+        s2 = np.zeros_like(C2)
+        for state in range(L):
+            s1[state] = self.s1(hbar, C1[state], C2[state])
+            s2[state] = self.s2(hbar, C1[state], C2[state])
+
+        # Build and diagonalize subspace Hamiltonian
+        G = np.zeros((L,L))
+
+        # Build and orthonormalize correction vectors
+
+        # 
+
+        print("\nCCEOM converged in %.3f seconds." % (time.time() - time_init))
+
+
+
+    def cis(self, nstates, o, v, H, contract):
+        """
+        Compute CIS excited states as guess vectors
+        """
+        no = o.stop - o.start
+        nv = v.stop - v.start
+        F = H.F
+        L = H.L
+
+        CIS_H = L[v,o,o,v].swapaxes(0,1).swapaxes(0,2).copy()
+        CIS_H += contract('ab,ij->iajb', F[v,v], np.eye(no))
+        CIS_H -= contract('ij,ab->iajb', F[o,o], np.eye(nv))
+
+        eps, c = np.linalg.eigh(np.reshape(CIS_H, (no*nv,no*nv)))
+
+        c = np.reshape(c.T[slice(0,nstates),:], (nstates, no, nv)).copy()
+
+        return eps, c
+
+
+    def s1(self, hbar, C1, C2):
+
+        contract = hbar.contract
+
+        s1 = contract('ie,ae->ia', C1, hbar.Hvv)
+        s1 -= contract('mi,ma->ia', hbar.Hoo, C1)
+        s1 += 2.0 * contract('maei,me->ia', hbar.Hovvo, C1)
+        s1 -= contract('maie,me->ia', hbar.Hovov, C1)
+        s1 += 2.0 * contract('miea,me->ia', C2, hbar.Hov)
+        s1 -= contract('imea,me->ia', C2, hbar.Hov)
+        s1 += 2.0 * contract('imef,amef->ia', C2, hbar.Hvovv)
+        s1 -= contract('imef,amfe->ia', C2, hbar.Hvovv)
+        s1 -= 2.0 * contract('mnie,mnae->ia', hbar.Hooov, C2)
+        s1 += contract('nmie,mnae->ia', hbar.Hooov, C2)
+
+        return s1
+
+
+    def s2(self, hbar, C1, C2):
+
+        contract = hbar.contract
+        L = hbar.ccwfn.H.L
+        t2 = hbar.ccwfn.t2
+        o = hbar.ccwfn.o
+        v = hbar.ccwfn.v
+
+        Zvv = 2.0 * contract('amef,mf->ae', hbar.Hvovv, C1)
+        Zvv -= contract('amfe,mf->ae', hbar.Hvovv, C1)
+        Zvv -= contract('nmaf,nmef->ae', C2, L[o,o,v,v])
+
+        Zoo = -2.0 * contract('mnie,ne->mi', hbar.Hooov, C1)
+        Zoo += contract('nmie,ne->mi', hbar.Hooov, C1)
+        Zoo -= contract('mnef,inef->mi', L[o,o,v,v], C2)
+
+        s2 = contract('ie,abej->ijab', C1, hbar.Hvvvo)
+        s2 -= contract('mbij,ma->ijab', hbar.Hovoo, C1)
+        s2 += contract('ijeb,ae->ijab', t2, Zvv)
+        s2 += contract('mi,mjab->ijab', Zoo, t2)
+        s2 += contract('ijeb,ae->ijab', C2, hbar.Hvv)
+        s2 -= contract('mi,mjab->ijab', hbar.Hoo, C2)
+        s2 += 0.5 * contract('mnij,mnab->ijab', hbar.Hoooo, C2)
+        s2 += 0.5 * contract('ijef,abef->ijab', C2, hbar.Hvvvv)
+        s2 -= contract('imeb,maje->ijab', C2, hbar.Hovov)
+        s2 -= contract('imea,mbej->ijab', C2, hbar.Hovvo)
+        s2 += 2.0 * contract('miea,mbej->ijab', C2, hbar.Hovvo)
+        s2 -= contract('miea,mbje->ijab', C2, hbar.Hovov)
+
+        return s2 + s2.swapaxes(0,1).swapaxes(2,3)

pycc/cchbar.py

@@ -56,16 +56,17 @@ def __init__(self, ccwfn):
         time_init = time.time()
 
         self.ccwfn = ccwfn
-
         self.contract = self.ccwfn.contract
 
-        o = ccwfn.o
-        v = ccwfn.v
         F = ccwfn.H.F
         ERI = ccwfn.H.ERI
         L = ccwfn.H.L
         t1 = ccwfn.t1
         t2 = ccwfn.t2
+        o = self.o = ccwfn.o
+        v = self.v = ccwfn.v
+        self.no = ccwfn.no
+        self.nv = ccwfn.nv
 
         self.Hov = self.build_Hov(o, v, F, L, t1)
         self.Hvv = self.build_Hvv(o, v, F, L, t1, t2)
@@ -79,20 +80,7 @@ def __init__(self, ccwfn):
         self.Hvvvo = self.build_Hvvvo(o, v, ERI, L, self.Hov, self.Hvvvv, t1, t2)
         self.Hovoo = self.build_Hovoo(o, v, ERI, L, self.Hov, self.Hoooo, t1, t2)
 
-        if isinstance(t1, torch.Tensor):
-            print("Hov norm = %20.15f" % torch.linalg.norm(self.Hov))
-            print("Hvv norm = %20.15f" % torch.linalg.norm(self.Hvv))
-            print("Hoo norm = %20.15f" % torch.linalg.norm(self.Hoo))
-            print("Hoooo norm = %20.15f" % torch.linalg.norm(self.Hoooo))
-            print("Hvvvv norm = %20.15f" % torch.linalg.norm(self.Hvvvv))
-        else:
-            print("Hov norm = %20.15f" % np.linalg.norm(self.Hov))
-            print("Hvv norm = %20.15f" % np.linalg.norm(self.Hvv))
-            print("Hoo norm = %20.15f" % np.linalg.norm(self.Hoo))
-            print("Hoooo norm = %20.15f" % np.linalg.norm(self.Hoooo))
-            print("Hvvvv norm = %20.15f" % np.linalg.norm(self.Hvvvv))
-
-        print("\nHBAR constructed in %.3f seconds.\n" % (time.time() - time_init))
+        print("\nHBAR constructed in %.3f seconds." % (time.time() - time_init))
 
     """
     For GPU implementation:

pycc/ccwfn.py

@@ -223,7 +223,7 @@ def __init__(self, scf_wfn, **kwargs):
                 self.H.ERI = torch.tensor(self.H.ERI, dtype=torch.complex64, device=self.device0)
                 self.H.L = torch.tensor(self.H.L, dtype=torch.complex64, device=self.device0)
 
-        print("CC object initialized in %.3f seconds." % (time.time() - time_init))
+        print("CCWFN object initialized in %.3f seconds." % (time.time() - time_init))
 
     def solve_cc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_diis=1):
         """
@@ -294,15 +294,15 @@ def solve_cc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_diis
             # check for convergence
             if isinstance(self.t1, torch.Tensor):
                 if ((torch.abs(ediff) < e_conv) and torch.abs(rms) < r_conv):
-                    print("\nCC has converged in %.3f seconds.\n" % (time.time() - ccsd_tstart))
+                    print("\nCCWFN converged in %.3f seconds.\n" % (time.time() - ccsd_tstart))
                     print("E(REF)  = %20.15f" % self.eref)
                     print("E(%s) = %20.15f" % (self.model, ecc))
                     print("E(TOT)  = %20.15f" % (ecc + self.eref))
                     self.ecc = ecc
                     return ecc
             else:
                 if ((abs(ediff) < e_conv) and abs(rms) < r_conv):
-                    print("\nCC has converged in %.3f seconds.\n" % (time.time() - ccsd_tstart))
+                    print("\nCCWFN converged in %.3f seconds.\n" % (time.time() - ccsd_tstart))
                     print("E(REF)  = %20.15f" % self.eref)
                     if (self.model == 'CCSD(T)'):
                         et = t_tjl(self)

pycc/tests/test_034_eomccsd.py

@@ -0,0 +1,41 @@
+"""
+Test EOM-CCSD solver.
+"""
+
+# Import package, test suite, and other packages as needed
+import psi4
+import pycc
+import pytest
+from ..data.molecules import *
+import numpy as np
+
+# H2O/cc-pVDZ
+def test_eomccsd_h2o():
+    # Psi4 Setup
+    psi4.set_memory('2 GB')
+    psi4.core.set_output_file('output.dat', False)
+    psi4.set_options({'basis': 'STO-3G',
+                      'scf_type': 'pk',
+                      'mp2_type': 'conv',
+                      'freeze_core': 'false',
+                      'e_convergence': 1e-12,
+                      'd_convergence': 1e-12,
+                      'r_convergence': 1e-12,
+                      'diis': 1})
+    mol = psi4.geometry(moldict["H2O_Teach"])
+    rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True)
+
+    maxiter = 75
+    e_conv = 1e-12
+    r_conv = 1e-12
+    cc = pycc.ccwfn(rhf_wfn)
+    ecc = cc.solve_cc(e_conv,r_conv,maxiter)
+#    epsi4 = -0.227888246840310
+#    ecfour = -0.2278882468404231
+#    assert (abs(epsi4 - ecc) < 1e-11)
+
+    hbar = pycc.cchbar(cc)
+
+    eom = pycc.cceom(hbar)
+
+    eom.solve_eom(nstates=3)