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-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 MODULE m_dyson_solver
28 
29  use defs_basis
30  use defs_abitypes
31  use m_xmpi
32  use m_errors
33  use m_abicore
34 
35  use m_time,          only : timab
36  use m_gwdefs,        only : sigparams_t
37  use m_numeric_tools, only : linfit, pade, dpade, newrap_step
38  use m_io_tools,      only : open_file
39  use m_fstrings,      only : int2char10
40  use m_hide_lapack,   only : xheev
41  use m_bz_mesh,       only : kmesh_t, get_BZ_item
42  use m_sigma,         only : sigma_t
43 
44  implicit none
45 
46  private
47 
48  public :: solve_dyson     ! Solve the Dyson equation for the QP energies.
49 
50  integer,private,parameter :: NR_MAX_NITER=1000
51   ! Max no of iterations in the Newton-Raphson method.
52 
53  real(dp),private,parameter :: NR_ABS_ROOT_ERR=0.0001/Ha_eV
54   ! Tolerance on the absolute error on the Newton-Raphson root.
55 
56 CONTAINS  !====================================================================

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

PARENTS

      m_dyson_solver

CHILDREN

      int2char10,wrtout

SOURCE

531 subroutine print_sigma_melems(ikcalc,ib1,ib2,nsp,htotal,hhartree,sigxme,sigcme,prefil)
532 
533 
534 !This section has been created automatically by the script Abilint (TD).
535 !Do not modify the following lines by hand.
536 #undef ABI_FUNC
537 #define ABI_FUNC 'print_sigma_melems'
538 !End of the abilint section
539 
540  implicit none
541 
542 ! Arguments ------------------------------------
543  !scalars
544  integer,intent(in) :: ikcalc,ib1,ib2,nsp
545  character(len=*),intent(in) :: prefil
546  !arrays
547  complex(dpc),intent(in) :: htotal(ib1:ib2,ib1:ib2,nsp),hhartree(ib1:ib2,ib1:ib2,nsp)
548  complex(dpc),intent(in) :: sigxme(ib1:ib2,ib1:ib2,nsp),sigcme(ib1:ib2,ib1:ib2,nsp)
549 
550 ! Local variables ------------------------------
551  integer,parameter :: MAX_NCOLS=14
552  integer :: isp,mc,mr,jj,ii,temp_unit,ount
553  character(len=10) :: sidx
554  character(len=500) :: msg
555  character(len=100) :: fmth,fmt1,fmt2,fmthh,kpt_index,fmtfile
556  character(len=fnlen) :: filename
557 ! *************************************************************************
558 
559  if (nsp==3.or.nsp>4) then
560    MSG_ERROR('nsp has wrong value in print_sigma_melems')
561  end if
562 
563  ount = std_out
564 
565  mc = ib2-ib1+1; if (mc>MAX_NCOLS) mc = MAX_NCOLS
566  mr = mc
567 
568  write(fmthh,*)'(2(a),2(I2,a))'
569  write(fmth,*)'(7x,',mc,'(i2,8x))'
570  write(fmt1,*)'(3x,i2,',mc,'f10.5)'
571  write(fmt2,*)'(5x   ,',mc,'f10.5,a)'
572 
573 ! First print to screen
574  do isp=1,nsp
575    write(msg,'(a)') ''
576    call wrtout(ount,msg,'COLL')
577    write(msg,fmthh) ch10,' Hermitianised matrix elements of Sigma (spin ',isp,' of ',nsp,'):'
578    call wrtout(ount,msg,'COLL')
579    write(msg,fmth)(jj,jj=1,mc)
580    call wrtout(ount,msg,'COLL') !header
581    do ii=ib1,ib1+mr-1
582      write(msg,fmt1)ii-ib1+1,DBLE(htotal(ii,ib1:(ib1+mc-1),isp))
583      call wrtout(ount,msg,'COLL') !real part
584      write(msg,fmt2)  AIMAG(htotal(ii,ib1:(ib1+mc-1),isp)),ch10
585      call wrtout(ount,msg,'COLL') !imag part
586    end do
587  end do !nsp
588 
589  write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output in _HTOTAL files"
590  call wrtout(ount,msg,'COLL')
591 
592  do isp=1,nsp
593    write(msg,fmthh) ch10,' H_Hartree matrix elements (spin ',isp,' of ',nsp,'):'
594    call wrtout(ount,msg,'COLL')
595    write(msg,fmth)(jj,jj=1,mc)
596    call wrtout(ount,msg,'COLL') !header
597    do ii=ib1,ib1+mr-1
598      write(msg,fmt1)ii-ib1+1,DBLE(hhartree(ii,ib1:(ib1+mc-1),isp))
599      call wrtout(ount,msg,'COLL') !real part
600      write(msg,fmt2)  AIMAG(hhartree(ii,ib1:(ib1+mc-1),isp)),ch10
601      call wrtout(ount,msg,'COLL') !imag part
602    end do
603  end do !nsp
604 
605  write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output in _HHARTREE files"
606  call wrtout(ount,msg,'COLL')
607 
608  do isp=1,nsp
609    write(msg,fmthh) ch10,' Sigma_x matrix elements (spin ',isp,' of ',nsp,'):'
610    call wrtout(ount,msg,'COLL')
611    write(msg,fmth)(jj,jj=1,mc)
612    call wrtout(ount,msg,'COLL') !header
613    do ii=ib1,ib1+mr-1
614      write(msg,fmt1)ii-ib1+1,DBLE(sigxme(ii,ib1:(ib1+mc-1),isp))
615      call wrtout(ount,msg,'COLL') !real part
616      write(msg,fmt2)  AIMAG(sigxme(ii,ib1:(ib1+mc-1),isp)),ch10
617      call wrtout(ount,msg,'COLL') !imag part
618    end do
619  end do !nsp
620 
621  write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output _SIGX files"
622  call wrtout(ount,msg,'COLL')
623 
624  do isp=1,nsp
625    write(msg,fmthh) ch10,' Sigma_c matrix elements (spin ',isp,' of ',nsp,'):'
626    call wrtout(ount,msg,'COLL')
627    write(msg,fmth)(jj,jj=1,mc)
628    call wrtout(ount,msg,'COLL') !header
629    do ii=ib1,ib1+mr-1
630      write(msg,fmt1)ii-ib1+1,DBLE(sigcme(ii,ib1:(ib1+mc-1),isp))
631      call wrtout(ount,msg,'COLL') !real part
632      write(msg,fmt2)  AIMAG(sigcme(ii,ib1:(ib1+mc-1),isp)),ch10
633      call wrtout(ount,msg,'COLL') !imag part
634    end do
635  end do !nsp
636 
637  write(msg,'(a,i2,a)')" Max ",MAX_NCOLS," elements printed. Full matrix output _SIGC files"
638  call wrtout(ount,msg,'COLL')
639 
640  ! Then print to file
641  ! Format is: row, column, value; with a blank space for each full
642  ! set of columns for easy plotting with the gnuplot splot command
643  write(fmtfile,*)'(3X,I6,2X,I6,',nsp,'(2(ES28.16E3,3x)))'
644 
645  call int2char10(ikcalc,sidx)
646  kpt_index = "_KPT"//TRIM(sidx)
647 
648  filename = TRIM(prefil)//'_HTOTAL'//TRIM(kpt_index)
649 
650  if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
651    MSG_ERROR(msg)
652  end if
653 
654  msg = '#   row    col.      Re(htotal(r,c)) Im(htotal(r,c))  for spin11   ... spin22 ... spin12 ... spin13'
655  call wrtout(temp_unit,msg,'COLL')
656  do ii=ib1,ib2
657    do jj=ib1,ib2
658      write(msg,fmtfile) ii,jj,(htotal(jj,ii,isp),isp=1,nsp)
659      call wrtout(temp_unit,msg,'COLL')
660    end do
661    call wrtout(temp_unit,"",'COLL')
662  end do
663  close(temp_unit)
664 
665  filename = TRIM(prefil)//'_HHARTREE'//TRIM(kpt_index)
666  if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
667    MSG_ERROR(msg)
668  end if
669 
670  msg = '#   row    col.      Re(hhartree(r,c))  Im(hhartree(r,c)  for spin11   ... spin22 ... spin12 ... spin13'
671  call wrtout(temp_unit,msg,'COLL')
672  do ii=ib1,ib2
673    do jj=ib1,ib2
674      write(msg,fmtfile) ii,jj,(hhartree(jj,ii,isp),isp=1,nsp)
675      call wrtout(temp_unit,msg,'COLL')
676    end do
677    call wrtout(temp_unit,"",'COLL')
678  end do
679  close(temp_unit)
680 
681  filename = TRIM(prefil)//'_SIGX'//TRIM(kpt_index)
682  if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
683    MSG_ERROR(msg)
684  end if
685 
686  write(msg,'(a)')'#   row    col.      Re(Sigx(r,c)) Im(Sigx(r,c) for spin11   ... spin22 ... spin12 ... spin13'
687  call wrtout(temp_unit,msg,'COLL')
688  do ii=ib1,ib2
689    do jj=ib1,ib2
690      write(msg,fmtfile) ii,jj,(sigxme(jj,ii,isp),isp=1,nsp)
691      call wrtout(temp_unit,msg,'COLL')
692    end do
693    call wrtout(temp_unit,"",'COLL')
694  end do
695  close(temp_unit)
696 
697  filename = TRIM(prefil)//'_SIGC'//TRIM(kpt_index)
698  if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
699    MSG_ERROR(msg)
700  end if
701 
702  write(msg,'(a)')'#   row    col.      Re(Sigc(r,c)) Im(Sigc(r,c) for spin11   ... spin22 ... spin12 ... spin21'
703  call wrtout(temp_unit,msg,'COLL')
704  do ii=ib1,ib2
705    do jj=ib1,ib2
706      write(msg,fmtfile) ii,jj,(sigcme(jj,ii,isp),isp=1,nsp)
707      call wrtout(temp_unit,msg,'COLL')
708    end do
709    call wrtout(temp_unit,"",'COLL')
710  end do
711 
712  close(temp_unit)
713 
714 end subroutine print_sigma_melems
715 
716 !----------------------------------------------------------------------
717 
718 END MODULE m_dyson_solver

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 analitic 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_tmp=(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.

PARENTS

      sigma

CHILDREN

      int2char10,wrtout

SOURCE

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