diff --git a/input/methods b/input/methods
index da5c368..4ef876d 100644
--- a/input/methods
+++ b/input/methods
@@ -11,9 +11,9 @@
 # phRPA* phRPAx* crRPA ppRPA 
   F      F       F     F 
 # G0F2* evGF2* qsGF2* G0F3 evGF3
-  T     F      F      F    F
+  F     F      F      F    F
 # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
   F     F     F     F        F      F
 # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
-  F      F      F      F      F      F
+  F      F      F      T      F      F
 # * unrestricted version available
diff --git a/input/options b/input/options
index 174b10c..f2119d9 100644
--- a/input/options
+++ b/input/options
@@ -9,10 +9,10 @@
 # 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 TDA_W reg
-      256    0.00001 T    5      T   0.0 F     F 
+      256    0.00001 T    5      F   0.0 F     F 
 # GT: maxSCF thresh  DIIS n_diis lin eta TDA_T reg
-      256    0.00001 T    5      T   0.0 F     F  
+      256    0.00001 T    5      F   0.0 F     F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: phBSE phBSE2 ppBSE dBSE dTDA
-       F     F      T     T    T     
+       F     F      F     F    T     
diff --git a/src/GT/G0T0eh.f90 b/src/GT/G0T0eh.f90
index 10354cf..bcd1a48 100644
--- a/src/GT/G0T0eh.f90
+++ b/src/GT/G0T0eh.f90
@@ -155,9 +155,11 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
 
   else 
 
-    write(*,*) ' *** Root search not yet implemented in G0T0eh *** '
+    write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
     write(*,*)
 
+    call GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT)
+
   end if
 
 ! Compute the RPA correlation energy based on the G0T0eh quasiparticle energies
diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90
index dfdce23..d4bb341 100644
--- a/src/GT/G0T0pp.f90
+++ b/src/GT/G0T0pp.f90
@@ -192,7 +192,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
     write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
     write(*,*)
      
-   call QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,eGTlin,eGT)
+   call GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,eGTlin,eGT)
 
   end if
 
diff --git a/src/GT/GTeh_QP_graph.f90 b/src/GT/GTeh_QP_graph.f90
new file mode 100644
index 0000000..01193ea
--- /dev/null
+++ b/src/GT/GTeh_QP_graph.f90
@@ -0,0 +1,75 @@
+subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT)
+
+  implicit none
+  include 'parameters.h'
+
+! Iput variables
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  
+  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rhoL(nBas,nBas,nS,2)
+  double precision,intent(in)   :: rhoR(nBas,nBas,nS,2)
+
+  double precision,intent(in)   :: eGTlin(nBas)
+  
+! Local variables
+  integer                       :: p
+  integer                       :: nIt
+  integer,parameter             :: maxIt = 64
+  double precision,parameter    :: thresh = 1d-6
+  double precision,external     :: GTeh_SigC,GTeh_dSigC
+  double precision              :: sigC,dsigC
+  double precision              :: f,df
+  double precision              :: w
+  
+! Output variables
+  double precision,intent(out)  :: eGT(nBas)
+
+  sigC = 0d0
+  dsigC = 0d0
+
+! Run Newton's algorithm to find the root
+  do p=nC+1,nBas-nR
+
+    write(*,*) '-----------------'
+    write(*,'(A10,I3)') 'Orbital ',p
+    write(*,*) '-----------------'
+
+    w = eGTlin(p)
+    nIt = 0
+    f = 1d0
+    write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f
+
+    do while (abs(f) > thresh .and. nIt < maxIt)
+
+       nIt = nIt + 1
+
+       sigC  = GTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR)
+       dsigC = GTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR)
+       f  = w - eHF(p)  - sigC 
+       df = 1d0 - dsigC
+    
+       w = w - f/df
+
+       write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sigC
+
+    end do
+
+    if(nIt == maxIt) then 
+      write(*,*) 'Newton root search has not converged!'
+    else
+      eGT(p) = w
+      write(*,'(A32,F16.10)')   'Quasiparticle energy (eV)   ',eGT(p)*HaToeV
+      write(*,*)
+    end if
+
+  end do
+  
+end subroutine 
diff --git a/src/GT/GTeh_SigC.f90 b/src/GT/GTeh_SigC.f90
new file mode 100644
index 0000000..b3d1228
--- /dev/null
+++ b/src/GT/GTeh_SigC.f90
@@ -0,0 +1,53 @@
+double precision function GTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
+
+! Compute diagonal of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: p
+  double precision,intent(in)   :: w
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  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)   :: Om(nS)
+  double precision,intent(in)   :: rhoL(nBas,nBas,nS,2)
+  double precision,intent(in)   :: rhoR(nBas,nBas,nS,2)
+
+! Local variables
+
+  integer                       :: i,a,m
+  double precision              :: num,eps
+
+! Initialize 
+
+  GTeh_SigC = 0d0
+
+! Occupied part of the correlation self-energy
+
+  do i=nC+1,nO
+     do m=1,nS
+        eps = w - e(i) + Om(m)
+        num = rhoL(i,p,m,1)*rhoR(i,p,m,2)
+        GTeh_SigC = GTeh_SigC + num*eps/(eps**2 + eta**2)
+     enddo
+  enddo
+
+! Virtual part of the correlation self-energy
+
+  do a=nO+1,nBas-nR
+     do m=1,nS
+        eps = w - e(a) - Om(m)
+        num = rhoL(p,a,m,1)*rhoR(p,a,m,2)
+        GTeh_SigC = GTeh_SigC + num*eps/(eps**2 + eta**2)
+     enddo
+  enddo
+
+end function 
diff --git a/src/GT/GTeh_dSigC.f90 b/src/GT/GTeh_dSigC.f90
new file mode 100644
index 0000000..3e91fb1
--- /dev/null
+++ b/src/GT/GTeh_dSigC.f90
@@ -0,0 +1,53 @@
+double precision function GTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
+
+! Compute the derivative of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: p
+  double precision,intent(in)   :: w
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  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)   :: Om(nS)
+  double precision,intent(in)   :: rhoL(nBas,nBas,nS,2)
+  double precision,intent(in)   :: rhoR(nBas,nBas,nS,2)
+
+! Local variables
+
+  integer                       :: i,a,m
+  double precision              :: num,eps
+
+! Initialize 
+
+  GTeh_dSigC = 0d0
+
+! Occupied part of the correlation self-energy
+
+  do i=nC+1,nO
+     do m=1,nS
+        eps = w - e(i) + Om(m)
+        num = rhoL(i,p,m,1)*rhoR(i,p,m,2)
+        GTeh_dSigC = GTeh_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+     enddo
+  enddo
+
+! Virtual part of the correlation self-energy
+
+  do a=nO+1,nBas-nR
+     do m=1,nS
+        eps = w - e(a) - Om(m)
+        num = rhoL(p,a,m,1)*rhoR(p,a,m,2)
+        GTeh_dSigC = GTeh_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+     enddo
+  enddo
+
+end function 
diff --git a/src/GT/GTeh_excitation_density.f90 b/src/GT/GTeh_excitation_density.f90
index 553048f..9735355 100644
--- a/src/GT/GTeh_excitation_density.f90
+++ b/src/GT/GTeh_excitation_density.f90
@@ -51,12 +51,12 @@ subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR)
               do m=1,nS
 
                  rhoL(p,q,m,1) = rhoL(p,q,m,1) + ERI(p,j,b,q)*X(jb,m) + ERI(p,b,j,q)*Y(jb,m)
-                 rhoL(p,q,m,2) = rhoL(p,q,m,2) + ERI(b,p,q,j)*Y(jb,m) + ERI(j,p,q,b)*X(jb,m)
+                 rhoL(p,q,m,2) = rhoL(p,q,m,2) + ERI(p,j,b,q)*Y(jb,m) + ERI(p,b,j,q)*X(jb,m)
 
                  rhoR(p,q,m,1) = rhoR(p,q,m,1) & 
-                               + (2d0*ERI(b,q,p,j) - ERI(b,q,j,p))*X(jb,m) + (2d0*ERI(j,q,p,b) - ERI(j,q,b,p))*Y(jb,m)
+                               + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*X(jb,m) + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*Y(jb,m)
                  rhoR(p,q,m,2) = rhoR(p,q,m,2) & 
-                               + (2d0*ERI(b,p,q,j) - ERI(b,p,j,q))*Y(jb,m) + (2d0*ERI(j,p,q,b) - ERI(j,p,b,q))*X(jb,m)
+                               + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*Y(jb,m) + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*X(jb,m)
 
               enddo
            enddo
diff --git a/src/GT/GTeh_self_energy_diag.f90 b/src/GT/GTeh_self_energy_diag.f90
index 58bb41e..ae5ba81 100644
--- a/src/GT/GTeh_self_energy_diag.f90
+++ b/src/GT/GTeh_self_energy_diag.f90
@@ -46,7 +46,7 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig
       do m=1,nS
 
         eps = e(p) - e(i) + Om(m)
-        num = rhoL(p,i,m,2)*rhoR(i,p,m,2)
+        num = rhoL(i,p,m,1)*rhoR(i,p,m,2)
         Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
         Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
@@ -61,7 +61,7 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig
       do m=1,nS
 
         eps = e(p) - e(a) - Om(m)
-        num = rhoL(a,p,m,1)*rhoR(p,a,m,1)
+        num = rhoL(p,a,m,1)*rhoR(p,a,m,2)
         Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
         Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
diff --git a/src/GT/QP_graph_GT.f90 b/src/GT/GTpp_QP_graph.f90
similarity index 80%
rename from src/GT/QP_graph_GT.f90
rename to src/GT/GTpp_QP_graph.f90
index 04fd870..d79291a 100644
--- a/src/GT/QP_graph_GT.f90
+++ b/src/GT/GTpp_QP_graph.f90
@@ -1,4 +1,4 @@
-subroutine QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2,eGTlin,eGT)
+subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,eGTlin,eGT)
 
   implicit none
   include 'parameters.h'
@@ -14,9 +14,9 @@ subroutine QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2,
   
   double precision,intent(in)   :: eta
   double precision,intent(in)   :: eHF(nBas)
-  double precision,intent(in)   :: Omega1(nVV)
+  double precision,intent(in)   :: Om1(nVV)
   double precision,intent(in)   :: rho1(nBas,nBas,nVV)
-  double precision,intent(in)   :: Omega2(nOO)
+  double precision,intent(in)   :: Om2(nOO)
   double precision,intent(in)   :: rho2(nBas,nBas,nOO)
 
   double precision,intent(in)   :: eGTlin(nBas)
@@ -26,7 +26,7 @@ subroutine QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2,
   integer                       :: nIt
   integer,parameter             :: maxIt = 64
   double precision,parameter    :: thresh = 1d-6
-  double precision,external     :: SigmaC_GT,dSigmaC_GT
+  double precision,external     :: GTpp_SigC,GTpp_dSigC
   double precision              :: sigC,dsigC
   double precision              :: f,df
   double precision              :: w
@@ -54,8 +54,8 @@ subroutine QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2,
 
        nIt = nIt + 1
 
-       sigC  =  SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2)
-       dsigC = dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2)
+       sigC  = GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2)
+       dsigC = GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2)
        write (*,*) sigC
        f  = w - eHF(p)  - sigC 
        df = 1d0 - dsigC
@@ -76,4 +76,4 @@ subroutine QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Omega1,rho1,Omega2,rho2,
 
   end do
   
-end subroutine QP_graph_GT
+end subroutine 
diff --git a/src/GT/SigmaC_GT.f90 b/src/GT/GTpp_SigC.f90
similarity index 74%
rename from src/GT/SigmaC_GT.f90
rename to src/GT/GTpp_SigC.f90
index 82a1cda..c151d65 100644
--- a/src/GT/SigmaC_GT.f90
+++ b/src/GT/GTpp_SigC.f90
@@ -1,4 +1,4 @@
-double precision function SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2)
+double precision function GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -17,9 +17,9 @@ double precision function SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rh
   integer,intent(in)            :: nR
   integer,intent(in)            :: nOO,nVV
   double precision,intent(in)   :: e(nBas)
-  double precision,intent(in)   :: Omega1(nVV)
+  double precision,intent(in)   :: Om1(nVV)
   double precision,intent(in)   :: rho1(nBas,nBas,nVV)
-  double precision,intent(in)   :: Omega2(nOO)
+  double precision,intent(in)   :: Om2(nOO)
   double precision,intent(in)   :: rho2(nBas,nBas,nOO)
 
 ! Local variables
@@ -29,7 +29,7 @@ double precision function SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rh
 
 ! Initialize 
 
-  SigmaC_GT = 0d0
+  GTpp_SigC = 0d0
 
 !----------------------------------------------
 ! Occupied part of the T-matrix self-energy 
@@ -37,8 +37,8 @@ double precision function SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rh
 
   do i=nC+1,nO
     do cd=1,nVV
-      eps = w + e(i) - Omega1(cd)
-      SigmaC_GT = SigmaC_GT + rho1(p,i,cd)**2*eps/(eps**2 + eta**2)
+      eps = w + e(i) - Om1(cd)
+      GTpp_SigC = GTpp_SigC + rho1(p,i,cd)**2*eps/(eps**2 + eta**2)
     enddo
   enddo
 
@@ -48,8 +48,8 @@ double precision function SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rh
 
   do a=nO+1,nBas-nR
     do kl=1,nOO
-      eps = w + e(a) - Omega2(kl)
-      SigmaC_GT = SigmaC_GT + rho2(p,a,kl)**2*eps/(eps**2 + eta**2)
+      eps = w + e(a) - Om2(kl)
+      GTpp_SigC = GTpp_SigC + rho2(p,a,kl)**2*eps/(eps**2 + eta**2)
     enddo
   enddo
      
diff --git a/src/GT/dSigmaC_GT.f90 b/src/GT/GTpp_dSigC.f90
similarity index 75%
rename from src/GT/dSigmaC_GT.f90
rename to src/GT/GTpp_dSigC.f90
index ee59d99..cced4dd 100644
--- a/src/GT/dSigmaC_GT.f90
+++ b/src/GT/GTpp_dSigC.f90
@@ -1,4 +1,4 @@
-double precision function dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2)
+double precision function GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -18,9 +18,9 @@ double precision function dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,r
   integer,intent(in)            :: nOO
   integer,intent(in)            :: nVV
   double precision,intent(in)   :: e(nBas)
-  double precision,intent(in)   :: Omega1(nVV)
+  double precision,intent(in)   :: Om1(nVV)
   double precision,intent(in)   :: rho1(nBas,nBas,nVV)
-  double precision,intent(in)   :: Omega2(nOO)
+  double precision,intent(in)   :: Om2(nOO)
   double precision,intent(in)   :: rho2(nBas,nBas,nOO)
 
 ! Local variables
@@ -30,7 +30,7 @@ double precision function dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,r
 
 ! Initialize 
 
-  dSigmaC_GT = 0d0
+  GTpp_dSigC = 0d0
 
 !----------------------------------------------
 ! Occupied part of the T-matrix self-energy 
@@ -39,8 +39,8 @@ double precision function dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,r
   do i=nC+1,nO
     do cd=1,nVV
 
-      eps = w + e(i) - Omega1(cd)
-      dSigmaC_GT = dSigmaC_GT - rho1(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+      eps = w + e(i) - Om1(cd)
+      GTpp_dSigC=  GTpp_dSigC - rho1(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
     enddo
   enddo
@@ -52,8 +52,8 @@ double precision function dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,r
 
   do a=nO+1,nBas-nR
     do kl=1,nOO
-      eps = w + e(a) - Omega2(kl)
-      dSigmaC_GT = dSigmaC_GT - rho2(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+      eps = w + e(a) - Om2(kl)
+      GTpp_dSigC = GTpp_dSigC - rho2(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
     enddo
   enddo
      
diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90
index c630726..ff784d7 100644
--- a/src/GW/G0W0.f90
+++ b/src/GW/G0W0.f90
@@ -149,7 +149,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
     write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
     write(*,*)
   
-    call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGWlin,eGW,regularize)
+    call GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGWlin,eGW,regularize)
 
   end if
 
diff --git a/src/GW/QP_graph.f90 b/src/GW/GW_QP_graph.f90
similarity index 85%
rename from src/GW/QP_graph.f90
rename to src/GW/GW_QP_graph.f90
index 7c5bed2..9c4ad13 100644
--- a/src/GW/QP_graph.f90
+++ b/src/GW/GW_QP_graph.f90
@@ -1,4 +1,4 @@
-subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize)
+subroutine GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize)
 
 ! Compute the graphical solution of the QP equation
 
@@ -27,7 +27,7 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize)
   integer                       :: nIt
   integer,parameter             :: maxIt = 64
   double precision,parameter    :: thresh = 1d-6
-  double precision,external     :: SigmaC,dSigmaC
+  double precision,external     :: GW_SigC,GW_dSigC
   double precision              :: sigC,dsigC
   double precision              :: f,df
   double precision              :: w
@@ -53,8 +53,8 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize)
     
       nIt = nIt + 1
 
-      sigC  =  SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho,regularize)
-      dsigC = dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho,regularize)
+      sigC  = GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho,regularize)
+      dsigC = GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho,regularize)
       f  = w - eHF(p) - SigC
       df = 1d0 - dsigC
     
@@ -80,4 +80,4 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize)
           
   end do
 
-end subroutine QP_graph
+end subroutine 
diff --git a/src/GW/SigmaC.f90 b/src/GW/GW_SigC.f90
similarity index 67%
rename from src/GW/SigmaC.f90
rename to src/GW/GW_SigC.f90
index 70b32fb..371086c 100644
--- a/src/GW/SigmaC.f90
+++ b/src/GW/GW_SigC.f90
@@ -1,4 +1,4 @@
-double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regularize)
+double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,regularize)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -17,7 +17,7 @@ double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regular
   integer,intent(in)            :: nR
   integer,intent(in)            :: nS
   double precision,intent(in)   :: e(nBas)
-  double precision,intent(in)   :: Omega(nS)
+  double precision,intent(in)   :: Om(nS)
   double precision,intent(in)   :: rho(nBas,nBas,nS)
   logical,intent(in)            :: regularize
 
@@ -29,37 +29,36 @@ double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regular
 
 ! Initialize 
 
-  SigmaC = 0d0
-
+  GW_SigC = 0d0
 
   if (regularize) then
   ! Occupied part of the correlation self-energy
      do i=nC+1,nO
         do jb=1,nS
-           eps = w - e(i) + Omega(jb)
-           Dpijb = e(p) - e(i) + Omega(jb)
-           SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*(1d0-exp(-2d0*eta*Dpijb*Dpijb))/eps
+           eps = w - e(i) + Om(jb)
+           Dpijb = e(p) - e(i) + Om(jb)
+           GW_SigC = GW_SigC + 2d0*rho(p,i,jb)**2*(1d0-exp(-2d0*eta*Dpijb*Dpijb))/eps
         enddo
      enddo
   ! Virtual part of the correlation self-energy
      do a=nO+1,nBas-nR
         do jb=1,nS
-           eps = w - e(a) - Omega(jb)
-           Dpajb = e(p) - e(a) - Omega(jb)
-           SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*(1d0-exp(-2d0*eta*Dpajb*Dpajb))/eps
+           eps = w - e(a) - Om(jb)
+           Dpajb = e(p) - e(a) - Om(jb)
+           GW_SigC = GW_SigC + 2d0*rho(p,a,jb)**2*(1d0-exp(-2d0*eta*Dpajb*Dpajb))/eps
         enddo
      enddo
 
      ! We add the static SRG term in the self-energy directly
      ! do i=nC+1,nO
      !    do jb=1,nS
-     !       Dpijb = e(p) - e(i) + Omega(jb)
+     !       Dpijb = e(p) - e(i) + Om(jb)
      !       SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*(exp(-2d0*eta*Dpijb*Dpijb)/Dpijb)
      !    enddo
      ! enddo
      ! do a=nO+1,nBas-nR
      !    do jb=1,nS
-     !       Dpajb = e(p) - e(a) - Omega(jb)
+     !       Dpajb = e(p) - e(a) - Om(jb)
      !       SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*(exp(-2d0*eta*Dpajb*Dpajb)/Dpajb)
      !    enddo
      ! enddo
@@ -68,17 +67,17 @@ double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regular
    ! Occupied part of the correlation self-energy
      do i=nC+1,nO
         do jb=1,nS
-           eps = w - e(i) + Omega(jb)
-           SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2)
+           eps = w - e(i) + Om(jb)
+           GW_SigC = GW_SigC + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2)
         enddo
      enddo
   ! Virtual part of the correlation self-energy
      do a=nO+1,nBas-nR
         do jb=1,nS
-           eps = w - e(a) - Omega(jb)
-           SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2)
+           eps = w - e(a) - Om(jb)
+           GW_SigC = GW_SigC + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2)
         enddo
      enddo
   end if 
 
-end function SigmaC
+end function 
diff --git a/src/GW/dSigmaC.f90 b/src/GW/GW_dSigC.f90
similarity index 65%
rename from src/GW/dSigmaC.f90
rename to src/GW/GW_dSigC.f90
index ef28b58..e5b18bf 100644
--- a/src/GW/dSigmaC.f90
+++ b/src/GW/GW_dSigC.f90
@@ -1,4 +1,4 @@
-double precision function dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regularize)
+double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,regularize)
 
 ! Compute the derivative of the correlation part of the self-energy
 
@@ -17,7 +17,7 @@ double precision function dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regula
   integer,intent(in)            :: nR
   integer,intent(in)            :: nS
   double precision,intent(in)   :: e(nBas)
-  double precision,intent(in)   :: Omega(nS)
+  double precision,intent(in)   :: Om(nS)
   double precision,intent(in)   :: rho(nBas,nBas,nS)
   logical,intent(in)            :: regularize
 
@@ -29,23 +29,23 @@ double precision function dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regula
 
 ! Initialize 
 
-  dSigmaC = 0d0
+  GW_dSigC = 0d0
 
   if (regularize) then
   ! Occupied part of the correlation self-energy
      do i=nC+1,nO
         do jb=1,nS
-           eps = w - e(i) + Omega(jb)
-           Dpijb = e(p) - e(i) + Omega(jb)
-           dSigmaC = dSigmaC - 2d0*rho(p,i,jb)**2*(1d0-exp(-2*eta*Dpijb*Dpijb))/(eps**2)
+           eps = w - e(i) + Om(jb)
+           Dpijb = e(p) - e(i) + Om(jb)
+           GW_dSigC = GW_dSigC - 2d0*rho(p,i,jb)**2*(1d0-exp(-2*eta*Dpijb*Dpijb))/(eps**2)
         enddo
      enddo
   ! Virtual part of the correlation self-energy
      do a=nO+1,nBas-nR
         do jb=1,nS
-           eps = w - e(a) - Omega(jb)
-           Dpajb = e(p) - e(a) - Omega(jb)
-           dSigmaC = dSigmaC - 2d0*rho(p,a,jb)**2*(1d0-exp(-2*eta*Dpajb*Dpajb))/(eps**2)
+           eps = w - e(a) - Om(jb)
+           Dpajb = e(p) - e(a) - Om(jb)
+           GW_dSigC = GW_dSigC - 2d0*rho(p,a,jb)**2*(1d0-exp(-2*eta*Dpajb*Dpajb))/(eps**2)
         enddo
      enddo
 
@@ -56,8 +56,8 @@ double precision function dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regula
         do j=nC+1,nO
            do b=nO+1,nBas-nR
               jb = jb + 1
-              eps = w - e(i) + Omega(jb)
-              dSigmaC = dSigmaC - 2d0*rho(p,i,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+              eps = w - e(i) + Om(jb)
+              GW_dSigC = GW_dSigC - 2d0*rho(p,i,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
            enddo
         enddo
      enddo
@@ -68,11 +68,11 @@ double precision function dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regula
         do j=nC+1,nO
            do b=nO+1,nBas-nR
               jb = jb + 1
-              eps = w - e(a) - Omega(jb)
-              dSigmaC = dSigmaC - 2d0*rho(p,a,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+              eps = w - e(a) - Om(jb)
+              GW_dSigC = GW_dSigC - 2d0*rho(p,a,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
            enddo
         enddo
      enddo
   end if
 
-end function dSigmaC
+end function 
diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90
index e333956..2cd9f26 100644
--- a/src/GW/evGW.f90
+++ b/src/GW/evGW.f90
@@ -161,7 +161,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
        write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
        write(*,*)
   
-       call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGW,eGW,regularize)
+       call GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGW,eGW,regularize)
  
     end if