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-2024 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
574 subroutine print_sigma_melems(ikcalc,ib1,ib2,nsp,htotal,hhartree,sigxme,sigcme,prefil) 575 576 ! Arguments ------------------------------------ 577 !scalars 578 integer,intent(in) :: ikcalc,ib1,ib2,nsp 579 character(len=*),intent(in) :: prefil 580 !arrays 581 complex(dpc),intent(in) :: htotal(ib1:ib2,ib1:ib2,nsp),hhartree(ib1:ib2,ib1:ib2,nsp) 582 complex(dpc),intent(in) :: sigxme(ib1:ib2,ib1:ib2,nsp),sigcme(ib1:ib2,ib1:ib2,nsp) 583 584 ! Local variables ------------------------------ 585 integer,parameter :: MAX_NCOLS = 14 586 integer :: isp,mc,mr,jj,ii,temp_unit,ount 587 character(len=10) :: sidx 588 character(len=500) :: msg 589 character(len=100) :: fmth,fmt1,fmt2,fmthh,kpt_index,fmtfile 590 character(len=fnlen) :: filename 591 592 ! ************************************************************************* 593 594 if (nsp==3.or.nsp>4) then 595 ABI_ERROR('nsp has wrong value in print_sigma_melems') 596 end if 597 598 ount = std_out 599 600 mc = ib2-ib1+1; if (mc>MAX_NCOLS) mc = MAX_NCOLS 601 mr = mc 602 603 write(fmthh,*)'(2(a),2(I2,a))' 604 write(fmth,*)'(7x,',mc,'(i2,8x))' 605 write(fmt1,*)'(3x,i2,',mc,'f10.5)' 606 write(fmt2,*)'(5x ,',mc,'f10.5,a)' 607 608 ! First print to screen 609 do isp=1,nsp 610 write(msg,'(a)') '' 611 call wrtout(ount,msg) 612 write(msg,fmthh) ch10,' Hermitianised matrix elements of Sigma (spin ',isp,' of ',nsp,'):' 613 call wrtout(ount,msg) 614 write(msg,fmth)(jj,jj=1,mc) 615 call wrtout(ount,msg) !header 616 do ii=ib1,ib1+mr-1 617 write(msg,fmt1)ii-ib1+1,DBLE(htotal(ii,ib1:(ib1+mc-1),isp)) 618 call wrtout(ount,msg) !real part 619 write(msg,fmt2) AIMAG(htotal(ii,ib1:(ib1+mc-1),isp)),ch10 620 call wrtout(ount,msg) !imag part 621 end do 622 end do !nsp 623 624 write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output in _HTOTAL files" 625 call wrtout(ount,msg) 626 627 do isp=1,nsp 628 write(msg,fmthh) ch10,' H_Hartree matrix elements (spin ',isp,' of ',nsp,'):' 629 call wrtout(ount,msg) 630 write(msg,fmth)(jj,jj=1,mc) 631 call wrtout(ount,msg) !header 632 do ii=ib1,ib1+mr-1 633 write(msg,fmt1)ii-ib1+1,DBLE(hhartree(ii,ib1:(ib1+mc-1),isp)) 634 call wrtout(ount,msg) !real part 635 write(msg,fmt2) AIMAG(hhartree(ii,ib1:(ib1+mc-1),isp)),ch10 636 call wrtout(ount,msg) !imag part 637 end do 638 end do !nsp 639 640 write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output in _HHARTREE files" 641 call wrtout(ount,msg) 642 643 do isp=1,nsp 644 write(msg,fmthh) ch10,' Sigma_x matrix elements (spin ',isp,' of ',nsp,'):' 645 call wrtout(ount,msg) 646 write(msg,fmth)(jj,jj=1,mc) 647 call wrtout(ount,msg) !header 648 do ii=ib1,ib1+mr-1 649 write(msg,fmt1)ii-ib1+1,DBLE(sigxme(ii,ib1:(ib1+mc-1),isp)) 650 call wrtout(ount,msg) !real part 651 write(msg,fmt2) AIMAG(sigxme(ii,ib1:(ib1+mc-1),isp)),ch10 652 call wrtout(ount,msg) !imag part 653 end do 654 end do !nsp 655 656 write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output _SIGX files" 657 call wrtout(ount,msg) 658 659 do isp=1,nsp 660 write(msg,fmthh) ch10,' Sigma_c matrix elements (spin ',isp,' of ',nsp,'):' 661 call wrtout(ount,msg) 662 write(msg,fmth)(jj,jj=1,mc) 663 call wrtout(ount,msg) !header 664 do ii=ib1,ib1+mr-1 665 write(msg,fmt1)ii-ib1+1,DBLE(sigcme(ii,ib1:(ib1+mc-1),isp)) 666 call wrtout(ount,msg) !real part 667 write(msg,fmt2) AIMAG(sigcme(ii,ib1:(ib1+mc-1),isp)),ch10 668 call wrtout(ount,msg) !imag part 669 end do 670 end do !nsp 671 672 write(msg,'(a,i2,a)')" Max ",MAX_NCOLS," elements printed. Full matrix output _SIGC files" 673 call wrtout(ount,msg) 674 675 ! Then print to file 676 ! Format is: row, column, value; with a blank space for each full 677 ! set of columns for easy plotting with the gnuplot splot command 678 write(fmtfile,*)'(3X,I6,2X,I6,',nsp,'(2(ES28.16E3,3x)))' 679 680 call int2char10(ikcalc,sidx) 681 kpt_index = "_KPT"//TRIM(sidx) 682 683 filename = TRIM(prefil)//'_HTOTAL'//TRIM(kpt_index) 684 685 if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then 686 ABI_ERROR(msg) 687 end if 688 689 msg = '# row col. Re(htotal(r,c)) Im(htotal(r,c)) for spin11 ... spin22 ... spin12 ... spin13' 690 call wrtout(temp_unit,msg) 691 do ii=ib1,ib2 692 do jj=ib1,ib2 693 write(msg,fmtfile) ii,jj,(htotal(jj,ii,isp),isp=1,nsp) 694 call wrtout(temp_unit,msg) 695 end do 696 call wrtout(temp_unit,"") 697 end do 698 close(temp_unit) 699 700 filename = TRIM(prefil)//'_HHARTREE'//TRIM(kpt_index) 701 if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then 702 ABI_ERROR(msg) 703 end if 704 705 msg = '# row col. Re(hhartree(r,c)) Im(hhartree(r,c) for spin11 ... spin22 ... spin12 ... spin13' 706 call wrtout(temp_unit,msg) 707 do ii=ib1,ib2 708 do jj=ib1,ib2 709 write(msg,fmtfile) ii,jj,(hhartree(jj,ii,isp),isp=1,nsp) 710 call wrtout(temp_unit,msg) 711 end do 712 call wrtout(temp_unit,"") 713 end do 714 close(temp_unit) 715 716 filename = TRIM(prefil)//'_SIGX'//TRIM(kpt_index) 717 if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then 718 ABI_ERROR(msg) 719 end if 720 721 write(msg,'(a)')'# row col. Re(Sigx(r,c)) Im(Sigx(r,c) for spin11 ... spin22 ... spin12 ... spin13' 722 call wrtout(temp_unit,msg) 723 do ii=ib1,ib2 724 do jj=ib1,ib2 725 write(msg,fmtfile) ii,jj,(sigxme(jj,ii,isp),isp=1,nsp) 726 call wrtout(temp_unit,msg) 727 end do 728 call wrtout(temp_unit,"") 729 end do 730 close(temp_unit) 731 732 filename = TRIM(prefil)//'_SIGC'//TRIM(kpt_index) 733 if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then 734 ABI_ERROR(msg) 735 end if 736 737 write(msg,'(a)')'# row col. Re(Sigc(r,c)) Im(Sigc(r,c) for spin11 ... spin22 ... spin12 ... spin21' 738 call wrtout(temp_unit,msg) 739 do ii=ib1,ib2 740 do jj=ib1,ib2 741 write(msg,fmtfile) ii,jj,(sigcme(jj,ii,isp),isp=1,nsp) 742 call wrtout(temp_unit,msg) 743 end do 744 call wrtout(temp_unit,"") 745 end do 746 747 close(temp_unit) 748 749 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
806 subroutine sigma_pade_eval(self, zz, val, dzdval) 807 808 !Arguments ------------------------------------ 809 class(sigma_pade_t),intent(in) :: self 810 complex(dp),intent(in) :: zz 811 complex(dp),intent(out) :: val 812 complex(dp),optional,intent(out) :: dzdval 813 814 ! ************************************************************************* 815 816 ! if zz in 2 or 3 quadrant, avoid branch cut in the complex plane using Sigma(-iw) = Sigma(iw)*. 817 if (real(zz) > zero) then 818 !if (real(zz) < zero) then 819 val = pade(self%npts, self%zmesh, self%sigc_cvals, zz) 820 if (present(dzdval)) dzdval = dpade(self%npts, self%zmesh, self%sigc_cvals, zz) 821 else 822 val = pade(self%npts, -self%zmesh, conjg(self%sigc_cvals), zz) 823 if (present(dzdval)) dzdval = dpade(self%npts, -self%zmesh, conjg(self%sigc_cvals), zz) 824 end if 825 826 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
763 subroutine sigma_pade_init(self, npts, zmesh, sigc_cvals, branch_cut) 764 765 !Arguments ------------------------------------ 766 class(sigma_pade_t),intent(out) :: self 767 integer,intent(in) :: npts 768 complex(dp),target,intent(in) :: zmesh(npts), sigc_cvals(npts) 769 character(len=*),intent(in) :: branch_cut 770 771 !Local variables------------------------------- 772 !scalars 773 ! integer :: ii 774 775 ! ************************************************************************* 776 777 self%npts = npts 778 779 ! Select first points according to wmax. 780 !if (wmax > zero) then 781 ! do ii=1,npts 782 ! if (real(zmesh(ii)) > wmax) exit 783 ! end do 784 ! self%npts = ii-1 785 !end if 786 787 self%branch_cut = branch_cut 788 self%zmesh => zmesh(1:self%npts) 789 self%sigc_cvals => sigc_cvals(1:self%npts) 790 791 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
843 subroutine sigma_pade_qp_solve(self, e0, v_meanf, sigx, z_guess, zsc, msg, ierr) 844 845 !Arguments ------------------------------------ 846 class(sigma_pade_t),intent(in) :: self 847 real(dp),intent(in) :: e0, v_meanf, sigx 848 complex(dp),intent(in) :: z_guess 849 complex(dp),intent(out) :: zsc 850 integer,intent(out) :: ierr 851 852 !Local variables------------------------------- 853 !scalars 854 integer :: iter 855 logical :: converged 856 complex(dpc) :: ctdpc, dct, dsigc, sigc 857 character(len=500) :: msg 858 859 ! ************************************************************************* 860 861 ! Use Newton-Rapson to find the root of: 862 ! f(z) = e0 - zz + Sigma_xc(z) - v_meanf 863 ! f'(z) = -1 + Sigma_c'(z) 864 865 iter = 0; converged = .FALSE.; ctdpc = cone 866 zsc = z_guess 867 do while (ABS(ctdpc) > NR_ABS_ROOT_ERR .or. iter < NR_MAX_NITER) 868 iter = iter + 1 869 870 call self%eval(zsc, sigc, dzdval=dsigc) 871 ctdpc = e0 - v_meanf + sigx + sigc - zsc 872 873 if (ABS(ctdpc) < NR_ABS_ROOT_ERR) then 874 converged=.TRUE.; EXIT 875 end if 876 dct = dsigc - one 877 zsc = newrap_step(zsc, ctdpc, dct) 878 end do 879 880 ierr = 0; msg = "" 881 if (.not. converged) then 882 write(msg,'(a,i0,3a,f8.4,a,f8.4)')& 883 'Newton-Raphson method not converged after ',NR_MAX_NITER,' iterations. ',ch10,& 884 'Absolute Error = ',ABS(ctdpc),' > ',NR_ABS_ROOT_ERR 885 ierr = 1 886 end if 887 888 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 ! Number of points 62 63 character(len=1) :: branch_cut 64 65 complex(dp),pointer :: zmesh(:) => null(), sigc_cvals(:) => null() 66 ! pointer to input mesh and values. 67 68 !real(d) :: wmax = -one 69 70 contains 71 72 procedure :: init => sigma_pade_init 73 ! Init object 74 75 procedure :: eval => sigma_pade_eval 76 ! Eval self-energy and derivative 77 78 procedure :: qp_solve => sigma_pade_qp_solve 79 ! Find the QP solution with Newton-Rapson method 80 81 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 indices included in the GW calculation for this k-point.
SOURCE
143 subroutine solve_dyson(ikcalc,minbnd,maxbnd,nomega_sigc,Sigp,Kmesh,sigcme,qp_ene,Sr,prtvol,Dtfil,comm) 144 145 !Arguments ------------------------------------ 146 !scalars 147 integer,intent(in) :: ikcalc,nomega_sigc,prtvol,minbnd,maxbnd,comm 148 type(kmesh_t),intent(in) :: Kmesh 149 type(Datafiles_type),intent(in) :: Dtfil 150 type(sigparams_t),intent(in) :: Sigp 151 type(sigma_t),intent(inout) :: Sr 152 !arrays 153 real(dp),intent(in) :: qp_ene(Sr%nbnds,Sr%nkibz,Sr%nsppol) 154 complex(dpc),intent(in) :: sigcme(nomega_sigc,minbnd:maxbnd,minbnd:maxbnd,Sigp%nsppol*Sigp%nsig_ab) 155 156 !Local variables------------------------------- 157 !scalars 158 integer,parameter :: master=0 159 integer :: iab,ib1,ib2,ikbz_gw,io,spin,is_idx,isym,iter,itim,jb, ie0 160 integer :: sk_ibz,kb,ld_matrix,mod10,nsploop,my_rank 161 real(dp) :: alpha,beta,smrt 162 complex(dpc) :: ctdpc,dct,dsigc,sigc,zz,phase 163 logical :: converged,ltest 164 character(len=500) :: msg 165 !type(sigma_pade_t) :: spade 166 !arrays 167 real(dp) :: kbz_gw(3),tsec(2) 168 real(dp),allocatable :: e0pde(:),eig(:),scme(:) 169 complex(dpc),allocatable :: hdp(:,:),tmpcdp(:),hhartree(:,:,:),htotal(:,:,:),h_tmp1(:,:),h_tmp2(:,:) 170 171 ! ************************************************************************* 172 173 DBG_ENTER("COLL") 174 175 call timab(490,1,tsec) ! csigme(Dyson) 176 177 my_rank = xmpi_comm_rank(comm) 178 179 mod10=MOD(Sigp%gwcalctyp,10) 180 181 ltest=(nomega_sigc==Sr%nomega_r+Sr%nomega4sd) 182 if (mod10==1) ltest=(nomega_sigc==Sr%nomega_i) 183 ABI_CHECK(ltest,'Wrong number of frequencies') 184 185 ! Index of the KS or QP energy. 186 !ioe0j=Sr%nomega4sd/2+1 187 188 ! min and Max band index for GW corrections (for this k-point). 189 ib1=MINVAL(Sigp%minbnd(ikcalc,:)) 190 ib2=MAXVAL(Sigp%maxbnd(ikcalc,:)) 191 192 ! Find the index of the k-point for sigma in the IBZ array. 193 ikbz_gw=Sigp%kptgw2bz(ikcalc) 194 call kmesh%get_BZ_item(ikbz_gw,kbz_gw,sk_ibz,isym,itim,phase) 195 196 sigc=czero; dsigc=czero 197 198 ! =========================================================== 199 ! ==== Solve the Dyson Equation and store results in Sr% ==== 200 ! =========================================================== 201 202 if (mod10 /= 1) then 203 ! =============================== 204 ! ==== Perturbative approach ==== 205 ! =============================== 206 207 ! Index of the KS or QP energy in sigme_tmp 208 ie0 = sr%nomega_r + Sr%nomega4sd/2+1 209 210 do spin=1,Sr%nsppol 211 do jb=ib1,ib2 212 ! === Get matrix elements of Sigma_c at energy E0 === 213 ! SigC(w) is linearly interpolated and the slope alpha is assumed as dSigC/dE 214 do iab=1,Sr%nsig_ab 215 is_idx=spin; if (Sr%nsig_ab>1) is_idx=iab 216 217 Sr%sigcmee0(jb,sk_ibz,is_idx) = sigcme(ie0,jb,jb,is_idx) 218 219 ABI_MALLOC(scme,(Sr%nomega4sd)) 220 ABI_MALLOC(e0pde,(Sr%nomega4sd)) 221 e0pde(:) = Sr%omega4sd(jb,sk_ibz,:,spin) 222 scme(:) = REAL(sigcme(Sr%nomega_r+1:Sr%nomega_r+Sr%nomega4sd,jb,jb,is_idx)) 223 224 if (Sr%nomega4sd==1) then 225 smrt = zero; alpha = zero 226 else 227 smrt = linfit(Sr%nomega4sd,e0pde(:),scme(:),alpha,beta) 228 end if 229 230 if (smrt>0.1/Ha_eV) then 231 write(msg,'(3a,i0,a,i0,2a,2(f22.15,2a))')& 232 'Values of Re Sig_c are not linear ',ch10,& 233 'band index = ',jb,' spin|component = ',is_idx,ch10,& 234 'root mean square= ',smrt,ch10,& 235 'estimated slope = ',alpha,ch10,& 236 'Omega [eV] SigC [eV]' 237 ABI_WARNING(msg) 238 do io=1,Sr%nomega4sd 239 write(msg,'(2f8.4)')e0pde(io)*Ha_eV,scme(io)*Ha_eV 240 call wrtout(std_out,msg) 241 end do 242 end if 243 244 ABI_FREE(scme) 245 ABI_FREE(e0pde) 246 ! 247 ! === Evaluate renormalization factor and QP correction === 248 ! * Z=(1-dSigma/domega(E0))^-1 249 ! * DeltaE_GW=E-E0= (Sigma(E0)-V_xc)/(1-dSigma/domega) 250 ! * If nspinor==2, this part is done at the end. 251 ! 252 Sr%dsigmee0(jb,sk_ibz,is_idx)=CMPLX(alpha,zero) 253 254 if (Sr%nsig_ab==1) then 255 Sr%ze0(jb,sk_ibz,spin)=one/(one-Sr%dsigmee0(jb,sk_ibz,spin)) 256 257 if (ABS(Sigp%mbpt_sciss) < tol6) then 258 Sr%degw(jb,sk_ibz,spin) = Sr%ze0(jb,sk_ibz,spin) * & 259 (Sr%sigxme(jb,sk_ibz,spin) + Sr%sigcmee0(jb,sk_ibz,spin) - Sr%e0(jb,sk_ibz,spin) + & 260 Sr%hhartree(jb,jb,sk_ibz,spin)) 261 262 Sr%egw(jb,sk_ibz,spin) = Sr%e0(jb,sk_ibz,spin) + Sr%degw(jb,sk_ibz,spin) 263 264 ! Estimate Sigma at the QP-energy: Sigma(E_qp)=Sigma(E0)+(E_qp-E0)*dSigma/dE 265 Sr%sigmee(jb,sk_ibz,spin) = & 266 Sr%sigxme(jb,sk_ibz,spin)+Sr%sigcmee0(jb,sk_ibz,spin)+Sr%degw(jb,sk_ibz,spin)*Sr%dsigmee0(jb,sk_ibz,spin) 267 268 else 269 ! If GW+scissor: e0 is replaced by qp_ene which contains the updated energy eigenvalue 270 Sr%degw(jb,sk_ibz,spin)= Sr%ze0(jb,sk_ibz,spin) * & 271 (Sr%sigxme(jb,sk_ibz,spin) + Sr%sigcmee0(jb,sk_ibz,spin) - qp_ene(jb,sk_ibz,spin) + & 272 Sr%hhartree(jb,jb,sk_ibz,spin)) 273 274 Sr%egw(jb,sk_ibz,spin) = qp_ene(jb,sk_ibz,spin) + Sr%degw(jb,sk_ibz,spin) 275 276 ! Estimate Sigma at the QP-energy: Sigma(E_qp)=Sigma(E0)+(E_qp-E0)*dSigma/dE 277 Sr%sigmee(jb,sk_ibz,spin)= & 278 Sr%sigxme(jb,sk_ibz,spin) + Sr%sigcmee0(jb,sk_ibz,spin) + & 279 Sr%degw(jb,sk_ibz,spin) * Sr%dsigmee0(jb,sk_ibz,spin) 280 281 ! RS: In the output, the gw corr with respect to e0 without mbpt_sciss is reported. 282 Sr%degw(jb,sk_ibz,spin) = Sr%egw(jb,sk_ibz,spin) - Sr%e0(jb,sk_ibz,spin) 283 end if 284 end if !Sigp%nsig_ab==1 285 286 ! Spectrum of Sigma 287 do io=1,Sr%nomega_r 288 Sr%sigcme (jb,sk_ibz,io,is_idx)= sigcme(io,jb,jb,is_idx) 289 Sr%sigxcme(jb,sk_ibz,io,is_idx)= Sr%sigxme(jb,sk_ibz,is_idx)+Sr%sigcme(jb,sk_ibz,io,is_idx) 290 end do 291 do io=1,Sr%nomega4sd 292 Sr%sigcme4sd (jb,sk_ibz,io,is_idx)= sigcme(Sr%nomega_r+io,jb,jb,is_idx) 293 Sr%sigxcme4sd(jb,sk_ibz,io,is_idx)= Sr%sigxme(jb,sk_ibz,is_idx)+Sr%sigcme4sd(jb,sk_ibz,io,is_idx) 294 end do 295 end do !iab 296 297 if (Sr%nsig_ab > 1) then 298 ABI_CHECK(ABS(Sigp%mbpt_sciss)<0.1d-4,'Scissor with spinor not coded') 299 !TODO this should be allocated with nsppol, recheck this part 300 301 ! Evaluate renormalization factor and QP correction. 302 ! Z=(1-dSigma/domega(E0))^-1 303 ! DeltaE_GW=E-E0= (Sigma(E0)-V_xc)/(1-dSigma/domega) 304 !write(std_out,'(a,i2,10f8.3)')' Correlation',jb,Sr%sigcmee0(jb,sk_ibz,:)*Ha_eV,SUM(Sr%sigcmee0(jb,sk_ibz,:))*Ha_eV 305 306 Sr%ze0 (jb,sk_ibz,1) = one/(one-SUM(Sr%dsigmee0(jb,sk_ibz,:))) 307 308 Sr%degw(jb,sk_ibz,1) = Sr%ze0(jb,sk_ibz,1) * & 309 (SUM(Sr%sigxme(jb,sk_ibz,:)+Sr%sigcmee0(jb,sk_ibz,:)+Sr%hhartree(jb,jb,sk_ibz,:))-Sr%e0(jb,sk_ibz,1)) 310 311 Sr%egw(jb,sk_ibz,1)=Sr%e0(jb,sk_ibz,1)+Sr%degw(jb,sk_ibz,1) 312 313 ! Estimate Sigma at the QP-energy. 314 do iab=1,Sr%nsig_ab 315 Sr%sigmee(jb,sk_ibz,iab)= & 316 Sr%sigxme(jb,sk_ibz,iab)+Sr%sigcmee0(jb,sk_ibz,iab)+Sr%degw(jb,sk_ibz,1)*Sr%dsigmee0(jb,sk_ibz,iab) 317 end do 318 end if 319 320 end do ! jb 321 end do ! spin 322 323 else 324 ! ============================= 325 ! === Analytic Continuation === 326 ! ============================= 327 ABI_CHECK(Sr%nsig_ab == 1, "AC with spinor not implemented") 328 329 ! Index of the KS or QP energy in sigme_tmp 330 !ie0 = sr%nomega_r + Sr%nomega4sd/2+1 331 332 do spin=1,Sr%nsppol 333 do jb=ib1,ib2 334 335 ABI_MALLOC(tmpcdp,(Sr%nomega_i)) 336 ! Calculate Sigc(E0), dSigc(E0) 337 zz = CMPLX(Sr%e0(jb,sk_ibz,spin), zero) 338 339 if (Sigp%mbpt_sciss > 0.1d-4) then 340 ! e0 is replaced by qp_ene which contains the updated energy eigenvalue 341 zz = CMPLX(qp_ene(jb,sk_ibz,spin), zero) 342 end if 343 344 !call spade%init(sr%nomega_i, sr%omega_i, tmpcdp, branch_cut=">") 345 !call spade%eval(zz, sigc_e0, dzdval=dsigc_de0) 346 347 ! Diagonal elements of sigcme 348 ! if zz in 2 or 3 quadrant, avoid branch cut in the complex plane using Sigma(-iw) = Sigma(iw)*. 349 do iab=1,Sr%nsig_ab 350 is_idx=spin; if (Sr%nsig_ab>1) is_idx=iab 351 if (real(zz) > zero) then 352 tmpcdp(:)=sigcme(:,jb,jb,is_idx) 353 Sr%sigcmee0(jb,sk_ibz,is_idx) = pade(Sr%nomega_i, Sr%omega_i, tmpcdp, zz) 354 Sr%dsigmee0(jb,sk_ibz,is_idx) = dpade(Sr%nomega_i, Sr%omega_i, tmpcdp, zz) 355 else 356 tmpcdp(:) = CONJG(sigcme(:,jb,jb,is_idx)) 357 Sr%sigcmee0(jb,sk_ibz,is_idx) = pade(Sr%nomega_i, CONJG(Sr%omega_i), tmpcdp, zz) 358 Sr%dsigmee0(jb,sk_ibz,is_idx) = dpade(Sr%nomega_i, CONJG(Sr%omega_i), tmpcdp, zz) 359 end if 360 end do !iab 361 362 ! Z = (1 - dSigma / domega(E0))^{-1} 363 if (Sr%nsig_ab == 1) then 364 Sr%ze0(jb,sk_ibz,spin) = one / (one - Sr%dsigmee0(jb,sk_ibz,spin)) 365 else 366 Sr%ze0(jb,sk_ibz,1) = one / (one - SUM(Sr%dsigmee0(jb,sk_ibz,:))) 367 end if 368 369 ! MG FIXME: Here we are solving the non-linear QP equation using the Pade' continuation + root finding 370 ! but this is very misleading because in the output file we are still reporting the Z factor 371 ! and there's no mention that the QP energies have been obtained from the non-linear equation!! 372 ! One should change the format used to print the results or at least warn the user! 373 374 ! Find roots of E^0-V_xc-V_U+Sig_x+Sig_c(z)-z, i.e E^qp. 375 ! using Newton-Raphson method and starting point E^0 376 zz = CMPLX(Sr%e0(jb,sk_ibz,spin), zero) 377 378 if (Sigp%mbpt_sciss>0.1d-4) then 379 ! e0 is replaced by qp_ene which contains the updated energy eigenvalue. 380 zz = CMPLX(qp_ene(jb,sk_ibz,spin),0.0) 381 end if 382 383 ! Solve the QP equation with Newton-Rapson starting from e0 384 !call spade%qp_solve(e0, v_meanf, sigx, zz, zsc, msg, ierr) 385 !qpe_pade_kcalc(ibc, ikcalc, spin) = zsc 386 !qp_solver_ierr(ibc, ikcalc, spin) = ierr 387 !if (ierr /= 0) then 388 ! ABI_WARNING(msg) 389 !end if 390 391 iter = 0; converged = .FALSE.; ctdpc = cone 392 do while (ABS(ctdpc) > NR_ABS_ROOT_ERR .or. iter < NR_MAX_NITER) 393 iter = iter + 1 394 sigc = czero; dsigc = czero 395 if (REAL(zz) > tol12) then 396 tmpcdp(:) = sigcme(:,jb,jb,spin) 397 sigc = pade(Sr%nomega_i, Sr%omega_i, tmpcdp, zz) 398 dsigc = dpade(Sr%nomega_i, Sr%omega_i, tmpcdp, zz) 399 else 400 tmpcdp(:) = CONJG(sigcme(:,jb,jb,spin)) 401 sigc = pade(Sr%nomega_i, CONJG(Sr%omega_i), tmpcdp, zz) 402 dsigc = dpade(Sr%nomega_i, CONJG(Sr%omega_i), tmpcdp, zz) 403 end if 404 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) & 405 + sigc - zz 406 if (ABS(ctdpc) < NR_ABS_ROOT_ERR) then 407 converged=.TRUE.; EXIT 408 end if 409 dct = dsigc - one 410 zz = newrap_step(zz, ctdpc, dct) 411 end do 412 413 if (.not. converged) then 414 write(msg,'(a,i0,3a,f8.4,a,f8.4)')& 415 'Newton-Raphson method not converged after ',NR_MAX_NITER,' iterations. ',ch10,& 416 'Absolute Error = ',ABS(ctdpc),' > ',NR_ABS_ROOT_ERR 417 ABI_WARNING(msg) 418 end if 419 420 ! Store the final result TODO re-shift everything according to efermi 421 Sr%egw(jb,sk_ibz,spin) = zz 422 Sr%degw(jb,sk_ibz,spin) = Sr%egw(jb,sk_ibz,spin) - Sr%e0(jb,sk_ibz,spin) 423 Sr%sigmee(jb,sk_ibz,spin) = Sr%sigxme(jb,sk_ibz,spin) + sigc 424 425 ! Spectra of Sigma, remember that Sr%nomega_r does not contains the frequencies 426 ! used to evaluate the derivative each frequency is obtained using the pade_expression 427 ! In sigma indeed we have: 428 ! nomega_sigc=Sr%nomega_r+Sr%nomega4sd; if (mod10==SIG_GW_AC) nomega_sigc=Sr%nomega_i 429 do io=1,Sr%nomega_r 430 zz=Sr%omega_r(io) 431 if (REAL(zz) > zero) then 432 tmpcdp(:) = sigcme(:,jb,jb,spin) 433 Sr%sigcme(jb,sk_ibz,io,spin) = pade(Sr%nomega_i, Sr%omega_i, tmpcdp, zz) 434 else 435 tmpcdp(:) = CONJG(sigcme(:,jb,jb,spin)) 436 Sr%sigcme(jb,sk_ibz,io,spin) = pade(Sr%nomega_i, CONJG(Sr%omega_i), tmpcdp, zz) 437 end if 438 Sr%sigxcme(jb,sk_ibz,io,spin) = Sr%sigxme(jb,sk_ibz,spin) + Sr%sigcme(jb,sk_ibz,io,spin) 439 end do 440 441 ! Save sigma values along the imaginary axis 442 do iab=1,Sr%nsig_ab 443 is_idx=spin; if (Sr%nsig_ab > 1) is_idx = iab 444 do io=1,Sr%nomega_i 445 Sr%sigcmesi (jb,sk_ibz,io,is_idx) = sigcme(io,jb,jb,is_idx) 446 Sr%sigxcmesi(jb,sk_ibz,io,is_idx) = Sr%sigxme(jb,sk_ibz,is_idx) + Sr%sigcmesi(jb,sk_ibz,io,is_idx) 447 end do 448 end do 449 450 ABI_FREE(tmpcdp) 451 452 end do !jb 453 end do !is 454 end if ! Analytic continuation. 455 456 ! === Diagonalize the QP Hamiltonian (forced to be Hermitian) === 457 ! Calculate Sr%en_qp_diago and Sr%eigvec_qp to be written in the QPS file. 458 ! TODO in case of AC results are wrong. 459 460 if (mod10 /= 1) then 461 462 ABI_MALLOC(hhartree, (ib1:ib2,ib1:ib2,Sr%nsppol*Sr%nsig_ab)) 463 hhartree = Sr%hhartree(ib1:ib2,ib1:ib2,sk_ibz,:) 464 465 ! If non self-consistent erase all off-diagonal elements 466 if (Sigp%gwcalctyp<20) then 467 do jb=ib1,ib2 468 do kb=ib1,ib2 469 if (jb == kb) CYCLE 470 hhartree(jb,kb,:) = czero 471 end do 472 end do 473 end if 474 475 ABI_MALLOC(htotal, (ib1:ib2,ib1:ib2,Sr%nsppol*Sr%nsig_ab)) 476 do spin=1,Sr%nsppol*Sr%nsig_ab 477 do jb=ib1,ib2 478 do kb=ib1,ib2 479 htotal(kb,jb,spin) = hhartree(kb,jb,spin) + Sr%x_mat(kb,jb,sk_ibz,spin) + sigcme(ie0,kb,jb,spin) 480 end do 481 end do 482 end do 483 484 ! Get the Hermitian part of htotal 485 ! In the noncollinear case A_{12}^{ab} = A_{21}^{ba}^* if A is Hermitian. 486 ABI_MALLOC(h_tmp1,(ib1:ib2,ib1:ib2)) 487 ABI_MALLOC(h_tmp2,(ib1:ib2,ib1:ib2)) 488 489 nsploop=Sr%nsppol; if (Sr%nsig_ab/=1) nsploop=2 490 do spin=1,nsploop 491 h_tmp1 = CONJG(htotal(:,:,spin)) 492 h_tmp2 = TRANSPOSE(h_tmp1) 493 h_tmp1 = htotal(:,:,spin) 494 htotal(:,:,spin)= half * (h_tmp1 + h_tmp2) 495 end do 496 497 ! Print the different matrix elements of sigma if QPSC and prtvol>9 498 if (Sigp%gwcalctyp >=20 .and. mod10 /= 1 .and. prtvol>9 .and. my_rank==master) then 499 call print_sigma_melems(ikcalc,ib1,ib2,Sr%nsppol*Sr%nsig_ab,htotal,hhartree,& 500 Sr%x_mat(ib1:ib2,ib1:ib2,sk_ibz,:),sigcme(ie0,:,:,:),Dtfil%filnam_ds(4)) 501 end if 502 503 if (Sr%nsig_ab==4) then 504 h_tmp1 = CONJG(htotal(:,:,4)) 505 h_tmp2 = TRANSPOSE(h_tmp1) 506 h_tmp1 = htotal(:,:,3) 507 htotal(:,:,3)= half * (h_tmp1 + h_tmp2) 508 509 h_tmp1 = CONJG(htotal(:,:,3)) 510 h_tmp2 = TRANSPOSE(h_tmp1) 511 htotal(:,:,4) = h_tmp2 512 end if 513 514 ! Solve Herm(htotal)*U = E*U 515 ld_matrix = ib2 - ib1 + 1 516 ABI_MALLOC(hdp, (ld_matrix,ld_matrix)) 517 ABI_MALLOC(eig, (ld_matrix)) 518 519 do spin=1,Sr%nsppol 520 if (Sr%nsig_ab==1) then 521 hdp=htotal(ib1:ib2,ib1:ib2,spin) 522 else 523 hdp = SUM(htotal(ib1:ib2,ib1:ib2,:), DIM=3) 524 end if 525 call xheev("Vectors","Upper", ld_matrix, hdp, eig) 526 527 Sr%eigvec_qp(ib1:ib2,ib1:ib2,sk_ibz,spin)=hdp(:,:) 528 Sr%en_qp_diago(ib1:ib2,sk_ibz,spin)=eig(:) 529 end do 530 531 ABI_FREE(hdp) 532 ABI_FREE(eig) 533 ABI_FREE(htotal) 534 ABI_FREE(hhartree) 535 ABI_FREE(h_tmp1) 536 ABI_FREE(h_tmp2) 537 538 end if 539 540 call timab(490,2,tsec) 541 542 DBG_EXIT("COLL") 543 544 end subroutine solve_dyson