TABLE OF CONTENTS
- ABINIT/m_dyson_solver
- m_dyson_solver/print_sigma_melems
- m_dyson_solver/sigma_pade_eval
- m_dyson_solver/sigma_pade_init
- m_dyson_solver/sigma_pade_qp_solve
- m_dyson_solver/sigma_pade_t
- m_dyson_solver/solve_dyson
ABINIT/m_dyson_solver [ 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