OK for ppLR

input/methods

@@ -9,7 +9,7 @@
 # CIS* CIS(D) CID CISD FCI
   F    F      F   F    F
 # phRPA* phRPAx* crRPA ppRPA 
-  F      F       F     F 
+  F      F       F     T 
 # G0F2* evGF2* qsGF2* G0F3 evGF3
   F     F      F      F    F
 # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW

input/options

@@ -5,7 +5,7 @@
 # CC: maxSCF thresh   DIIS n_diis
       64     0.0000001  T    5
 # spin: TDA singlet triplet spin_conserved spin_flip 
-        F   T       F       T              T 
+        F   T       T       T              T 
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm reg
       256    0.00001 T    5      T   0.0 0      F
 # GW: maxSCF thresh  DIIS n_diis lin eta COHSEX TDA_W reg
@@ -13,6 +13,6 @@
 # GT: maxSCF thresh  DIIS n_diis lin eta TDA_T reg
       256    0.00001  T    5      T   0.0 F     F  
 # ACFDT: AC Kx  XBS
-         F  T   T
+         T  F   F
 # BSE: phBSE phBSE2 ppBSE dBSE dTDA evDyn 
        F     F      F     T   F     F

src/GT/G0T0pp.f90

@@ -108,8 +108,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
   allocate(Bpp(nVVab,nOOab),Cpp(nVVab,nVVab),Dpp(nOOab,nOOab))
 
   if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI_MO,Bpp)
-                 call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eHF,ERI_MO,Cpp)
+                 call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVab,1d0,eHF,ERI_MO,Cpp)
-                 call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eHF,ERI_MO,Dpp)
+                 call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,1d0,eHF,ERI_MO,Dpp)
 
   call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
 
@@ -130,8 +130,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
   allocate(Bpp(nVVaa,nOOaa),Cpp(nVVaa,nVVaa),Dpp(nOOaa,nOOaa))
 
   if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI_MO,Bpp)
-                 call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eHF,ERI_MO,Cpp)
+                 call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVaa,1d0,eHF,ERI_MO,Cpp)
-                 call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eHF,ERI_MO,Dpp)
+                 call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,1d0,eHF,ERI_MO,Dpp)
 
   call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
 
@@ -220,8 +220,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
   allocate(Bpp(nVVab,nOOab),Cpp(nVVab,nVVab),Dpp(nOOab,nOOab))
 
   if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI_MO,Bpp)
-                 call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eGT,ERI_MO,Cpp)
+                 call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVab,1d0,eGT,ERI_MO,Cpp)
-                 call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eGT,ERI_MO,Dpp)
+                 call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,1d0,eGT,ERI_MO,Dpp)
 
   call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
 
@@ -233,8 +233,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
   allocate(Bpp(nVVaa,nOOaa),Cpp(nVVaa,nVVaa),Dpp(nOOaa,nOOaa))
 
   if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI_MO,Bpp)
-                 call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eGT,ERI_MO,Cpp)
+                 call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVaa,1d0,eGT,ERI_MO,Cpp)
-                 call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eGT,ERI_MO,Dpp)
+                 call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,1d0,eGT,ERI_MO,Dpp)
 
   call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
 

src/GT/evGTpp.f90

@@ -138,8 +138,8 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
     allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
 
     if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp)
-                   call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eHF,ERI_MO,Cpp)
+                   call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI_MO,Cpp)
-                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eHF,ERI_MO,Dpp)
+                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI_MO,Dpp)
 
     call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
 
@@ -157,8 +157,8 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
     allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
 
     if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp)
-                   call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eHF,ERI_MO,Cpp)
+                   call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI_MO,Cpp)
-                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eHF,ERI_MO,Dpp)
+                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI_MO,Dpp)
 
     call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
 

src/GT/qsGTpp.f90

@@ -194,8 +194,8 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
     allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
 
     if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp)
-                   call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO,Cpp)
+                   call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp)
-                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO,Dpp)
+                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp)
 
     call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
 
@@ -207,8 +207,8 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
     allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
 
     if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp)
-                   call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO,Cpp)
+                   call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp)
-                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO,Dpp)
+                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp)
 
     call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
 

src/LR/ppLR.f90

@@ -1,6 +1,6 @@
 subroutine ppLR(TDA,nOO,nVV,B,C,D,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
 
-! Compute the p-p channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013)
+! Compute the pp channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013)
 
   implicit none
   include 'parameters.h'
@@ -16,8 +16,6 @@ subroutine ppLR(TDA,nOO,nVV,B,C,D,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
   
 ! Local variables
 
-  integer                       :: ab,cd,ij,kl
-  integer                       :: p,q,r,s
   double precision              :: trace_matrix
   double precision              :: EcRPA1
   double precision              :: EcRPA2
@@ -77,15 +75,15 @@ subroutine ppLR(TDA,nOO,nVV,B,C,D,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
 
   ! Split the various quantities in p-p and h-h parts
 
-    call sort_ppRPA(nOO,nVV,Om(:),Z(:,:),Om1(:),X1(:,:),Y1(:,:),Om2(:),X2(:,:),Y2(:,:))
+    call sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
 
   end if
 
 ! Compute the RPA correlation energy
 
-  EcRPA = 0.5d0*( sum(Om1(:)) - sum(Om2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )
+  EcRPA = 0.5d0*( sum(Om1) - sum(Om2) - trace_matrix(nVV,C) - trace_matrix(nOO,D) )
-  EcRPA1 = +sum(Om1(:)) - trace_matrix(nVV,C(:,:))
+  EcRPA1 = +sum(Om1) - trace_matrix(nVV,C)
-  EcRPA2 = -sum(Om2(:)) - trace_matrix(nOO,D(:,:))
+  EcRPA2 = -sum(Om2) - trace_matrix(nOO,D)
 
   if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) & 
     print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!'

src/RPA/ppACFDT.f90

@@ -1,4 +1,4 @@
-subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
+subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcAC)
 
 ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for pp sector
 
@@ -17,7 +17,6 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
   integer,intent(in)            :: nO
   integer,intent(in)            :: nV
   integer,intent(in)            :: nR
-  integer,intent(in)            :: nS
 
   double precision,intent(in)   :: e(nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
@@ -32,6 +31,9 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
   integer                       :: nOO
   integer                       :: nVV
 
+  double precision,allocatable  :: Bpp(:,:)
+  double precision,allocatable  :: Cpp(:,:)
+  double precision,allocatable  :: Dpp(:,:)
   double precision,allocatable  :: Om1(:)
   double precision,allocatable  :: X1(:,:)
   double precision,allocatable  :: Y1(:,:)
@@ -64,7 +66,8 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
     nVV = nV*(nV+1)/2
 
     allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),rho1(nBas,nBas,nVV), & 
-             Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO),rho2(nBas,nBas,nOO))
+             Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO),rho2(nBas,nBas,nOO), &
+             Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
 
     write(*,*) '--------------'
     write(*,*) 'Singlet states'
@@ -79,9 +82,13 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
  
       lambda = rAC(iAC)
 
-      call ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om1,X1,Y1,Om2,X2,Y2,EcAC(ispin))
+      if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
+                   call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
+                   call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
 
-      call ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y1,X2,Y2,Ec(iAC,ispin))
+      call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcAC(ispin))
+
+      call ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,ERI,nOO,nVV,X1,Y1,X2,Y2,Ec(iAC,ispin))
 
       write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
 
@@ -94,7 +101,7 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
     write(*,*) '-----------------------------------------------------------------------------------'
     write(*,*)
 
-    deallocate(Om1,X1,Y1,rho1,Om2,X2,Y2,rho2)
+    deallocate(Om1,X1,Y1,rho1,Om2,X2,Y2,rho2,Bpp,Cpp,Dpp)
 
   end if
  
@@ -108,7 +115,8 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
     nVV = nV*(nV-1)/2
 
     allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),rho1(nBas,nBas,nVV), & 
-             Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO),rho2(nBas,nBas,nOO))
+             Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO),rho2(nBas,nBas,nOO), &
+             Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
 
     write(*,*) '--------------'
     write(*,*) 'Triplet states'
@@ -123,11 +131,13 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
  
       lambda = rAC(iAC)
 
-      ! Initialize T matrix
+      if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
+                   call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
+                   call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
 
-      call ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om1,X1,Y1,Om2,X2,Y2,EcAC(ispin))
+      call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcAC(ispin))
 
-      call ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y1,X2,Y2,Ec(iAC,ispin))
+      call ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,ERI,nOO,nVV,X1,Y1,X2,Y2,Ec(iAC,ispin))
 
       write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
 
@@ -140,7 +150,7 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
     write(*,*) '-----------------------------------------------------------------------------------'
     write(*,*)
 
-    deallocate(Om1,X1,Y1,rho1,Om2,X2,Y2,rho2)
+    deallocate(Om1,X1,Y1,rho1,Om2,X2,Y2,rho2,Bpp,Cpp,Dpp)
 
   end if
 

src/RPA/ppACFDT_correlation_energy.f90

@@ -1,4 +1,4 @@
-subroutine ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y1,X2,Y2,EcAC)
+subroutine ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,ERI,nOO,nVV,X1,Y1,X2,Y2,EcAC)
 
 ! Compute the correlation energy via the adiabatic connection formula for the pp sector
 
@@ -13,7 +13,6 @@ subroutine ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y
   integer,intent(in)            :: nO
   integer,intent(in)            :: nV
   integer,intent(in)            :: nR
-  integer,intent(in)            :: nS
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
 
   integer,intent(in)            :: nOO
@@ -53,4 +52,3 @@ subroutine ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y
        - trace_matrix(nVV,C) - trace_matrix(nOO,D)
 
 end subroutine 
-

src/RPA/ppRPA.f90

@@ -1,6 +1,6 @@
-subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipole_int,e)
+subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
 
-! Perform pp-RPA calculation
+! Perform ppRPA calculation
 
   implicit none
   include 'parameters.h'
@@ -17,7 +17,7 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo
   integer,intent(in)            :: nV
   integer,intent(in)            :: nR
   double precision,intent(in)   :: ENuc
-  double precision,intent(in)   :: ERHF
+  double precision,intent(in)   :: EHF
   double precision,intent(in)   :: e(nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
@@ -25,9 +25,11 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo
 ! Local variables
 
   integer                       :: ispin
-  integer                       :: nS
   integer                       :: nOO
   integer                       :: nVV
+  double precision,allocatable  :: Bpp(:,:)
+  double precision,allocatable  :: Cpp(:,:)
+  double precision,allocatable  :: Dpp(:,:)
   double precision,allocatable  :: Om1(:)
   double precision,allocatable  :: X1(:,:)
   double precision,allocatable  :: Y1(:,:)
@@ -35,7 +37,7 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo
   double precision,allocatable  :: X2(:,:)
   double precision,allocatable  :: Y2(:,:)
 
-  double precision              :: Ec_ppRPA(nspin)
+  double precision              :: EcRPA(nspin)
   double precision              :: EcAC(nspin)
 
 ! Hello world
@@ -48,12 +50,8 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo
 
 ! Initialization
 
-  Ec_ppRPA(:) = 0d0
+  EcRPA(:) = 0d0
-  EcAC(:)     = 0d0
+  EcAC(:)  = 0d0
-
-! Useful quantities
-
-  nS = nO*nV
 
 ! Singlet manifold
 
@@ -69,16 +67,21 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo
     nOO = nO*(nO+1)/2
     nVV = nV*(nV+1)/2
 
-    allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO))
+    allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
+             Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
 
-    call ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI,Om1,X1,Y1,Om2,X2,Y2,Ec_ppRPA(ispin))
+    if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
+                 call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,e,ERI,Cpp)
+                 call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,e,ERI,Dpp)
+
+    call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA(ispin))
 
 !   call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
 
-    call print_excitation('pp-BSE (N+2)',ispin,nVV,Om1)
+    call print_excitation('ppRPA (N+2) ',ispin,nVV,Om1)
-    call print_excitation('pp-BSE (N-2)',ispin,nOO,Om2)
+    call print_excitation('ppRPA (N-2) ',ispin,nOO,Om2)
 
-    deallocate(Om1,X1,Y1,Om2,X2,Y2)
+    deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp)
 
   endif
 
@@ -96,25 +99,30 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo
     nOO = nO*(nO-1)/2
     nVV = nV*(nV-1)/2
 
-    allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO))
+    allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
+             Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
 
-    call ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI,Om1,X1,Y1,Om2,X2,Y2,Ec_ppRPA(ispin))
+    if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
+                 call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,e,ERI,Cpp)
+                 call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,e,ERI,Dpp)
+
+    call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA(ispin))
 
 !   call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
 
-    call print_excitation('pp-BSE (N+2)',ispin,nVV,Om1)
+    call print_excitation('ppRPA (N+2) ',ispin,nVV,Om1)
-    call print_excitation('pp-BSE (N-2)',ispin,nOO,Om2)
+    call print_excitation('ppRPA (N-2) ',ispin,nOO,Om2)
 
-    deallocate(Om1,X1,Y1,Om2,X2,Y2)
+    deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp)
 
   endif
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (singlet) =',Ec_ppRPA(1)
+  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (singlet) =',EcRPA(1)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (triplet) =',3d0*Ec_ppRPA(2)
+  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (triplet) =',3d0*EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy           =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2)
+  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy           =',EcRPA(1) + 3d0*EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy                 =',ENuc + ERHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2)
+  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy                 =',ENuc + EHF + EcRPA(1) + 3d0*EcRPA(2)
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
@@ -122,22 +130,22 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo
 
   if(doACFDT) then
 
-    write(*,*) '---------------------------------------------------------'
+    write(*,*) '--------------------------------------------------------'
-    write(*,*) 'Adiabatic connection version of pp-RPA correlation energy'
+    write(*,*) 'Adiabatic connection version of ppRPA correlation energy'
-    write(*,*) '---------------------------------------------------------'
+    write(*,*) '--------------------------------------------------------'
     write(*,*)
 
-    call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
+    call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcAC)
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcAC(1),' au'
     write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcAC(2),' au'
     write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy           =',EcAC(1) + EcAC(2),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy                 =',ENuc + ERHF + EcAC(1) + EcAC(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy                 =',ENuc + EHF + EcAC(1) + EcAC(2),' au'
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
   end if
 
-end subroutine ppRPA
+end subroutine