diff --git a/input/methods b/input/methods
index 7530736..9e9a79e 100644
--- a/input/methods
+++ b/input/methods
@@ -2,14 +2,14 @@
   T   F   F
 # MP2 MP3 MP2-F12
   F   F   F
-# CCD CCSD CCSD(T)
-  F   F    T
+# CCD CCSD CCSD(T) ringCCD ladderCCD
+  F   F    F       T       T
 # CIS RPA TDHF ppRPA ADC
-  F   T   T    F     F
+  F   T   T    T     F
 # GF2 GF3
   F   F
 # G0W0 evGW qsGW 
-  F    F    F
+  T    F    F
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2
diff --git a/src/QuAcK/CCSD.f90 b/src/QuAcK/CCSD.f90
index bbe6d46..e005062 100644
--- a/src/QuAcK/CCSD.f90
+++ b/src/QuAcK/CCSD.f90
@@ -223,14 +223,14 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
 
   end if
 
-    write(*,*)
-    write(*,*)'----------------------------------------------------'
-    write(*,*)'                 CCSD energy                        '
-    write(*,*)'----------------------------------------------------'
-    write(*,'(1X,A20,1X,F15.10)')' E(CCSD)  = ',ECCSD  
-    write(*,'(1X,A20,1X,F10.6)')' Ec(CCSD) = ',EcCCSD 
-    write(*,*)'----------------------------------------------------'
-    write(*,*)
+  write(*,*)
+  write(*,*)'----------------------------------------------------'
+  write(*,*)'                 CCSD energy                        '
+  write(*,*)'----------------------------------------------------'
+  write(*,'(1X,A20,1X,F15.10)')' E(CCSD)  = ',ECCSD  
+  write(*,'(1X,A20,1X,F10.6)')' Ec(CCSD) = ',EcCCSD 
+  write(*,*)'----------------------------------------------------'
+  write(*,*)
 
 ! Deallocate memory
 
diff --git a/src/QuAcK/Ec_AC.f90 b/src/QuAcK/Ec_AC.f90
index a4aa102..b7f2b7c 100644
--- a/src/QuAcK/Ec_AC.f90
+++ b/src/QuAcK/Ec_AC.f90
@@ -1,4 +1,4 @@
-subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC)
+subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,EcAC)
 
 ! Compute the correlation energy via the adiabatic connection formula
 
@@ -12,6 +12,7 @@ subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC)
   integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: XpY(nS,nS)
+  double precision,intent(in)   :: XmY(nS,nS)
 
 ! Local variables
 
@@ -20,7 +21,10 @@ subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC)
   double precision              :: delta_spin
   double precision              :: delta_dRPA
   double precision,allocatable  :: P(:,:)
-  double precision,allocatable  :: V(:,:)
+  double precision,allocatable  :: Ap(:,:)
+  double precision,allocatable  :: Bp(:,:)
+  double precision,allocatable  :: X(:,:)
+  double precision,allocatable  :: Y(:,:)
   double precision,external     :: trace_matrix
 
 ! Output variables
@@ -40,7 +44,7 @@ subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC)
 
 ! Memory allocation
 
-  allocate(P(nS,nS),V(nS,nS))
+  allocate(P(nS,nS),Ap(nS,nS),Bp(nS,nS),X(nS,nS),Y(nS,nS))
 
 ! Compute P = (X+Y)(X+Y) - 1
 
@@ -50,7 +54,7 @@ subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC)
     P(ia,ia) = P(ia,ia) - 1d0
   enddo
 
-! Compute Viajb = (ia|bj)
+! Compute Aiajb = (ia|bj) and Biajb = (ia|jb)
 
   ia = 0
   do i=nC+1,nO
@@ -61,16 +65,28 @@ subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC)
         do b=nO+1,nBas-nR
           jb = jb + 1
 
-            V(ia,jb) = (1d0 + delta_spin)*ERI(i,b,a,j) 
+            Ap(ia,jb) = (1d0 + delta_spin)*ERI(i,b,a,j) 
+            Bp(ia,jb) = (1d0 + delta_spin)*ERI(i,j,b,a) 
 
         enddo
       enddo
     enddo
   enddo
 
-! Compute Tr(VP)
+! Compute Tr(A x P)
 
-  EcAC = trace_matrix(nS,matmul(V,P))
+! EcAC = trace_matrix(nS,matmul(Ap,P))
+
+! print*,'EcAC =',EcAC
+
+  X(:,:) = 0.5d0*(XpY(:,:) + XmY(:,:))
+  Y(:,:) = 0.5d0*(XpY(:,:) - XmY(:,:))
+
+  EcAC = trace_matrix(nS,matmul(X,matmul(Bp,transpose(Y))) + matmul(Y,matmul(Bp,transpose(X)))) &
+       + trace_matrix(nS,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) &
+       - trace_matrix(nS,Ap)
+
+! print*,'EcAC =',EcAC
 
 end subroutine Ec_AC
 
diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90
index 2d6df3a..0343877 100644
--- a/src/QuAcK/G0W0.f90
+++ b/src/QuAcK/G0W0.f90
@@ -38,6 +38,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
   double precision,allocatable  :: Z(:)
   double precision,allocatable  :: Omega(:,:)
   double precision,allocatable  :: XpY(:,:,:)
+  double precision,allocatable  :: XmY(:,:,:)
   double precision,allocatable  :: rho(:,:,:,:)
   double precision,allocatable  :: rhox(:,:,:,:)
 
@@ -79,7 +80,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
 
 ! Memory allocation
 
-  allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),  & 
+  allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),  & 
            rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin))
 
   AC = .true.
@@ -88,7 +89,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
 ! Compute linear response
 
   call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & 
-                       rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
+                       rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
 
 ! Compute correlation part of the self-energy 
 
@@ -129,11 +130,11 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
       EcBSE(ispin) = 0d0
 
       call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
-                           rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
+                           rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
       call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
 
       call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, &
-                           rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
+                           rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
       call print_excitation('BSE  ',ispin,nS,Omega(:,ispin))
 
     end if
@@ -146,11 +147,11 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
       EcBSE(ispin) = 0d0
 
       call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
-                           rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
+                           rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
       call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
 
       call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, &
-                           rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
+                           rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
       call print_excitation('BSE  ',ispin,nS,Omega(:,ispin))
 
     end if
@@ -192,13 +193,13 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
           lambda = rAC(iAC)
 
           call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, &
-                               rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
+                               rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
           call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
   
           call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,eG0W0,ERI, &
-                           rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin))
+                           rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
 
-          call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin))
+          call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
 
           write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin)
 
@@ -230,13 +231,13 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
           lambda = rAC(iAC)
 
           call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, &
-                               rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
+                               rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
           call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
          
           call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,eG0W0,ERI, &
-                               rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin))
+                               rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
 
-          call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin))
+          call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
 
           write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin)
 
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index 2df031c..ef2e2c4 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -7,6 +7,7 @@ program QuAcK
   logical                       :: doRHF,doUHF,doMOM 
   logical                       :: doMP2,doMP3,doMP2F12
   logical                       :: doCCD,doCCSD,doCCSDT
+  logical                       :: do_ring_CCD,do_ladder_CCD
   logical                       :: doCIS,doRPA,doTDHF
   logical                       :: doppRPA,doADC
   logical                       :: doGF2,doGF3
@@ -108,14 +109,15 @@ program QuAcK
 
 ! Which calculations do you want to do?
 
-  call read_methods(doRHF,doUHF,doMOM,    &
-                    doMP2,doMP3,doMP2F12, &
-                    doCCD,doCCSD,doCCSDT, &
-                    doCIS,doRPA,doTDHF,   & 
-                    doppRPA,doADC,        &
-                    doGF2,doGF3,          &
-                    doG0W0,doevGW,doqsGW, &
-                    doG0T0,doevGT,doqsGT, &
+  call read_methods(doRHF,doUHF,doMOM,         &
+                    doMP2,doMP3,doMP2F12,      &
+                    doCCD,doCCSD,doCCSDT,      &
+                    do_ring_CCD,do_ladder_CCD, &
+                    doCIS,doRPA,doTDHF,        & 
+                    doppRPA,doADC,             &
+                    doGF2,doGF3,               &
+                    doG0W0,doevGW,doqsGW,      &
+                    doG0T0,doevGT,doqsGT,      &
                     doMCMP2)
 
 ! Read options for methods
@@ -351,10 +353,7 @@ program QuAcK
   if(doCCD) then
 
     call cpu_time(start_CCD)
-!   call ring_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
-!   call ladder_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
     call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
-    call cpu_time(end_CCD)
 
     t_CCD = end_CCD - start_CCD
     write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CCD,' seconds'
@@ -380,6 +379,38 @@ program QuAcK
 
   end if
 
+!------------------------------------------------------------------------
+! Perform ring CCD calculation
+!------------------------------------------------------------------------
+
+  if(do_ring_CCD) then
+
+    call cpu_time(start_CCD)
+    call ring_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
+    call cpu_time(end_CCD)
+
+    t_CCD = end_CCD - start_CCD
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ring CCD = ',t_CCD,' seconds'
+    write(*,*)
+
+  end if
+
+!------------------------------------------------------------------------
+! Perform ladder CCD calculation
+!------------------------------------------------------------------------
+
+  if(do_ladder_CCD) then
+
+    call cpu_time(start_CCD)
+    call ladder_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
+    call cpu_time(end_CCD)
+
+    t_CCD = end_CCD - start_CCD
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CCD,' seconds'
+    write(*,*)
+
+  end if
+
 !------------------------------------------------------------------------
 ! Compute CIS excitations
 !------------------------------------------------------------------------
diff --git a/src/QuAcK/RPA.f90 b/src/QuAcK/RPA.f90
index 4b1c561..2a5765b 100644
--- a/src/QuAcK/RPA.f90
+++ b/src/QuAcK/RPA.f90
@@ -29,6 +29,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
   integer                       :: ispin
   double precision,allocatable  :: Omega(:,:)
   double precision,allocatable  :: XpY(:,:,:)
+  double precision,allocatable  :: XmY(:,:,:)
 
   double precision              :: rho
   double precision              :: EcRPA(nspin)
@@ -65,7 +66,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
 
 ! Memory allocation
 
-  allocate(Omega(nS,nspin),XpY(nS,nS,nspin))
+  allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
 
   AC = .true.
   allocate(EcACRPA(nAC,nspin),EcAC(nAC,nspin))
@@ -77,7 +78,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
     ispin = 1
 
     call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
-                         EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
+                         EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
     call print_excitation('RPA  ',ispin,nS,Omega(:,ispin))
 
   endif
@@ -89,7 +90,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
     ispin = 2
 
     call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
-                         EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
+                         EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
     call print_excitation('RPA  ',ispin,nS,Omega(:,ispin))
 
   endif
@@ -131,9 +132,9 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
           lambda = rAC(iAC)
 
           call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho, &
-                               EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin))
+                               EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
 
-          call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin))
+          call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
 
           write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin)
 
@@ -165,9 +166,9 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
           lambda = rAC(iAC)
 
           call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho, &
-                               EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin))
+                               EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
 
-          call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin))
+          call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
 
           write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin)
 
diff --git a/src/QuAcK/TDHF.f90 b/src/QuAcK/TDHF.f90
index 50f22e2..c52ecee 100644
--- a/src/QuAcK/TDHF.f90
+++ b/src/QuAcK/TDHF.f90
@@ -29,6 +29,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
   integer                       :: ispin
   double precision,allocatable  :: Omega(:,:)
   double precision,allocatable  :: XpY(:,:,:)
+  double precision,allocatable  :: XmY(:,:,:)
 
   double precision              :: rho
   double precision              :: EcRPA(nspin)
@@ -65,7 +66,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
 
 ! Memory allocation
 
-  allocate(Omega(nS,nspin),XpY(nS,nS,nspin))
+  allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
 
   AC = .true.
   allocate(EcACRPA(nAC,nspin),EcAC(nAC,nspin))
@@ -77,7 +78,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
     ispin = 1
 
     call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
-                         EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
+                         EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
     call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
 
   endif
@@ -89,7 +90,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
     ispin = 2
 
     call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
-                         EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
+                         EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
     call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
 
   endif
@@ -129,15 +130,11 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
         do iAC=1,nAC
  
           lambda = rAC(iAC)
-
-!         call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, &
-!                              rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
-!         call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
   
           call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
-                           rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin))
+                           rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
 
-          call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin))
+          call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
 
           write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin)
 
@@ -167,15 +164,11 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
         do iAC=1,nAC
  
           lambda = rAC(iAC)
-
-!         call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, &
-!                              rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
-!         call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
          
           call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
-                               rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin))
+                               rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
 
-          call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin))
+          call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
 
           write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin)
 
diff --git a/src/QuAcK/ladder_CCD.f90 b/src/QuAcK/ladder_CCD.f90
index fd6da4a..fd3adc6 100644
--- a/src/QuAcK/ladder_CCD.f90
+++ b/src/QuAcK/ladder_CCD.f90
@@ -182,12 +182,13 @@ subroutine ladder_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
 
   endif
 
-! Moller-Plesset energies
-
   write(*,*)
-  write(*,'(1X,A15,1X,F10.6)') 'Ec(MP2)     = ',EcMP2
-  write(*,'(1X,A15,1X,F10.6)') 'Ec(MP3)     = ',EcMP3
-  write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4
+  write(*,*)'----------------------------------------------------'
+  write(*,*)'              ladder-CCD energy                     '
+  write(*,*)'----------------------------------------------------'
+  write(*,'(1X,A30,1X,F15.10)')' E(ladder-CCD)  = ',ECCD
+  write(*,'(1X,A30,1X,F10.6)')' Ec(ladder-CCSD) = ',EcCCD
+  write(*,*)'----------------------------------------------------'
   write(*,*)
 
 end subroutine ladder_CCD
diff --git a/src/QuAcK/linear_response.f90 b/src/QuAcK/linear_response.f90
index cda7170..4a9e66c 100644
--- a/src/QuAcK/linear_response.f90
+++ b/src/QuAcK/linear_response.f90
@@ -1,4 +1,4 @@
-subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho,EcRPA,Omega,XpY)
+subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho,EcRPA,Omega,XpY,XmY)
 
 ! Compute linear response
 
@@ -22,7 +22,9 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
 ! Output variables
 
   double precision,intent(out)  :: EcRPA
-  double precision,intent(out)  :: Omega(nS),XpY(nS,nS)
+  double precision,intent(out)  :: Omega(nS)
+  double precision,intent(out)  :: XpY(nS,nS)
+  double precision,intent(out)  :: XmY(nS,nS)
 
 
 ! Memory allocation
@@ -83,6 +85,9 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
   XpY = matmul(transpose(Z),AmBSq)
   call DA(nS,1d0/sqrt(abs(Omega)),XpY)
 
+  XmY = matmul(transpose(Z),AmBSq)
+  call DA(nS,sqrt(abs(Omega)),XmY)
+
 ! print*,'X+Y'
 ! call matout(nS,nS,XpY)
 
diff --git a/src/QuAcK/print_excitation.f90 b/src/QuAcK/print_excitation.f90
index 2f76bff..a32831e 100644
--- a/src/QuAcK/print_excitation.f90
+++ b/src/QuAcK/print_excitation.f90
@@ -7,7 +7,7 @@ subroutine print_excitation(method,ispin,nS,Omega)
 
 ! Input variables
 
-  character*4,intent(in)             :: method
+  character*6,intent(in)             :: method
   integer,intent(in)                 :: ispin,nS
   double precision,intent(in)        :: Omega(nS)
 
diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90
index 602f0b1..638dbf6 100644
--- a/src/QuAcK/read_methods.f90
+++ b/src/QuAcK/read_methods.f90
@@ -1,11 +1,12 @@
-subroutine read_methods(doRHF,doUHF,doMOM,    & 
-                        doMP2,doMP3,doMP2F12, & 
-                        doCCD,doCCSD,doCCSDT, & 
-                        doCIS,doRPA,doTDHF,   & 
-                        doppRPA,doADC,        & 
-                        doGF2,doGF3,          & 
-                        doG0W0,doevGW,doqsGW, & 
-                        doG0T0,doevGT,doqsGT, & 
+subroutine read_methods(doRHF,doUHF,doMOM,         & 
+                        doMP2,doMP3,doMP2F12,      & 
+                        doCCD,doCCSD,doCCSDT,      & 
+                        do_ring_CCD,do_ladder_CCD, &
+                        doCIS,doRPA,doTDHF,        & 
+                        doppRPA,doADC,             & 
+                        doGF2,doGF3,               & 
+                        doG0W0,doevGW,doqsGW,      & 
+                        doG0T0,doevGT,doqsGT,      & 
                         doMCMP2)
 
 ! Read desired methods 
@@ -17,6 +18,7 @@ subroutine read_methods(doRHF,doUHF,doMOM,    &
   logical,intent(out)           :: doRHF,doUHF,doMOM
   logical,intent(out)           :: doMP2,doMP3,doMP2F12
   logical,intent(out)           :: doCCD,doCCSD,doCCSDT
+  logical,intent(out)           :: do_ring_CCD,do_ladder_CCD
   logical,intent(out)           :: doCIS,doRPA,doTDHF,doppRPA,doADC
   logical,intent(out)           :: doGF2,doGF3  
   logical,intent(out)           :: doG0W0,doevGW,doqsGW
@@ -45,6 +47,9 @@ subroutine read_methods(doRHF,doUHF,doMOM,    &
   doCCSD  = .false.
   doCCSDT = .false.
 
+  do_ring_CCD   = .false.
+  do_ladder_CCD = .false.
+
   doCIS   = .false.
   doRPA   = .false.
   doTDHF  = .false.
@@ -83,10 +88,12 @@ subroutine read_methods(doRHF,doUHF,doMOM,    &
 ! Read CC methods
 
   read(1,*) 
-  read(1,*) answer1,answer2,answer3
-  if(answer1 == 'T') doCCD   = .true.
-  if(answer2 == 'T') doCCSD  = .true.
-  if(answer3 == 'T') doCCSDT = .true.
+  read(1,*) answer1,answer2,answer3,answer4,answer5
+  if(answer1 == 'T') doCCD         = .true.
+  if(answer2 == 'T') doCCSD        = .true.
+  if(answer3 == 'T') doCCSDT       = .true.
+  if(answer4 == 'T') do_ring_CCD   = .true.
+  if(answer5 == 'T') do_ladder_CCD = .true.
 
 ! Read excited state methods
 
diff --git a/src/QuAcK/ring_CCD.f90 b/src/QuAcK/ring_CCD.f90
index 4cef290..717ceda 100644
--- a/src/QuAcK/ring_CCD.f90
+++ b/src/QuAcK/ring_CCD.f90
@@ -182,12 +182,13 @@ subroutine ring_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
 
   endif
 
-! Moller-Plesset energies
-
   write(*,*)
-  write(*,'(1X,A15,1X,F10.6)') 'Ec(MP2)     = ',EcMP2
-  write(*,'(1X,A15,1X,F10.6)') 'Ec(MP3)     = ',EcMP3
-  write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4
+  write(*,*)'----------------------------------------------------'
+  write(*,*)'              ring-CCD energy                       '
+  write(*,*)'----------------------------------------------------'
+  write(*,'(1X,A30,1X,F15.10)')' E(ring-CCD) = ',ECCD  
+  write(*,'(1X,A30,1X,F15.10)')' Ec(ring-CCD) = ',EcCCD 
+  write(*,*)'----------------------------------------------------'
   write(*,*)
 
 end subroutine ring_CCD