diff --git a/input/options b/input/options
index f2119d9..0b1a49b 100644
--- a/input/options
+++ b/input/options
@@ -9,9 +9,9 @@
 # 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      F   0.0 F     F 
+      256    0.00001 T    5      T   0.0 F     F 
 # GT: maxSCF thresh  DIIS n_diis lin eta TDA_T reg
-      256    0.00001 T    5      F   0.0 F     F  
+      256    0.00001 T    5      T   0.0 F     F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: phBSE phBSE2 ppBSE dBSE dTDA
diff --git a/src/GF/G0F2.f90 b/src/GF/G0F2.f90
index e32a9ca..f223287 100644
--- a/src/GF/G0F2.f90
+++ b/src/GF/G0F2.f90
@@ -35,7 +35,6 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
   double precision              :: Ec
   double precision              :: EcBSE(nspin)
   double precision,allocatable  :: eGF(:)
-  double precision,allocatable  :: eGFlin(:)
   double precision,allocatable  :: SigC(:)
   double precision,allocatable  :: Z(:)
 
@@ -49,7 +48,7 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
 
 ! Memory allocation
 
-  allocate(SigC(nBas),Z(nBas),eGF(nBas),eGFlin(nBas))
+  allocate(SigC(nBas),Z(nBas),eGF(nBas))
 
   if(linearize) then 
   
@@ -69,20 +68,17 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
     call GF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,eHF,eHF,ERI,SigC,Z)
 
   end if
-
-  eGFlin(:) = eHF(:) + Z(:)*SigC(:)
   
   if(linearize) then
 
-    eGF(:) = eGFlin(:)
+    eGF(:) = eHF(:) + Z(:)*SigC(:)
 
- else
+  else
 
     write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
     write(*,*)
 
-    call QP_graph_GF2(eta,nBas,nC,nO,nV,nR,eHF,eGFlin,ERI,eGF)
-
+    call GF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGF)
 
   end if
 
diff --git a/src/GF/QP_graph_GF2.f90 b/src/GF/GF2_QP_graph.f90
similarity index 77%
rename from src/GF/QP_graph_GF2.f90
rename to src/GF/GF2_QP_graph.f90
index 2e273cf..9899a46 100644
--- a/src/GF/QP_graph_GF2.f90
+++ b/src/GF/GF2_QP_graph.f90
@@ -1,4 +1,4 @@
-subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
+subroutine GF2_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,ERI,eGF)
 
 ! Compute the graphical solution of the GF2 QP equation
 
@@ -10,7 +10,6 @@ subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
   double precision,intent(in)   :: eta
   integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
   double precision,intent(in)   :: eHF(nBas)
-  double precision,intent(in)   :: eGF2lin(nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
 
 ! Local variables
@@ -19,14 +18,14 @@ subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
   integer                       :: nIt
   integer,parameter             :: maxIt = 64
   double precision,parameter    :: thresh = 1d-6
-  double precision,external     :: SigmaC_GF2,dSigmaC_GF2
+  double precision,external     :: GF2_SigC,GF2_dSigC
   double precision              :: sigC,dsigC
   double precision              :: f,df
   double precision              :: w
   
 ! Output variables
 
-  double precision,intent(out)  :: eGF2(nBas)
+  double precision,intent(out)  :: eGF(nBas)
 
 
 ! Run Newton's algorithm to find the root
@@ -37,7 +36,7 @@ subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
     write(*,'(A10,I3)') 'Orbital ',p
     write(*,*) '-----------------'
 
-    w = eGF2lin(p)
+    w = eHF(p)
     nIt = 0
     f = 1d0
     write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f
@@ -46,8 +45,8 @@ subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
     
       nIt = nIt + 1
 
-      sigC  =  SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
-      dsigC = dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
+      sigC  = GF2_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
+      dsigC = GF2_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
       f  = w - eHF(p) - sigC
       df = 1d0 - dsigC
     
@@ -64,9 +63,9 @@ subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
 
     else
 
-      eGF2(p) = w
+      eGF(p) = w
 
-      write(*,'(A32,F16.10)')   'Quasiparticle energy (eV)   ',eGF2(p)*HaToeV
+      write(*,'(A32,F16.10)')   'Quasiparticle energy (eV)   ',eGF(p)*HaToeV
       write(*,*)
 
    end if
diff --git a/src/GF/SigmaC_GF2.f90 b/src/GF/GF2_SigC.f90
similarity index 74%
rename from src/GF/SigmaC_GF2.f90
rename to src/GF/GF2_SigC.f90
index 0fa97aa..193f6fc 100644
--- a/src/GF/SigmaC_GF2.f90
+++ b/src/GF/GF2_SigC.f90
@@ -1,4 +1,4 @@
-double precision function SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
+double precision function GF2_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -19,7 +19,7 @@ double precision function SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
   integer                       :: i,j,a,b
   double precision              :: eps
 
-  SigmaC_GF2 = 0d0
+  GF2_SigC = 0d0
 
   ! Occupied part of the correlation self-energy
 
@@ -28,7 +28,7 @@ double precision function SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
       do a=nO+1,nBas-nR
 
         eps = w + eHF(a) - eHF(i) - eHF(j)
-        SigmaC_GF2 = SigmaC_GF2 + (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*eps/(eps**2 + eta**2)
+        GF2_SigC = GF2_SigC + (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*eps/(eps**2 + eta**2)
 
         end do
      end do
@@ -41,7 +41,7 @@ double precision function SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
       do b=nO+1,nBas-nR
 
         eps = w + eHF(i) - eHF(a) - eHF(b)
-        SigmaC_GF2 = SigmaC_GF2 + (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*eps/(eps**2 + eta**2)
+        GF2_SigC = GF2_SigC + (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*eps/(eps**2 + eta**2)
 
       end do
     end do
diff --git a/src/GF/dSigmaC_GF2.f90 b/src/GF/GF2_dSigC.f90
similarity index 72%
rename from src/GF/dSigmaC_GF2.f90
rename to src/GF/GF2_dSigC.f90
index e0c3109..6d24dfb 100644
--- a/src/GF/dSigmaC_GF2.f90
+++ b/src/GF/GF2_dSigC.f90
@@ -1,4 +1,4 @@
-double precision function dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
+double precision function GF2_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -21,7 +21,7 @@ double precision function dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
 
 ! Initialize 
 
-  dSigmaC_GF2 = 0d0
+  GF2_dSigC = 0d0
 
   ! Occupied part of the correlation self-energy
 
@@ -30,7 +30,7 @@ double precision function dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
       do a=nO+1,nBas-nR
 
         eps = w + eHF(a) - eHF(i) - eHF(j)
-        dSigmaC_GF2 = dSigmaC_GF2 - (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+        GF2_dSigC = GF2_dSigC - (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
       end do
     end do
@@ -43,7 +43,7 @@ double precision function dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
       do b=nO+1,nBas-nR
 
         eps = w + eHF(i) - eHF(a) - eHF(b)
-        dSigmaC_GF2 = dSigmaC_GF2 - (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+        GF2_dSigC = GF2_dSigC - (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
       end do
     end do
diff --git a/src/GF/evGF2.f90 b/src/GF/evGF2.f90
index 75feffa..bf85965 100644
--- a/src/GF/evGF2.f90
+++ b/src/GF/evGF2.f90
@@ -90,14 +90,19 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr
 
     end if
 
+    ! Solve the quasi-particle equation
+ 
     if(linearize) then
-
-      eGF(:) = eHF(:) + Z(:)*SigC(:)
-
-    else
-
+ 
       eGF(:) = eHF(:) + SigC(:)
-
+ 
+    else
+ 
+      write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
+      write(*,*)
+ 
+      call GF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGF)
+ 
     end if
 
     Conv = maxval(abs(eGF - eOld))
diff --git a/src/GT/G0T0eh.f90 b/src/GT/G0T0eh.f90
index bcd1a48..11b22f1 100644
--- a/src/GT/G0T0eh.f90
+++ b/src/GT/G0T0eh.f90
@@ -62,7 +62,6 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
   double precision,allocatable  :: rhoR(:,:,:,:)
 
   double precision,allocatable  :: eGT(:)
-  double precision,allocatable  :: eGTlin(:)
 
   double precision,allocatable  :: KA_sta(:,:)
   double precision,allocatable  :: KB_sta(:,:)
@@ -102,7 +101,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
 ! Memory allocation
 
   allocate(Aph(nS,nS),Bph(nS,nS),Sig(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & 
-           rhoL(nBas,nBas,nS,2),rhoR(nBas,nBas,nS,2),eGT(nBas),eGTlin(nBas))
+           rhoL(nBas,nBas,nS,2),rhoR(nBas,nBas,nS,2),eGT(nBas))
 
 !---------------------------------
 ! Compute (triplet) RPA screening 
@@ -142,8 +141,6 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
 ! Solve the quasi-particle equation !
 !-----------------------------------!
 
-  eGTlin(:) = eHF(:) + Z(:)*Sig(:)
-
   ! Linearized or graphical solution?
 
   if(linearize) then 
@@ -151,14 +148,14 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
     write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
     write(*,*)
 
-    eGT(:) = eGTlin(:)
+    eGT(:) = eHF(:) + Z(:)*Sig(:)
 
   else 
 
     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)
+    call GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eHF,eGT)
 
   end if
 
diff --git a/src/GT/GTeh_QP_graph.f90 b/src/GT/GTeh_QP_graph.f90
index 01193ea..f525298 100644
--- a/src/GT/GTeh_QP_graph.f90
+++ b/src/GT/GTeh_QP_graph.f90
@@ -51,8 +51,8 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT)
 
        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)
+       sigC  = GTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGTlin,Om,rhoL,rhoR)
+       dsigC = GTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGTlin,Om,rhoL,rhoR)
        f  = w - eHF(p)  - sigC 
        df = 1d0 - dsigC
     
diff --git a/src/GT/GTeh_SigC.f90 b/src/GT/GTeh_SigC.f90
index b3d1228..cf50115 100644
--- a/src/GT/GTeh_SigC.f90
+++ b/src/GT/GTeh_SigC.f90
@@ -35,7 +35,7 @@ double precision function GTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
   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)
+        num = rhoL(i,p,m,1)*rhoR(i,p,m,1)
         GTeh_SigC = GTeh_SigC + num*eps/(eps**2 + eta**2)
      enddo
   enddo
@@ -45,7 +45,7 @@ double precision function GTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
   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)
+        num = rhoL(p,a,m,1)*rhoR(p,a,m,1)
         GTeh_SigC = GTeh_SigC + num*eps/(eps**2 + eta**2)
      enddo
   enddo
diff --git a/src/GT/GTeh_dSigC.f90 b/src/GT/GTeh_dSigC.f90
index 3e91fb1..01412d5 100644
--- a/src/GT/GTeh_dSigC.f90
+++ b/src/GT/GTeh_dSigC.f90
@@ -35,7 +35,7 @@ double precision function GTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
   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)
+        num = rhoL(i,p,m,1)*rhoR(i,p,m,1)
         GTeh_dSigC = GTeh_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
      enddo
   enddo
@@ -45,7 +45,7 @@ double precision function GTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR)
   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)
+        num = rhoL(p,a,m,1)*rhoR(p,a,m,1)
         GTeh_dSigC = GTeh_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
      enddo
   enddo
diff --git a/src/GT/GTeh_excitation_density.f90 b/src/GT/GTeh_excitation_density.f90
index 9735355..e08062f 100644
--- a/src/GT/GTeh_excitation_density.f90
+++ b/src/GT/GTeh_excitation_density.f90
@@ -54,9 +54,9 @@ subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR)
                  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(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)
+                               + (2d0*ERI(b,p,q,j) - ERI(b,p,j,q))*X(jb,m) + (2d0*ERI(j,p,q,b) - ERI(j,p,b,q))*Y(jb,m)
                  rhoR(p,q,m,2) = rhoR(p,q,m,2) & 
-                               + (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)
+                               + (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)
 
               enddo
            enddo
diff --git a/src/GT/GTeh_self_energy_diag.f90 b/src/GT/GTeh_self_energy_diag.f90
index ae5ba81..753f5b3 100644
--- a/src/GT/GTeh_self_energy_diag.f90
+++ b/src/GT/GTeh_self_energy_diag.f90
@@ -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(p,a,m,1)*rhoR(p,a,m,2)
+        num = rhoL(p,a,m,2)*rhoR(p,a,m,1)
         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/evGTeh.f90 b/src/GT/evGTeh.f90
index 554a9c3..9e85485 100644
--- a/src/GT/evGTeh.f90
+++ b/src/GT/evGTeh.f90
@@ -148,6 +148,22 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
 
     eGT(:) = eHF(:) + Sig(:)
 
+    if(linearize) then
+
+       write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
+       write(*,*)
+
+       eGT(:) = eGT(:)
+
+    else
+
+       write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
+       write(*,*)
+
+       call GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eOld,eGT)
+
+    end if
+
     ! Convergence criteria
 
     Conv = maxval(abs(eGT - eOld))
diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90
index ff784d7..ba3e7b1 100644
--- a/src/GW/G0W0.f90
+++ b/src/GW/G0W0.f90
@@ -57,7 +57,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
   double precision,allocatable  :: rho(:,:,:)
 
   double precision,allocatable  :: eGW(:)
-  double precision,allocatable  :: eGWlin(:)
 
 ! Output variables
 
@@ -94,8 +93,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
 
 ! Memory allocation
 
-  allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & 
-           eGW(nBas),eGWlin(nBas))
+  allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),eGW(nBas))
 
 !-------------------!
 ! Compute screening !
@@ -133,8 +131,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
 ! Solve the quasi-particle equation !
 !-----------------------------------!
 
-  eGWlin(:) = eHF(:) + Z(:)*SigC(:)
-
   ! Linearized or graphical solution?
 
   if(linearize) then 
@@ -142,14 +138,14 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
     write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
     write(*,*)
 
-    eGW(:) = eGWlin(:)
+    eGW(:) = eHF(:) + Z(:)*SigC(:)
 
   else 
 
     write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
     write(*,*)
   
-    call GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGWlin,eGW,regularize)
+    call GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eHF,eGW)
 
   end if
 
@@ -168,7 +164,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
 
 ! Deallocate memory
 
-  deallocate(SigC,Z,Om,XpY,XmY,rho,eGWlin)
+  deallocate(SigC,Z,Om,XpY,XmY,rho)
 
 ! Perform BSE calculation
 
diff --git a/src/GW/GW_QP_graph.f90 b/src/GW/GW_QP_graph.f90
index 9c4ad13..d862284 100644
--- a/src/GW/GW_QP_graph.f90
+++ b/src/GW/GW_QP_graph.f90
@@ -1,4 +1,4 @@
-subroutine GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regularize)
+subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eGW)
 
 ! Compute the graphical solution of the QP equation
 
@@ -15,11 +15,10 @@ subroutine GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regulari
   integer,intent(in)            :: nS
   double precision,intent(in)   :: eta
   double precision,intent(in)   :: eHF(nBas)
-  double precision,intent(in)   :: Omega(nS)
+  double precision,intent(in)   :: Om(nS)
   double precision,intent(in)   :: rho(nBas,nBas,nS)
 
   double precision,intent(in)   :: eGWlin(nBas)
-  logical,intent(in)            :: regularize
 
 ! Local variables
 
@@ -53,8 +52,8 @@ subroutine GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW,regulari
     
       nIt = nIt + 1
 
-      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)
+      sigC  = GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho)
+      dsigC = GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eGWlin,Om,rho)
       f  = w - eHF(p) - SigC
       df = 1d0 - dsigC
     
diff --git a/src/GW/GW_SigC.f90 b/src/GW/GW_SigC.f90
index 371086c..0067c06 100644
--- a/src/GW/GW_SigC.f90
+++ b/src/GW/GW_SigC.f90
@@ -1,4 +1,4 @@
-double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,regularize)
+double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -19,65 +19,32 @@ double precision function GW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,regulariz
   double precision,intent(in)   :: e(nBas)
   double precision,intent(in)   :: Om(nS)
   double precision,intent(in)   :: rho(nBas,nBas,nS)
-  logical,intent(in)            :: regularize
 
 ! Local variables
 
-  integer                       :: i,a,jb
+  integer                       :: i,a,m
   double precision              :: eps
-  double precision              :: Dpijb,Dpajb
 
 ! Initialize 
 
   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) + 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) - 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
+! Occupied part of the correlation self-energy
 
-     ! 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) + 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) - Om(jb)
-     !       SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*(exp(-2d0*eta*Dpajb*Dpajb)/Dpajb)
-     !    enddo
-     ! enddo
-  
-  else
-   ! Occupied part of the correlation self-energy
-     do i=nC+1,nO
-        do jb=1,nS
-           eps = w - e(i) + Om(jb)
-           GW_SigC = GW_SigC + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2)
-        enddo
+  do i=nC+1,nO
+     do m=1,nS
+        eps = w - e(i) + Om(m)
+        GW_SigC = GW_SigC + 2d0*rho(p,i,m)**2*eps/(eps**2 + eta**2)
      enddo
-  ! Virtual part of the correlation self-energy
-     do a=nO+1,nBas-nR
-        do jb=1,nS
-           eps = w - e(a) - Om(jb)
-           GW_SigC = GW_SigC + 2d0*rho(p,a,jb)**2*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)
+        GW_SigC = GW_SigC + 2d0*rho(p,a,m)**2*eps/(eps**2 + eta**2)
      enddo
-  end if 
+  enddo
 
 end function 
diff --git a/src/GW/GW_dSigC.f90 b/src/GW/GW_dSigC.f90
index e5b18bf..e31d9fc 100644
--- a/src/GW/GW_dSigC.f90
+++ b/src/GW/GW_dSigC.f90
@@ -1,4 +1,4 @@
-double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,regularize)
+double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
 
 ! Compute the derivative of the correlation part of the self-energy
 
@@ -19,60 +19,32 @@ double precision function GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,regulari
   double precision,intent(in)   :: e(nBas)
   double precision,intent(in)   :: Om(nS)
   double precision,intent(in)   :: rho(nBas,nBas,nS)
-  logical,intent(in)            :: regularize
 
 ! Local variables
 
-  integer                       :: i,j,a,b,jb
+  integer                       :: i,a,m
   double precision              :: eps
-  double precision              :: Dpijb,Dpajb
 
 ! Initialize 
 
   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) + 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) - 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
+! Occupied part of the correlation self-energy
 
-  else
-   ! Occupied part of the correlation self-energy
-     do i=nC+1,nO
-        jb = 0
-        do j=nC+1,nO
-           do b=nO+1,nBas-nR
-              jb = jb + 1
-              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
+  do i=nC+1,nO
+     do m=1,nS
+       eps = w - e(i) + Om(m)
+       GW_dSigC = GW_dSigC - 2d0*rho(p,i,m)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
      enddo
+  enddo
 
 ! Virtual part of the correlation self-energy
-     do a=nO+1,nBas-nR
-        jb = 0
-        do j=nC+1,nO
-           do b=nO+1,nBas-nR
-              jb = jb + 1
-              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
+
+  do a=nO+1,nBas-nR
+     do m=1,nS
+       eps = w - e(a) - Om(m)
+       GW_dSigC = GW_dSigC - 2d0*rho(p,a,m)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
      enddo
-  end if
+  enddo
 
 end function 
diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90
index 2cd9f26..3df7fe7 100644
--- a/src/GW/evGW.f90
+++ b/src/GW/evGW.f90
@@ -145,23 +145,19 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
 
     ! Solve the quasi-particle equation
 
-    eGW(:) = eHF(:) + SigC(:)
-
-    ! Linearized or graphical solution?
-
     if(linearize) then 
  
        write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
        write(*,*)
 
-       eGW(:) = eGW(:)
+       eGW(:) = eHF(:) + SigC(:)
 
     else 
 
        write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
        write(*,*)
   
-       call GW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGW,eGW,regularize)
+       call GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eGW)
  
     end if