diff --git a/src/CI/RCIS.f90 b/src/CI/RCIS.f90
index aed3196..6353191 100644
--- a/src/CI/RCIS.f90
+++ b/src/CI/RCIS.f90
@@ -48,7 +48,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in
   if(singlet) then
 
     ispin = 1
-    call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
+    call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
  
     if(dump_matrix) then
       print*,'CIS matrix (singlet state)'
@@ -84,7 +84,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in
   if(triplet) then
 
     ispin = 2
-    call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
+    call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
  
     if(dump_matrix) then
       print*,'CIS matrix (triplet state)'
diff --git a/src/GF/GGF2_phBSE.f90 b/src/GF/GGF2_phBSE.f90
index a705ecf..ef635a6 100644
--- a/src/GF/GGF2_phBSE.f90
+++ b/src/GF/GGF2_phBSE.f90
@@ -57,7 +57,7 @@ subroutine GGF2_phBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,E
 
   ! Compute phBSE@GF2 excitation energies
 
-  call phLR(TDA,nS,A_sta,B_sta,EcBSE,OmBSE,XpY,XmY)
+  call phGLR(TDA,nS,A_sta,B_sta,EcBSE,OmBSE,XpY,XmY)
   call print_excitation_energies('phBSE@GGF2','generalized',nS,OmBSE)
   call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
 
diff --git a/src/GF/RGF2_phBSE.f90 b/src/GF/RGF2_phBSE.f90
index 0e424b5..4d3ea5b 100644
--- a/src/GF/RGF2_phBSE.f90
+++ b/src/GF/RGF2_phBSE.f90
@@ -54,8 +54,8 @@ subroutine RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
     ispin = 1
     EcBSE(ispin) = 0d0
 
-                 call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
+                 call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta)
+    if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
 
     ! Compute static kernel
 
@@ -67,7 +67,7 @@ subroutine RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
 
     ! Compute phBSE@GF2 excitation energies
 
-    call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
+    call phRLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
     call print_excitation_energies('phBSE@GF2','singlet',nS,OmBSE)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
 
@@ -87,8 +87,8 @@ subroutine RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
     ispin = 2
     EcBSE(ispin) = 0d0
 
-                 call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
+                 call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta)
+    if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
 
     ! Compute static kernel
 
@@ -100,7 +100,7 @@ subroutine RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
 
     ! Compute phBSE@GF2 excitation energies
 
-    call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
+    call phRLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
     call print_excitation_energies('phBSE@GF2','triplet',nS,OmBSE)
     call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
 
diff --git a/src/GT/RG0T0eh.f90 b/src/GT/RG0T0eh.f90
index 4d7710c..ebc1c14 100644
--- a/src/GT/RG0T0eh.f90
+++ b/src/GT/RG0T0eh.f90
@@ -109,10 +109,10 @@ subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
 
   ispin = 2
 
-  call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
-  if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                 call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
+  if(.not.TDA_T) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-  call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+  call phRLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
   if(print_T) call print_excitation_energies('phRPA@RHF','triplet',nS,Om)
 
@@ -158,10 +158,10 @@ subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
 
 ! Compute the RPA correlation energy based on the G0T0eh quasiparticle energies
 
-  call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
-  if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                 call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
+  if(.not.TDA_T) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-  call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+  call phRLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
 !--------------!
 ! Dump results !
diff --git a/src/GT/RGTpp_phACFTD.f90 b/src/GT/RGTpp_phACFTD.f90
index 23425c6..0db7693 100644
--- a/src/GT/RGTpp_phACFTD.f90
+++ b/src/GT/RGTpp_phACFTD.f90
@@ -173,13 +173,13 @@ subroutine RGTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
 
       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)
+                   call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,Aph)
+      if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
  
                    Aph(:,:) = Aph(:,:) - TAt(:,:) - TAs(:,:)
       if(.not.TDA) Bph(:,:) = Bph(:,:) - TBt(:,:) - TBs(:,:)
  
-      call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
+      call phRLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
 
       call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
 
@@ -257,13 +257,13 @@ subroutine RGTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
 
       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)
+                   call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,Aph)
+      if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
                    Aph(:,:) = Aph(:,:) + TAs(:,:) - TAt(:,:)
       if(.not.TDA) Bph(:,:) = Bph(:,:) + TBs(:,:) - TBt(:,:)
 
-      call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
+      call phRLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
 
       call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
 
diff --git a/src/GT/RGTpp_phBSE.f90 b/src/GT/RGTpp_phBSE.f90
index dcdb6ff..62a03b7 100644
--- a/src/GT/RGTpp_phBSE.f90
+++ b/src/GT/RGTpp_phBSE.f90
@@ -137,13 +137,13 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
 
     ! Compute BSE singlet excitation energies
 
-                 call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                 call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
+    if(.not.TDA) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
                  Aph(:,:) = Aph(:,:) + TAt(:,:) ! TAt(:,:)
     if(.not.TDA) Bph(:,:) = Bph(:,:) + TBt(:,:) ! TBt(:,:)
 
-    call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
+    call phRLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
 
     call print_excitation_energies('phBSE@GTpp','singlet',nS,OmBSE)
     call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
@@ -167,13 +167,13 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n
 
     ! Compute BSE triplet excitation energies
 
-                 call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                 call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
+    if(.not.TDA) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
                  Aph(:,:) = Aph(:,:) + 1d0*TAt(:,:) - TAs(:,:)
     if(.not.TDA) Bph(:,:) = Bph(:,:) + 1d0*TBt(:,:) - TBs(:,:)
 
-    call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
+    call phRLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
 
     call print_excitation_energies('phBSE@GTpp','triplet',nS,OmBSE)
     call phLR_transition_vectors(.false.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
diff --git a/src/GT/evRGTeh.f90 b/src/GT/evRGTeh.f90
index 107e915..ff86982 100644
--- a/src/GT/evRGTeh.f90
+++ b/src/GT/evRGTeh.f90
@@ -123,10 +123,10 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
 
    ! Compute screening
 
-    call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
-    if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                   call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
+    if(.not.TDA_T) call phrLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-    call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+    call phRLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
     if(print_T) call print_excitation_energies('phRPA@RHF','triplet',nS,Om)
 
diff --git a/src/GT/qsRGTeh.f90 b/src/GT/qsRGTeh.f90
index 5e0263e..63374b6 100644
--- a/src/GT/qsRGTeh.f90
+++ b/src/GT/qsRGTeh.f90
@@ -197,10 +197,10 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
 
     ! Compute linear response
 
-    call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph)
-    if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
+                   call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph)
+    if(.not.TDA_T) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
 
-    call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+    call phRLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
     if(print_T) call print_excitation_energies('phRPA@RHF','triplet',nS,Om)
 
@@ -221,10 +221,6 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
     ! Solve the quasi-particle equation
 
     F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + Sigp(:,:)
-    if(nBas .ne. nOrb) then
-      call AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1))
-      call MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1))
-    endif
 
     ! Compute commutator and convergence criteria
 
@@ -241,17 +237,10 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
 
     ! Diagonalize Hamiltonian in AO basis
 
-    if(nBas .eq. nOrb) then
-      Fp = matmul(transpose(X), matmul(F, X))
-      cp(:,:) = Fp(:,:)
-      call diagonalize_matrix(nOrb, cp, eGT)
-      c = matmul(X, cp)
-    else
-      Fp = matmul(transpose(c), matmul(F, c))
-      cp(:,:) = Fp(:,:)
-      call diagonalize_matrix(nOrb, cp, eGT)
-      c = matmul(c, cp)
-    endif
+    Fp = matmul(transpose(X), matmul(F, X))
+    cp(:,:) = Fp(:,:)
+    call diagonalize_matrix(nOrb, cp, eGT)
+    c = matmul(X, cp)
 
     ! Compute new density matrix in the AO basis
 
diff --git a/src/GW/RG0W0.f90 b/src/GW/RG0W0.f90
index 6efcd42..5362af0 100644
--- a/src/GW/RG0W0.f90
+++ b/src/GW/RG0W0.f90
@@ -103,10 +103,10 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 ! Compute screening !
 !-------------------!
 
-                 call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
-  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                 call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
+  if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-  call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+  call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
   if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
 
@@ -160,10 +160,10 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
 ! Compute the RPA correlation energy
 
-                 call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
-  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                 call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
+  if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-  call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+  call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
 !--------------!
 ! Dump results !
diff --git a/src/GW/RGW_phACFDT.f90 b/src/GW/RGW_phACFDT.f90
index 154fa6f..4b2286a 100644
--- a/src/GW/RGW_phACFDT.f90
+++ b/src/GW/RGW_phACFDT.f90
@@ -93,10 +93,10 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
 
   isp_W = 1
 
-                 call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
-  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                 call phRLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
+  if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-  call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+  call phRLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
   call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
                call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA)
@@ -123,10 +123,10 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
 
       if(doXBS) then
 
-                       call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Aph)
-        if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+                       call phRLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Aph)
+        if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
    
-        call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+        call phRLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
         call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
    
                      call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
@@ -134,13 +134,13 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
 
       end if
 
-                   call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph)
-      if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+                   call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph)
+      if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
                    Aph(:,:) = Aph(:,:) + KA(:,:)
       if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
 
-      call phLR(TDA,nS,Aph,Bph,EcAC(ispin),Om,XpY,XmY)
+      call phRLR(TDA,nS,Aph,Bph,EcAC(ispin),Om,XpY,XmY)
 
       call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
 
@@ -180,10 +180,10 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
 
       if(doXBS) then
 
-                       call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Aph)
-        if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+                       call phRLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Aph)
+        if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
   
-        call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+        call phRLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
         call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
   
                      call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
@@ -191,13 +191,13 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
 
       end if  
 
-                   call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph)
-      if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+                   call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph)
+      if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
                    Aph(:,:) = Aph(:,:) + KA(:,:)
       if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
 
-      call phLR(TDA,nS,Aph,Bph,EcAC(ispin),Om,XpY,XmY)
+      call phRLR(TDA,nS,Aph,Bph,EcAC(ispin),Om,XpY,XmY)
 
       call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
 
diff --git a/src/GW/RGW_phBSE.f90 b/src/GW/RGW_phBSE.f90
index 4e1a926..6cf1d09 100644
--- a/src/GW/RGW_phBSE.f90
+++ b/src/GW/RGW_phBSE.f90
@@ -83,10 +83,10 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
   isp_W = 1
   EcRPA = 0d0
 
-                 call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
-  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                 call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
+  if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-  call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+  call phRLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
   call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
                call RGW_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta)
@@ -102,8 +102,8 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
 
     ! Compute BSE excitation energies
 
-                 call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                 call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
+    if(.not.TDA) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
     ! Second-order BSE static kernel
   
@@ -127,7 +127,7 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
                  Aph(:,:) = Aph(:,:) + KA_sta(:,:)
     if(.not.TDA) Bph(:,:) = Bph(:,:) + KB_sta(:,:)
 
-    call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
+    call phRLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
 
     call print_excitation_energies('phBSE@GW@RHF','singlet',nS,OmBSE)
     call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
@@ -160,13 +160,13 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
 
     ! Compute BSE excitation energies
 
-                 call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                 call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
+    if(.not.TDA) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
                  Aph(:,:) = Aph(:,:) + KA_sta(:,:)
     if(.not.TDA) Bph(:,:) = Bph(:,:) + KB_sta(:,:)
 
-    call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
+    call phRLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
 
     call print_excitation_energies('phBSE@GW@RHF','triplet',nS,OmBSE)
     call phLR_transition_vectors(.false.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
diff --git a/src/GW/RGW_ppBSE.f90 b/src/GW/RGW_ppBSE.f90
index c535657..2bd1767 100644
--- a/src/GW/RGW_ppBSE.f90
+++ b/src/GW/RGW_ppBSE.f90
@@ -98,11 +98,10 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
   allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nOrb,nOrb,nS), &
            Aph(nS,nS),Bph(nS,nS))
  
-  call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
+                 call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
+  if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
-
-  call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+  call phRLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
 !  call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmRPA,XpY_RPA,XmY_RPA)
 
   call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
diff --git a/src/GW/ccRG0W0.f90 b/src/GW/ccRG0W0.f90
index f5ea24e..1ee417d 100644
--- a/src/GW/ccRG0W0.f90
+++ b/src/GW/ccRG0W0.f90
@@ -93,10 +93,10 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,nS,ERI,ENuc,ERHF
 
   allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS))
 
-  call phLR_A(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
-  call phLR_B(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+  call phRLR_A(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
+  call phRLR_B(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-  call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+  call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
   call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
 
diff --git a/src/GW/evRGW.f90 b/src/GW/evRGW.f90
index 3d89073..90677ce 100644
--- a/src/GW/evRGW.f90
+++ b/src/GW/evRGW.f90
@@ -120,10 +120,10 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
    ! Compute screening
 
-    call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
-    if(.not.TDA_W) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                   call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
+    if(.not.TDA_W) call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-    call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+    call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
     ! Compute spectral weights
 
diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90
index 43f5d9f..7adf22b 100644
--- a/src/GW/qsRGW.f90
+++ b/src/GW/qsRGW.f90
@@ -202,10 +202,10 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
     ! Compute linear response
 
-    call phLR_A(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
-    if(.not.TDA_W) call phLR_B(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
+                   call phRLR_A(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
+    if(.not.TDA_W) call phRLR_B(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
 
-    call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+    call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
     if(print_W) call print_excitation_energies('phRPA@GW@RHF','singlet',nS,Om)
 
     call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho)
diff --git a/src/GW/ufRG0W0.f90 b/src/GW/ufRG0W0.f90
index 010dd45..d4b5304 100644
--- a/src/GW/ufRG0W0.f90
+++ b/src/GW/ufRG0W0.f90
@@ -96,10 +96,10 @@ subroutine ufRG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
 
     allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS))
  
-    call phLR_A(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
-    call phLR_B(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+    call phRLR_A(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
+    call phRLR_B(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
  
-    call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+    call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
     if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
  
diff --git a/src/GW/ufRGW.f90 b/src/GW/ufRGW.f90
index 970ef14..dc5194a 100644
--- a/src/GW/ufRGW.f90
+++ b/src/GW/ufRGW.f90
@@ -96,10 +96,10 @@ subroutine ufRGW(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
 
     ispin = 1
 
-    call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
-    call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+    call phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
+    call phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-    call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+    call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
     if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
 
diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90
index 7b7801f..780a001 100644
--- a/src/HF/RHF_search.f90
+++ b/src/HF/RHF_search.f90
@@ -126,8 +126,8 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
 
     ispin = 1
  
-    call phLR_A(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph)
-    call phLR_B(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
+    call phRLR_A(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph)
+    call phRLR_B(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
  
     AB(:,:) = Aph(:,:) + Bph(:,:)
  
diff --git a/src/HF/RHF_stability.f90 b/src/HF/RHF_stability.f90
index 154a4f9..b8b033b 100644
--- a/src/HF/RHF_stability.f90
+++ b/src/HF/RHF_stability.f90
@@ -38,8 +38,8 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
  
   ispin = 1
 
-  call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
-  call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
+  call phRLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
+  call phRLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
 
   AB(:,:) = A(:,:) + B(:,:)
 
@@ -115,8 +115,8 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
  
   ispin = 2
 
-  call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
-  call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
+  call phRLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
+  call phRLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
 
   AB(:,:) = A(:,:) + B(:,:)
 
diff --git a/src/LR/phLR.f90 b/src/LR/phRLR.f90
similarity index 97%
rename from src/LR/phLR.f90
rename to src/LR/phRLR.f90
index 6040ad7..4eb74bc 100644
--- a/src/LR/phLR.f90
+++ b/src/LR/phRLR.f90
@@ -1,4 +1,4 @@
-subroutine phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+subroutine phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
 ! Compute linear response
 
diff --git a/src/LR/phLR_A.f90 b/src/LR/phRLR_A.f90
similarity index 98%
rename from src/LR/phLR_A.f90
rename to src/LR/phRLR_A.f90
index 708799a..75e0352 100644
--- a/src/LR/phLR_A.f90
+++ b/src/LR/phRLR_A.f90
@@ -1,4 +1,4 @@
-subroutine phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
+subroutine phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
 
 ! Compute resonant block of the ph channel
 
diff --git a/src/LR/phLR_B.f90 b/src/LR/phRLR_B.f90
similarity index 97%
rename from src/LR/phLR_B.f90
rename to src/LR/phRLR_B.f90
index 9645952..19cc7cb 100644
--- a/src/LR/phLR_B.f90
+++ b/src/LR/phRLR_B.f90
@@ -1,4 +1,4 @@
-subroutine phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+subroutine phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
 ! Compute the coupling block of the ph channel
 
diff --git a/src/RPA/crACFDT.f90 b/src/RPA/crACFDT.f90
index 296d9f5..73127bc 100644
--- a/src/RPA/crACFDT.f90
+++ b/src/RPA/crACFDT.f90
@@ -87,10 +87,10 @@ subroutine crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,
  
       lambda = -rAC(iAC)
 
-      call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
-      if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+                   call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
+      if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
-      call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
+      call phRLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
 
       call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
 
@@ -128,10 +128,10 @@ subroutine crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,
  
       lambda = -rAC(iAC)
 
-      call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
-      if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+                   call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
+      if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
-      call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
+      call phRLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
 
       call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
 
diff --git a/src/RPA/crGRPA.f90 b/src/RPA/crGRPA.f90
index 33db3a8..c84166a 100644
--- a/src/RPA/crGRPA.f90
+++ b/src/RPA/crGRPA.f90
@@ -58,8 +58,8 @@ subroutine crGRPA(dotest,TDA,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
 
   allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS),Bph(nS,nS))
 
-               call phLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
-  if(.not.TDA) call phLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
+               call phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
+  if(.not.TDA) call phGLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
 
   call phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
   call print_excitation_energies('crRPA@GHF','spinorbital',nS,Om)
diff --git a/src/RPA/crRRPA.f90 b/src/RPA/crRRPA.f90
index e2c64df..537ba7f 100644
--- a/src/RPA/crRRPA.f90
+++ b/src/RPA/crRRPA.f90
@@ -69,10 +69,10 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
 
     ispin = 1
 
-    call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
+                 call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
+    if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
 
-    call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
+    call phRLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('crRPA@RHF','singlet',nS,Om)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
@@ -84,10 +84,10 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
 
     ispin = 2
 
-    call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
+                 call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
+    if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
 
-    call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
+    call phRLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('crRPA@RHF','triplet',nS,Om)
     call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
diff --git a/src/RPA/phACFDT.f90 b/src/RPA/phACFDT.f90
index c6e9e2d..d3e1e90 100644
--- a/src/RPA/phACFDT.f90
+++ b/src/RPA/phACFDT.f90
@@ -97,10 +97,10 @@ subroutine phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,
  
       lambda = rAC(iAC)
 
-      call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
-      if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+                   call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
+      if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
-      call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
+      call phRLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
 
       call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
 
@@ -138,10 +138,10 @@ subroutine phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,
  
       lambda = rAC(iAC)
 
-      call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
-      if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+                   call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
+      if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
-      call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
+      call phRLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY)
 
       call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
 
diff --git a/src/RPA/phGRPAx.f90 b/src/RPA/phGRPAx.f90
index acf4b0c..71ae6d9 100644
--- a/src/RPA/phGRPAx.f90
+++ b/src/RPA/phGRPAx.f90
@@ -59,7 +59,7 @@ subroutine phGRPAx(dotest,TDA,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
   allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS),Bph(nS,nS))
 
                call phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
-  if(.not.TDA) call phLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+  if(.not.TDA) call phGLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
   call phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
   call print_excitation_energies('phRPAx@GHF','spinorbital',nS,Om)
diff --git a/src/RPA/phRRPA.f90 b/src/RPA/phRRPA.f90
index 48f3c5c..dec2f79 100644
--- a/src/RPA/phRRPA.f90
+++ b/src/RPA/phRRPA.f90
@@ -73,17 +73,10 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
 
     ispin = 1
 
-    !call wall_time(t1)
-    call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+                 call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
+    if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
-    call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
-    !call wall_time(t2)
-    !print *, "wall time for dRPA on CPU (sec) = ", t2 - t1
-    !do ia = 1, nS
-    !  write(112, *) Om(ia) 
-    !enddo
-    !stop
+    call phRLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
@@ -95,10 +88,10 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
 
     ispin = 2
 
-    call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+                 call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
+    if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
-    call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
+    call phRLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('phRPA@RHF','triplet',nS,Om)
     call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
diff --git a/src/RPA/phRRPAx.f90 b/src/RPA/phRRPAx.f90
index 38a6243..b0a32ed 100644
--- a/src/RPA/phRRPAx.f90
+++ b/src/RPA/phRRPAx.f90
@@ -70,10 +70,10 @@ subroutine phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO
 
     ispin = 1
 
-    call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                 call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
+    if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-    call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
+    call phRLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('phRPAx@RHF','singlet',nS,Om)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
@@ -85,10 +85,10 @@ subroutine phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO
 
     ispin = 2
 
-    call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+                 call phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
+    if(.not.TDA) call phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
   
-    call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
+    call phRLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('phRPAx@RHF','triplet',nS,Om)
     call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)