diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml
index f327c25..b437672 100644
--- a/.github/workflows/CI.yaml
+++ b/.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
diff --git a/devtools/conda-envs/psi.yaml b/devtools/conda-envs/psi.yaml
index 9ca2bd1..3395a05 100644
--- a/devtools/conda-envs/psi.yaml
+++ b/devtools/conda-envs/psi.yaml
@@ -1,6 +1,5 @@
 name: p4dev
 channels:
-  - psi4/label/dev
   - conda-forge
 dependencies:
   - psi4
diff --git a/pycc/ccdensity.py b/pycc/ccdensity.py
index a23155e..f911474 100644
--- a/pycc/ccdensity.py
+++ b/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
 
diff --git a/pycc/cclambda.py b/pycc/cclambda.py
index f4162c8..aad49ec 100644
--- a/pycc/cclambda.py
+++ b/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)
diff --git a/pycc/cctriples.py b/pycc/cctriples.py
index 48754d1..82f9526 100644
--- a/pycc/cctriples.py
+++ b/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,7 +231,6 @@ 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])
@@ -220,8 +238,8 @@ def l3_ijk(i, j, k, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDenom=T
     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,7 +295,7 @@ 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,:])
@@ -285,8 +303,8 @@ def l3_abc(a, b, c, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDenom=T
     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,7 +492,7 @@ 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,:])
@@ -482,8 +500,8 @@ def l3_bc(b, c, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDenom=True)
     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):
diff --git a/pycc/ccwfn.py b/pycc/ccwfn.py
index 2795329..c90e710 100644
--- a/pycc/ccwfn.py
+++ b/pycc/ccwfn.py
@@ -13,7 +13,7 @@
 from .utils import helper_diis, cc_contract
 from .hamiltonian import Hamiltonian
 from .local import Local
-from .cctriples import t_tjl, t3c_ijk
+from .cctriples import t_tjl, t3c_ijk, t3d_ijk, t3c_abc, t3d_abc
 from .lccwfn import lccwfn
 
 class ccwfn(object):
@@ -52,7 +52,7 @@ class ccwfn(object):
         the final CC correlation energy
 
     Methods
-    ------- 
+    -------
     solve_cc()
         Solves the CC T amplitude equations
     residuals()
@@ -65,26 +65,28 @@ def __init__(self, scf_wfn, **kwargs):
         ----------
         scf_wfn : Psi4 Wavefunction Object
             computed by Psi4 energy() method
-        
+
         Returns
         -------
         None
         """
 
         time_init = time.time()
-             
+
         valid_cc_models = ['CCD', 'CC2', 'CCSD', 'CCSD(T)', 'CC3']
         model = kwargs.pop('model','CCSD').upper()
         if model not in valid_cc_models:
             raise Exception("%s is not an allowed CC model." % (model))
         self.model = model
-        
+
         # models requiring singles
         self.need_singles = ['CCSD', 'CCSD(T)', 'CC2', 'CC3']
 
         # models requiring T1-transformed integrals
         self.need_t1_transform = ['CC3']
 
+        self.make_t3_density = kwargs.pop('make_t3_density', False)
+
         valid_local_models = [None, 'PNO', 'PAO','PNO++']
         local = kwargs.pop('local', None)
         # TODO: case-protect this kwarg
@@ -154,7 +156,7 @@ def __init__(self, scf_wfn, **kwargs):
                 self.Local.trans_integrals(self.o, self.v)
                 self.Local.overlaps(self.Local.QL)
                 self.lccwfn = lccwfn(self.o, self.v,self.no, self.nv, self.H, self.local, self.model, self.eref, self.Local)
-        
+
         # denominators
         eps_occ = np.diag(self.H.F)[o]
         eps_vir = np.diag(self.H.F)[v]
@@ -174,13 +176,13 @@ def __init__(self, scf_wfn, **kwargs):
         if precision.upper() not in valid_precision:
             raise Exception('%s is not an allowed precision arithmetic.' % (precision))
         self.precision = precision.upper()
-         
+
         valid_device = ['CPU', 'GPU']
         device = kwargs.pop('device', 'CPU')
         if device.upper() not in valid_device:
             raise Exception("%s is not an allowed device." % (device))
         self.device = device.upper()
-        
+
         if self.precision == 'SP':
             self.H.F = np.float32(self.H.F)
             self.t1 = np.float32(self.t1)
@@ -188,12 +190,12 @@ def __init__(self, scf_wfn, **kwargs):
             self.Dia = np.float32(self.Dia)
             self.Dijab = np.float32(self.Dijab)
             self.H.ERI = np.float32(self.H.ERI)
-            self.H.L = np.float32(self.H.L)    
+            self.H.L = np.float32(self.H.L)
 
         # Initiate the object for a generalized contraction function 
         # for GPU or CPU.  
         self.contract = cc_contract(device=self.device)
-       
+
         # Convert the arrays to torch.Tensors if the calculation is on GPU.
         # Send the copy of F, t1, t2 to GPU.
         # ERI will be kept on GPU
@@ -224,7 +226,8 @@ def __init__(self, scf_wfn, **kwargs):
                 self.H.L = torch.tensor(self.H.L, dtype=torch.complex64, device=self.device0)
                 
         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):
         """
         Parameters
@@ -253,7 +256,7 @@ def solve_cc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_diis
         L = self.H.L
         Dia = self.Dia
         Dijab = self.Dijab
-        
+
         contract = self.contract
 
         ecc = self.cc_energy(o, v, F, L, self.t1, self.t2)
@@ -283,8 +286,8 @@ def solve_cc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_diis
                 rms = contract('ia,ia->', r1/Dia, r1/Dia)
                 rms += contract('ijab,ijab->', r2/Dijab, r2/Dijab)
                 if isinstance(r1, torch.Tensor):
-                    rms = torch.sqrt(rms) 
-                else:               
+                    rms = torch.sqrt(rms)
+                else:
                     rms = np.sqrt(rms)
 
             ecc = self.cc_energy(o, v, F, L, self.t1, self.t2)
@@ -305,10 +308,13 @@ def solve_cc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_diis
                     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)
                         print("E(CCSD) = %20.15f" % ecc)
-                        print("E(T)  = %20.15f" % et)
-                        ecc = ecc + et   
+                        if self.make_t3_density is True:
+                            et = self.t3_density()
+                        else:
+                            et = t_tjl(self)
+                        print("E(T)    = %20.15f" % et)
+                        ecc = ecc + et
                     else:
                         print("E(%s) = %20.15f" % (self.model, ecc))
                     self.ecc = ecc
@@ -335,7 +341,7 @@ def residuals(self, F, t1, t2):
         r1, r2: NumPy arrays
             New T1 and T2 residuals: r_mu = <mu|HBAR|0>
         """
-    
+
         contract = self.contract
 
         o = self.o
@@ -355,41 +361,41 @@ def residuals(self, F, t1, t2):
 
         r1 = self.r_T1(o, v, F, ERI, L, t1, t2, Fae, Fme, Fmi)
         r2 = self.r_T2(o, v, F, ERI, L, t1, t2, Fae, Fme, Fmi, Wmnij, Wmbej, Wmbje, Zmbij)
-             
+
         if isinstance(Fae, torch.Tensor):
-            del Fae, Fmi, Wmnij, Wmbej, Wmbje, Zmbij 
-        
-        if self.model == 'CC3':      
+            del Fae, Fmi, Wmnij, Wmbej, Wmbje, Zmbij
+
+        if self.model == 'CC3':
             Wmnij_cc3 = self.build_cc3_Wmnij(o, v, ERI, t1)
             Wmbij_cc3 = self.build_cc3_Wmbij(o, v, ERI, t1, Wmnij_cc3)
             Wmnie_cc3 = self.build_cc3_Wmnie(o, v, ERI, t1)
             Wamef_cc3 = self.build_cc3_Wamef(o, v, ERI, t1)
             Wabei_cc3 = self.build_cc3_Wabei(o, v, ERI, t1)
-            
+
             if isinstance(t1, torch.Tensor):
                 X1 = torch.zeros_like(t1)
                 X2 = torch.zeros_like(t2)
             else:
                 X1 = np.zeros_like(t1)
                 X2 = np.zeros_like(t2)
-              
+
             for i in range(no):
                 for j in range(no):
-                    for k in range(no):                       
+                    for k in range(no):
                         t3 = t3c_ijk(o, v, i, j, k, t2, Wabei_cc3, Wmbij_cc3, F,
 contract, WithDenom=True)
-                        
-                        X1[i] += contract('abc,bc->a', t3 - t3.swapaxes(0,2), L[j,k,v,v])                       
+
+                        X1[i] += contract('abc,bc->a', t3 - t3.swapaxes(0,2), L[j,k,v,v])
                         X2[i,j] += contract('abc,c->ab', t3 - t3.swapaxes(0,2), Fme[k])
                         X2[i,j] += contract('abc,dbc->ad', 2 * t3 - t3.swapaxes(1,2) - t3.swapaxes(0,2), Wamef_cc3.swapaxes(0,1)[k])
                         X2[i] -= contract('abc,lc->lab', 2 * t3 - t3.swapaxes(1,2) - t3.swapaxes(0,2), Wmnie_cc3[j,k])
-            
+
             r1 += X1
             r2 += X2 + X2.swapaxes(0,1).swapaxes(2,3)
-                       
+
             if isinstance(t3, torch.Tensor):
-                del Fme, Wmnij_cc3, Wmbij_cc3, Wmnie_cc3, Wamef_cc3, Wabei_cc3     
-             
+                del Fme, Wmnij_cc3, Wmbij_cc3, Wmnie_cc3, Wamef_cc3, Wabei_cc3
+
         return r1, r2
 
     def build_tau(self, t1, t2, fact1=1.0, fact2=1.0):
@@ -529,7 +535,7 @@ def r_T1(self, o, v, F, ERI, L, t1, t2, Fae, Fme, Fmi):
         if self.model == 'CCD':
             if isinstance(t1, torch.Tensor):
                 r_T1 = torch.zero_like(t1)
-            else: 
+            else:
                 r_T1 = np.zeros_like(t1)
         else:
             if isinstance(t1, torch.Tensor):
@@ -573,7 +579,7 @@ def r_T2(self, o, v, F, ERI, L, t1, t2, Fae, Fme, Fmi, Wmnij, Wmbej, Wmbje, Zmbi
             r_T2 = r_T2 + 0.5 * contract('ma,mbij->ijab', t1, contract('nb,mnij->mbij', t1, Wmnij))
             r_T2 = r_T2 + 0.5 * contract('jf,abif->ijab', t1, contract('ie,abef->abif', t1, ERI[v,v,v,v]))
             r_T2 = r_T2 - contract('ma,mbij->ijab', t1, Zmbij)
-            r_T2 = r_T2 - contract('ma,mbij->ijab', t1, contract('ie,mbej->mbij', t1, ERI[o,v,v,o])) 
+            r_T2 = r_T2 - contract('ma,mbij->ijab', t1, contract('ie,mbej->mbij', t1, ERI[o,v,v,o]))
             r_T2 = r_T2 - contract('mb,maji->ijab', t1, contract('ie,maje->maji', t1, ERI[o,v,o,v]))
             r_T2 = r_T2 + contract('ie,abej->ijab', t1, ERI[v,v,v,o])
             r_T2 = r_T2 - contract('ma,mbij->ijab', t1, ERI[o,v,o,o])
@@ -601,7 +607,7 @@ def r_T2(self, o, v, F, ERI, L, t1, t2, Fae, Fme, Fmi, Wmnij, Wmbej, Wmbje, Zmbi
             r_T2 = r_T2 - contract('imeb,maje->ijab', tmp, ERI[o,v,o,v])
             r_T2 = r_T2 + contract('ie,abej->ijab', t1, ERI[v,v,v,o])
             r_T2 = r_T2 - contract('ma,mbij->ijab', t1, ERI[o,v,o,o])
-            
+
             if isinstance(tmp, torch.Tensor):
                 del tmp
 
@@ -629,7 +635,7 @@ def build_cc3_Wmbij(self, o, v, ERI, t1, Wmnij):
             W = ERI[o,v,o,o].copy()
         W = W - contract('mnij,nb->mbij', Wmnij, t1)
         W = W + contract('mbie,je->mbij', ERI[o,v,o,v], t1)
-        if isinstance(t1, torch.Tensor):            
+        if isinstance(t1, torch.Tensor):
             tmp = ERI[o,v,v,o].clone().to(self.device1) + contract('mbef,jf->mbej', ERI[o,v,v,v], t1)
         else:
             tmp = ERI[o,v,v,o].copy() + contract('mbef,jf->mbej', ERI[o,v,v,v], t1)
@@ -652,7 +658,7 @@ def build_cc3_Wamef(self, o, v, ERI, t1):
         else:
             W = ERI[v,o,v,v].copy()
         W = W - contract('na,nmef->amef', t1, ERI[o,o,v,v])
-        return W   
+        return W
 
     def build_cc3_Wabei(self, o, v, ERI, t1):
         contract =self.contract
@@ -695,7 +701,7 @@ def build_cc3_Wabei(self, o, v, ERI, t1):
         # Wabei
         W = Z_abei + Z_eiab.swapaxes(0,2).swapaxes(1,3)
         return W
- 
+
     def cc_energy(self, o, v, F, L, t1, t2):
         contract = self.contract
         if self.model == 'CCD':
@@ -704,3 +710,139 @@ def cc_energy(self, o, v, F, L, t1, t2):
             ecc = 2.0 * contract('ia,ia->', F[o,v], t1)
             ecc = ecc + contract('ijab,ijab->', self.build_tau(t1, t2), L[o,o,v,v])
         return ecc
+
+    def t3_density(self):
+        """
+        Computes (T) contributions to Lambda equations and one-/two-electron densities
+        """
+
+        contract = self.contract
+
+        o = self.o
+        v = self.v
+        no = self.no
+        nv = self.nv
+        t1 = self.t1
+        t2 = self.t2
+        F = self.H.F
+        ERI = self.H.ERI
+        L = self.H.L
+
+        Dvv = np.zeros((nv,nv))
+        Doo = np.zeros((no,no))
+        Dov = np.zeros((no,nv))
+        Goovv = np.zeros_like(t2)
+        Gooov = np.zeros((no,no,no,nv))
+        Gvvvo = np.zeros((nv,nv,nv,no))
+        S1 = np.zeros_like(t1)
+        S2 = np.zeros_like(t2)
+        Z3 = np.zeros((nv,nv,nv))
+        X1 = np.zeros_like(t1)
+        X2 = np.zeros_like(t2)
+
+        for i in range(no):
+            for j in range(no):
+                for k in range(no):
+                    M3 = t3c_ijk(o, v, i, j, k, t2, ERI[v,v,v,o], ERI[o,v,o,o], F, contract, True)
+                    N3 = t3d_ijk(o, v, i, j, k, t1, t2, ERI[o,o,v,v], F, contract, True)
+                    X3 = 8*M3 - 4*M3.swapaxes(0,1) - 4*M3.swapaxes(1,2) - 4*M3.swapaxes(0,2) + 2*np.moveaxis(M3, 0, 2) + 2*np.moveaxis(M3, 2, 0)
+                    Y3 = 8*N3 - 4*N3.swapaxes(0,1) - 4*N3.swapaxes(1,2) - 4*N3.swapaxes(0,2) + 2*np.moveaxis(N3, 0, 2) + 2*np.moveaxis(N3, 2, 0)
+
+                    # Doubles contribution (T) correction (Viking's formulation)
+                    X2[i,j] += contract('abc,c->ab',(M3 - M3.swapaxes(0,2)), F[k,v])
+                    X2[i,j] += contract('abc,dbc->ad', (2*M3 - M3.swapaxes(1,2) - M3.swapaxes(0,2)),ERI[v,k,v,v])
+                    X2[i] -= contract('abc,lc->lab', (2*M3 - M3.swapaxes(1,2) - M3.swapaxes(0,2)),ERI[j,k,o,v])
+
+                    # (T) contribution to vir-vir block of one-electron density
+                    Dvv += 0.5 * contract('acd,bcd->ab', M3, (X3 + Y3))
+
+                    # (T) contribution to occ-vir block of one-electron density
+                    Dov[i] += contract('abc,bc->a', (M3 - M3.swapaxes(0,2)), (4*t2[j,k] - 2*t2[j,k].T))
+
+                    # (T) contributions to two-electron density
+                    Z3 = 2*(M3 - M3.swapaxes(1,2)) - (M3.swapaxes(0,1) - np.moveaxis(M3, 2, 0))
+                    Goovv[i,j,:,:] += 4*contract('c,abc->ab', t1[k,:], Z3)
+                    Gooov[j,i] -= contract('abc,lbc->la', (2*X3 + Y3), t2[:,k])
+                    Gvvvo[:,:,:,j] += contract('abc,cd->abd', (2*X3 + Y3), t2[k,i,:,:])
+
+                    # (T) contribution to Lambda_1 residual
+                    S1[i] += contract('abc,bc->a', 2*(M3 - M3.swapaxes(0,1)), L[j,k,v,v])
+                    # (T) contribution to Lambda_2 residual
+                    S2[i] -= contract('abc,lc->lab', (2*X3 + Y3), ERI[j,k,o,v])
+                    S2[i,j] += contract('abc,dcb->ad', (2*X3 + Y3), ERI[k,v,v,v])
+
+        S2 = S2 + S2.swapaxes(0,1).swapaxes(2,3)
+
+        # (T) contribution to occ-occ block of one-electron density
+        for a in range(nv):
+            for b in range(nv):
+                for c in range(nv):
+                    M3 = t3c_abc(o, v, a, b, c, t2, ERI[v,v,v,o], ERI[o,v,o,o], F, contract, True)
+                    N3 = t3d_abc(o, v, a, b, c, t1, t2, ERI[o,o,v,v], F, contract, True)
+                    X3 = 8*M3 - 4*M3.swapaxes(0,1) - 4*M3.swapaxes(1,2) - 4*M3.swapaxes(0,2) + 2*np.moveaxis(M3, 0, 2) + 2*np.moveaxis(M3, 2, 0)
+                    Y3 = 8*N3 - 4*N3.swapaxes(0,1) - 4*N3.swapaxes(1,2) - 4*N3.swapaxes(0,2) + 2*np.moveaxis(N3, 0, 2) + 2*np.moveaxis(N3, 2, 0)
+                    Doo -= 0.5 * contract('ikl,jkl->ij', M3, (X3 + Y3))
+
+        self.Dvv = Dvv
+        self.Doo = Doo
+        self.Dov = Dov # Need to add this even though it doesn't contribute to the energy for RHF references
+
+        self.Goovv = Goovv
+        self.Gooov = Gooov
+        self.Gvvvo = Gvvvo
+        self.S1 = S1
+        self.S2 = S2
+
+        # (T) correction
+        ET = contract('ia,ia->', t1, S1) # NB: Factor of two is already included in S1 definition
+        ET += contract('ijab,ijab->', (4.0*t2 - 2.0*t2.swapaxes(2,3)), X2)
+
+#        print("Dvv:")
+#        it = np.nditer(self.Dvv, flags=['multi_index'])
+#        for val in it:
+#            if np.abs(val) > 1e-12:
+#                print("%s %20.15f" % (it.multi_index, val))
+#
+#        print("Doo:")
+#        it = np.nditer(self.Doo, flags=['multi_index'])
+#        for val in it:
+#            if np.abs(val) > 1e-12:
+#                print("%s %20.15f" % (it.multi_index, val))
+#
+#        print("Dov:")
+#        it = np.nditer(self.Dov, flags=['multi_index'])
+#        for val in it:
+#            if np.abs(val) > 1e-12:
+#                print("%s %20.15f" % (it.multi_index, val))
+#
+#        print("S1 Amplitudes:")
+#        it = np.nditer(self.S1, flags=['multi_index'])
+#        for val in it:
+#            if np.abs(val) > 1e-12:
+#                print("%s %20.15f" % (it.multi_index, val))
+#
+#        print("S2 Amplitudes:")
+#        it = np.nditer(self.S2, flags=['multi_index'])
+#        for val in it:
+#            if np.abs(val) > 1e-12:
+#                print("%s %20.15f" % (it.multi_index, val))
+#
+#        print("Goovv Density:")
+#        it = np.nditer(self.Goovv, flags=['multi_index'])
+#        for val in it:
+#            if np.abs(val) > 1e-12:
+#                print("%s %20.15f" % (it.multi_index, val))
+#
+#        print("Gooov Density:")
+#        it = np.nditer(self.Gooov, flags=['multi_index'])
+#        for val in it:
+#            if np.abs(val) > 1e-12:
+#                print("%s %20.15f" % (it.multi_index, val))
+#
+#        print("Gvvvo Density:")
+#        it = np.nditer(self.Gvvvo, flags=['multi_index'])
+#        for val in it:
+#            if np.abs(val) > 1e-12:
+#                print("%s %20.15f" % (it.multi_index, val))
+
+        return ET
diff --git a/pycc/rt/rtcc.py b/pycc/rt/rtcc.py
index 69d3007..91fae0d 100644
--- a/pycc/rt/rtcc.py
+++ b/pycc/rt/rtcc.py
@@ -207,26 +207,24 @@ def extract_amps(self, y):
 
         return t1, t2, l1, l2, phase
 
-    def dipole(self, t1, t2, l1, l2, withref = True, magnetic = False):
+    def dipole(self, t1, t2, l1, l2, magnetic = False):
         """
         Parameters
         ----------
         t1, t2, l1, l2 : NumPy arrays
             current cluster amplitudes
-        withref        : Bool (default = True)
-            include reference contribution to the OPDM
         magnetic       : Bool (default = False)
             compute magnetic dipole rather than electric
 
         Returns
         -------
         x, y, z : scalars
-            Cartesian components of the dipole moment
+            Cartesian components of the correlated dipole moment
         """
         if self.ccwfn.model == 'CC3':
-            (opdm, opdm_cc3) = self.ccdensity.compute_onepdm(t1, t2, l1, l2, withref=withref)
+            (opdm, opdm_cc3) = self.ccdensity.compute_onepdm(t1, t2, l1, l2)
         else:
-            opdm = self.ccdensity.compute_onepdm(t1, t2, l1, l2, withref=withref)
+            opdm = self.ccdensity.compute_onepdm(t1, t2, l1, l2)
 
         if magnetic:
             ints = self.m
@@ -423,12 +421,12 @@ def step(self,ODE,yi,t,ref=False):
         ret = {}
         t1, t2, l1, l2, phase = self.extract_amps(y)
         ret['ecc'] = self.lagrangian(t,t1,t2,l1,l2)
-        mu_x, mu_y, mu_z = self.dipole(t1,t2,l1,l2,withref=ref,magnetic=False)
+        mu_x, mu_y, mu_z = self.dipole(t1,t2,l1,l2,magnetic=False)
         ret['mu_x'] = mu_x
         ret['mu_y'] = mu_y
         ret['mu_z'] = mu_z
         if self.magnetic:
-            m_x, m_y, m_z = self.dipole(t1,t2,l1,l2,withref=ref,magnetic=True)
+            m_x, m_y, m_z = self.dipole(t1,t2,l1,l2,magnetic=True)
             ret['m_x'] = m_x
             ret['m_y'] = m_y
             ret['m_z'] = m_z
@@ -509,12 +507,12 @@ def propagate(self, ODE, yi, tf, ti=0, ref=False, chk=False, tchk=False,
         # initial properties
         t1, t2, l1, l2, phase = self.extract_amps(yi)
         ret[key]['ecc'] = self.lagrangian(ti,t1,t2,l1,l2)
-        mu_x, mu_y, mu_z = self.dipole(t1,t2,l1,l2,withref=ref,magnetic=False)
+        mu_x, mu_y, mu_z = self.dipole(t1,t2,l1,l2,magnetic=False)
         ret[key]['mu_x'] = mu_x
         ret[key]['mu_y'] = mu_y
         ret[key]['mu_z'] = mu_z
         if self.magnetic:
-            m_x, m_y, m_z = self.dipole(t1,t2,l1,l2,withref=ref,magnetic=True)
+            m_x, m_y, m_z = self.dipole(t1,t2,l1,l2,magnetic=True)
             ret[key]['m_x'] = m_x
             ret[key]['m_y'] = m_y
             ret[key]['m_z'] = m_z
diff --git a/pycc/tests/test_007_dipole.py b/pycc/tests/test_007_dipole.py
index 8ddefdb..e805667 100644
--- a/pycc/tests/test_007_dipole.py
+++ b/pycc/tests/test_007_dipole.py
@@ -6,10 +6,11 @@
 import psi4
 import pycc
 import pytest
+import numpy as np
 from ..data.molecules import *
 
 def test_dipole_h2_2_cc_pvdz():
-    """He cc-pVDZ"""
+    """H4 cc-pVDZ"""
     psi4.set_memory('2 GiB')
     psi4.core.set_output_file('output.dat', False)
     psi4.set_options({'basis': 'cc-pVDZ',
@@ -41,9 +42,9 @@ def test_dipole_h2_2_cc_pvdz():
     y0 = rtcc.collect_amps(cc.t1, cc.t2, cclambda.l1, cclambda.l2, ecc).astype('complex128')
     t1, t2, l1, l2, phase = rtcc.extract_amps(y0)
 
-    ref = [0, 0, 0.005371586416860086] # au
+    ref = np.array([0, 0, -0.0007395036977002]) # computed by removing SCF from original ref
+
     mu_x, mu_y, mu_z = rtcc.dipole(t1, t2, l1, l2)
-    opdm = rtcc.ccdensity.compute_onepdm(t1, t2, l1, l2, withref = True)
 
     assert (abs(ref[0] - mu_x) < 1E-10)
     assert (abs(ref[1] - mu_y) < 1E-10)
diff --git a/pycc/tests/test_021_rk4.py b/pycc/tests/test_021_rk4.py
index 0cf92ed..1e297fe 100644
--- a/pycc/tests/test_021_rk4.py
+++ b/pycc/tests/test_021_rk4.py
@@ -87,7 +87,7 @@ def test_rtcc_water_cc_pvdz():
         """
         
     print(mu_z)
-    mu_z_ref = -0.34894577
+    mu_z_ref = -0.0780067603267549 # computed by removing SCF from original ref
     assert (abs(mu_z_ref - mu_z.real) < 1e-4)
     
     #return (dip_x, dip_y, dip_z, time_points)
diff --git a/pycc/tests/test_022_adap_int.py b/pycc/tests/test_022_adap_int.py
index 2ee3576..4cf520a 100644
--- a/pycc/tests/test_022_adap_int.py
+++ b/pycc/tests/test_022_adap_int.py
@@ -90,7 +90,7 @@ def test_rtcc_water_cc_pvdz():
         """
         
     print(mu_z)
-    mu_z_ref = -0.34894577
+    mu_z_ref = -0.0780067603267549 # computed by removing SCF from original ref
     assert (abs(mu_z_ref - mu_z.real) < 1e-3)
     
     #return (dip_x, dip_y, dip_z, time_points)
diff --git a/pycc/tests/test_023_ms_int.py b/pycc/tests/test_023_ms_int.py
index 2bf5ece..b48f477 100644
--- a/pycc/tests/test_023_ms_int.py
+++ b/pycc/tests/test_023_ms_int.py
@@ -102,7 +102,7 @@ def test_rtcc_water_cc_pvdz():
         """
         
     print(mu_z)
-    mu_z_ref = -0.34894577
+    mu_z_ref = -0.0780067603267549 # computed by removing SCF from original ref
     assert (abs(mu_z_ref - mu_z.real) < 1e-1)
     
     #return (dip_x, dip_y, dip_z, time_points)
diff --git a/pycc/tests/test_024_contract_cpu.py b/pycc/tests/test_024_contract_cpu.py
index 71f255c..94922cd 100644
--- a/pycc/tests/test_024_contract_cpu.py
+++ b/pycc/tests/test_024_contract_cpu.py
@@ -90,7 +90,7 @@ def test_rtcc_water_cc_pvdz():
         """
         
     print(mu_z)
-    mu_z_ref = -0.34894577
+    mu_z_ref = -0.0780067603267549 # computed by removing SCF from original ref
     assert (abs(mu_z_ref - mu_z.real) < 1e-4)
     
     #return (dip_x, dip_y, dip_z, time_points)
diff --git a/pycc/tests/test_030_sp.py b/pycc/tests/test_030_sp.py
index d242a4c..1d52141 100644
--- a/pycc/tests/test_030_sp.py
+++ b/pycc/tests/test_030_sp.py
@@ -74,7 +74,7 @@ def test_rtcc_water_cc_pvdz():
     ecc0 = rtcc.lagrangian(t0, t1, t2, l1, l2)
    
     # Check dipole moment 
-    mu0_z_ref = -0.3489459218340155 
+    mu0_z_ref = -0.0780069121607703 # computed by removing SCF from original ref
     assert(abs(mu0_z_ref - mu0_z) < 1e-6)
     
     while t < tf:
@@ -87,6 +87,6 @@ def test_rtcc_water_cc_pvdz():
     print(mu_z)
    
     # Check the dipole value at time step 1
-    mu_z_ref = -0.3489459204749751
+    mu_z_ref = -0.0780069121607703 # computed by removing SCF from original ref
     assert (abs(mu_z_ref - mu_z.real) < 1e-6)    
 
diff --git a/pycc/tests/test_031_cc3.py b/pycc/tests/test_031_cc3.py
index eaaa4a9..cb0f0a3 100644
--- a/pycc/tests/test_031_cc3.py
+++ b/pycc/tests/test_031_cc3.py
@@ -45,14 +45,12 @@ def test_cc3_h2o():
     # no laser
     rtcc = pycc.rtcc(cc, cclambda, ccdensity, None, magnetic = False)
     
-    # Nuclear dipole from Psi4: 1.12729111248 au
-    #nuc_dipole = mol.nuclear_dipole()
-    # Total dipole from CFOUR: [0, 0, 0.7703875967] au
-    ref = [0, 0, -0.3569035158] # au
+    CFOUR = [0, 0, 0.7703875967] # CFOUR total dipole (CC3 + SCF + nuclear)
+    scf = rhf_wfn.variable('SCF DIPOLE') # PSI4 reference dipole (SCF + nuclear)
+    ref = CFOUR - scf # Final reference: CC3 only
 
-    mu_x, mu_y, mu_z = rtcc.dipole(cc.t1, cc.t2, cclambda.l1, cclambda.l2, withref=True)
+    mu_x, mu_y, mu_z = rtcc.dipole(cc.t1, cc.t2, cclambda.l1, cclambda.l2)
 
-    assert (abs(ref[0] - np.real(mu_x)) < 1E-10)
     assert (abs(ref[1] - np.real(mu_y)) < 1E-10)
     assert (abs(ref[2] - np.real(mu_z)) < 1E-10)
 
diff --git a/pycc/tests/test_034_ccsd_t_density.py b/pycc/tests/test_034_ccsd_t_density.py
new file mode 100644
index 0000000..e924b7b
--- /dev/null
+++ b/pycc/tests/test_034_ccsd_t_density.py
@@ -0,0 +1,85 @@
+"""
+Test CCSD equation solution using various molecule test cases.
+"""
+
+# Import package, test suite, and other packages as needed
+import psi4
+import pycc
+import pytest
+from ..data.molecules import *
+from ..cctriples import t_vikings, t_vikings_inverted, t_tjl
+
+def test_ccsd_t_h2o():
+    """H2O cc-pVDZ"""
+    # 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(
+"""
+O 0.000000000000000   0.000000000000000   0.143225857166674
+H 0.000000000000000  -1.638037301628121  -1.136549142277225
+H 0.000000000000000   1.638037301628121  -1.136549142277225
+symmetry c1
+units bohr
+"""
+)
+    rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True)
+
+    maxiter = 75
+    e_conv = 1e-12
+    r_conv = 1e-12
+
+    etot = psi4.energy('CCSD(T)')
+    ecc_psi4 = psi4.variable('CCSD(T) CORRELATION ENERGY')
+
+    cc = pycc.ccwfn(rhf_wfn, model='ccsd(t)', make_t3_density=True)
+    ecc = cc.solve_cc(e_conv, r_conv, maxiter, max_diis=0)
+    assert (abs(ecc_psi4 - ecc) < 1e-11)
+    hbar = pycc.cchbar(cc)
+    cclambda = pycc.cclambda(cc, hbar)
+    lcc = cclambda.solve_lambda(e_conv, r_conv, maxiter, max_diis=0)
+    ccdensity = pycc.ccdensity(cc, cclambda)
+    ecc_density = ccdensity.compute_energy()
+    eone = ccdensity.eone
+    etwo = ccdensity.etwo
+
+    lambda_psi4 = -0.069084521221746
+    eone_psi4 = 0.104463374777302
+    etwo_psi4 = -0.175243393781829
+    assert (abs(lambda_psi4 - lcc) < 1e-11)
+    assert (abs(eone_psi4 - eone) < 1e-11)
+    assert (abs(etwo_psi4 - etwo) < 1e-11)
+
+    psi4.core.clean()
+
+    psi4.set_options({'basis': 'cc-pVDZ'})
+    rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True)
+
+    etot = psi4.energy('CCSD(T)')
+    ecc_psi4 = psi4.variable('CCSD(T) CORRELATION ENERGY')
+
+    cc = pycc.ccwfn(rhf_wfn, model='ccsd(t)', make_t3_density=True)
+    ecc = cc.solve_cc(e_conv, r_conv, maxiter)
+    assert (abs(ecc_psi4 - ecc) < 1e-11)
+    hbar = pycc.cchbar(cc)
+    cclambda = pycc.cclambda(cc, hbar)
+    lcc = cclambda.solve_lambda(e_conv, r_conv)
+    ccdensity = pycc.ccdensity(cc, cclambda)
+    ecc_density = ccdensity.compute_energy()
+    eone = ccdensity.eone
+    etwo = ccdensity.etwo
+
+    lambda_psi4 = -0.227199866607450
+    eone_psi4 = 0.251210862963227
+    etwo_psi4 = -0.479006477929931
+    assert (abs(lambda_psi4 - lcc) < 1e-11)
+    assert (abs(eone_psi4 - eone) < 1e-11)
+    assert (abs(etwo_psi4 - etwo) < 1e-11)
diff --git a/versioneer.py b/versioneer.py
index 64fea1c..3aa5da3 100644
--- a/versioneer.py
+++ b/versioneer.py
@@ -339,9 +339,9 @@ def get_config_from_root(root):
     # configparser.NoOptionError (if it lacks "VCS="). See the docstring at
     # the top of versioneer.py for instructions on writing your setup.cfg .
     setup_cfg = os.path.join(root, "setup.cfg")
-    parser = configparser.SafeConfigParser()
+    parser = configparser.ConfigParser()
     with open(setup_cfg, "r") as f:
-        parser.readfp(f)
+        parser.read_file(f)
     VCS = parser.get("versioneer", "VCS")  # mandatory
 
     def get(parser, name):