diff --git a/src/GW/complex_RGW_plot_self_energy.f90 b/src/GW/complex_RGW_plot_self_energy.f90
index b6ae2e9..9271eff 100644
--- a/src/GW/complex_RGW_plot_self_energy.f90
+++ b/src/GW/complex_RGW_plot_self_energy.f90
@@ -1,4 +1,4 @@
-subroutine complex_RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
+subroutine complex_RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,Re_eHF,Im_eHF,Re_eGW,Im_eGW,Om,rho)
 
 ! Dump several GW quantities for external plotting
 
@@ -14,10 +14,10 @@ subroutine complex_RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
   integer,intent(in)            :: nV
   integer,intent(in)            :: nR
   integer,intent(in)            :: nS
-  double precision,intent(in)   :: eHF(nBas)
-  double precision,intent(in)   :: eGW(nBas)
-  double precision,intent(in)   :: Om(nS)
-  double precision,intent(in)   :: rho(nBas,nBas,nS)
+  double precision,intent(in)   :: Re_eHF(nBas),Im_eHF(nBas)
+  double precision,intent(in)   :: Re_eGW(nBas),Im_eGW(nBas)
+  complex*16,intent(in)         :: Om(nS)
+  complex*16,intent(in)         :: rho(nBas,nBas,nS)
 
 ! Local variables
 
@@ -28,23 +28,27 @@ subroutine complex_RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
   double precision,allocatable  :: w(:)
   double precision,allocatable  :: ReSigC(:,:),ImSigC(:,:)
   double precision,allocatable  :: Z(:,:)
+  double precision,allocatable  :: ReDS(:,:),ImDS(:,:)
   double precision,allocatable  :: A(:,:)
 
 ! Construct grid
 
   nGrid = 5000
-  allocate(w(nGrid),ReSigC(nBas,nGrid),ImSigC(nBas,nGrid),Z(nBas,nGrid),A(nBas,nGrid))
+  allocate(w(nGrid),ReSigC(nBas,nGrid),ImSigC(nBas,nGrid),Z(nBas,nGrid),&
+          ReDS(nBas,nGrid),ImDS(nBas,nGrid),A(nBas,nGrid))
 
 ! Initialize 
 
   ReSigC(:,:) = 0d0
   ImSigC(:,:) = 0d0
   Z(:,:)      = 0d0
+  ReDS(:,:)      = 0d0
+  ImDS(:,:)      = 0d0
 
 ! Minimum and maximum frequency values
 
-  wmin = -5d0
-  wmax = +5d0
+  wmin = 0d0
+  wmax = +1d0
   dw = (wmax - wmin)/dble(ngrid)
 
   do g=1,nGrid
@@ -56,20 +60,19 @@ subroutine complex_RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
   do g=1,nGrid
     do p=nC+1,nBas-nR
 
-      ReSigC(p,g) = RGW_Re_SigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
-      ImSigC(p,g) = RGW_Im_SigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
-      Z(p,g)      = RGW_Re_dSigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
-
+      call complex_RGW_SigC_dSigC(p,eta,nBas,nC,nO,nV,nR,nS,w(g),0d0,&
+              Re_eGW,Im_eGW,Om,rho,ReSigC(p,g),ImSigC(p,g),ReDS(p,g),ImDS(p,g))
     end do
   end do
- 
-  Z(:,:) = 1d0/(1d0 + Z(:,:))
+  
+  Z(:,:) = (1d0-ReDS(:,:))/((1d0 - ReDS(:,:))**2 + ImDS(:,:)**2)
 
 ! Compute spectral function
 
   do g=1,nGrid
     do p=nC+1,nBas-nR
-      A(p,g) = abs(ImSigC(p,g))/((w(g) - eHF(p) - ReSigC(p,g))**2 + ImSigC(p,g)**2)
+      A(p,g) = abs(0d0 - Im_eHF(p) - ImSigC(p,g))&
+              /((w(g) - Re_eHF(p) - ReSigC(p,g))**2 + (0d0 - Im_eHF(p) - ImSigC(p,g))**2)
     end do
   end do
 
@@ -84,16 +87,18 @@ subroutine complex_RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
 
   do g=1,nGrid
     write(8 ,*) w(g)*HaToeV,(ReSigC(p,g)*HaToeV,p=nC+1,nBas-nR)
-    write(9 ,*) w(g)*HaToeV,((w(g)-eHF(p))*HaToeV,p=nC+1,nBas-nR)
+    write(9 ,*) w(g)*HaToeV,((w(g)-Re_eHF(p))*HaToeV,p=nC+1,nBas-nR)
     write(10,*) w(g)*HaToeV,(Z(p,g),p=nC+1,nBas-nR)
     write(11,*) w(g)*HaToeV,(A(p,g),p=nC+1,nBas-nR)
   end do
 
-! Closing files
+! Closing files and deallocation
 
   close(unit=8)
   close(unit=9)
   close(unit=10)
   close(unit=11)
+  deallocate(w,ReSigC,ImSigC,Z,&
+          ReDS,ImDS,A)
 
 end subroutine 
diff --git a/src/GW/complex_cRG0W0.f90 b/src/GW/complex_cRG0W0.f90
index df5e28a..d4dd612 100644
--- a/src/GW/complex_cRG0W0.f90
+++ b/src/GW/complex_cRG0W0.f90
@@ -44,7 +44,7 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
 ! Local variables
 
   logical                       :: print_W   = .false.
-  logical                       :: plot_self = .false.
+  logical                       :: plot_self = .true.
   logical                       :: dRPA_W
   integer                       :: isp_W
   integer                       :: p
@@ -156,7 +156,7 @@ end if
 
 ! Plot self-energy, renormalization factor, and spectral function
 !
-  if(plot_self) call complex_RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
+  if(plot_self) call complex_RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,Re_eHF,Im_eHF,Re_eGW,Im_eGW,Om,rho)
 !
 !! Cumulant expansion 
 !