TABLE OF CONTENTS


ABINIT/m_dyson_solver [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dyson_solver

FUNCTION

  This module contains procedures to solve the Dyson equation to find QP energies.

COPYRIGHT

 Copyright (C) 2008-2022 ABINIT group (MG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 MODULE m_dyson_solver
23 
24  use defs_basis
25  use m_xmpi
26  use m_errors
27  use m_abicore
28  use m_dtfil
29 
30  use m_time,          only : timab
31  use m_gwdefs,        only : sigparams_t
32  use m_numeric_tools, only : linfit, pade, dpade, newrap_step
33  use m_io_tools,      only : open_file
34  use m_fstrings,      only : int2char10
35  use m_hide_lapack,   only : xheev
36  use m_bz_mesh,       only : kmesh_t
37  use m_sigma,         only : sigma_t
38 
39  implicit none
40 
41  private

m_dyson_solver/print_sigma_melems [ Functions ]

[ Top ] [ m_dyson_solver ] [ Functions ]

NAME

  print_sigma_melems

FUNCTION

  This routine prints the Hermitian and the non-hermitian part of the matrix
  elements of Sigma, as well as the individual contributions.
  The first 14x14 are printed to screen, and the full matrices are printed
  to files: sigma_melems_, sigma_nonH_melems_, sigma_Hart_melems_,
            sigma_x_melems, and sigma_c_melems

INPUTS

  ikcalc  : index of k-point
  ib1,ib2 : starting and ending band indices
  nsp     : no. of spin elements
  htotal  : Hermitianised matrix elements of Sigma
  hhartree : Hartree contribution to matrix elements
  sigxme  : Sigma_x contribution to matrix elements
  sigcme  : Sigma_c contribution to matrix elements
  prefil : prefix for output files.

OUTPUT

SOURCE

568 subroutine print_sigma_melems(ikcalc,ib1,ib2,nsp,htotal,hhartree,sigxme,sigcme,prefil)
569 
570 ! Arguments ------------------------------------
571  !scalars
572  integer,intent(in) :: ikcalc,ib1,ib2,nsp
573  character(len=*),intent(in) :: prefil
574  !arrays
575  complex(dpc),intent(in) :: htotal(ib1:ib2,ib1:ib2,nsp),hhartree(ib1:ib2,ib1:ib2,nsp)
576  complex(dpc),intent(in) :: sigxme(ib1:ib2,ib1:ib2,nsp),sigcme(ib1:ib2,ib1:ib2,nsp)
577 
578 ! Local variables ------------------------------
579  integer,parameter :: MAX_NCOLS = 14
580  integer :: isp,mc,mr,jj,ii,temp_unit,ount
581  character(len=10) :: sidx
582  character(len=500) :: msg
583  character(len=100) :: fmth,fmt1,fmt2,fmthh,kpt_index,fmtfile
584  character(len=fnlen) :: filename
585 
586 ! *************************************************************************
587 
588  if (nsp==3.or.nsp>4) then
589    ABI_ERROR('nsp has wrong value in print_sigma_melems')
590  end if
591 
592  ount = std_out
593 
594  mc = ib2-ib1+1; if (mc>MAX_NCOLS) mc = MAX_NCOLS
595  mr = mc
596 
597  write(fmthh,*)'(2(a),2(I2,a))'
598  write(fmth,*)'(7x,',mc,'(i2,8x))'
599  write(fmt1,*)'(3x,i2,',mc,'f10.5)'
600  write(fmt2,*)'(5x   ,',mc,'f10.5,a)'
601 
602  ! First print to screen
603  do isp=1,nsp
604    write(msg,'(a)') ''
605    call wrtout(ount,msg)
606    write(msg,fmthh) ch10,' Hermitianised matrix elements of Sigma (spin ',isp,' of ',nsp,'):'
607    call wrtout(ount,msg)
608    write(msg,fmth)(jj,jj=1,mc)
609    call wrtout(ount,msg) !header
610    do ii=ib1,ib1+mr-1
611      write(msg,fmt1)ii-ib1+1,DBLE(htotal(ii,ib1:(ib1+mc-1),isp))
612      call wrtout(ount,msg) !real part
613      write(msg,fmt2)  AIMAG(htotal(ii,ib1:(ib1+mc-1),isp)),ch10
614      call wrtout(ount,msg) !imag part
615    end do
616  end do !nsp
617 
618  write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output in _HTOTAL files"
619  call wrtout(ount,msg)
620 
621  do isp=1,nsp
622    write(msg,fmthh) ch10,' H_Hartree matrix elements (spin ',isp,' of ',nsp,'):'
623    call wrtout(ount,msg)
624    write(msg,fmth)(jj,jj=1,mc)
625    call wrtout(ount,msg) !header
626    do ii=ib1,ib1+mr-1
627      write(msg,fmt1)ii-ib1+1,DBLE(hhartree(ii,ib1:(ib1+mc-1),isp))
628      call wrtout(ount,msg) !real part
629      write(msg,fmt2)  AIMAG(hhartree(ii,ib1:(ib1+mc-1),isp)),ch10
630      call wrtout(ount,msg) !imag part
631    end do
632  end do !nsp
633 
634  write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output in _HHARTREE files"
635  call wrtout(ount,msg)
636 
637  do isp=1,nsp
638    write(msg,fmthh) ch10,' Sigma_x matrix elements (spin ',isp,' of ',nsp,'):'
639    call wrtout(ount,msg)
640    write(msg,fmth)(jj,jj=1,mc)
641    call wrtout(ount,msg) !header
642    do ii=ib1,ib1+mr-1
643      write(msg,fmt1)ii-ib1+1,DBLE(sigxme(ii,ib1:(ib1+mc-1),isp))
644      call wrtout(ount,msg) !real part
645      write(msg,fmt2)  AIMAG(sigxme(ii,ib1:(ib1+mc-1),isp)),ch10
646      call wrtout(ount,msg) !imag part
647    end do
648  end do !nsp
649 
650  write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output _SIGX files"
651  call wrtout(ount,msg)
652 
653  do isp=1,nsp
654    write(msg,fmthh) ch10,' Sigma_c matrix elements (spin ',isp,' of ',nsp,'):'
655    call wrtout(ount,msg)
656    write(msg,fmth)(jj,jj=1,mc)
657    call wrtout(ount,msg) !header
658    do ii=ib1,ib1+mr-1
659      write(msg,fmt1)ii-ib1+1,DBLE(sigcme(ii,ib1:(ib1+mc-1),isp))
660      call wrtout(ount,msg) !real part
661      write(msg,fmt2)  AIMAG(sigcme(ii,ib1:(ib1+mc-1),isp)),ch10
662      call wrtout(ount,msg) !imag part
663    end do
664  end do !nsp
665 
666  write(msg,'(a,i2,a)')" Max ",MAX_NCOLS," elements printed. Full matrix output _SIGC files"
667  call wrtout(ount,msg)
668 
669  ! Then print to file
670  ! Format is: row, column, value; with a blank space for each full
671  ! set of columns for easy plotting with the gnuplot splot command
672  write(fmtfile,*)'(3X,I6,2X,I6,',nsp,'(2(ES28.16E3,3x)))'
673 
674  call int2char10(ikcalc,sidx)
675  kpt_index = "_KPT"//TRIM(sidx)
676 
677  filename = TRIM(prefil)//'_HTOTAL'//TRIM(kpt_index)
678 
679  if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
680    ABI_ERROR(msg)
681  end if
682 
683  msg = '#   row    col.      Re(htotal(r,c)) Im(htotal(r,c))  for spin11   ... spin22 ... spin12 ... spin13'
684  call wrtout(temp_unit,msg)
685  do ii=ib1,ib2
686    do jj=ib1,ib2
687      write(msg,fmtfile) ii,jj,(htotal(jj,ii,isp),isp=1,nsp)
688      call wrtout(temp_unit,msg)
689    end do
690    call wrtout(temp_unit,"")
691  end do
692  close(temp_unit)
693 
694  filename = TRIM(prefil)//'_HHARTREE'//TRIM(kpt_index)
695  if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
696    ABI_ERROR(msg)
697  end if
698 
699  msg = '#   row    col.      Re(hhartree(r,c))  Im(hhartree(r,c)  for spin11   ... spin22 ... spin12 ... spin13'
700  call wrtout(temp_unit,msg)
701  do ii=ib1,ib2
702    do jj=ib1,ib2
703      write(msg,fmtfile) ii,jj,(hhartree(jj,ii,isp),isp=1,nsp)
704      call wrtout(temp_unit,msg)
705    end do
706    call wrtout(temp_unit,"")
707  end do
708  close(temp_unit)
709 
710  filename = TRIM(prefil)//'_SIGX'//TRIM(kpt_index)
711  if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
712    ABI_ERROR(msg)
713  end if
714 
715  write(msg,'(a)')'#   row    col.      Re(Sigx(r,c)) Im(Sigx(r,c) for spin11   ... spin22 ... spin12 ... spin13'
716  call wrtout(temp_unit,msg)
717  do ii=ib1,ib2
718    do jj=ib1,ib2
719      write(msg,fmtfile) ii,jj,(sigxme(jj,ii,isp),isp=1,nsp)
720      call wrtout(temp_unit,msg)
721    end do
722    call wrtout(temp_unit,"")
723  end do
724  close(temp_unit)
725 
726  filename = TRIM(prefil)//'_SIGC'//TRIM(kpt_index)
727  if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
728    ABI_ERROR(msg)
729  end if
730 
731  write(msg,'(a)')'#   row    col.      Re(Sigc(r,c)) Im(Sigc(r,c) for spin11   ... spin22 ... spin12 ... spin21'
732  call wrtout(temp_unit,msg)
733  do ii=ib1,ib2
734    do jj=ib1,ib2
735      write(msg,fmtfile) ii,jj,(sigcme(jj,ii,isp),isp=1,nsp)
736      call wrtout(temp_unit,msg)
737    end do
738    call wrtout(temp_unit,"")
739  end do
740 
741  close(temp_unit)
742 
743 end subroutine print_sigma_melems

m_dyson_solver/sigma_pade_eval [ Functions ]

[ Top ] [ m_dyson_solver ] [ Functions ]

NAME

  sigma_pade_eval

FUNCTION

  Evaluate the Pade' at the complex point `zz`. Return result in val and, optional,
  the derivative at zz in `dzdval`

SOURCE

788 subroutine sigma_pade_eval(self, zz, val, dzdval)
789 
790 !Arguments ------------------------------------
791  class(sigma_pade_t),intent(in) :: self
792  complex(dp),intent(in) :: zz
793  complex(dp),intent(out) :: val
794  complex(dp),optional,intent(out) :: dzdval
795 
796 ! *************************************************************************
797 
798  ! if zz in 2 or 3 quadrant, avoid branch cut in the complex plane using Sigma(-iw) = Sigma(iw)*.
799  if (real(zz) > zero) then
800  !if (real(zz) < zero) then
801    val = pade(self%npts, self%zmesh, self%sigc_cvals, zz)
802    if (present(dzdval)) dzdval = dpade(self%npts, self%zmesh, self%sigc_cvals, zz)
803  else
804    val = pade(self%npts, -self%zmesh, conjg(self%sigc_cvals), zz)
805    if (present(dzdval)) dzdval = dpade(self%npts, -self%zmesh, conjg(self%sigc_cvals), zz)
806  end if
807 
808 end subroutine sigma_pade_eval

m_dyson_solver/sigma_pade_init [ Functions ]

[ Top ] [ m_dyson_solver ] [ Functions ]

NAME

  sigma_pade_init

FUNCTION

  Initialize the Pade' from the `npts` values of Sigma_c(iw) given on the mesh `zmesh`.

SOURCE

757 subroutine sigma_pade_init(self, npts, zmesh, sigc_cvals, branch_cut)
758 
759 !Arguments ------------------------------------
760  class(sigma_pade_t),intent(out) :: self
761  integer,intent(in) :: npts
762  complex(dp),target,intent(in) :: zmesh(npts), sigc_cvals(npts)
763  character(len=*),intent(in) :: branch_cut
764 
765 ! *************************************************************************
766 
767  self%npts = npts
768  self%branch_cut = branch_cut
769 
770  self%zmesh => zmesh
771  self%sigc_cvals => sigc_cvals
772 
773 end subroutine sigma_pade_init

m_dyson_solver/sigma_pade_qp_solve [ Functions ]

[ Top ] [ m_dyson_solver ] [ Functions ]

NAME

  sigma_pade_qp_solve

FUNCTION

  Use the Pade' approximant and Newton-Rapson method  to solve the QP equation
  in the complex pane starting from the initial guess `z_guess`.

INPUTS

SOURCE

825 subroutine sigma_pade_qp_solve(self, e0, v_meanf, sigx, z_guess, zsc, msg, ierr)
826 
827 !Arguments ------------------------------------
828  class(sigma_pade_t),intent(in) :: self
829  real(dp),intent(in) :: e0, v_meanf, sigx
830  complex(dp),intent(in) :: z_guess
831  complex(dp),intent(out) :: zsc
832  integer,intent(out) :: ierr
833 
834 !Local variables-------------------------------
835 !scalars
836  integer :: iter
837  logical :: converged
838  complex(dpc) :: ctdpc, dct, dsigc, sigc
839  character(len=500) :: msg
840 
841 ! *************************************************************************
842 
843  ! Use Newton-Rapson to find the root of:
844  ! f(z) = e0 - zz + Sigma_xc(z) - v_meanf
845  ! f'(z) = -1 + Sigma_c'(z)
846 
847  iter = 0; converged = .FALSE.; ctdpc = cone
848  zsc = z_guess
849  do while (ABS(ctdpc) > NR_ABS_ROOT_ERR .or. iter < NR_MAX_NITER)
850    iter = iter + 1
851 
852    call self%eval(zsc, sigc, dzdval=dsigc)
853    ctdpc = e0 - v_meanf + sigx + sigc - zsc
854 
855    if (ABS(ctdpc) < NR_ABS_ROOT_ERR) then
856      converged=.TRUE.; EXIT
857    end if
858    dct = dsigc - one
859    zsc = newrap_step(zsc, ctdpc, dct)
860  end do
861 
862  ierr = 0; msg = ""
863  if (.not. converged) then
864    write(msg,'(a,i0,3a,f8.4,a,f8.4)')&
865      'Newton-Raphson method not converged after ',NR_MAX_NITER,' iterations. ',ch10,&
866      'Absolute Error = ',ABS(ctdpc),' > ',NR_ABS_ROOT_ERR
867    ierr = 1
868  end if
869 
870 end subroutine sigma_pade_qp_solve

m_dyson_solver/sigma_pade_t [ Types ]

[ Top ] [ m_dyson_solver ] [ Types ]

NAME

 sigma_pade_t

FUNCTION

  Small object to perform the analytic continuation with Pade' and
  find the QP solution with Newton-Rapson method

SOURCE

58  type, public :: sigma_pade_t
59 
60     integer :: npts
61     character(len=1) :: branch_cut
62     complex(dp),pointer :: zmesh(:) => null(), sigc_cvals(:) => null()
63 
64  contains
65 
66    procedure :: init => sigma_pade_init
67    ! Init object
68 
69    procedure :: eval => sigma_pade_eval
70    ! Eval self-energy and derivative
71 
72    procedure :: qp_solve => sigma_pade_qp_solve
73    ! Find the QP solution with Newton-Rapson method
74 
75  end type sigma_pade_t

m_dyson_solver/solve_dyson [ Functions ]

[ Top ] [ m_dyson_solver ] [ Functions ]

NAME

 solve_dyson

FUNCTION

  Solve the Dyson equation for the QP energies. Two different methods are coded:
  The first one is based on the standard perturbative approach in which the self-energy
  is linearly expanded around the previous single-particle energy (KS energy if one-shot)
  and the derivative is evaluated by finite differences.
  In the second method (AC), the values of the self-energy operator on the real axis are obtained
  by means of an analytic continuation based on the Pade extrapolation.

INPUTS

  ikcalc=Index of the considered k-point in the Sigp%kptgw2bz array.
  nomega_sigc=Number of frequencies used to evaluate the correlation part of Sigma.
  Sigp<sigparams_t>=Structure gathering parameters on the calculation of Sigma.
     %minbnd and %maxbnd= min and Max band index for GW correction (for this k-point)
     %gwcalctyp=Type of the GW calculation.
     %mbpt_sciss=Scissor energy
  Sr<sigma_t>=Structure containing the matrix elements of the self-energy INOUT
     %nbnds=Number of bands in G0.
     %nsppol=Number of independent spin polarizations.
     %nsig_ab=Numner of components in the self-energy operator.
     %nomega_r=Number of real frequencies for spectral function.
     %nomega4sd=Number of real frequencies used to evalute the derivative of Sigma.
     %nomega_i=Number of imaginary frequencies for AC.
     %omega_i=Purely imaginary frequencies for AC.
  Kmesh<kmesh_t>=Info on the K-mesh for the wavefunctions.
     %nkibz=Number of points in the IBZ
  sigcme=(nomega_sigc,ib1:ib2,ib1:ib2,nsppol)=Matrix elements of Sigma_c.
  qp_ene(nbnds,nkibz,nsppol)= KS or QP energies, only used in case of calculation with scissor operator.
  comm=MPI communicator.

OUTPUT

  Sr<sigma_t>=Structure containing the matrix elements of the self-energy:
     %sigxme(ib1:ib2,jkibz,nsspol)=Diagonal elements of Sigma_x
     %sigcmee0(ib1:ib2,jkibz,nsppol)=Matrix elements of Sigma_c at the initial energy E0.
     %dsigmee0(jb,ib1:ib2,nsppol)=Derivate of sigma at the energy E0.
     %ze0(ib1:ib2,jkibz,is)=Renormalization factor at the energy E0.
     %degw(ib1:ib2,jkibz,is)= QP correction  i.e DeltaE_GW=E-E0
     %egw(ib1:ib2,jkibz,is)=QP energy
     %sigmee(ib1:ib2,jkibz,is)=Self-energy evaluated at the QP energy.
     %sigcme (ib1:ib2,jkibz,io,is)= Sigma_c as a function of frequency.
     %sigxcme(ib1:ib2,jkibz,io,is)= Sigma_xc as a function of frequency.
     %sigcme4sd (ib1:ib2,jkibz,io,is)= Diagonal matrix elements of \Sigma_c  at frequencies around the KS eigenvalue
     %sigxcme4sd(ib1:ib2,jkibz,io,is)= Diagonal matrix elements of \Sigma_xc at frequencies around the KS eigenvalue
    where ib1 and ib2 are the band indeces included in the GW calculation for this k-point.

SOURCE

137 subroutine solve_dyson(ikcalc,minbnd,maxbnd,nomega_sigc,Sigp,Kmesh,sigcme,qp_ene,Sr,prtvol,Dtfil,comm)
138 
139 !Arguments ------------------------------------
140 !scalars
141  integer,intent(in) :: ikcalc,nomega_sigc,prtvol,minbnd,maxbnd,comm
142  type(kmesh_t),intent(in) :: Kmesh
143  type(Datafiles_type),intent(in) :: Dtfil
144  type(sigparams_t),intent(in) :: Sigp
145  type(sigma_t),intent(inout) :: Sr
146 !arrays
147  real(dp),intent(in) :: qp_ene(Sr%nbnds,Sr%nkibz,Sr%nsppol)
148  complex(dpc),intent(in) :: sigcme(nomega_sigc,minbnd:maxbnd,minbnd:maxbnd,Sigp%nsppol*Sigp%nsig_ab)
149 
150 !Local variables-------------------------------
151 !scalars
152  integer,parameter :: master=0
153  integer :: iab,ib1,ib2,ikbz_gw,io,spin,is_idx,isym,iter,itim,jb, ie0
154  integer :: sk_ibz,kb,ld_matrix,mod10,nsploop,my_rank
155  real(dp) :: alpha,beta,smrt
156  complex(dpc) :: ctdpc,dct,dsigc,sigc,zz,phase
157  logical :: converged,ltest
158  character(len=500) :: msg
159  !type(sigma_pade_t) :: spade
160 !arrays
161  real(dp) :: kbz_gw(3),tsec(2)
162  real(dp),allocatable :: e0pde(:),eig(:),scme(:)
163  complex(dpc),allocatable :: hdp(:,:),tmpcdp(:),hhartree(:,:,:),htotal(:,:,:),h_tmp1(:,:),h_tmp2(:,:)
164 
165 ! *************************************************************************
166 
167  DBG_ENTER("COLL")
168 
169  call timab(490,1,tsec) ! csigme(Dyson)
170 
171  my_rank = xmpi_comm_rank(comm)
172 
173  mod10=MOD(Sigp%gwcalctyp,10)
174 
175  ltest=(nomega_sigc==Sr%nomega_r+Sr%nomega4sd)
176  if (mod10==1) ltest=(nomega_sigc==Sr%nomega_i)
177  ABI_CHECK(ltest,'Wrong number of frequencies')
178 
179  ! Index of the KS or QP energy.
180  !ioe0j=Sr%nomega4sd/2+1
181 
182  ! min and Max band index for GW corrections (for this k-point).
183  ib1=MINVAL(Sigp%minbnd(ikcalc,:))
184  ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
185 
186  ! Find the index of the k-point for sigma in the IBZ array.
187  ikbz_gw=Sigp%kptgw2bz(ikcalc)
188  call kmesh%get_BZ_item(ikbz_gw,kbz_gw,sk_ibz,isym,itim,phase)
189 
190  sigc=czero; dsigc=czero
191 
192  ! ===========================================================
193  ! ==== Solve the Dyson Equation and store results in Sr% ====
194  ! ===========================================================
195 
196  if (mod10 /= 1) then
197    ! ===============================
198    ! ==== Perturbative approach ====
199    ! ===============================
200 
201    ! Index of the KS or QP energy in sigme_tmp
202    ie0 = sr%nomega_r + Sr%nomega4sd/2+1
203 
204    do spin=1,Sr%nsppol
205      do jb=ib1,ib2
206        ! === Get matrix elements of Sigma_c at energy E0 ===
207        ! SigC(w) is linearly interpolated and the slope alpha is assumed as dSigC/dE
208        do iab=1,Sr%nsig_ab
209          is_idx=spin; if (Sr%nsig_ab>1) is_idx=iab
210 
211          Sr%sigcmee0(jb,sk_ibz,is_idx) = sigcme(ie0,jb,jb,is_idx)
212 
213          ABI_MALLOC(scme,(Sr%nomega4sd))
214          ABI_MALLOC(e0pde,(Sr%nomega4sd))
215          e0pde(:) = Sr%omega4sd(jb,sk_ibz,:,spin)
216          scme(:)  = REAL(sigcme(Sr%nomega_r+1:Sr%nomega_r+Sr%nomega4sd,jb,jb,is_idx))
217 
218          if (Sr%nomega4sd==1) then
219            smrt = zero; alpha = zero
220          else
221            smrt = linfit(Sr%nomega4sd,e0pde(:),scme(:),alpha,beta)
222          end if
223 
224          if (smrt>0.1/Ha_eV) then
225            write(msg,'(3a,i0,a,i0,2a,2(f22.15,2a))')&
226              'Values of Re Sig_c are not linear ',ch10,&
227              'band index = ',jb,' spin|component = ',is_idx,ch10,&
228              'root mean square= ',smrt,ch10,&
229              'estimated slope = ',alpha,ch10,&
230              'Omega [eV] SigC [eV]'
231            ABI_WARNING(msg)
232            do io=1,Sr%nomega4sd
233              write(msg,'(2f8.4)')e0pde(io)*Ha_eV,scme(io)*Ha_eV
234              call wrtout(std_out,msg)
235            end do
236          end if
237 
238          ABI_FREE(scme)
239          ABI_FREE(e0pde)
240          !
241          ! === Evaluate renormalization factor and QP correction ===
242          ! * Z=(1-dSigma/domega(E0))^-1
243          ! * DeltaE_GW=E-E0= (Sigma(E0)-V_xc)/(1-dSigma/domega)
244          ! * If nspinor==2, this part is done at the end.
245          !
246          Sr%dsigmee0(jb,sk_ibz,is_idx)=CMPLX(alpha,zero)
247 
248          if (Sr%nsig_ab==1) then
249            Sr%ze0(jb,sk_ibz,spin)=one/(one-Sr%dsigmee0(jb,sk_ibz,spin))
250 
251            if (ABS(Sigp%mbpt_sciss) < tol6) then
252              Sr%degw(jb,sk_ibz,spin) = Sr%ze0(jb,sk_ibz,spin) * &
253                (Sr%sigxme(jb,sk_ibz,spin) + Sr%sigcmee0(jb,sk_ibz,spin) - Sr%e0(jb,sk_ibz,spin) + &
254                 Sr%hhartree(jb,jb,sk_ibz,spin))
255 
256              Sr%egw(jb,sk_ibz,spin) = Sr%e0(jb,sk_ibz,spin) + Sr%degw(jb,sk_ibz,spin)
257 
258              ! Estimate Sigma at the QP-energy: Sigma(E_qp)=Sigma(E0)+(E_qp-E0)*dSigma/dE
259              Sr%sigmee(jb,sk_ibz,spin) = &
260                Sr%sigxme(jb,sk_ibz,spin)+Sr%sigcmee0(jb,sk_ibz,spin)+Sr%degw(jb,sk_ibz,spin)*Sr%dsigmee0(jb,sk_ibz,spin)
261 
262            else
263              ! If GW+scissor: e0 is replaced by qp_ene which contains the updated energy eigenvalue
264              Sr%degw(jb,sk_ibz,spin)= Sr%ze0(jb,sk_ibz,spin) * &
265                (Sr%sigxme(jb,sk_ibz,spin) + Sr%sigcmee0(jb,sk_ibz,spin) - qp_ene(jb,sk_ibz,spin) + &
266                 Sr%hhartree(jb,jb,sk_ibz,spin))
267 
268              Sr%egw(jb,sk_ibz,spin) = qp_ene(jb,sk_ibz,spin) + Sr%degw(jb,sk_ibz,spin)
269 
270              ! Estimate Sigma at the QP-energy: Sigma(E_qp)=Sigma(E0)+(E_qp-E0)*dSigma/dE
271              Sr%sigmee(jb,sk_ibz,spin)= &
272                Sr%sigxme(jb,sk_ibz,spin) + Sr%sigcmee0(jb,sk_ibz,spin) + &
273                Sr%degw(jb,sk_ibz,spin) * Sr%dsigmee0(jb,sk_ibz,spin)
274 
275              ! RS: In the output, the gw corr with respect to e0 without mbpt_sciss is reported.
276              Sr%degw(jb,sk_ibz,spin) = Sr%egw(jb,sk_ibz,spin) - Sr%e0(jb,sk_ibz,spin)
277            end if
278          end if !Sigp%nsig_ab==1
279 
280          ! Spectrum of Sigma
281          do io=1,Sr%nomega_r
282            Sr%sigcme (jb,sk_ibz,io,is_idx)= sigcme(io,jb,jb,is_idx)
283            Sr%sigxcme(jb,sk_ibz,io,is_idx)= Sr%sigxme(jb,sk_ibz,is_idx)+Sr%sigcme(jb,sk_ibz,io,is_idx)
284          end do
285          do io=1,Sr%nomega4sd
286            Sr%sigcme4sd (jb,sk_ibz,io,is_idx)= sigcme(Sr%nomega_r+io,jb,jb,is_idx)
287            Sr%sigxcme4sd(jb,sk_ibz,io,is_idx)= Sr%sigxme(jb,sk_ibz,is_idx)+Sr%sigcme4sd(jb,sk_ibz,io,is_idx)
288          end do
289        end do !iab
290 
291        if (Sr%nsig_ab > 1) then
292          ABI_CHECK(ABS(Sigp%mbpt_sciss)<0.1d-4,'Scissor with spinor not coded')
293          !TODO this should be allocated with nsppol, recheck this part
294 
295          ! Evaluate renormalization factor and QP correction.
296          ! Z=(1-dSigma/domega(E0))^-1
297          ! DeltaE_GW=E-E0= (Sigma(E0)-V_xc)/(1-dSigma/domega)
298          !write(std_out,'(a,i2,10f8.3)')' Correlation',jb,Sr%sigcmee0(jb,sk_ibz,:)*Ha_eV,SUM(Sr%sigcmee0(jb,sk_ibz,:))*Ha_eV
299 
300          Sr%ze0 (jb,sk_ibz,1) = one/(one-SUM(Sr%dsigmee0(jb,sk_ibz,:)))
301 
302          Sr%degw(jb,sk_ibz,1) = Sr%ze0(jb,sk_ibz,1) * &
303            (SUM(Sr%sigxme(jb,sk_ibz,:)+Sr%sigcmee0(jb,sk_ibz,:)+Sr%hhartree(jb,jb,sk_ibz,:))-Sr%e0(jb,sk_ibz,1))
304 
305          Sr%egw(jb,sk_ibz,1)=Sr%e0(jb,sk_ibz,1)+Sr%degw(jb,sk_ibz,1)
306 
307          ! Estimate Sigma at the QP-energy.
308          do iab=1,Sr%nsig_ab
309           Sr%sigmee(jb,sk_ibz,iab)= &
310             Sr%sigxme(jb,sk_ibz,iab)+Sr%sigcmee0(jb,sk_ibz,iab)+Sr%degw(jb,sk_ibz,1)*Sr%dsigmee0(jb,sk_ibz,iab)
311          end do
312        end if
313 
314      end do ! jb
315    end do ! spin
316 
317  else
318    ! =============================
319    ! === Analytic Continuation ===
320    ! =============================
321    ABI_CHECK(Sr%nsig_ab == 1, "AC with spinor not implemented")
322 
323    ! Index of the KS or QP energy in sigme_tmp
324    !ie0 = sr%nomega_r + Sr%nomega4sd/2+1
325 
326    do spin=1,Sr%nsppol
327      do jb=ib1,ib2
328 
329       ABI_MALLOC(tmpcdp,(Sr%nomega_i))
330       ! Calculate Sigc(E0), dSigc(E0)
331       zz = CMPLX(Sr%e0(jb,sk_ibz,spin), zero)
332 
333       if (Sigp%mbpt_sciss > 0.1d-4) then
334         ! e0 is replaced by qp_ene which contains the updated energy eigenvalue
335         zz = CMPLX(qp_ene(jb,sk_ibz,spin), zero)
336       end if
337 
338       !call spade%init(sr%nomega_i, sr%omega_i, tmpcdp, branch_cut=">")
339       !call spade%eval(zz, sigc_e0, dzdval=dsigc_de0)
340 
341       ! Diagonal elements of sigcme
342       ! if zz in 2 or 3 quadrant, avoid branch cut in the complex plane using Sigma(-iw) = Sigma(iw)*.
343       do iab=1,Sr%nsig_ab
344         is_idx=spin; if (Sr%nsig_ab>1) is_idx=iab
345         if (real(zz) > zero) then
346           tmpcdp(:)=sigcme(:,jb,jb,is_idx)
347           Sr%sigcmee0(jb,sk_ibz,is_idx) =  pade(Sr%nomega_i, Sr%omega_i, tmpcdp, zz)
348           Sr%dsigmee0(jb,sk_ibz,is_idx) = dpade(Sr%nomega_i, Sr%omega_i, tmpcdp, zz)
349         else
350           tmpcdp(:) = CONJG(sigcme(:,jb,jb,is_idx))
351           Sr%sigcmee0(jb,sk_ibz,is_idx) =  pade(Sr%nomega_i, CONJG(Sr%omega_i), tmpcdp, zz)
352           Sr%dsigmee0(jb,sk_ibz,is_idx) = dpade(Sr%nomega_i, CONJG(Sr%omega_i), tmpcdp, zz)
353         end if
354       end do !iab
355 
356       ! Z = (1 - dSigma / domega(E0))^{-1}
357       if (Sr%nsig_ab == 1) then
358         Sr%ze0(jb,sk_ibz,spin) = one / (one - Sr%dsigmee0(jb,sk_ibz,spin))
359       else
360         Sr%ze0(jb,sk_ibz,1) = one / (one - SUM(Sr%dsigmee0(jb,sk_ibz,:)))
361       end if
362 
363       ! MG FIXME: Here we are solving the non-linear QP equation using the Pade' continuation + root finding
364       ! but this is very misleading because in the output file we are still reporting the Z factor
365       ! and there's no mention that the QP energies have been obtained from the non-linear equation!!
366       ! One should change the format used to print the results or at least warn the user!
367 
368       ! Find roots of E^0-V_xc-V_U+Sig_x+Sig_c(z)-z, i.e E^qp.
369       ! using Newton-Raphson method and starting point E^0
370       zz = CMPLX(Sr%e0(jb,sk_ibz,spin), zero)
371 
372       if (Sigp%mbpt_sciss>0.1d-4) then
373         ! e0 is replaced by qp_ene which contains the updated energy eigenvalue.
374         zz = CMPLX(qp_ene(jb,sk_ibz,spin),0.0)
375       end if
376 
377       ! Solve the QP equation with Newton-Rapson starting from e0
378       !call spade%qp_solve(e0, v_meanf, sigx, zz, zsc, msg, ierr)
379       !qpe_pade_kcalc(ibc, ikcalc, spin) = zsc
380       !qp_solver_ierr(ibc, ikcalc, spin) = ierr
381       !if (ierr /= 0) then
382       !  ABI_WARNING(msg)
383       !end if
384 
385       iter = 0; converged = .FALSE.; ctdpc = cone
386       do while (ABS(ctdpc) > NR_ABS_ROOT_ERR .or. iter < NR_MAX_NITER)
387         iter = iter + 1
388         sigc = czero; dsigc = czero
389         if (REAL(zz) > tol12) then
390           tmpcdp(:) = sigcme(:,jb,jb,spin)
391           sigc  =  pade(Sr%nomega_i, Sr%omega_i, tmpcdp, zz)
392           dsigc = dpade(Sr%nomega_i, Sr%omega_i, tmpcdp, zz)
393         else
394           tmpcdp(:) = CONJG(sigcme(:,jb,jb,spin))
395           sigc  =  pade(Sr%nomega_i, CONJG(Sr%omega_i), tmpcdp, zz)
396           dsigc = dpade(Sr%nomega_i, CONJG(Sr%omega_i), tmpcdp, zz)
397         end if
398         ctdpc = Sr%e0(jb,sk_ibz,spin) - Sr%vxcme(jb,sk_ibz,spin) - Sr%vUme(jb,sk_ibz,spin) + Sr%sigxme(jb,sk_ibz,spin) &
399                 + sigc - zz
400         if (ABS(ctdpc) < NR_ABS_ROOT_ERR) then
401           converged=.TRUE.; EXIT
402         end if
403         dct = dsigc - one
404         zz = newrap_step(zz, ctdpc, dct)
405       end do
406 
407       if (.not. converged) then
408         write(msg,'(a,i0,3a,f8.4,a,f8.4)')&
409           'Newton-Raphson method not converged after ',NR_MAX_NITER,' iterations. ',ch10,&
410           'Absolute Error = ',ABS(ctdpc),' > ',NR_ABS_ROOT_ERR
411         ABI_WARNING(msg)
412       end if
413 
414       ! Store the final result TODO re-shift everything according to efermi
415       Sr%egw(jb,sk_ibz,spin) = zz
416       Sr%degw(jb,sk_ibz,spin) = Sr%egw(jb,sk_ibz,spin) - Sr%e0(jb,sk_ibz,spin)
417       Sr%sigmee(jb,sk_ibz,spin) = Sr%sigxme(jb,sk_ibz,spin) + sigc
418 
419       ! Spectra of Sigma, remember that Sr%nomega_r does not contains the frequencies
420       ! used to evaluate the derivative each frequency is obtained using the pade_expression
421       ! In sigma indeed we have:
422       !     nomega_sigc=Sr%nomega_r+Sr%nomega4sd; if (mod10==SIG_GW_AC) nomega_sigc=Sr%nomega_i
423       do io=1,Sr%nomega_r
424         zz=Sr%omega_r(io)
425         if (REAL(zz) > zero) then
426           tmpcdp(:) = sigcme(:,jb,jb,spin)
427           Sr%sigcme(jb,sk_ibz,io,spin) = pade(Sr%nomega_i, Sr%omega_i, tmpcdp, zz)
428         else
429           tmpcdp(:) = CONJG(sigcme(:,jb,jb,spin))
430           Sr%sigcme(jb,sk_ibz,io,spin) = pade(Sr%nomega_i, CONJG(Sr%omega_i), tmpcdp, zz)
431         end if
432         Sr%sigxcme(jb,sk_ibz,io,spin) = Sr%sigxme(jb,sk_ibz,spin) + Sr%sigcme(jb,sk_ibz,io,spin)
433       end do
434 
435       ! Save sigma values along the imaginary axis
436       do iab=1,Sr%nsig_ab
437         is_idx=spin; if (Sr%nsig_ab > 1) is_idx = iab
438         do io=1,Sr%nomega_i
439           Sr%sigcmesi (jb,sk_ibz,io,is_idx) = sigcme(io,jb,jb,is_idx)
440           Sr%sigxcmesi(jb,sk_ibz,io,is_idx) = Sr%sigxme(jb,sk_ibz,is_idx) + Sr%sigcmesi(jb,sk_ibz,io,is_idx)
441         end do
442       end do
443 
444       ABI_FREE(tmpcdp)
445 
446      end do !jb
447    end do !is
448  end if ! Analytic continuation.
449 
450  ! === Diagonalize the QP Hamiltonian (forced to be Hermitian) ===
451  ! Calculate Sr%en_qp_diago and Sr%eigvec_qp to be written in the QPS file.
452  ! TODO in case of AC results are wrong.
453 
454 if (mod10 /= 1) then
455 
456  ABI_MALLOC(hhartree, (ib1:ib2,ib1:ib2,Sr%nsppol*Sr%nsig_ab))
457  hhartree = Sr%hhartree(ib1:ib2,ib1:ib2,sk_ibz,:)
458 
459  ! If non self-consistent erase all off-diagonal elements
460  if (Sigp%gwcalctyp<20) then
461    do jb=ib1,ib2
462      do kb=ib1,ib2
463       if (jb == kb) CYCLE
464       hhartree(jb,kb,:) = czero
465      end do
466    end do
467  end if
468 
469  ABI_MALLOC(htotal, (ib1:ib2,ib1:ib2,Sr%nsppol*Sr%nsig_ab))
470  do spin=1,Sr%nsppol*Sr%nsig_ab
471    do jb=ib1,ib2
472      do kb=ib1,ib2
473       htotal(kb,jb,spin) = hhartree(kb,jb,spin) + Sr%x_mat(kb,jb,sk_ibz,spin) + sigcme(ie0,kb,jb,spin)
474      end do
475    end do
476  end do
477 
478  ! Get the Hermitian part of htotal
479  ! In the noncollinear case A_{12}^{ab} = A_{21}^{ba}^* if A is Hermitian.
480  ABI_MALLOC(h_tmp1,(ib1:ib2,ib1:ib2))
481  ABI_MALLOC(h_tmp2,(ib1:ib2,ib1:ib2))
482 
483  nsploop=Sr%nsppol; if (Sr%nsig_ab/=1) nsploop=2
484  do spin=1,nsploop
485    h_tmp1 = CONJG(htotal(:,:,spin))
486    h_tmp2 = TRANSPOSE(h_tmp1)
487    h_tmp1 = htotal(:,:,spin)
488    htotal(:,:,spin)= half * (h_tmp1 + h_tmp2)
489  end do
490 
491  ! Print the different matrix elements of sigma if QPSC and prtvol>9
492  if (Sigp%gwcalctyp >=20 .and. mod10 /= 1 .and. prtvol>9 .and. my_rank==master) then
493    call print_sigma_melems(ikcalc,ib1,ib2,Sr%nsppol*Sr%nsig_ab,htotal,hhartree,&
494                            Sr%x_mat(ib1:ib2,ib1:ib2,sk_ibz,:),sigcme(ie0,:,:,:),Dtfil%filnam_ds(4))
495  end if
496 
497  if (Sr%nsig_ab==4) then
498    h_tmp1 = CONJG(htotal(:,:,4))
499    h_tmp2 = TRANSPOSE(h_tmp1)
500    h_tmp1 = htotal(:,:,3)
501    htotal(:,:,3)= half * (h_tmp1 + h_tmp2)
502 
503    h_tmp1 = CONJG(htotal(:,:,3))
504    h_tmp2 = TRANSPOSE(h_tmp1)
505    htotal(:,:,4) = h_tmp2
506  end if
507 
508  ! Solve Herm(htotal)*U = E*U
509  ld_matrix = ib2 - ib1 + 1
510  ABI_MALLOC(hdp, (ld_matrix,ld_matrix))
511  ABI_MALLOC(eig, (ld_matrix))
512 
513  do spin=1,Sr%nsppol
514    if (Sr%nsig_ab==1) then
515      hdp=htotal(ib1:ib2,ib1:ib2,spin)
516    else
517      hdp = SUM(htotal(ib1:ib2,ib1:ib2,:), DIM=3)
518    end if
519    call xheev("Vectors","Upper", ld_matrix, hdp, eig)
520 
521    Sr%eigvec_qp(ib1:ib2,ib1:ib2,sk_ibz,spin)=hdp(:,:)
522    Sr%en_qp_diago(ib1:ib2,sk_ibz,spin)=eig(:)
523  end do
524 
525  ABI_FREE(hdp)
526  ABI_FREE(eig)
527  ABI_FREE(htotal)
528  ABI_FREE(hhartree)
529  ABI_FREE(h_tmp1)
530  ABI_FREE(h_tmp2)
531 
532 end if
533 
534  call timab(490,2,tsec)
535 
536  DBG_EXIT("COLL")
537 
538 end subroutine solve_dyson