diff --git a/src/GT/dynamic_Tmatrix_A.f90 b/src/GT/GTpp_dynamic_kernel_Aph.f90
similarity index 94%
rename from src/GT/dynamic_Tmatrix_A.f90
rename to src/GT/GTpp_dynamic_kernel_Aph.f90
index a5c8ec9..0ab8a6f 100644
--- a/src/GT/dynamic_Tmatrix_A.f90
+++ b/src/GT/GTpp_dynamic_kernel_Aph.f90
@@ -1,4 +1,4 @@
-subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Om1,Om2,rho1,rho2,OmBSE,TA,ZA)
+subroutine GTpp_dynamic_kernel_Aph(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Om1,Om2,rho1,rho2,OmBSE,TA,ZA)
 
 ! Compute the dynamic part of the Bethe-Salpeter equation matrices for GT
 
@@ -92,4 +92,4 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Om1,Om2,
 
 !$omp end parallel do
 
-end subroutine dynamic_Tmatrix_A
+end subroutine 
diff --git a/src/GT/dynamic_Tmatrix_B.f90 b/src/GT/GTpp_dynamic_kernel_Bph.f90
similarity index 94%
rename from src/GT/dynamic_Tmatrix_B.f90
rename to src/GT/GTpp_dynamic_kernel_Bph.f90
index 6d67b78..41097ba 100644
--- a/src/GT/dynamic_Tmatrix_B.f90
+++ b/src/GT/GTpp_dynamic_kernel_Bph.f90
@@ -1,4 +1,4 @@
-subroutine dynamic_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,TB,ZB)
+subroutine GTpp_dynamic_kernel_Bph(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,TB,ZB)
 
 ! Compute the off-diagonal dynamic part of the Bethe-Salpeter equation matrices for GT
 
@@ -88,4 +88,4 @@ subroutine dynamic_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
     end do
   end do
 
-end subroutine dynamic_Tmatrix_B
+end subroutine 
diff --git a/src/GT/GTpp_phACFTD.f90 b/src/GT/GTpp_phACFTD.f90
index e15d00c..45c88f9 100644
--- a/src/GT/GTpp_phACFTD.f90
+++ b/src/GT/GTpp_phACFTD.f90
@@ -134,8 +134,8 @@ subroutine GTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet
 
         call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
 
-                     call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TAs)
-        if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TBs)
+                     call GTpp_static_kernel_Aph(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TAs)
+        if(.not.TDA) call GTpp_static_kernel_Bph(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TBs)
 
         isp_T  = 2
         iblock = 4
@@ -152,16 +152,16 @@ subroutine GTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet
 
         call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
-                     call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TAt)
-        if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TBt)
+                     call GTpp_static_kernel_Aph(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TAt)
+        if(.not.TDA) call Gtpp_static_kernel_Bph(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TBt)
 
       end if
 
                    call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,Aph)
       if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
  
-                   Aph(:,:) = Aph(:,:) + TAt(:,:) + TAs(:,:)
-      if(.not.TDA) Bph(:,:) = Bph(:,:) + TBt(:,:) + TBs(:,:)
+                   Aph(:,:) = Aph(:,:) - TAt(:,:) - TAs(:,:)
+      if(.not.TDA) Bph(:,:) = Bph(:,:) - TBt(:,:) - TBs(:,:)
  
       call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
 
@@ -218,8 +218,8 @@ subroutine GTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet
 
         call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
 
-                     call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TAs)
-        if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TBs)
+                     call GTpp_static_kernel_Aph(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TAs)
+        if(.not.TDA) call GTpp_static_kernel_Bph(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TBs)
 
         isp_T  = 2
         iblock = 4
@@ -236,16 +236,16 @@ subroutine GTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet
 
         call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
-                     call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TAt)
-        if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TBt)
+                     call GTpp_static_kernel_Aph(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TAt)
+        if(.not.TDA) call GTpp_static_kernel_Bph(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TBt)
 
       end if
 
                    call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,Aph)
       if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
-                   Aph(:,:) = Aph(:,:) + TAt(:,:) - TAs(:,:)
-      if(.not.TDA) Bph(:,:) = Bph(:,:) + TBt(:,:) - TBs(:,:)
+                   Aph(:,:) = Aph(:,:) + TAs(:,:) - TAt(:,:)
+      if(.not.TDA) Bph(:,:) = Bph(:,:) + TBs(:,:) - TBt(:,:)
 
       call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
 
diff --git a/src/GT/GTpp_phBSE.f90 b/src/GT/GTpp_phBSE.f90
index e0f5137..7f62806 100644
--- a/src/GT/GTpp_phBSE.f90
+++ b/src/GT/GTpp_phBSE.f90
@@ -99,8 +99,8 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n
 
   deallocate(Bpp,Cpp,Dpp)
 
-               call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TAab)
-  if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TBab)
+               call GTpp_static_kernel_Aph(eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TAab)
+  if(.not.TDA) call GTpp_static_kernel_Bph(eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TBab)
 
 !----------------------------------------!
 ! Compute T-matrix for alpha-alpha block !
@@ -119,8 +119,8 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n
 
   deallocate(Bpp,Cpp,Dpp)
 
-               call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TAaa)
-  if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TBaa)
+               call GTpp_static_kernel_Aph(eta,nBas,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TAaa)
+  if(.not.TDA) call GTpp_static_kernel_Bph(eta,nBas,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TBaa)
 
 !------------------!
 ! Singlet manifold !
@@ -153,9 +153,9 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n
 
       else
  
-        call Bethe_Salpeter_Tmatrix_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, &
-                                                         Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, & 
-                                                         dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa)
+        call GTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, &
+                                             Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, & 
+                                             dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa)
       end if
  
     end if
@@ -193,9 +193,9 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n
 
       else
  
-        call Bethe_Salpeter_Tmatrix_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, &
-                                                         Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, & 
-                                                         dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa)
+        call GTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, &
+                                             Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, & 
+                                             dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa)
       end if
  
     end if
diff --git a/src/GT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 b/src/GT/GTpp_phBSE_dynamic_perturbation.f90
similarity index 86%
rename from src/GT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90
rename to src/GT/GTpp_phBSE_dynamic_perturbation.f90
index 273cc1f..b0f6bbe 100644
--- a/src/GT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90
+++ b/src/GT/GTpp_phBSE_dynamic_perturbation.f90
@@ -1,6 +1,5 @@
-subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa,         &
-                                                       Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, &
-                                                       dipole_int,OmBSE,XpY,XmY,TAab,TAaa)
+subroutine GTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa,Om1ab,Om2ab,Om1aa,Om2aa, & 
+                                           rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT,dipole_int,OmBSE,XpY,XmY,TAab,TAaa)
 
 ! Compute dynamical effects via perturbation theory for BSE@GT
 
@@ -82,11 +81,11 @@ subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO
 
     ! Compute dynamical T-matrix for alpha-beta block  
 
-    call dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,eGT,Om1ab,Om2ab,rho1ab,rho2ab,OmBSE(ia),dTAab,ZAab)
+    call GTpp_dynamic_kernel_Aph(eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,eGT,Om1ab,Om2ab,rho1ab,rho2ab,OmBSE(ia),dTAab,ZAab)
  
     ! Compute dynamical T-matrix for alpha-beta block  
 
-    call dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,eGT,Om1aa,Om2aa,rho1aa,rho2aa,OmBSE(ia),dTAaa,ZAaa)
+    call GTpp_dynamic_kernel_Aph(eta,nBas,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,eGT,Om1aa,Om2aa,rho1aa,rho2aa,OmBSE(ia),dTAaa,ZAaa)
  
     X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
     Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
@@ -130,4 +129,4 @@ subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO
   write(*,*) '---------------------------------------------------------------------------------------------------'
   write(*,*) 
 
-end subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation
+end subroutine 
diff --git a/src/GT/GTpp_ppBSE.f90 b/src/GT/GTpp_ppBSE.f90
index 5207aee..3ece5ea 100644
--- a/src/GT/GTpp_ppBSE.f90
+++ b/src/GT/GTpp_ppBSE.f90
@@ -62,6 +62,7 @@ subroutine GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n
   integer                       :: nVVt
 
   double precision              :: EcRPA(nspin)
+  double precision,allocatable  :: Bpp(:,:),Cpp(:,:),Dpp(:,:)
   double precision,allocatable  :: TBab(:,:),TCab(:,:),TDab(:,:)
   double precision,allocatable  :: TBaa(:,:),TCaa(:,:),TDaa(:,:)
 
@@ -94,30 +95,41 @@ subroutine GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n
 
     iblock = 3
  
-    EcRPA(ispin) = 0d0
+    allocate(Bpp(nVVab,nOOab),Cpp(nVVab,nVVab),Dpp(nOOab,nOOab))
  
-    call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eT,ERI,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
+    if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI,Bpp)
+                   call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVab,1d0,eT,ERI,Cpp)
+                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,1d0,eT,ERI,Dpp)
  
+    call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
+ 
+    deallocate(Bpp,Cpp,Dpp)
     allocate(TBab(nVVs,nOOs),TCab(nVVs,nVVs),TDab(nOOs,nOOs))
  
-    if(.not.TDA) call static_Tmatrix_B_pp(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOs,nVVs,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TBab)
-                 call static_Tmatrix_C_pp(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOs,nVVs,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TCab)
-                 call static_Tmatrix_D_pp(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOs,nVVs,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TDab)
+    if(.not.TDA_T) call GTpp_static_kernel_Bpp(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOs,nVVs,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TBab)
+                   call GTpp_static_kernel_Cpp(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOs,nVVs,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TCab)
+                   call GTpp_static_kernel_Dpp(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOs,nVVs,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TDab)
 
   !----------------------------------------!
   ! Compute T-matrix for alpha-alpha block !
   !----------------------------------------!
 
     iblock = 4
-    EcRPA(ispin) = 0d0
  
-    call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eT,ERI,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
-
+    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,Bpp)
+                   call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVaa,1d0,eT,ERI,Cpp)
+                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,1d0,eT,ERI,Dpp)
+ 
+    call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
+ 
+    deallocate(Bpp,Cpp,Dpp)
     allocate(TBaa(nVVs,nOOs),TCaa(nVVs,nVVs),TDaa(nOOs,nOOs))
  
-    if(.not.TDA) call static_Tmatrix_B_pp(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOs,nVVs,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TBaa)
-                 call static_Tmatrix_C_pp(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOs,nVVs,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TCaa)
-                 call static_Tmatrix_D_pp(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOs,nVVs,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TDaa)
+    if(.not.TDA_T) call GTpp_static_kernel_Bpp(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOs,nVVs,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TBaa)
+                   call GTpp_static_kernel_Cpp(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOs,nVVs,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TCaa)
+                   call GTpp_static_kernel_Dpp(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOs,nVVs,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TDaa)
 
   !----------------------------------!
   ! pp/hh sectors for singlet states !
@@ -125,14 +137,22 @@ subroutine GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n
 
     EcBSE(ispin) = 0d0
 
-    allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs))
-  
-    call linear_response_pp_BSE(ispin,TDA,.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI, &
-                                TBaa+TBab,TCaa+TCab,TDaa+TDab,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcBSE(ispin))
+    allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
+             Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
+
+    if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
+                 call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
+                 call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
+
+    Bpp(:,:) = Bpp(:,:) - TBab(:,:) - TBaa(:,:) 
+    Cpp(:,:) = Cpp(:,:) - TCab(:,:) - TCaa(:,:) 
+    Dpp(:,:) = Dpp(:,:) - TDab(:,:) - TDaa(:,:) 
+
+    call ppLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcBSE(ispin))
 
     call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,dipole_int,Om1s,X1s,Y1s,Om2s,X2s,Y2s)
 
-    deallocate(Om1s,X1s,Y1s,Om2s,X2s,Y2s,TBab,TCab,TDab,TBaa,TCaa,TDaa)
+    deallocate(Om1s,X1s,Y1s,Om2s,X2s,Y2s,TBab,TCab,TDab,TBaa,TCaa,TDaa,Bpp,Cpp,Dpp) 
 
   end if
 
@@ -152,32 +172,42 @@ subroutine GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n
   !---------------------------------------!
 
     iblock = 3
-    EcRPA(ispin) = 0d0
- 
-    call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eT,ERI,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
- 
+
+    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,Bpp)
+                   call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVab,1d0,eT,ERI,Cpp)
+                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,1d0,eT,ERI,Dpp)
+
+    call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
+
+    deallocate(Bpp,Cpp,Dpp)
     allocate(TBab(nVVt,nOOt),TCab(nVVt,nVVt),TDab(nOOt,nOOt))
  
-    if(.not.TDA) call static_Tmatrix_B_pp(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOt,nVVt,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TBab)
-                 call static_Tmatrix_C_pp(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOt,nVVt,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TCab)
-                 call static_Tmatrix_D_pp(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOt,nVVt,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TDab)
+    if(.not.TDA_T) call GTpp_static_kernel_Bpp(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOt,nVVt,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TBab)
+                   call GTpp_static_kernel_Cpp(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOt,nVVt,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TCab)
+                   call GTpp_static_kernel_Dpp(ispin,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOt,nVVt,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TDab)
  
   !----------------------------------------!
   ! Compute T-matrix for alpha-alpha block !
   !----------------------------------------!
 
-    ispin  = 2
     iblock = 4
  
-    EcRPA(ispin) = 0d0
- 
-    call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eT,ERI,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
+    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,Bpp)
+                   call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVaa,1d0,eT,ERI,Cpp)
+                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,1d0,eT,ERI,Dpp)
+
+    call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
+
+    deallocate(Bpp,Cpp,Dpp)
     allocate(TBaa(nVVt,nOOt),TCaa(nVVt,nVVt),TDaa(nOOt,nOOt))
  
-    if(.not.TDA) call static_Tmatrix_B_pp(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOt,nVVt,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TBaa)
-                 call static_Tmatrix_C_pp(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOt,nVVt,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TCaa)
-                 call static_Tmatrix_D_pp(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOt,nVVt,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TDaa)
+    if(.not.TDA_T) call GTpp_static_kernel_Bpp(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOt,nVVt,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TBaa)
+                   call GTpp_static_kernel_Cpp(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOt,nVVt,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TCaa)
+                   call GTpp_static_kernel_Dpp(ispin,eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,nOOt,nVVt,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TDaa)
 
   !----------------------------------!
   ! pp/hh sectors for triplet states !
@@ -185,14 +215,23 @@ subroutine GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n
 
     EcBSE(ispin) = 0d0
 
-    allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt))
+    allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
+             Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
+
   
-    call linear_response_pp_BSE(ispin,TDA,.true.,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI, &
-                                TBaa-TBab,TCaa-TCab,TDaa-TDab,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcBSE(ispin))
+    if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
+                 call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
+                 call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
+
+    Bpp(:,:) = Bpp(:,:) + TBab(:,:) - TBaa(:,:)
+    Cpp(:,:) = Cpp(:,:) + TCab(:,:) - TCaa(:,:)
+    Dpp(:,:) = Dpp(:,:) + TDab(:,:) - TDaa(:,:)
+
+    call ppLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcBSE(ispin))
 
     call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,dipole_int,Om1t,X1t,Y1t,Om2t,X2t,Y2t)
 
-    deallocate(Om1t,X1t,Y1t,Om2t,X2t,Y2t,TBab,TCab,TDab,TBaa,TCaa,TDaa)
+    deallocate(Om1t,X1t,Y1t,Om2t,X2t,Y2t,TBab,TCab,TDab,TBaa,TCaa,TDaa,Bpp,Cpp,Dpp)
 
   end if
 
diff --git a/src/GT/static_Tmatrix_A.f90 b/src/GT/GTpp_static_kernel_Aph.f90
similarity index 83%
rename from src/GT/static_Tmatrix_A.f90
rename to src/GT/GTpp_static_kernel_Aph.f90
index f4c42b2..a80eee1 100644
--- a/src/GT/static_Tmatrix_A.f90
+++ b/src/GT/GTpp_static_kernel_Aph.f90
@@ -1,4 +1,4 @@
-subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,Omega2,rho2,TA)
+subroutine GTpp_static_kernel_Aph(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,Omega2,rho2,KA)
 
 ! Compute the OOVV block of the static T-matrix
 
@@ -30,12 +30,12 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,O
 
 ! Output variables
 
-  double precision,intent(out)  :: TA(nS,nS)
+  double precision,intent(out)  :: KA(nS,nS)
 
-  TA(:,:) = 0d0
+  KA(:,:) = 0d0
 
   jb = 0
-!$omp parallel do default(private) shared(TA,Omega1,Omega2,rho1,rho2,nO,nBas,nVV,nOO,chi,eps,eta,nC,nR,lambda)
+!$omp parallel do default(private) shared(KA,Omega1,Omega2,rho1,rho2,nO,nBas,nVV,nOO,chi,eps,eta,nC,nR,lambda)
   do j=nC+1,nO
     do b=nO+1,nBas-nR
       jb = (b-nO) + (j-1)*(nBas-nO) 
@@ -58,7 +58,7 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,O
             chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2)
           enddo
 
-          TA(ia,jb) = lambda*chi
+          KA(ia,jb) = lambda*chi
 
         enddo
       enddo
@@ -67,4 +67,4 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,O
 
 !$omp end parallel do
 
-end subroutine static_Tmatrix_A
+end subroutine 
diff --git a/src/GT/static_Tmatrix_B.f90 b/src/GT/GTpp_static_kernel_Bph.f90
similarity index 86%
rename from src/GT/static_Tmatrix_B.f90
rename to src/GT/GTpp_static_kernel_Bph.f90
index cbe2f47..34eafa6 100644
--- a/src/GT/static_Tmatrix_B.f90
+++ b/src/GT/GTpp_static_kernel_Bph.f90
@@ -1,4 +1,4 @@
-subroutine static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,Omega2,rho2,TB)
+subroutine GTpp_static_kernel_Bph(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,Omega2,rho2,KB)
 
 ! Compute the OVVO block of the static T-matrix
 
@@ -30,9 +30,9 @@ subroutine static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,O
 
 ! Output variables
 
-  double precision,intent(out)  :: TB(nS,nS)
+  double precision,intent(out)  :: KB(nS,nS)
 
-  TB(:,:) = 0d0
+  KB(:,:) = 0d0
 
   ia = 0
   do i=nC+1,nO
@@ -55,11 +55,11 @@ subroutine static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,O
             chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
           enddo
 
-          TB(ia,jb) = lambda*chi
+          KB(ia,jb) = lambda*chi
 
         enddo
       enddo
     enddo
   enddo
 
-end subroutine static_Tmatrix_B
+end subroutine 
diff --git a/src/GT/static_Tmatrix_B_pp.f90 b/src/GT/GTpp_static_kernel_Bpp.f90
similarity index 95%
rename from src/GT/static_Tmatrix_B_pp.f90
rename to src/GT/GTpp_static_kernel_Bpp.f90
index 4a4eeec..8865260 100644
--- a/src/GT/static_Tmatrix_B_pp.f90
+++ b/src/GT/GTpp_static_kernel_Bpp.f90
@@ -1,4 +1,4 @@
-subroutine static_Tmatrix_B_pp(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TB)
+subroutine GTpp_static_kernel_Bpp(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TB)
 
 ! Compute the VVOO block of the static T-matrix
 
@@ -148,4 +148,4 @@ subroutine static_Tmatrix_B_pp(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lamb
 
   end if
 
-end subroutine static_Tmatrix_B_pp
+end subroutine 
diff --git a/src/GT/static_Tmatrix_C_pp.f90 b/src/GT/GTpp_static_kernel_Cpp.f90
similarity index 95%
rename from src/GT/static_Tmatrix_C_pp.f90
rename to src/GT/GTpp_static_kernel_Cpp.f90
index 0151b9b..9f2ac33 100644
--- a/src/GT/static_Tmatrix_C_pp.f90
+++ b/src/GT/GTpp_static_kernel_Cpp.f90
@@ -1,4 +1,4 @@
-subroutine static_Tmatrix_C_pp(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TC)
+subroutine GTpp_static_kernel_Cpp(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TC)
 
 ! Compute the VVVV block of the static T-matrix
 
@@ -151,4 +151,4 @@ subroutine static_Tmatrix_C_pp(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lamb
 
   end if
 
-end subroutine static_Tmatrix_C_pp
+end subroutine 
diff --git a/src/GT/static_Tmatrix_D_pp.f90 b/src/GT/GTpp_static_kernel_Dpp.f90
similarity index 95%
rename from src/GT/static_Tmatrix_D_pp.f90
rename to src/GT/GTpp_static_kernel_Dpp.f90
index 30a59eb..84081aa 100644
--- a/src/GT/static_Tmatrix_D_pp.f90
+++ b/src/GT/GTpp_static_kernel_Dpp.f90
@@ -1,4 +1,4 @@
-subroutine static_Tmatrix_D_pp(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TD)
+subroutine GTpp_static_kernel_Dpp(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TD)
 
 ! Compute the OOOO block of the static T-matrix 
 
@@ -149,4 +149,4 @@ subroutine static_Tmatrix_D_pp(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lamb
 
   end if
 
-end subroutine static_Tmatrix_D_pp
+end subroutine 
diff --git a/src/GT/evGTpp.f90 b/src/GT/evGTpp.f90
index 1412d78..a58a4ef 100644
--- a/src/GT/evGTpp.f90
+++ b/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,nVVs,1d0,eHF,ERI_MO,Cpp)
-                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI_MO,Dpp)
+                   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,1d0,eGT,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,nVVt,1d0,eHF,ERI_MO,Cpp)
-                   call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI_MO,Dpp)
+                   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,1d0,eGT,ERI_MO,Dpp)
 
     call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
 
diff --git a/src/GW/GW_ppBSE.f90 b/src/GW/GW_ppBSE.f90
index 83a9dd2..f611d32 100644
--- a/src/GW/GW_ppBSE.f90
+++ b/src/GW/GW_ppBSE.f90
@@ -47,6 +47,10 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole
   double precision,allocatable  :: XmY_RPA(:,:)
   double precision,allocatable  :: rho_RPA(:,:,:)
 
+  double precision,allocatable  :: Bpp(:,:)
+  double precision,allocatable  :: Cpp(:,:)
+  double precision,allocatable  :: Dpp(:,:)
+
   double precision,allocatable  :: Om1(:)
   double precision,allocatable  :: X1(:,:)
   double precision,allocatable  :: Y1(:,:)
@@ -91,30 +95,35 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole
     write(*,*) 
 
     ispin = 1
-    EcBSE(ispin) = 0d0
 
     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), &
              WB(nVV,nOO),WC(nVV,nVV),WD(nOO,nOO))
 
+    ! Compute BSE excitation energies
+
     if(.not.TDA) call static_screening_WB_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WB)
     call static_screening_WC_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WC)
     call static_screening_WD_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WD)
 
-    ! Compute BSE excitation energies
 
-    call linear_response_pp_BSE(ispin,TDA,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, &
-                                Om1,X1,Y1,Om2,X2,Y2,EcBSE(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,eGW,ERI,Cpp)
+                 call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
 
-!   call print_excitation('pp-BSE (N+2)',ispin,nVV,Om1)
-!   call print_excitation('pp-BSE (N-2)',ispin,nOO,Om2)
+    Bpp(:,:) = Bpp(:,:) + WB(:,:)
+    Cpp(:,:) = Cpp(:,:) + WC(:,:)
+    Dpp(:,:) = Dpp(:,:) + WD(:,:)
+
+    call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
 
     call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
 
-    deallocate(Om1,X1,Y1,Om2,X2,Y2,WB,WC,WD)
+    deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,WB,WC,WD)
 
   end if
 
@@ -135,25 +144,31 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole
     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), &
              WB(nVV,nOO),WC(nVV,nVV),WD(nOO,nOO))
 
+    ! Compute BSE excitation energies
+
     if(.not.TDA) call static_screening_WB_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WB)
     call static_screening_WC_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WC)
     call static_screening_WD_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WD)
 
-    ! Compute BSE excitation energies
 
-    call linear_response_pp_BSE(ispin,TDA,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, &
-                                Om1,X1,Y1,Om2,X2,Y2,EcBSE(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,eGW,ERI,Cpp)
+                 call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
 
-!   call print_excitation('pp-BSE (N+2)',ispin,nVV,Om1)
-!   call print_excitation('pp-BSE (N-2)',ispin,nOO,Om2)
+    Bpp(:,:) = Bpp(:,:) + WB(:,:)
+    Cpp(:,:) = Cpp(:,:) + WC(:,:)
+    Dpp(:,:) = Dpp(:,:) + WD(:,:)
+
+    call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
 
     call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
 
-    deallocate(Om1,X1,Y1,Om2,X2,Y2,WB,WC,WD)
+    deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,WB,WC,WD)
 
   end if
 
diff --git a/src/LR/linear_response_pp_BSE.f90 b/src/LR/linear_response_pp_BSE.f90
deleted file mode 100644
index f336a4a..0000000
--- a/src/LR/linear_response_pp_BSE.f90
+++ /dev/null
@@ -1,121 +0,0 @@
-subroutine linear_response_pp_BSE(ispin,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,WB,WC,WD,Omega1,X1,Y1,Omega2,X2,Y2,EcBSE)
-
-! Compute the p-p channel of BSE
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  integer,intent(in)            :: ispin
-  logical,intent(in)            :: TDA
-  logical,intent(in)            :: BSE
-  integer,intent(in)            :: nBas
-  integer,intent(in)            :: nC
-  integer,intent(in)            :: nO
-  integer,intent(in)            :: nV
-  integer,intent(in)            :: nR
-  integer,intent(in)            :: nOO
-  integer,intent(in)            :: nVV
-  double precision,intent(in)   :: lambda
-  double precision,intent(in)   :: e(nBas)
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: WB(nVV,nOO)
-  double precision,intent(in)   :: WC(nVV,nVV)
-  double precision,intent(in)   :: WD(nOO,nOO)
-  
-! Local variables
-
-  integer                       :: ab,cd,ij,kl
-  integer                       :: p,q,r,s
-  double precision              :: trace_matrix
-  double precision              :: EcBSE1
-  double precision              :: EcBSE2
-  double precision,allocatable  :: B(:,:)
-  double precision,allocatable  :: C(:,:)
-  double precision,allocatable  :: D(:,:)
-  double precision,allocatable  :: M(:,:)
-  double precision,allocatable  :: Z(:,:)
-  double precision,allocatable  :: Omega(:)
-
-! Output variables
-
-  double precision,intent(out)  :: Omega1(nVV)
-  double precision,intent(out)  :: X1(nVV,nVV)
-  double precision,intent(out)  :: Y1(nOO,nVV)
-  double precision,intent(out)  :: Omega2(nOO)
-  double precision,intent(out)  :: X2(nVV,nOO)
-  double precision,intent(out)  :: Y2(nOO,nOO)
-  double precision,intent(out)  :: EcBSE
-
-! Memory allocation
-
-  allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO),M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),Omega(nOO+nVV))
-
-!-------------------------------------------------!
-! Solve the p-p eigenproblem                      !
-!-------------------------------------------------!
-!                                                 !
-!  |  C   B | | X1  X2 |   | w1  0  | | X1  X2 |  !
-!  |        | |        | = |        | |        |  !
-!  | -Bt -D | | Y1  Y2 |   | 0   w2 | | Y1  Y2 |  !
-!                                                 !
-!-------------------------------------------------!
-
-! Build B, C and D matrices for the pp channel
-
-  call ppLR_C(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C)
-  call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D)
-
-  if(BSE) then
-
-    C(:,:) = C(:,:) - WC(:,:)
-    D(:,:) = D(:,:) - WD(:,:)
-
-  end if
-
-  if(TDA) then
-
-    X1(:,:) = +C(:,:)
-    Y1(:,:) = 0d0
-    if(nVV > 0) call diagonalize_matrix(nVV,X1,Omega1)
- 
-    X2(:,:) = 0d0
-    Y2(:,:) = -D(:,:)
-    if(nOO > 0) call diagonalize_matrix(nOO,Y2,Omega2)
-
-  else
-
-    call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B)
-    if(BSE) B(:,:) = B(:,:) - WB(:,:)
-
-  ! Diagonal blocks 
-
-    M(    1:nVV    ,    1:nVV)     = + C(1:nVV,1:nVV)
-    M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - D(1:nOO,1:nOO)
-
-  ! Off-diagonal blocks
-
-    M(    1:nVV    ,nVV+1:nOO+nVV) = -           B(1:nVV,1:nOO)
-    M(nVV+1:nOO+nVV,    1:nVV)     = + transpose(B(1:nVV,1:nOO))
-
-  ! Diagonalize the p-p matrix
-
-    if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Omega,Z)
-
-  ! Split the various quantities in p-p and h-h parts
-
-    call sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
-
-  end if
-
-! Compute the BSE correlation energy
-
-  EcBSE = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )
-  EcBSE1 = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:))
-  EcBSE2 = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:))
-
-  if(abs(EcBSE - EcBSE1) > 1d-6 .or. abs(EcBSE - EcBSE2) > 1d-6) & 
-    print*,'!!! Issue in pp-BSE linear reponse calculation BSE1 != BSE2 !!!'
-
-end subroutine