fix root search in GTpp and plot routine

input/methods

@@ -13,7 +13,7 @@
 # G0F2* evGF2* qsGF2* G0F3 evGF3
   F     F      F      F    F
 # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
-  T     F     F     F        F      F
+  F     F     F     F        F      F
 # G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh
-  F       F       F       T      F      F
+  T       F       F       F      F      F
 # * unrestricted version available

src/GT/G0T0pp.f90

@@ -197,6 +197,9 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
 
   end if
 
+  call GTpp_plot_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, &
+                             Om1t,rho1t,Om2t,rho2t)
+
 !----------------------------------------------
 ! Dump results
 !----------------------------------------------

src/GT/GTeh_QP_graph.f90

@@ -1,5 +1,6 @@
 subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z)
 
+
 ! Compute the graphical solution of the QP equation
 
   implicit none
@@ -67,8 +68,8 @@ subroutine GTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eGT,Z)
 
     if(nIt == maxIt) then 
 
-      write(*,*) 'Newton root search has not converged!'
       eGT(p) = eGTlin(p)
+      write(*,*) 'Newton root search has not converged!'
 
     else
 

src/GT/GTpp_QP_graph.f90

@@ -1,10 +1,13 @@
 subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & 
                          Om1t,rho1t,Om2t,rho2t,eGTlin,eGT,Z)
 
+! Compute the graphical solution of the QP equation
+
   implicit none
   include 'parameters.h'
 
-! Iput variables
+! Input variables
+
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
@@ -23,6 +26,7 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s
   double precision,intent(in)   :: eGTlin(nBas)
   
 ! Local variables
+
   integer                       :: p
   integer                       :: nIt
   integer,parameter             :: maxIt = 64
@@ -41,6 +45,7 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s
   dSigC = 0d0
 
 ! Run Newton's algorithm to find the root
+
   do p=nC+1,nBas-nR
 
     write(*,*) '-----------------'
@@ -48,7 +53,6 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s
     write(*,*) '-----------------'
 
     w = eGTlin(p)
-    write(*,*) 'HERE', eGTlin(p), eHF(p)
     nIt = 0
     f = 1d0
     write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f
@@ -57,12 +61,11 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s
 
       nIt = nIt + 1
 
-       SigC  = GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
+      SigC  = GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGTlin,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
-       dSigC = GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
+      dSigC = GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGTlin,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
       f  = w - eHF(p) - SigC 
       df = 1d0/(1d0 - dSigC)
-    
+      w = w - df*f
-       w = w - f/df
 
       write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,df,f
 
@@ -80,6 +83,7 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s
 
       write(*,'(A32,F16.10)')   'Quasiparticle energy (eV)   ',eGT(p)*HaToeV
       write(*,*)
+
     end if
 
 

src/GT/GTpp_plot_self_energy.f90

@@ -0,0 +1,102 @@
+subroutine GTpp_plot_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, &
+                                 Om1t,rho1t,Om2t,rho2t)
+
+! Dump several GTpp quantities for external plotting
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nOOs,nOOt
+  integer,intent(in)            :: nVVs,nVVt
+
+  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: eGT(nBas)
+  double precision,intent(in)   :: Om1s(nVVs),Om1t(nVVt)
+  double precision,intent(in)   :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt)
+  double precision,intent(in)   :: Om2s(nOOs),Om2t(nOOt)
+  double precision,intent(in)   :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt)
+
+! Local variables
+
+  integer                       :: p,g
+  integer                       :: nGrid
+  double precision              :: wmin,wmax,dw
+  double precision,external     :: GTpp_SigC,GTpp_dSigC
+  double precision,allocatable  :: w(:)
+  double precision,allocatable  :: SigC(:,:)
+  double precision,allocatable  :: Z(:,:)
+  double precision,allocatable  :: S(:,:)
+
+! Construct grid
+
+  nGrid = 1000
+  allocate(w(nGrid),SigC(nBas,nGrid),Z(nBas,nGrid),S(nBas,nGrid))
+
+! Initialize 
+
+  SigC(:,:) = 0d0
+  Z(:,:)    = 0d0
+
+! Minimum and maximum frequency values
+
+  wmin = -5d0
+  wmax = +5d0
+  dw = (wmax - wmin)/dble(ngrid)
+
+  do g=1,nGrid
+    w(g) = wmin + dble(g)*dw
+  enddo
+
+! Occupied part of the self-energy and renormalization factor
+ 
+  do g=1,nGrid
+    do p=nC+1,nBas-nR
+
+      SigC(p,g) = GTpp_SigC(p,w(g),eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
+      Z(p,g)    = GTpp_dSigC(p,w(g),eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
+
+    end do
+  end do
+ 
+  Z(:,:) = 1d0/(1d0 + Z(:,:))
+
+! Compute spectral function
+
+  do g=1,nGrid
+    do p=nC+1,nBas-nR
+      S(p,g) = eta/((w(g) - eHF(p) - SigC(p,g))**2 + eta**2)
+    enddo
+  enddo
+
+  S(:,:) = S(:,:)/pi
+
+! Dump quantities in files as a function of w
+
+  open(unit=8 ,file='GTpp_SigC.dat')
+  open(unit=9 ,file='GTpp_freq.dat')
+  open(unit=10 ,file='GTpp_Z.dat')
+  open(unit=11 ,file='GTpp_A.dat')
+
+  do g=1,nGrid
+    write(8 ,*) w(g)*HaToeV,(SigC(p,g)*HaToeV,p=nC+1,nBas-nR)
+    write(9 ,*) w(g)*HaToeV,((w(g)-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,(S(p,g),p=nC+1,nBas-nR)
+  enddo
+
+! Closing files
+
+  close(unit=8)
+  close(unit=9)
+  close(unit=10)
+  close(unit=11)
+
+end subroutine