Draft of PNO Lambda CCD and CCSD

pycc/data/molecules.py

@@ -217,7 +217,37 @@
 H 3 0.75 2 90.0 1 60.0
 symmetry c1
 """
+# H2_4 chain
+h2_4 = """
+H 0.000000 0.000000 0.000000
+H 0.750000 0.000000 0.000000
+H 0.000000 1.500000 0.000000
+H 0.375000 1.500000 -0.649520
+H 0.000000 3.000000 0.000000
+H -0.375000 3.000000 -0.649520
+H 0.000000 4.500000 -0.000000
+H -0.750000 4.500000 -0.000000
+symmetry c1
+"""
 
+# H2_7 chain
+h2_7 = """
+H 0.000000 0.000000 0.000000
+H 0.750000 0.000000 0.000000
+H 0.000000 1.500000 0.000000
+H 0.375000 1.500000 -0.649520
+H 0.000000 3.000000 0.000000
+H -0.375000 3.000000 -0.649520
+H 0.000000 4.500000 -0.000000
+H -0.750000 4.500000 -0.000000
+H 0.000000 6.000000 -0.000000
+H -0.375000 6.000000 0.649520
+H 0.000000 7.500000 -0.000000
+H 0.375000 7.500000 -0.649520
+H 0.000000 9.000000 -0.000000
+H 0.750000 9.000000 0.000000
+symmetry c1
+"""
 # (S)-1,3-dimethylallene
 sdma = """
 C  0.000000  0.000000  0.414669
@@ -282,6 +312,8 @@
 moldict["uracil"] = uracil
 moldict["benzene"] = benzene
 moldict["(H2)_2"] = h2_2
+moldict["(H2)_4"] = h2_4 
+moldict["(H2)_7"] = h2_7
 moldict["(S)-dimethylallene"] = sdma
 moldict["(S)-2-chloropropionitrile"] = s2cpn
 moldict["(R)-methylthiirane"] = rmt

pycc/lccwfn.py

@@ -112,6 +112,7 @@ def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8,start_diis
         #ldiis = helper_ldiis(self.t1, self.t2, max_diis)
 
         elcc = self.lcc_energy(self.Local.Fov,self.Local.Loovv,self.t1, self.t2)
+
         print("CC Iter %3d: lCC Ecorr = %.15f dE = % .5E MP2" % (0,elcc,-elcc))
 
         for niter in range(1, maxiter+1):
@@ -125,11 +126,10 @@ def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8,start_diis
             rms_t2 = 0
 
             for i in range(self.no):
+
                 ii = i*self.no + i
 
-                #need to change to reshape
-                for a in range(self.Local.dim[ii]):
-                    self.t1[i][a] += r1[i][a]/(self.H.F[i,i] - self.Local.eps[ii][a])
+                self.t1[i] -= r1[i]/(self.Local.eps[ii].reshape(-1,) - self.H.F[i,i])
 
                 rms_t1 += contract('Z,Z->',r1[i],r1[i])
 
@@ -141,7 +141,7 @@ def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8,start_diis
 
                     rms_t2 += contract('ZY,ZY->',r2[ij],r2[ij])
 
-            rms = np.sqrt(rms_t2)
+            rms = np.sqrt(rms_t1 + rms_t2)
             elcc = self.lcc_energy(self.Local.Fov,self.Local.Loovv,self.t1, self.t2)
             ediff = elcc - elcc_last
             print("lCC Iter %3d: lCC Ecorr = %.15f  dE = % .5E  rms = % .5E" % (niter, elcc, ediff, rms))
@@ -155,7 +155,7 @@ def solve_lcc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8,start_diis
                 self.elcc = elcc
                 #print(Timer.timers)
                 return elcc
-
+      
             #ldiis.add_error_vector(self.t1,self.t2)
             #if niter >= start_diis:
                 #self.t1, self.t2 = ldiis.extrapolate(self.t1, self.t2)
@@ -203,7 +203,7 @@ def build_Fae(self, Fae_ij, L, Fvv, Fov, Sijmm, Sijmn, t1, t2):
         o = self.o
         v = self.v
         QL = self.QL
-
+        
         if self.model == 'CCD':
             for ij in range(self.no*self.no):
                 i = ij // self.no
@@ -414,6 +414,7 @@ def build_Wmbej(self, Wmbej_ijim, ERI, L, ERIoovo, Sijnn, Sijnj, Sijjn, t1, t2):
                 i = ij // self.no
                 j = ij % self.no
                 jj = j*self.no + j
+
                 for m in range(self.no):
                     im = i*self.no + m
 
@@ -601,7 +602,6 @@ def r_T1(self, r1_ii, Fov , ERI, L, Loovo, Sijmm, Sijim, Sijmn, t1, t2, Fae, Fmi
                 for mn in range(self.no*self.no):
                     m = mn // self.no
                     n = mn % self.no
-                    imn = i*(self.no**2) + mn
                     iimn =ii*(self.no**2) + mn 
 
                     tmp4 = Sijmn[iimn] @ t2[mn]
@@ -737,7 +737,7 @@ def r_T2(self,r2_ij, ERI, ERIoovv, ERIvvvv, ERIovoo, Sijmm, Sijim, Sijmj, Sijnn,
                         r2 += 0.5 * contract('a,b->ab',tmp2_0, tmp13) * Wmnij[m,n,i,j]
 
                 nr2.append(r2) 
-
+    
         for i in range(self.no):
             for j in range(self.no):
                 ij = i*self.no + j
@@ -755,7 +755,7 @@ def lcc_energy(self, Fov, Loovv, t1, t2):
         ecc_ii = 0
         ecc_ij = 0
         ecc = 0
-
+        
         if self.model == 'CCD':
             for i in range(self.no):
                 for j in range(self.no):

pycc/local.py

@@ -76,7 +76,6 @@ def __init__(self, local, C, nfzc, no, nv, H, cutoff, it2_opt,
         self.lindep_cut = lindep_cut
         self.e_conv = e_conv
         self.r_conv = r_conv
-
         self._build()
 
     def _build(self):
@@ -722,7 +721,7 @@ def _local_MP2_loop(self,local,S):
             ediff = emp2 - elast
 
             print("MP2 Iter %3d: MP2 Ecorr = %.15f dE = % .5E rmsd = % .5E" % (niter, emp2, ediff, rmsd))
-    
+
     def filter_t2amps(self,r2):
         no = self.no
         nv = self.nv
@@ -826,12 +825,15 @@ def trans_integrals(self, o, v):
         ERIoovo = []
         ERIooov = []
         ERIovvv = []
+        ERIvovv = []
         ERIvvvv = []
         ERIoovv = []
         ERIovvo = []
         ERIvvvo = []
         ERIovov = []
         ERIovoo = []
+        ERIvoov = []
+        ERIvovo = []
         Loovv = []
         Lovvv = []
         Looov = []
@@ -858,47 +860,55 @@ def trans_integrals(self, o, v):
             ERIoovv.append(contract('ijAb,bB->ijAB',tmp,QL[ij]))
 
             ERIovvo.append(ERIoovv[ij].swapaxes(1,3))
+
+            ERIvoov.append(ERIovvo[ij].swapaxes(0,1).swapaxes(2,3))
 
             tmp1 = contract('iajb,aA->iAjb',self.H.ERI[o,v,o,v], QL[ij])
             ERIovov.append(contract('iAjb,bB->iAjB',tmp1, QL[ij]))
 
+            ERIvovo.append(ERIovov[ij].swapaxes(0,1).swapaxes(2,3))
+
             tmp2 = contract('iabc,aA->iAbc',self.H.ERI[o,v,v,v], QL[ij])
             tmp2 = contract('iAbc,bB->iABc',tmp2, QL[ij])
             ERIovvv.append(contract('iABc,cC->iABC',tmp2, QL[ij]))
 
-            tmp3 = ERIovvv[ij].swapaxes(0,1).swapaxes(2,3)
-            ERIvvvo.append(tmp3.swapaxes(1,3))
+            ERIvovv.append(ERIovvv[ij].swapaxes(0,1).swapaxes(2,3))
+
+            ERIvvvo.append(ERIvovv[ij].swapaxes(1,3))
 
-            tmp4 = contract('abcd,aA->Abcd',self.H.ERI[v,v,v,v], QL[ij])
-            tmp4 = contract('Abcd,bB->ABcd',tmp4, QL[ij])
-            tmp4 = contract('ABcd,cC->ABCd',tmp4, QL[ij])
-            ERIvvvv.append(contract('ABCd,dD->ABCD',tmp4, QL[ij]))
+            tmp3 = contract('abcd,aA->Abcd',self.H.ERI[v,v,v,v], QL[ij])
+            tmp3 = contract('Abcd,bB->ABcd',tmp3, QL[ij])
+            tmp3 = contract('ABcd,cC->ABCd',tmp3, QL[ij])
+            ERIvvvv.append(contract('ABCd,dD->ABCD',tmp3, QL[ij]))
 
             Loovo.append(contract('ijak,aA->ijAk', self.H.L[o,o,v,o],QL[ij]))
 
             Looov.append(Loovo[ij].swapaxes(0,1).swapaxes(2,3))
 
-            tmp5 = contract('ijab,aA->ijAb',self.H.L[o,o,v,v], QL[ij])
-            Loovv.append(contract('ijAb,bB->ijAB',tmp5,QL[ij]))
+            tmp4 = contract('ijab,aA->ijAb',self.H.L[o,o,v,v], QL[ij])
+            Loovv.append(contract('ijAb,bB->ijAB',tmp4,QL[ij]))
 
             Lovvo.append(Loovv[ij].swapaxes(1,3))
 
-            tmp6 = contract('iabc,aA->iAbc',self.H.L[o,v,v,v], QL[ij])
-            tmp6 = contract('iAbc,bB->iABc',tmp6, QL[ij])
-            Lovvv.append(contract('iABc,cC->iABC',tmp6, QL[ij]))
+            tmp5 = contract('iabc,aA->iAbc',self.H.L[o,v,v,v], QL[ij])
+            tmp5 = contract('iAbc,bB->iABc',tmp5, QL[ij])
+            Lovvv.append(contract('iABc,cC->iABC',tmp5, QL[ij]))
 
             self.QL = QL
             self.Fov = Fov
             self.Fvv = Fvv
             self.ERIoovo = ERIoovo
             self.ERIooov = ERIooov
             self.ERIovvv = ERIovvv
+            self.ERIvovv = ERIvovv
             self.ERIvvvv = ERIvvvv
             self.ERIoovv = ERIoovv
             self.ERIovvo = ERIovvo
             self.ERIvvvo = ERIvvvo
             self.ERIovov = ERIovov
             self.ERIovoo = ERIovoo
+            self.ERIvoov = ERIvoov
+            self.ERIvovo = ERIvovo
             self.Loovv = Loovv
             self.Lovvv = Lovvv
             self.Looov = Looov
@@ -920,11 +930,17 @@ def overlaps(self, QL):
         To do
         -----
         Compare the timings for the use of stored overlap terms versus "on the fly" overlap terms 
+
+        Redundant overlap terms generated - another pull request will be made to update this with changes to
+        other object class that utilizes this (lccwfn object) 
         """
         no = self.no 
 
+        Sijii = []
+        Sijjj = []
         Sijmm = []
         Sijim = []
+        Sijmi = []
         Sijmj = []
         Sijnn = []
         Sijin = []
@@ -933,17 +949,26 @@ def overlaps(self, QL):
         Sijmn = []
 
         for i in range(no):
+            ii = i*no + i
             for j in range(no):
                 ij = i*no + j
+                jj = j*no + j
+                ji = j*no + i 
+
+                Sijii.append(QL[ij].T @ QL[ii])    
+                Sijjj.append(QL[ij].T @ QL[jj])
+
                 for m in range(no):
                     mm = m*no + m
                     im = i*no + m
                     mj = m*no + j
+                    mi = m*no + i
 
                     Sijmm.append(QL[ij].T @ QL[mm])
                     Sijim.append(QL[ij].T @ QL[im]) 
                     Sijmj.append(QL[ij].T @ QL[mj])  
-
+                    Sijmi.append(QL[ij].T @ QL[mi])
+
                 for n in range(no):
                     nn = n*no + n
                     _in = i*no + n 
@@ -958,9 +983,12 @@ def overlaps(self, QL):
                 for mn in range(no*no):
                     Sijmn.append(QL[ij].T @ QL[mn])
 
+        self.Sijii = Sijii
+        self.Sijjj = Sijjj
         self.Sijmj = Sijmj
         self.Sijmm = Sijmm 
         self.Sijim = Sijim
+        self.Sijmi = Sijmi
         self.Sijnn = Sijnn
         self.Sijin = Sijin
         self.Sijnj = Sijnj

pycc/tests/test_013_pnocc.py

@@ -29,7 +29,7 @@ def test_pno_ccsd():
     r_conv = 1e-12
     max_diis = 8
 
-    ccsd = pycc.ccwfn(rhf_wfn, local='PNO', local_cutoff=1e-5, it2_opt=False)
+    ccsd = pycc.ccwfn(rhf_wfn, local='PNO', local_cutoff=1e-5, it2_opt=False, filter=True)
     eccsd = ccsd.solve_cc(e_conv, r_conv, maxiter)
 
     hbar = pycc.cchbar(ccsd)
@@ -63,7 +63,7 @@ def test_pno_ccsd_opt():
     r_conv = 1e-12
     max_diis = 8
 
-    ccsd = pycc.ccwfn(rhf_wfn, local='PNO', local_cutoff=1e-5)
+    ccsd = pycc.ccwfn(rhf_wfn, local='PNO', local_cutoff=1e-5, filter=True)
     eccsd = ccsd.solve_cc(e_conv, r_conv, maxiter)
 
     hbar = pycc.cchbar(ccsd)

pycc/tests/test_018_paocc.py

@@ -37,7 +37,7 @@ def test_pao_H8():
     r_conv = 1e-12
     max_diis = 8
 
-    ccsd = pycc.ccwfn(rhf_wfn, local='PAO', local_cutoff=2e-2)
+    ccsd = pycc.ccwfn(rhf_wfn, local='PAO', local_cutoff=2e-2, filter=True)
 
     eccsd = ccsd.solve_cc(e_conv, r_conv, maxiter, max_diis)
 
@@ -68,7 +68,7 @@ def test_pao_h2o():
     r_conv = 1e-7
     max_diis = 8
 
-    ccsd = pycc.ccwfn(rhf_wfn, local='PAO', local_cutoff=2e-2)
+    ccsd = pycc.ccwfn(rhf_wfn, local='PAO', local_cutoff=2e-2, filter=True)
 
     eccsd = ccsd.solve_cc(e_conv, r_conv, maxiter, max_diis)
 
@@ -99,7 +99,7 @@ def test_pao_h2o_frzc():
     r_conv = 1e-7
     max_diis = 8
 
-    ccsd = pycc.ccwfn(rhf_wfn, local='PAO', local_cutoff=2e-2)
+    ccsd = pycc.ccwfn(rhf_wfn, local='PAO', local_cutoff=2e-2, filter=True)
 
     eccsd = ccsd.solve_cc(e_conv, r_conv, maxiter, max_diis)
 

pycc/tests/test_019_localrt.py

@@ -53,7 +53,7 @@ def test_rtpno(datadir):
     r_conv = 1e-13
     max_diis = 8
 
-    cc = pycc.ccwfn(rhf_wfn, local='PNO', local_cutoff=1e-5)
+    cc = pycc.ccwfn(rhf_wfn, local='PNO', local_cutoff=1e-5, filter=True)
     ecc = cc.solve_cc(e_conv, r_conv, maxiter, max_diis)
     hbar = pycc.cchbar(cc)
     cclambda = pycc.cclambda(cc, hbar)
@@ -109,7 +109,7 @@ def test_rtpao(datadir):
     r_conv = 1e-13
     max_diis = 8
 
-    cc = pycc.ccwfn(rhf_wfn, local='PAO', local_cutoff=1e-2)
+    cc = pycc.ccwfn(rhf_wfn, local='PAO', local_cutoff=1e-2, filter=True)
     ecc = cc.solve_cc(e_conv, r_conv, maxiter, max_diis)
     hbar = pycc.cchbar(cc)
     cclambda = pycc.cclambda(cc, hbar)

pycc/tests/test_028_pnoppcc.py

@@ -30,7 +30,7 @@ def test_pnopp_ccsd():
     r_conv = 1e-12
     max_diis = 8
 
-    ccsd = pycc.ccwfn(rhf_wfn, local='PNO++', local_cutoff=1e-7, it2_opt=False)
+    ccsd = pycc.ccwfn(rhf_wfn, local='PNO++', local_cutoff=1e-7, it2_opt=False, filter=True)
     eccsd = ccsd.solve_cc(e_conv, r_conv, maxiter)
 
     hbar = pycc.cchbar(ccsd)
@@ -64,7 +64,7 @@ def test_pnopp_ccsd_opt():
     r_conv = 1e-12
     max_diis = 8
 
-    ccsd = pycc.ccwfn(rhf_wfn, local='PNO++', local_cutoff=1e-9)
+    ccsd = pycc.ccwfn(rhf_wfn, local='PNO++', local_cutoff=1e-9, filter=True)
     eccsd = ccsd.solve_cc(e_conv, r_conv, maxiter)
 
     hbar = pycc.cchbar(ccsd)