diff --git a/src/CC/static_screening.f90 b/src/CC/static_screening.f90
index 9e7f299..cf9f184 100644
--- a/src/CC/static_screening.f90
+++ b/src/CC/static_screening.f90
@@ -53,7 +53,7 @@ subroutine static_screening(nBas,nC,nO,nV,nR,eW,ERI,dbERI)
   TDA = .false.
 
   call linear_response(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,Om,XpY,XmY)
-  call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
+  call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
 
   do p=1,nBas
     do q=1,nBas
diff --git a/src/GW/Bethe_Salpeter.f90 b/src/GW/Bethe_Salpeter.f90
index 77c9360..2687f6d 100644
--- a/src/GW/Bethe_Salpeter.f90
+++ b/src/GW/Bethe_Salpeter.f90
@@ -68,7 +68,7 @@ subroutine Bethe_Salpeter(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBa
 
   call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
                        EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-  call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+  call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
   call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta)
   call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta)
@@ -167,4 +167,4 @@ subroutine Bethe_Salpeter(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBa
 
   end if
 
-end subroutine Bethe_Salpeter
+end subroutine 
diff --git a/src/GW/Bethe_Salpeter_pp.f90 b/src/GW/Bethe_Salpeter_pp.f90
index 0e7a69e..cb5ae31 100644
--- a/src/GW/Bethe_Salpeter_pp.f90
+++ b/src/GW/Bethe_Salpeter_pp.f90
@@ -67,7 +67,7 @@ subroutine Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,E
 
   call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
                        EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-  call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+  call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
 !-------------------
 ! Singlet manifold
diff --git a/src/GW/Bethe_Salpeter_pp_so.f90 b/src/GW/Bethe_Salpeter_pp_so.f90
index 16333ea..ab5f671 100644
--- a/src/GW/Bethe_Salpeter_pp_so.f90
+++ b/src/GW/Bethe_Salpeter_pp_so.f90
@@ -67,7 +67,7 @@ subroutine Bethe_Salpeter_pp_so(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
 
   call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
                        EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-  call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+  call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
   write(*,*) '****************'
   write(*,*) '*** SpinOrbs ***'
diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90
index b79d77f..f448df7 100644
--- a/src/GW/G0W0.f90
+++ b/src/GW/G0W0.f90
@@ -130,7 +130,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
 ! Compute spectral weights !
 !--------------------------!
 
-  call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
+  call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
 
 !------------------------!
 ! Compute GW self-energy !
diff --git a/src/GW/excitation_density.f90 b/src/GW/GW_excitation_density.f90
similarity index 90%
rename from src/GW/excitation_density.f90
rename to src/GW/GW_excitation_density.f90
index 0baa5fb..b00f3ca 100644
--- a/src/GW/excitation_density.f90
+++ b/src/GW/GW_excitation_density.f90
@@ -1,4 +1,4 @@
-subroutine excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
+subroutine GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
 
 ! Compute excitation densities
 
@@ -40,4 +40,4 @@ subroutine excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
   !$OMP END DO
   !$OMP END PARALLEL
 
-end subroutine excitation_density
+end subroutine 
diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90
index 7745817..03acb89 100644
--- a/src/GW/SRG_qsGW.f90
+++ b/src/GW/SRG_qsGW.f90
@@ -192,7 +192,7 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE
 
     call wall_time(tex1)
     
-    call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
+    call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
 
     call wall_time(tex2)
     tex=tex+tex2-tex1
diff --git a/src/GW/XBSE.f90 b/src/GW/XBSE.f90
index 9e8eb92..6b4bae0 100644
--- a/src/GW/XBSE.f90
+++ b/src/GW/XBSE.f90
@@ -64,7 +64,7 @@ subroutine XBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
 
   call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
                        EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-  call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+  call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
   call XBSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta,eW,eGW)
   call XBSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta,eW)
diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90
index d5bfbce..bbbf4fb 100644
--- a/src/GW/evGW.f90
+++ b/src/GW/evGW.f90
@@ -157,7 +157,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
 
    ! Compute spectral weights
 
-    call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
+    call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
 
     ! Compute correlation part of the self-energy 
 
diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90
index dd2e82a..d399f4c 100644
--- a/src/GW/qsGW.f90
+++ b/src/GW/qsGW.f90
@@ -184,7 +184,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
 
     ! Compute correlation part of the self-energy 
 
-    call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
+    call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
 
      if(regularize) then
 
diff --git a/src/GW/soBSE.f90 b/src/GW/soBSE.f90
index 28266df..f04a5b3 100644
--- a/src/GW/soBSE.f90
+++ b/src/GW/soBSE.f90
@@ -62,7 +62,7 @@ subroutine soBSE(TDA_W,TDA,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,eGW)
 
   call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
                        EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-  call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+  call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
   call static_screening_WA_so(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA_sta)
   call static_screening_WB_so(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB_sta)
@@ -119,4 +119,4 @@ subroutine soBSE(TDA_W,TDA,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,eGW)
 
 ! call matout(nO*nO,nV*nV,t)
 
-end subroutine soBSE
+end subroutine 
diff --git a/src/GW/soG0W0.f90 b/src/GW/soG0W0.f90
index 8d7f815..2d39b12 100644
--- a/src/GW/soG0W0.f90
+++ b/src/GW/soG0W0.f90
@@ -143,7 +143,7 @@ subroutine soG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,e
 ! Compute spectral weights !
 !--------------------------!
 
-  call excitation_density(nBas2,nC2,nO2,nR2,nS2,sERI,XpY_RPA,rho_RPA)
+  call GW_excitation_density(nBas2,nC2,nO2,nR2,nS2,sERI,XpY_RPA,rho_RPA)
 
 !------------------------!
 ! Compute GW self-energy !
diff --git a/src/GW/ufG0W0.f90 b/src/GW/ufG0W0.f90
index 9dc0b98..72887d2 100644
--- a/src/GW/ufG0W0.f90
+++ b/src/GW/ufG0W0.f90
@@ -222,7 +222,7 @@ subroutine ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,TDA_W)
      !--------------------------!
      ! Compute spectral weights !
      !--------------------------!
-     call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+     call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
      !---------!
      ! Block F !
diff --git a/src/RPA/ACFDT.f90 b/src/RPA/ACFDT.f90
index 3a55c37..9cc066e 100644
--- a/src/RPA/ACFDT.f90
+++ b/src/RPA/ACFDT.f90
@@ -73,7 +73,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB
 
   call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
                        EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-  call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+  call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
   call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA)
   call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB)
@@ -101,7 +101,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB
 
         call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, &
                              EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-        call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+        call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 !       call print_excitation('W^lambda:   ',isp_W,nS,OmRPA)
 
         call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA)
@@ -153,7 +153,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB
 
         call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, &
                              EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-        call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+        call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
         call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA)
         call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB)
diff --git a/src/RPA/ACFDT_cr.f90 b/src/RPA/ACFDT_cr.f90
index 4df5b02..ab08fd2 100644
--- a/src/RPA/ACFDT_cr.f90
+++ b/src/RPA/ACFDT_cr.f90
@@ -74,7 +74,7 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta
 
   call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
                        EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-  call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+  call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
   call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA)
   call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB)
@@ -102,7 +102,7 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta
 
         call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, &
                              EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-        call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+        call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 !       call print_excitation('W^lambda:   ',isp_W,nS,OmRPA)
 
         call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA)
@@ -154,7 +154,7 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta
 
         call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, &
                              EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-        call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+        call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
         call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA)
         call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB)