Merge branch 'main' into eomcc

.github/workflows/CI.yaml

@@ -53,7 +53,7 @@ jobs:
       run: |
         conda info
         conda list --show-channel-urls
-        conda install psi4 python=${{ matrix.python-version }} -c psi4/label/dev
+        conda install psi4 python=${{ matrix.python-version }} -c conda-forge 
         ls -l $CONDA
 
     - name: Test Psi4 Python Loading

devtools/conda-envs/psi.yaml

@@ -1,6 +1,5 @@
 name: p4dev
 channels:
-  - psi4/label/dev
   - conda-forge
 dependencies:
   - psi4

pycc/ccdensity.py

@@ -66,7 +66,7 @@ def __init__(self, ccwfn, cclambda, onlyone=False):
         self.ccwfn = ccwfn
         self.cclambda = cclambda
         self.contract = self.ccwfn.contract
-       
+
         o = ccwfn.o
         v = ccwfn.v
         no = ccwfn.no
@@ -83,7 +83,7 @@ def __init__(self, ccwfn, cclambda, onlyone=False):
         self.Dvo = self.build_Dvo(l1)
         self.Dvv = self.build_Dvv(t1, t2, l1, l2)
         self.Doo = self.build_Doo(t1, t2, l1, l2)
-        
+
         self.onlyone = onlyone
 
         if onlyone is False:
@@ -93,7 +93,7 @@ def __init__(self, ccwfn, cclambda, onlyone=False):
             self.Dvvvo = self.build_Dvvvo(t1, t2, l1, l2)
             self.Dovov = self.build_Dovov(t1, t2, l1, l2)
             self.Doovv = self.build_Doovv(t1, t2, l1, l2)
-        
+
         print("\nCCDENSITY constructed in %.3f seconds.\n" % (time.time() - time_init))
 
     def compute_energy(self):
@@ -114,9 +114,10 @@ def compute_energy(self):
         v = self.ccwfn.v
         F = self.ccwfn.H.F
         ERI = self.ccwfn.H.ERI
-             
+
         contract = self.contract 
 
+        # We assume here that the Brillouin condition holds
         oo_energy = contract('ij,ij->', F[o,o], self.Doo)
         vv_energy = contract('ab,ab->', F[v,v], self.Dvv)
         eone = oo_energy + vv_energy
@@ -134,8 +135,8 @@ def compute_energy(self):
             oovv_energy = 0.5 * contract('ijab,ijab->', ERI[o,o,v,v], self.Doovv)
             etwo = oooo_energy + vvvv_energy + ooov_energy + vvvo_energy + ovov_energy + oovv_energy
 
-            print("OOOV Energy = %20.15f" % oooo_energy)
-            print("OOOV Energy = %20.15f" % vvvv_energy)
+            print("OOOO Energy = %20.15f" % oooo_energy)
+            print("VVVV Energy = %20.15f" % vvvv_energy)
             print("OOOV Energy = %20.15f" % ooov_energy)
             print("VVVO Energy = %20.15f" % vvvo_energy)
             print("OVOV Energy = %20.15f" % ovov_energy)
@@ -146,22 +147,22 @@ def compute_energy(self):
         print("CC Correlation Energy  = %20.15f" % ecc)
 
         self.ecc = ecc
+        self.eone = eone
+        self.etwo = etwo
 
         return ecc
 
-    def compute_onepdm(self, t1, t2, l1, l2, withref=False):
+    def compute_onepdm(self, t1, t2, l1, l2):
         """
         Parameters
         ----------
         t1, t2, l1, l2 : NumPy arrays
             current cluster amplitudes
-        withref : Boolean (default: False)
-            include the reference contribution if True
 
         Returns
         -------
         onepdm : NumPy array
-            the CC one-electron density as a single, full matrix
+            the CC one-electron density as a single, full matrix (only the correlated contribution)
         """
         o = self.ccwfn.o
         v = self.ccwfn.v
@@ -173,7 +174,7 @@ def compute_onepdm(self, t1, t2, l1, l2, withref=False):
         L = self.ccwfn.H.L
 
         if isinstance(t1, torch.Tensor):
-            if self.ccwfn.precision == 'DP':      
+            if self.ccwfn.precision == 'DP':
                 opdm = torch.zeros((nt, nt), dtype=torch.complex128, device=self.ccwfn.device1)
             elif self.ccwfn.precision == 'SP':
                 opdm = torch.zeros((nt, nt), dtype=torch.complex64, device=self.ccwfn.device1)
@@ -183,18 +184,6 @@ def compute_onepdm(self, t1, t2, l1, l2, withref=False):
             elif self.ccwfn.precision == 'SP':
                 opdm = np.zeros((nt, nt), dtype='complex64')
         opdm[o,o] = self.build_Doo(t1, t2, l1, l2)
-        if withref is True:
-            if isinstance(t1, torch.Tensor):
-                if self.ccwfn.precision == 'DP':
-                    opdm[o,o] += 2.0 * torch.eye(no, dtype=torch.complex128, device=self.ccwfn.device1)  # Reference contribution
-                elif self.ccwfn.precision == 'SP':
-                    opdm[o,o] += 2.0 * torch.eye(no, dtype=torch.complex64, device=self.ccwfn.device1) 
-            else:
-                if self.ccwfn.precision == 'DP':
-                    opdm[o,o] += 2.0 * np.eye(no)  # Reference contribution
-                elif self.ccwfn.precision == 'SP':
-                    opdm[o,o] += 2.0 * np.eye(no, dtype=np.complex64) 
-
         opdm[v,v] = self.build_Dvv(t1, t2, l1, l2)
         opdm[o,v] = self.build_Dov(t1, t2, l1, l2)
         opdm[v,o] = self.build_Dvo(l1)
@@ -206,20 +195,20 @@ def compute_onepdm(self, t1, t2, l1, l2, withref=False):
             Wovoo = self.ccwfn.build_cc3_Wmbij(o, v, ERI, t1, Woooo)
             Wvovv = self.ccwfn.build_cc3_Wamef(o, v, ERI, t1)
             Wooov = self.ccwfn.build_cc3_Wmnie(o, v, ERI, t1)
-           
+
             opdm[o,v] += self.build_cc3_Dov(o, v, no, nv, F, L, t1, t2, l1, l2, Wvvvo, Wovoo, Fov, Wvovv, Wooov)
 
             # Density matrix blocks in contractions with T1-transformed dipole integrals
             if isinstance(t1, torch.Tensor):
                 opdm_cc3 = torch.zeros_like(opdm)
-            else: 
+            else:
                 opdm_cc3 = np.zeros_like(opdm)
             opdm_cc3[o,o] += self.build_cc3_Doo(o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov)
             opdm_cc3[v,v] += self.build_cc3_Dvv(o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov)
 
             return (opdm, opdm_cc3)
-                        
-        else:       
+
+        else:
             return opdm
 
     def build_Doo(self, t1, t2, l1, l2):  # complete
@@ -229,6 +218,9 @@ def build_Doo(self, t1, t2, l1, l2):  # complete
         else:
             Doo = -1.0 * contract('ie,je->ij', t1, l1)
             Doo -= contract('imef,jmef->ij', t2, l2)
+            # (T) contributions computed in ccwfn.t3_density()
+            if self.ccwfn.model == 'CCSD(T)':
+                Doo += self.ccwfn.Doo
 
         return Doo
 
@@ -240,6 +232,9 @@ def build_Dvv(self, t1, t2, l1, l2):  # complete
         else:
             Dvv = contract('mb,ma->ab', t1, l1)
             Dvv += contract('mnbe,mnae->ab', t2, l2)
+            # (T) contributions computed in ccwfn.t3_density()
+            if self.ccwfn.model == 'CCSD(T)':
+                Dvv += self.ccwfn.Dvv
 
         return Dvv
 
@@ -270,11 +265,14 @@ def build_Dov(self, t1, t2, l1, l2):  # complete
             tmp = contract('mnef,mnaf->ea', l2, t2)
             Dov -= contract('ea,ie->ia', tmp, t1)
 
+            if self.ccwfn.model == 'CCSD(T)':
+                Dov += self.ccwfn.Dov
+
             if isinstance(tmp, torch.Tensor):
                 del tmp
 
         return Dov
-    
+
     # CC3 contributions to the one electron densities
     def build_cc3_Dov(self, o, v, no, nv, F, L, t1, t2, l1, l2, Wvvvo, Wovoo, Fov, Wvovv, Wooov):
         contract = self.contract  
@@ -386,6 +384,10 @@ def build_Dooov(self, t1, t2, l1, l2):  # complete
             tmp = contract('kmef,jf->kmej', l2, t1)
             tmp = contract('kmej,ie->kmij', tmp, t1)
             Dooov += contract('kmij,ma->ijka', tmp, t1)
+
+            # (T) contributions to twopdm computed in ccwfn.t3_density()
+            if self.ccwfn.model == 'CCSD(T)':
+                Dooov += self.ccwfn.Gooov
 
             if isinstance(tmp, torch.Tensor):
                 del tmp
@@ -434,6 +436,10 @@ def build_Dvvvo(self, t1, t2, l1, l2):  # complete
             tmp = contract('nmci,na->amci', tmp, t1)
             Dvvvo -= contract('amci,mb->abci', tmp, t1)
 
+            # (T) contributions to twopdm computed in ccwfn.t3_density()
+            if self.ccwfn.model == 'CCSD(T)':
+                Dvvvo += self.ccwfn.Gvvvo
+
             if isinstance(tmp, torch.Tensor):
                 del tmp
                 del Gvv
@@ -550,6 +556,10 @@ def build_Doovv(self, t1, t2, l1, l2):
             tmp = contract('nb,mnij->mbij', t1, tmp)
             Doovv += contract('ma,mbij->ijab', t1, tmp)
 
+            # (T) contributions to twopdm computed in ccwfn.t3_density()
+            if self.ccwfn.model == 'CCSD(T)':
+                Doovv += self.ccwfn.Goovv
+
             if isinstance(tmp, torch.Tensor):
                 del tmp, tmp1, tmp2, Goo, Gvv
 

pycc/cclambda.py

@@ -458,6 +458,10 @@ def r_L1(self, o, v, l1, l2, Hov, Hvv, Hoo, Hovvo, Hovov, Hvvvo, Hovoo, Hvovv, H
             else: 
                 r_l1 = 2.0 * Hov.copy()
 
+            # Add (T) contributions to L1
+            if self.ccwfn.model == 'CCSD(T)':
+                r_l1 = r_l1 + self.ccwfn.S1
+
             r_l1 = r_l1 + contract('ie,ea->ia', l1, Hvv)
             r_l1 = r_l1 - contract('ma,im->ia', l1, Hoo)          
             r_l1 = r_l1 + contract('imef,efam->ia', l2, Hvvvo)
@@ -483,6 +487,7 @@ def r_L2(self, o, v, l1, l2, L, Hov, Hvv, Hoo, Hoooo, Hvvvv, Hovvo, Hovov, Hvvvo
                 r_l2 = L[o,o,v,v].clone().to(self.ccwfn.device1)
             else:
                 r_l2 = L[o,o,v,v].copy()
+
             r_l2 = r_l2 + contract('ijeb,ea->ijab', l2, Hvv)
             r_l2 = r_l2 - contract('mjab,im->ijab', l2, Hoo)
             r_l2 = r_l2 + 0.5 * contract('mnab,ijmn->ijab', l2, Hoooo)
@@ -497,6 +502,11 @@ def r_L2(self, o, v, l1, l2, L, Hov, Hvv, Hoo, Hoooo, Hvvvv, Hovvo, Hovov, Hvvvo
                 r_l2 = L[o,o,v,v].clone().to(self.ccwfn.device1)
             else:
                 r_l2 = L[o,o,v,v].copy()
+
+            # Add (T) contributions to L2
+            if self.ccwfn.model == 'CCSD(T)':
+                r_l2 = r_l2 + 0.5 * self.ccwfn.S2
+
             r_l2 = r_l2 + 2.0 * contract('ia,jb->ijab', l1, Hov)
             r_l2 = r_l2 - contract('ja,ib->ijab', l1, Hov)
             r_l2 = r_l2 + 2.0 * contract('ie,ejab->ijab', l1, Hvovv)

pycc/cctriples.py

@@ -82,6 +82,25 @@ def t3d_ijk(o, v, i, j, k, t1, t2, Woovv, F, contract, WithDenom=True):
     else:
         return t3
 
+def t3d_abc(o, v, a, b, c, t1, t2, Woovv, F, contract, WithDenom=True):
+    t3 = contract('ij,k->ijk', Woovv[:,:,a,b], t1[:,c])
+    t3 += contract('ik,j->ijk', Woovv[:,:,a,c], t1[:,b])
+    t3 += contract('jk,i->ijk', Woovv[:,:,b,c], t1[:,a])
+    Fov = F[o,v]
+    t3 += contract('ij,k->ijk', t2[:,:,a,b], Fov[:,c])
+    t3 += contract('ik,j->ijk', t2[:,:,a,c], Fov[:,b])
+    t3 += contract('jk,i->ijk', t2[:,:,b,c], Fov[:,a])
+
+    if WithDenom is True:
+        no = o.stop
+        Fo = np.diag(F)[o]
+        denom = np.zeros_like(t3)
+        denom += Fo.reshape(-1,1,1) + Fo.reshape(-1,1) + Fo
+        denom -= F[a+no,a+no] + F[b+no,b+no] + F[c+no,c+no]
+        return t3/denom
+    else:
+        return t3
+
 
 # Lee and Rendell's formulation
 def t_tjl(ccwfn):
@@ -212,16 +231,15 @@ def l3_ijk(i, j, k, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDenom=T
     l3 += contract('b,ca->abc', Fov[j], l2[k,i]) - contract('c,ba->abc', Fov[j], l2[k,i])
     l3 += contract('c,ba->abc', Fov[k], l2[j,i]) - contract('b,ca->abc', Fov[k], l2[j,i])
 
-
     tmp_W = 2 * Wvovv - Wvovv.swapaxes(2,3)
     W = contract('eab,ce->abc',  tmp_W[:,j,:,:], l2[k,i])
     W += contract('eac,be->abc', tmp_W[:,k,:,:], l2[j,i])
     W += contract('eba,ce->abc', tmp_W[:,i,:,:], l2[k,j])
     W += contract('eca,be->abc', tmp_W[:,i,:,:], l2[j,k])
     W += contract('ebc,ae->abc', tmp_W[:,k,:,:], l2[i,j])
     W += contract('ecb,ae->abc', tmp_W[:,j,:,:], l2[i,k])
-    
-    W -= contract('ebc,ea->abc', Wvovv[:,i,:,:], l2[j,k,:,:])    
+
+    W -= contract('ebc,ea->abc', Wvovv[:,i,:,:], l2[j,k,:,:])
     W -= contract('ecb,ea->abc', Wvovv[:,i,:,:], l2[k,j,:,:])
     W -= contract('eba,ec->abc', Wvovv[:,k,:,:], l2[j,i,:,:])
     W -= contract('eac,eb->abc', Wvovv[:,j,:,:], l2[i,k,:,:])
@@ -242,9 +260,9 @@ def l3_ijk(i, j, k, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDenom=T
     W += contract('mc,mab->abc', Wooov[j,i,:,:], l2[k])
     W += contract('ma,mcb->abc', Wooov[j,k,:,:], l2[i])
     W += contract('mb,mac->abc', Wooov[k,i,:,:], l2[j])
-  
+
     l3 += W
-   
+
     if WithDenom is True:
         if isinstance(l2, torch.Tensor):
             Fv = torch.diag(F)[v]
@@ -277,16 +295,16 @@ def l3_abc(a, b, c, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDenom=T
     l3 += contract('j,ki->ijk', Fov[:,b], l2[:,:,c,a]) - contract('j,ki->ijk', Fov[:,c], l2[:,:,b,a])
     l3 += contract('k,ji->ijk', Fov[:,c], l2[:,:,b,a]) - contract('k,ji->ijk', Fov[:,b], l2[:,:,c,a])
 
- 
+
     tmp_W = 2 * Wvovv - Wvovv.swapaxes(2,3)
     W = contract('ej,kie->ijk',  tmp_W[:,:,a,b], l2[:,:,c,:])
     W += contract('ek,jie->ijk', tmp_W[:,:,a,c], l2[:,:,b,:])
     W += contract('ei,kje->ijk', tmp_W[:,:,b,a], l2[:,:,c,:])
     W += contract('ei,jke->ijk', tmp_W[:,:,c,a], l2[:,:,b,:])
     W += contract('ek,ije->ijk', tmp_W[:,:,b,c], l2[:,:,a,:])
     W += contract('ej,ike->ijk', tmp_W[:,:,c,b], l2[:,:,a,:])
-    
-    W -= contract('ei,jke->ijk', Wvovv[:,:,b,c], l2[:,:,:,a])    
+
+    W -= contract('ei,jke->ijk', Wvovv[:,:,b,c], l2[:,:,:,a])
     W -= contract('ei,kje->ijk', Wvovv[:,:,c,b], l2[:,:,:,a])
     W -= contract('ek,jie->ijk', Wvovv[:,:,b,a], l2[:,:,:,c])
     W -= contract('ej,ike->ijk', Wvovv[:,:,a,c], l2[:,:,:,b])
@@ -307,9 +325,9 @@ def l3_abc(a, b, c, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDenom=T
     W += contract('jim,km->ijk', Wooov[:,:,:,c], l2[:,:,a,b])
     W += contract('jkm,im->ijk', Wooov[:,:,:,a], l2[:,:,c,b])
     W += contract('kim,jm->ijk', Wooov[:,:,:,b], l2[:,:,a,c])
-  
+
     l3 += W
-   
+
     if WithDenom is True:
         no = o.stop
         if isinstance(l2, torch.Tensor):
@@ -355,9 +373,9 @@ def l3_ijk_alt(i, j, k, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDen
     W -= contract('mc,mba->abc', Wooov[i,k,:,:], l2[j])
     W -= contract('mb,mac->abc', Wooov[k,j,:,:], l2[i])
     W -= contract('mc,mab->abc', Wooov[j,k,:,:], l2[i])
-    
-    l3 += 2 * W - W.swapaxes(0,1) - W.swapaxes(0,2)    
-    
+
+    l3 += 2 * W - W.swapaxes(0,1) - W.swapaxes(0,2)
+
     if WithDenom is True:
         if isinstance(l2, torch.Tensor):
             Fv = torch.diag(F)[v]
@@ -405,7 +423,7 @@ def l3_abc_alt(a, b, c, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDen
     W -= contract('jkm,im->ijk', Wooov[:,:,:,c], l2[:,:,a,b])
 
     l3 += 2 * W - W.swapaxes(0,1) - W.swapaxes(0,2)
-        
+
     if WithDenom is True:
         no = o.stop
         if isinstance(l2, torch.Tensor):
@@ -423,7 +441,7 @@ def l3_abc_alt(a, b, c, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDen
 # Triples drivers that are useful for density matrix calculation
 # W_bc(ijka)
 def t3c_bc(o, v, b, c, t2, Wvvvo, Wovoo, F, contract, WithDenom=True):
-    
+
     t3 = contract('aei,kje->ijka', Wvvvo[b], t2[:,:,c])
     t3 += contract('aei,jke->ijka', Wvvvo[c], t2[:,:,b])
     t3 += contract('aek,jie->ijka', Wvvvo[:,c], t2[:,:,b])
@@ -474,16 +492,16 @@ def l3_bc(b, c, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDenom=True)
     l3 += contract('j,kia->ijka', Fov[:,b], l2[:,:,c]) - contract('j,kia->ijka', Fov[:,c], l2[:,:,b])
     l3 += contract('k,jia->ijka', Fov[:,c], l2[:,:,b]) - contract('k,jia->ijka', Fov[:,b], l2[:,:,c])
 
- 
+
     tmp_W = 2 * Wvovv - Wvovv.swapaxes(2,3)
     W = contract('eja,kie->ijka',  tmp_W[:,:,:,b], l2[:,:,c,:])
     W += contract('eka,jie->ijka', tmp_W[:,:,:,c], l2[:,:,b,:])
     W += contract('eia,kje->ijka', tmp_W[:,:,b], l2[:,:,c,:])
     W += contract('eia,jke->ijka', tmp_W[:,:,c], l2[:,:,b,:])
     W += contract('ek,ijae->ijka', tmp_W[:,:,b,c], l2)
     W += contract('ej,ikae->ijka', tmp_W[:,:,c,b], l2)
-    
-    W -= contract('ei,jkea->ijka', Wvovv[:,:,b,c], l2)    
+
+    W -= contract('ei,jkea->ijka', Wvovv[:,:,b,c], l2)
     W -= contract('ei,kjea->ijka', Wvovv[:,:,c,b], l2)
     W -= contract('eka,jie->ijka', Wvovv[:,:,b], l2[:,:,:,c])
     W -= contract('eja,ike->ijka', Wvovv[:,:,:,c], l2[:,:,:,b])
@@ -504,9 +522,9 @@ def l3_bc(b, c, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDenom=True)
     W += contract('jim,kma->ijka', Wooov[:,:,:,c], l2[:,:,:,b])
     W += contract('jkma,im->ijka', Wooov, l2[:,:,c,b])
     W += contract('kim,jma->ijka', Wooov[:,:,:,b], l2[:,:,:,c])
-  
+
     l3 += W
-   
+
     if WithDenom is True:
         no = o.stop
         if isinstance(l2, torch.Tensor):