TABLE OF CONTENTS
ABINIT/setup_sigma [ Functions ]
NAME
setup_sigma
FUNCTION
Initialize the data type containing parameters for a sigma calculation.
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 .
INPUTS
acell(3)=length scales of primitive translations (bohr) wfk_fname=Name of the WFK file. Dtset<type(dataset_type)>=all input variables for this dataset Dtfil<type(datafiles_type)>=variables related to files rprim(3,3)=dimensionless real space primitive translations ngfft(18)=information on the (fine) FFT grid used for the density. Psps <Pseudopotential_type)>=Info on pseudopotential, only for consistency check of the WFK file
OUTPUT
Sigp<sigparams_t>=Parameters governing the self-energy calculation. Kmesh <kmesh_t>=Structure describing the k-point sampling. Qmesh <kmesh_t>=Structure describing the q-point sampling. Cryst<crystal_t>=Info on unit cell and symmetries. Gsph_Max<gsphere_t>=Info on the G-sphere Gsph_c<gsphere_t>=Info on the G-sphere for W and Sigma_c Gsph_x<gsphere_t>=Info on the G-sphere for and Sigma_x Hdr_wfk<hdr_type>=The header of the WFK file Hdr_out<hdr_type>=The header to be used for the results of sigma calculations. Vcp<vcoul_t>= Datatype gathering information on the coulombian interaction and the cutoff technique. Er<Epsilonm1_results>=Datatype storing data used to construct the screening (partially Initialized in OUTPUT) KS_BSt<ebands_t>=The KS energies and occupation factors. gwc_ngfft(18), gwx_ngfft(18)= FFT meshes for the oscillator strengths used for the correlated and the exchange part of the self-energy, respectively. comm=MPI communicator.
PARENTS
sigma
CHILDREN
SOURCE
48 #if defined HAVE_CONFIG_H 49 #include "config.h" 50 #endif 51 52 #include "abi_common.h" 53 54 subroutine setup_sigma(codvsn,wfk_fname,acell,rprim,ngfftf,Dtset,Dtfil,Psps,Pawtab,& 55 & gwx_ngfft,gwc_ngfft,Hdr_wfk,Hdr_out,Cryst,Kmesh,Qmesh,KS_BSt,Gsph_Max,Gsph_x,Gsph_c,Vcp,Er,Sigp,comm) 56 57 use defs_basis 58 use defs_datatypes 59 use defs_abitypes 60 use defs_wvltypes 61 use m_profiling_abi 62 use m_errors 63 use m_xmpi 64 use m_nctk 65 use m_hdr 66 67 use m_gwdefs, only : GW_Q0_DEFAULT, SIG_GW_AC, sigparams_t, sigma_is_herm, sigma_needs_w 68 use m_io_tools, only : file_exists 69 use m_fstrings, only : basename, sjoin, ktoa, ltoa 70 use m_geometry, only : mkrdim, metric 71 use m_crystal, only : crystal_print, idx_spatial_inversion, crystal_t 72 use m_crystal_io, only : crystal_from_hdr 73 use m_bz_mesh, only : kmesh_t, kmesh_init, has_BZ_item, isamek, get_ng0sh, kmesh_print,& 74 & get_bz_item, has_IBZ_item, find_qmesh 75 use m_ebands, only : ebands_init, enclose_degbands, get_valence_idx, ebands_update_occ, ebands_report_gap, & 76 & get_gaps, gaps_free, gaps_t, gaps_print 77 use m_vcoul, only : vcoul_t, vcoul_init, vcoul_free 78 use m_fft_mesh, only : setmesh 79 use m_gsphere, only : gsphere_t, gsph_init, merge_and_sort_kg, gsph_extend, setshells 80 use m_screening, only : init_er_from_file, epsilonm1_results 81 use m_pawtab, only : pawtab_type 82 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free 83 use m_io_kss, only : make_gvec_kss 84 use m_wfk, only : wfk_read_eigenvalues 85 use m_xcdata, only : get_xclevel 86 87 !This section has been created automatically by the script Abilint (TD). 88 !Do not modify the following lines by hand. 89 #undef ABI_FUNC 90 #define ABI_FUNC 'setup_sigma' 91 use interfaces_14_hidewrite 92 use interfaces_56_io_mpi 93 use interfaces_70_gw, except_this_one => setup_sigma 94 !End of the abilint section 95 96 implicit none 97 98 !Arguments ------------------------------------ 99 !scalars 100 integer,intent(in) :: comm 101 character(len=6),intent(in) :: codvsn 102 character(len=*),intent(in) :: wfk_fname 103 type(Datafiles_type),intent(in) :: Dtfil 104 type(Dataset_type),intent(inout) :: Dtset 105 type(Pseudopotential_type),intent(in) :: Psps 106 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw) 107 type(sigparams_t),intent(out) :: Sigp 108 type(Epsilonm1_results),intent(out) :: Er 109 type(ebands_t),intent(out) :: KS_BSt 110 type(kmesh_t),intent(out) :: Kmesh,Qmesh 111 type(crystal_t),intent(out) :: Cryst 112 type(gsphere_t),intent(out) :: Gsph_Max,Gsph_x,Gsph_c 113 type(Hdr_type),intent(out) :: Hdr_wfk,Hdr_out 114 type(vcoul_t),intent(out) :: Vcp 115 !arrays 116 integer,intent(in) :: ngfftf(18) 117 integer,intent(out) :: gwc_ngfft(18),gwx_ngfft(18) 118 real(dp),intent(in) :: acell(3),rprim(3,3) 119 120 !Local variables------------------------------- 121 !scalars 122 integer,parameter :: pertcase0=0,master=0 123 integer :: bantot,enforce_sym,gwcalctyp,ib,ibtot,icutcoul_eff,ii,ikcalc,ikibz,io,isppol,itypat,jj,method 124 integer :: mod10,mqmem,mband,ng_kss,nsheps,ikcalc2bz,ierr,gap_err,ng 125 integer :: gwc_nfftot,gwx_nfftot,nqlwl,test_npwkss,my_rank,nprocs,ik,nk_found,ifo,timrev,usefock_ixc 126 integer :: iqbz,isym,iq_ibz,itim,ic,pinv,ig1,ng_sigx,spin,gw_qprange,ivcoul_init,nvcoul_init,xclevel_ixc 127 real(dp),parameter :: OMEGASIMIN=0.01d0,tol_enediff=0.001_dp*eV_Ha 128 real(dp) :: domegas,domegasi,ucvol,rcut 129 logical,parameter :: linear_imag_mesh=.TRUE. 130 logical :: ltest,remove_inv,changed,found 131 character(len=500) :: msg 132 character(len=fnlen) :: fname,fcore,string 133 type(wvl_internal_type) :: wvl 134 type(gaps_t) :: gaps 135 !arrays 136 integer :: ng0sh_opt(3),G0(3),q_umklp(3),kpos(6) 137 integer,allocatable :: npwarr(:),val_indeces(:,:) 138 integer,pointer :: gvec_kss(:,:),gsphere_sigx_p(:,:) 139 integer,pointer :: test_gvec_kss(:,:) 140 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),sq(3),q_bz(3),gamma_point(3,1) 141 real(dp),pointer :: energies_p(:,:,:) 142 real(dp),allocatable :: doccde(:),eigen(:),occfact(:),qlwl(:,:) 143 type(Pawrhoij_type),allocatable :: Pawrhoij(:) 144 type(vcoul_t) :: Vcp_ks 145 146 ! ************************************************************************* 147 148 DBG_ENTER('COLL') 149 150 ! Check for calculations that are not implemented 151 ltest=ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)==Dtset%nband(1)) 152 ABI_CHECK(ltest,'Dtset%nband(:) must be constant') 153 154 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 155 156 ! Basic parameters 157 Sigp%ppmodel = Dtset%ppmodel 158 Sigp%gwcalctyp = Dtset%gwcalctyp 159 Sigp%nbnds = Dtset%nband(1) 160 Sigp%symsigma = Dtset%symsigma 161 Sigp%zcut = Dtset%zcut 162 Sigp%mbpt_sciss = Dtset%mbpt_sciss 163 164 timrev= 2 ! This information is not reported in the header 165 ! 1 => do not use time-reversal symmetry 166 ! 2 => take advantage of time-reversal symmetry 167 if (any(dtset%kptopt == [3, 4])) timrev = 1 168 169 ! === For HF, SEX or COHSEX use Hybertsen-Louie PPM (only $\omega=0$) === 170 ! * Use fake screening for HF. 171 ! FIXME Why, we should not redefine Sigp%ppmodel 172 gwcalctyp=Sigp%gwcalctyp 173 mod10 =MOD(Sigp%gwcalctyp,10) 174 if (mod10==5.or.mod10==6.or.mod10==7) Sigp%ppmodel=2 175 if (mod10<5.and.MOD(Sigp%gwcalctyp,1)/=1) then ! * One shot GW (PPM or contour deformation). 176 if (Dtset%nomegasrd==1) then ! avoid division by zero! 177 Sigp%nomegasrd =1 178 Sigp%maxomega4sd=zero 179 Sigp%deltae =zero 180 else 181 Sigp%nomegasrd = Dtset%nomegasrd 182 Sigp%maxomega4sd = Dtset%omegasrdmax 183 Sigp%deltae = (2*Sigp%maxomega4sd)/(Sigp%nomegasrd-1) 184 endif 185 else 186 ! For AC no need to evaluate derivative by finite differences. 187 Sigp%nomegasrd =1 188 Sigp%maxomega4sd=zero 189 Sigp%deltae =zero 190 end if 191 192 ! For analytic continuation define the number of imaginary frequencies for Sigma 193 ! Tests show than more than 12 freqs in the Pade approximant worsen the results! 194 Sigp%nomegasi=0 195 196 if (mod10==1) then 197 Sigp%nomegasi =Dtset%nomegasi 198 Sigp%omegasimax=Dtset%omegasimax 199 Sigp%omegasimin=OMEGASIMIN 200 write(msg,'(4a,i3,2(2a,f8.3),a)')ch10,& 201 & ' Parameters for analytic continuation : ',ch10,& 202 & ' number of imaginary frequencies for sigma = ',Sigp%nomegasi,ch10,& 203 & ' min frequency for sigma on imag axis [eV] = ',Sigp%omegasimin*Ha_eV,ch10,& 204 & ' max frequency for sigma on imag axis [eV] = ',Sigp%omegasimax*Ha_eV,ch10 205 call wrtout(std_out,msg,'COLL') 206 207 !TODO this should not be done here but in init_sigma_t 208 ABI_MALLOC(Sigp%omegasi,(Sigp%nomegasi)) 209 210 if (linear_imag_mesh) then ! * Linear mesh along the imaginary axis. 211 domegasi=Sigp%omegasimax/(Sigp%nomegasi-1) 212 do io=1,Sigp%nomegasi 213 Sigp%omegasi(io)=CMPLX(zero,(io-1)*domegasi) 214 end do 215 else ! * Logarithmic mesh along the imaginary axis. 216 MSG_ERROR("AC + log mesh not implemented") 217 !domegasi=(Sigp%omegasimax/Sigp%omegasimin)**(one/(Sigp%nomegasi-1)) 218 !Sigp%omegasi(1)=czero; ldi=domegasi 219 !do io=2,Sigp%nomegasi 220 ! omega(io)=CMPLX(zero,ldi*Sigp%omegasimin) 221 ! Sigp%omegasi(io)=ldi*domegasi 222 !end do 223 end if 224 225 write(msg,'(4a)')ch10,& 226 & ' setup_sigma : calculating Sigma(iw)',& 227 & ' at imaginary frequencies [eV] (Fermi Level set to 0) ',ch10 228 call wrtout(std_out,msg,'COLL') 229 call wrtout(ab_out,msg,'COLL') 230 do io=1,Sigp%nomegasi 231 write(msg,'(2(f10.3,2x))')Sigp%omegasi(io)*Ha_eV 232 call wrtout(std_out,msg,'COLL') 233 call wrtout(ab_out,msg,'COLL') 234 end do 235 236 ltest=(Sigp%omegasimax>0.1d-4.and.Sigp%nomegasi>0) 237 ABI_CHECK(ltest,'Wrong value of omegasimax or nomegasi') 238 if (Sigp%gwcalctyp/=1) then ! only one shot GW is allowed for AC. 239 MSG_ERROR("SC-GW with analytic continuation is not coded") 240 end if 241 end if 242 243 if (Sigp%symsigma/=0.and.gwcalctyp>=20) then 244 MSG_WARNING("SC-GW with symmetries is still under development. Use at your own risk!") 245 end if 246 247 ! Setup parameters for Spectral function. 248 if (Dtset%gw_customnfreqsp/=0) then 249 Sigp%nomegasr = Dtset%gw_customnfreqsp 250 MSG_WARNING('Custom grid for spectral function specified. Assuming experienced user.') 251 if (Dtset%gw_customnfreqsp/=0) then 252 Dtset%nfreqsp = Dtset%gw_customnfreqsp 253 MSG_WARNING('nfreqsp has been set to the same number as gw_customnfreqsp') 254 end if 255 else 256 Sigp%nomegasr =Dtset%nfreqsp 257 Sigp%minomega_r=Dtset%freqspmin 258 Sigp%maxomega_r=Dtset%freqspmax 259 end if 260 261 if (Sigp%nomegasr>0) then 262 if (Dtset%gw_customnfreqsp==0) then 263 ! Check 264 if (Sigp%minomega_r >= Sigp%maxomega_r) then 265 MSG_ERROR('freqspmin must be smaller than freqspmax!') 266 end if 267 if(Sigp%nomegasr==1) then 268 domegas=0.d0 269 else 270 domegas=(Sigp%maxomega_r-Sigp%minomega_r)/(Sigp%nomegasr-1) 271 endif 272 !TODO this should be moved to Sr% and done in init_sigma_t 273 ABI_MALLOC(Sigp%omega_r,(Sigp%nomegasr)) 274 do io=1,Sigp%nomegasr 275 Sigp%omega_r(io) = CMPLX(Sigp%minomega_r + domegas*(io-1),zero) 276 end do 277 write(msg,'(4a,i8,3(2a,f8.3),a)')ch10,& 278 & ' Parameters for the calculation of the spectral function : ',ch10,& 279 & ' Number of points = ',Sigp%nomegasr,ch10,& 280 & ' Min frequency [eV] = ',Sigp%minomega_r*Ha_eV,ch10,& 281 & ' Max frequency [eV] = ',Sigp%maxomega_r*Ha_eV,ch10,& 282 & ' Frequency step [eV] = ',domegas*Ha_eV,ch10 283 call wrtout(std_out,msg,'COLL') 284 else 285 Sigp%minomega_r = MINVAL(Dtset%gw_freqsp(:)) 286 Sigp%maxomega_r = MAXVAL(Dtset%gw_freqsp(:)) 287 !TODO this should be moved to Sr% and done in init_sigma_t 288 ABI_MALLOC(Sigp%omega_r,(Sigp%nomegasr)) 289 do io=1,Sigp%nomegasr 290 Sigp%omega_r(io) = CMPLX(Dtset%gw_freqsp(io),zero) 291 end do 292 write(msg,'(4a,i8,2(2a,f8.3),3a)')ch10,& 293 & ' Parameters for the calculation of the spectral function : ',ch10,& 294 & ' Number of points = ',Sigp%nomegasr,ch10,& 295 & ' Min frequency [eV] = ',Sigp%minomega_r*Ha_eV,ch10,& 296 & ' Max frequency [eV] = ',Sigp%maxomega_r*Ha_eV,ch10,& 297 & ' A custom set of frequencies is used! See the input file for values.',ch10 298 call wrtout(std_out,msg,'COLL') 299 end if 300 else 301 !In indefo all these quantities are set to zero 302 !Sigp%nomegasr=1 303 !allocate(Sigp%omega_r(Sigp%nomegasr)) 304 !Sigp%omega_r(1)=0 305 end if 306 307 ! Dimensional primitive translations rprimd (from input), gprimd, metrics and unit cell volume 308 call mkrdim(acell,rprim,rprimd) 309 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 310 311 Sigp%npwwfn=Dtset%npwwfn 312 Sigp%npwx =Dtset%npwsigx 313 314 ! Read parameters of the WFK, verifify them and retrieve all G-vectors. 315 call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm) 316 mband = MAXVAL(Hdr_wfk%nband) 317 318 remove_inv = .FALSE. 319 call hdr_vs_dtset(Hdr_wfk,Dtset) 320 321 test_npwkss = 0 322 call make_gvec_kss(Dtset%nkpt,Dtset%kptns,Hdr_wfk%ecut_eff,Dtset%symmorphi,Dtset%nsym,Dtset%symrel,Dtset%tnons,& 323 & gprimd,Dtset%prtvol,test_npwkss,test_gvec_kss,ierr) 324 ABI_CHECK(ierr==0,"Fatal error in make_gvec_kss") 325 326 ABI_MALLOC(gvec_kss,(3,test_npwkss)) 327 gvec_kss = test_gvec_kss 328 ng_kss = test_npwkss 329 330 ng = MIN(SIZE(gvec_kss,DIM=2),SIZE(test_gvec_kss,DIM=2)) 331 ierr = 0 332 do ig1=1,ng 333 if (ANY(gvec_kss(:,ig1)/=test_gvec_kss(:,ig1))) then 334 ierr=ierr+1 335 write(std_out,*)" gvec_kss ",ig1,"/",ng,gvec_kss(:,ig1),test_gvec_kss(:,ig1) 336 end if 337 end do 338 ABI_CHECK(ierr == 0, "Mismatch between gvec_kss and test_gvec_kss") 339 340 ABI_FREE(test_gvec_kss) 341 342 ! Get important dimensions from the WFK header 343 Sigp%nsppol =Hdr_wfk%nsppol 344 Sigp%nspinor=Hdr_wfk%nspinor 345 Sigp%nsig_ab=Hdr_wfk%nspinor**2 ! TODO Is it useful calculating only diagonal terms? 346 347 if (Sigp%nbnds>mband) then 348 Sigp%nbnds =mband 349 Dtset%nband(:)=mband 350 Dtset%mband =MAXVAL(Dtset%nband) 351 write(msg,'(3a,i4,a)')& 352 & 'Number of bands found less then required',ch10,& 353 & 'calculation will proceed with nbnds = ',mband,ch10 354 MSG_WARNING(msg) 355 end if 356 357 ! Check input 358 if (Sigp%ppmodel==3.or.Sigp%ppmodel==4) then 359 if (gwcalctyp>=10) then 360 write(msg,'(a,i3,a)')' The ppmodel chosen and gwcalctyp ',Dtset%gwcalctyp,' are not compatible. ' 361 MSG_ERROR(msg) 362 end if 363 if (Sigp%nspinor==2) then 364 write(msg,'(a,i3,a)')' The ppmodel chosen and nspinor ',Sigp%nspinor,' are not compatible. ' 365 MSG_ERROR(msg) 366 end if 367 end if 368 369 ! Create crystal_t data type 370 call crystal_from_hdr(Cryst,Hdr_wfk,timrev,remove_inv) 371 call crystal_print(Cryst) 372 373 if (Sigp%npwwfn>ng_kss) then ! cannot use more G"s for the wfs than those stored on file 374 Sigp%npwwfn =ng_kss 375 Dtset%npwwfn =ng_kss 376 write(msg,'(2a,(a,i8,a))')& 377 & 'Number of G-vectors for WFS found in the KSS file is less than required',ch10,& 378 & 'calculation will proceed with npwwfn = ',Sigp%npwwfn,ch10 379 MSG_WARNING(msg) 380 end if 381 382 if (Sigp%npwx>ng_kss) then 383 ! Have to recalcuate the (large) sphere for Sigma_x. 384 pinv=1; if (remove_inv.and.Cryst%timrev==2) pinv=-1 385 gamma_point(:,1) = (/zero,zero,zero/); nullify(gsphere_sigx_p) 386 387 call merge_and_sort_kg(1,gamma_point,Dtset%ecutsigx,Cryst%nsym,pinv,Cryst%symrel,& 388 & Cryst%gprimd,gsphere_sigx_p,Dtset%prtvol) 389 390 ng_sigx=SIZE(gsphere_sigx_p,DIM=2) 391 Sigp%npwx = ng_sigx 392 Dtset%npwsigx = ng_sigx 393 394 write(msg,'(2a,(a,i8,a))')& 395 & 'Number of G-vectors for Sigma_x found in the KSS file is less than required',ch10,& 396 & 'calculation will proceed with npwsigx = ',Sigp%npwx,ch10 397 MSG_WARNING(msg) 398 399 ltest = (Sigp%npwx >= ng_kss) 400 ABI_CHECK(ltest,"Sigp%npwx<ng_kss!") 401 402 ! * Fill gvec_kss with larger sphere. 403 ABI_FREE(gvec_kss) 404 ABI_MALLOC(gvec_kss,(3,Sigp%npwx)) 405 gvec_kss = gsphere_sigx_p 406 ABI_FREE(gsphere_sigx_p) 407 end if 408 409 ! Set up of the k-points and tables in the whole BZ === 410 ! TODO Recheck symmorphy and inversion 411 call kmesh_init(Kmesh,Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.FALSE.) 412 !call kmesh_init(Kmesh,Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.TRUE.) 413 414 ! Some required information are not filled up inside kmesh_init 415 ! So doing it here, even though it is not clean 416 Kmesh%kptrlatt(:,:) =Dtset%kptrlatt(:,:) 417 Kmesh%nshift =Dtset%nshiftk 418 ABI_ALLOCATE(Kmesh%shift,(3,Kmesh%nshift)) 419 Kmesh%shift(:,:) =Dtset%shiftk(:,1:Dtset%nshiftk) 420 421 call kmesh_print(Kmesh,"K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL") 422 call kmesh_print(Kmesh,"K-mesh for the wavefunctions",ab_out, 0, "COLL") 423 424 ! === Initialize the band structure datatype === 425 ! * Copy WFK energies and occupations up to Sigp%nbnds==Dtset%nband(:) 426 bantot=SUM(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)) 427 ABI_MALLOC(doccde,(bantot)) 428 ABI_MALLOC(eigen,(bantot)) 429 ABI_MALLOC(occfact,(bantot)) 430 doccde(:)=zero; eigen(:)=zero; occfact(:)=zero 431 432 jj=0; ibtot=0 433 do isppol=1,Dtset%nsppol 434 do ikibz=1,Dtset%nkpt 435 do ib=1,Hdr_wfk%nband(ikibz+(isppol-1)*Dtset%nkpt) 436 ibtot=ibtot+1 437 if (ib<=Sigp%nbnds) then 438 jj=jj+1 439 occfact(jj)=Hdr_wfk%occ(ibtot) 440 eigen (jj)=energies_p(ib,ikibz,isppol) 441 end if 442 end do 443 end do 444 end do 445 ABI_FREE(energies_p) 446 447 ! Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of 448 ! symmetry operations in the old GW code (symmorphy and inversion) 449 ltest=(ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz))<tol6)) 450 ABI_CHECK(ltest,'Mismatch between Dtset%wtk and Kmesh%wt') 451 452 ABI_MALLOC(npwarr,(Dtset%nkpt)) 453 npwarr(:)=Sigp%npwwfn 454 455 call ebands_init(bantot,KS_BSt,Dtset%nelect,doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,& 456 & Kmesh%nibz,npwarr,Dtset%nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,& 457 & dtset%charge, dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig,& 458 & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 459 460 ABI_FREE(doccde) 461 ABI_FREE(eigen) 462 ABI_FREE(npwarr) 463 464 ! Calculate KS occupation numbers and ks_vbk(nkibz,nsppol) ==== 465 ! ks_vbk gives the (valence|last Fermi band) index for each k and spin. 466 ! spinmagntarget is passed to fermi.F90 to fix the problem with newocc in case of magnetic metals 467 call ebands_update_occ(KS_BSt,Dtset%spinmagntarget,prtvol=0) 468 469 gap_err = get_gaps(KS_BSt, gaps) 470 call gaps_print(gaps, unit=std_out) 471 call ebands_report_gap(KS_BSt, unit=std_out) 472 473 ABI_MALLOC(val_indeces,(KS_BSt%nkpt,KS_BSt%nsppol)) 474 val_indeces = get_valence_idx(KS_BSt) 475 476 ! Create Sigma header 477 ! TODO Fix problems with symmorphy and k-points 478 call hdr_init(KS_BSt,codvsn,Dtset,Hdr_out,Pawtab,pertcase0,Psps,wvl) 479 480 ! Get Pawrhoij from the header of the WFK file 481 ABI_DT_MALLOC(Pawrhoij,(Cryst%natom*Dtset%usepaw)) 482 if (Dtset%usepaw==1) then 483 call pawrhoij_alloc(Pawrhoij,1,Dtset%nspden,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab) 484 call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij) 485 end if 486 call hdr_update(hdr_out,bantot,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1)) 487 488 ABI_FREE(occfact) 489 call pawrhoij_free(Pawrhoij) 490 ABI_DT_FREE(Pawrhoij) 491 492 ! =========================================================== 493 ! ==== Setup of k-points and bands for the GW corrections ==== 494 ! =========================================================== 495 ! * maxbdgw and minbdgw are the Max and min band index for GW corrections over k-points. 496 ! They are used to dimension the wavefunctions and to calculate the matrix elements. 497 ! 498 if (dtset%nkptgw==0) then 499 ! Use qp_range to select the interesting k-points and the corresponing bands. 500 ! 501 ! 0 --> Compute the QP corrections only for the fundamental and the optical gap. 502 ! +num --> Compute the QP corrections for all the k-points in the irreducible zone and include `num` 503 ! bands above and below the Fermi level. 504 ! -num --> Compute the QP corrections for all the k-points in the irreducible zone. 505 ! Include all occupied states and `num` empty states. 506 507 call wrtout(std_out, "nkptgw == 0 ==> Automatic selection of k-points and bands for the corrections.") 508 if (gap_err /=0 .and. dtset%gw_qprange==0) then 509 msg = "Problem while computing the fundamental and optical gap (likely metal). Will replace gw_qprange=0 with gw_qprange=1" 510 MSG_WARNING(msg) 511 dtset%gw_qprange = 1 512 end if 513 gw_qprange = dtset%gw_qprange 514 515 if (dtset%ucrpa>0) then 516 dtset%nkptgw=Kmesh%nbz 517 Sigp%nkptgw =dtset%nkptgw 518 ABI_MALLOC(Sigp%kptgw,(3,Sigp%nkptgw)) 519 ABI_MALLOC(Sigp%minbnd,(Sigp%nkptgw,Sigp%nsppol)) 520 ABI_MALLOC(Sigp%maxbnd,(Sigp%nkptgw,Sigp%nsppol)) 521 Sigp%kptgw(:,:)=Kmesh%bz(:,:) 522 Sigp%minbnd=1 523 Sigp%maxbnd=Sigp%nbnds 524 525 else if (gw_qprange/=0) then 526 ! Include all the k-points in the IBZ. 527 ! Note that kptgw == ebands%kptns so we can use a single ik index in the loop over k-points 528 ! No need to map kptgw onto ebands%kptns. 529 dtset%nkptgw=Kmesh%nibz 530 Sigp%nkptgw =dtset%nkptgw 531 ABI_MALLOC(Sigp%kptgw,(3,Sigp%nkptgw)) 532 ABI_MALLOC(Sigp%minbnd,(Sigp%nkptgw,Sigp%nsppol)) 533 ABI_MALLOC(Sigp%maxbnd,(Sigp%nkptgw,Sigp%nsppol)) 534 Sigp%kptgw(:,:)=Kmesh%ibz(:,:) 535 Sigp%minbnd=1 536 Sigp%maxbnd=Sigp%nbnds 537 538 if (gw_qprange>0) then 539 ! All k-points: Add buffer of bands above and below the Fermi level. 540 do spin=1,Sigp%nsppol 541 do ik=1,Sigp%nkptgw 542 Sigp%minbnd(ik,spin) = MAX(val_indeces(ik,spin) - gw_qprange, 1) 543 Sigp%maxbnd(ik,spin) = MIN(val_indeces(ik,spin) + gw_qprange + 1, Sigp%nbnds) 544 end do 545 end do 546 547 else 548 ! All k-points: include all occupied states and -gw_qprange empty states. 549 Sigp%minbnd = 1 550 do spin=1,Sigp%nsppol 551 do ik=1,Sigp%nkptgw 552 Sigp%maxbnd(ik,spin) = MIN(val_indeces(ik,spin) - gw_qprange, Sigp%nbnds) 553 end do 554 end do 555 end if 556 557 else 558 ! gw_qprange is not specified in the input. 559 ! Include the optical and the fundamental KS gap. 560 ! The main problem here is that kptgw and nkptgw do not depend on the spin and therefore 561 ! we have compute the union of the k-points where the fundamental and the optical gaps are located. 562 ! 563 ! Find the list of `interesting` kpoints. 564 ABI_CHECK(gap_err == 0, "gw_qprange 0 cannot be used because I cannot find the gap (gap_err !=0)") 565 nk_found = 1; kpos(1) = gaps%fo_kpos(1,1) 566 567 do spin=1,Sigp%nsppol 568 do ifo=1,3 569 ik = gaps%fo_kpos(ifo, spin) 570 found = .FALSE.; jj = 0 571 do while (.not. found .and. jj < nk_found) 572 jj = jj + 1; found = (kpos(jj) == ik) 573 end do 574 if (.not. found) then 575 nk_found = nk_found + 1; kpos(nk_found) = ik 576 end if 577 end do 578 end do 579 580 ! Now we can define the list of k-points and the bands range. 581 dtset%nkptgw=nk_found 582 Sigp%nkptgw =dtset%nkptgw 583 584 ABI_MALLOC(Sigp%kptgw,(3,Sigp%nkptgw)) 585 ABI_MALLOC(Sigp%minbnd,(Sigp%nkptgw,Sigp%nsppol)) 586 ABI_MALLOC(Sigp%maxbnd,(Sigp%nkptgw,Sigp%nsppol)) 587 588 do ii=1,Sigp%nkptgw 589 ik = kpos(ii) 590 Sigp%kptgw(:,ii)=Kmesh%ibz(:,ik) 591 do spin=1,Sigp%nsppol 592 Sigp%minbnd(ii,spin) = val_indeces(ik,spin) 593 Sigp%maxbnd(ii,spin) = val_indeces(ik,spin) + 1 594 end do 595 end do 596 end if 597 598 else 599 ! Treat only the k-points and bands specified in the input file. 600 Sigp%nkptgw=dtset%nkptgw 601 ABI_MALLOC(Sigp%kptgw,(3,Sigp%nkptgw)) 602 ABI_MALLOC(Sigp%minbnd,(Sigp%nkptgw,Sigp%nsppol)) 603 ABI_MALLOC(Sigp%maxbnd,(Sigp%nkptgw,Sigp%nsppol)) 604 605 do spin=1,Sigp%nsppol 606 Sigp%minbnd(:,spin)=dtset%bdgw(1,:,spin) 607 Sigp%maxbnd(:,spin)=dtset%bdgw(2,:,spin) 608 end do 609 610 do ii=1,3 611 do ikcalc=1,Sigp%nkptgw 612 Sigp%kptgw(ii,ikcalc)=Dtset%kptgw(ii,ikcalc) 613 end do 614 end do 615 616 do spin=1,Sigp%nsppol 617 do ikcalc=1,Sigp%nkptgw 618 if (Dtset%bdgw(2,ikcalc,spin)>Sigp%nbnds) then 619 write(msg,'(a,2i0,2(a,i0),2a,i0)')& 620 & "For (k,s) ",ikcalc,spin," bdgw= ",Dtset%bdgw(2,ikcalc,spin), " > nbnds=",Sigp%nbnds,ch10,& 621 & "Calculation will continue with bdgw =",Sigp%nbnds 622 MSG_COMMENT(msg) 623 Dtset%bdgw(2,ikcalc,spin)=Sigp%nbnds 624 end if 625 end do 626 end do 627 628 end if 629 630 ! Make sure that all the degenerate states are included. 631 ! * We will have to average the GW corrections over degenerate states if symsigma=1 is used. 632 ! * KS states belonging to the same irreducible representation should be included in the basis set used for SCGW. 633 if (Sigp%symsigma/=0 .or. gwcalctyp>=10) then 634 do isppol=1,Sigp%nsppol 635 do ikcalc=1,Sigp%nkptgw 636 637 if (has_IBZ_item(Kmesh,Sigp%kptgw(:,ikcalc),ikibz,G0)) then 638 call enclose_degbands(KS_BSt,ikibz,isppol,Sigp%minbnd(ikcalc,isppol),Sigp%maxbnd(ikcalc,isppol),changed,tol_enediff) 639 if (changed) then 640 write(msg,'(2(a,i0),2a,2(1x,i0))')& 641 & "Not all the degenerate states at ikcalc= ",ikcalc,", spin= ",isppol,ch10,& 642 & "were included in the bdgw set. bdgw has been changed to: ",Sigp%minbnd(ikcalc,isppol),Sigp%maxbnd(ikcalc,isppol) 643 MSG_COMMENT(msg) 644 end if 645 else 646 MSG_ERROR(sjoin('k-point', ktoa(Sigp%kptgw(:,ikcalc)), 'not in IBZ')) 647 end if 648 649 end do 650 end do 651 end if 652 653 !if (.not. associated(Dtset%bdgw)) then 654 ! ABI_MALLOC(Dtset%bdgw, (2,Sigp%nkptgw,Sigp%nsppol)) 655 !end if 656 !do spin=1,Sigp%nsppol 657 ! Dtset%bdgw(1,:,spin) = Sigp%minbnd(:,spin) 658 ! Dtset%bdgw(2,:,spin) = Sigp%maxbnd(:,spin) 659 !end do 660 661 Sigp%minbdgw=MINVAL(Sigp%minbnd) 662 Sigp%maxbdgw=MAXVAL(Sigp%maxbnd) 663 664 ABI_MALLOC(Sigp%kptgw2bz,(Sigp%nkptgw)) 665 ! 666 !=== Check if the k-points are in the BZ === 667 !FB TODO Honestly the code is not able to treat k-points, which are not in the IBZ. 668 !This extension should require to change the code in different places. 669 !Therefore, one should by now prevent the user from calculating sigma for a k-point not in the IBZ. 670 671 do ikcalc=1,Sigp%nkptgw 672 if (has_BZ_item(Kmesh,Sigp%kptgw(:,ikcalc),ikcalc2bz,G0)) then 673 !found = has_IBZ_item(Kmesh,Sigp%kptgw(:,ikcalc),ikcalc2bz,G0) 674 Sigp%kptgw2bz(ikcalc) = ikcalc2bz 675 else 676 MSG_ERROR(sjoin('k-point:', ktoa(Sigp%kptgw(:,ikcalc)), 'not in the kbz set')) 677 end if 678 end do 679 680 ! Check if there are duplicated k-point in Sigp% 681 do ii=1,Sigp%nkptgw 682 do jj=ii+1,Sigp%nkptgw 683 if (isamek(Sigp%kptgw(:,ii),Sigp%kptgw(:,jj),G0)) then 684 write(msg,'(5a)')& 685 & 'kptgw contains duplicated k-points. This is not allowed since ',ch10,& 686 & 'the QP corrections for this k-point will be calculated more than once. ',ch10,& 687 & 'Check your input file. ' 688 MSG_ERROR(msg) 689 end if 690 end do 691 end do 692 ! 693 ! Warn the user if SCGW run and not all the k-points are included. 694 if (gwcalctyp>=10 .and. Sigp%nkptgw/=Hdr_wfk%nkpt) then 695 write(msg,'(3a,2(a,i0),2a)')ch10,& 696 & " COMMENT: In a self-consistent GW run, the QP corrections should be calculated for all the k-points of the KSS file ",ch10,& 697 & " but nkptgw= ",Sigp%nkptgw," and WFK nkpt= ",Hdr_wfk%nkpt,ch10,& 698 & " Assuming expert user. Execution will continue. " 699 call wrtout(ab_out,msg,"COLL") 700 end if 701 702 ! Setup of the table used in the case of SCGW on wavefunctions to reduce the number 703 ! of elements <i,kgw,s|\Sigma|j,kgw,s> that have to be calculated. No use of symmetries, except for Hermiticity. 704 call sigma_tables(Sigp,Kmesh) 705 706 ! === Read external file and initialize basic dimension of Er% === 707 ! TODO use mqmem as input variable instead of gwmem 708 709 ! === If required, use a matrix for $\Sigma_c$ which is smaller than that stored on file === 710 ! * By default the entire matrix is read and used, 711 ! * Define consistently npweps and ecuteps for \Sigma_c according the input 712 if (Dtset%npweps>0.or.Dtset%ecuteps>0) then 713 ! This should not happen: the Dtset array should not be modified after having been initialized. 714 if (Dtset%npweps>0) Dtset%ecuteps=zero 715 nsheps=0 716 call setshells(Dtset%ecuteps,Dtset%npweps,nsheps,Dtset%nsym,gmet,gprimd,Dtset%symrel,'eps',ucvol) 717 end if 718 719 mqmem=0; if (Dtset%gwmem/10==1) mqmem=1 720 721 if (Dtset%getscr/=0.or.Dtset%irdscr/=0) then 722 fname=Dtfil%fnameabi_scr 723 else if (Dtset%getsuscep/=0.or.Dtset%irdsuscep/=0) then 724 fname=Dtfil%fnameabi_sus 725 else 726 fname=Dtfil%fnameabi_scr 727 !FIXME this has to be cleaned, in tgw2_3 Dtset%get* and Dtset%ird* are not defined 728 !MSG_ERROR("getsuscep or irdsuscep are not defined") 729 end if 730 ! 731 ! === Setup of q-mesh in the whole BZ === 732 ! * Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed, 733 ! epsilon^-1(Sq) should be symmetrized in csigme using a different expression (G-G_o is needed) 734 ! 735 if (sigma_needs_w(Sigp)) then 736 if (.not. file_exists(fname)) then 737 fname = nctk_ncify(fname) 738 MSG_COMMENT(sjoin("File not found. Will try netcdf file:", fname)) 739 end if 740 741 call init_Er_from_file(Er,fname,mqmem,Dtset%npweps,comm) 742 743 Sigp%npwc=Er%npwe 744 if (Sigp%npwc>Sigp%npwx) then 745 Sigp%npwc=Sigp%npwx 746 MSG_COMMENT("Found npw_correlation > npw_exchange, Imposing npwc=npwx") 747 ! There is a good reason for doing so, see csigme.F90 and the size of the arrays 748 ! rhotwgp and rhotwgp: we need to define a max size and we opt for Sigp%npwx. 749 end if 750 Er%npwe=Sigp%npwc 751 Dtset%npweps=Er%npwe 752 call kmesh_init(Qmesh,Cryst,Er%nqibz,Er%qibz,Dtset%kptopt) 753 754 else 755 Er%npwe =1 756 Sigp%npwc =1 757 Dtset%npweps=1 758 call find_qmesh(Qmesh,Cryst,Kmesh) 759 ABI_MALLOC(Er%gvec,(3,1)) 760 Er%gvec(:,1) = (/0,0,0/) 761 end if 762 763 call kmesh_print(Qmesh,"Q-mesh for screening function",std_out,Dtset%prtvol,"COLL") 764 call kmesh_print(Qmesh,"Q-mesh for screening function",ab_out ,0 ,"COLL") 765 766 do iqbz=1,Qmesh%nbz 767 call get_BZ_item(Qmesh,iqbz,q_bz,iq_ibz,isym,itim,umklp=q_umklp) 768 769 if (ANY(q_umklp/=0)) then 770 sq = (3-2*itim)*MATMUL(Cryst%symrec(:,:,isym),Qmesh%ibz(:,iq_ibz)) 771 write(std_out,*) sq,Qmesh%bz(:,iqbz) 772 write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')& 773 & 'qpoint ',Qmesh%bz(:,iqbz),' is the symmetric of ',Qmesh%ibz(:,iq_ibz),ch10,& 774 & 'through operation ',Cryst%symrec(:,:,isym),' and itim ',itim,ch10,& 775 & 'however a non zero umklapp G_o vector is required and this is not yet allowed' 776 MSG_ERROR(msg) 777 end if 778 end do 779 ! 780 ! === Find optimal value for G-sphere enlargment due to oscillator matrix elements === 781 ! * Here I have to be sure that Qmesh%bz is always inside the BZ, not always true size bz is buggy 782 ! * -one is used because we loop over all the possibile differences, unlike screening 783 784 call get_ng0sh(Sigp%nkptgw,Sigp%kptgw,Kmesh%nbz,Kmesh%bz,Qmesh%nbz,Qmesh%bz,-one,ng0sh_opt) 785 call wrtout(std_out, sjoin(' Optimal value for ng0sh ', ltoa(ng0sh_opt)), "COLL") 786 Sigp%mG0=ng0sh_opt 787 788 ! G-sphere for W and Sigma_c is initialized from the SCR file. 789 call gsph_init(Gsph_c,Cryst,Er%npwe,gvec=Er%gvec) 790 call gsph_init(Gsph_x,Cryst,Sigp%npwx,gvec=gvec_kss) 791 792 ! === Make biggest G-sphere of Sigp%npwvec vectors === 793 Sigp%npwvec=MAX(Sigp%npwwfn,Sigp%npwx) 794 call gsph_init(Gsph_Max,Cryst,Sigp%npwvec,gvec=gvec_kss) 795 796 !BEGINDEBUG 797 ! Make sure that the two G-spheres are equivalent. 798 ierr=0 799 if (sigma_needs_w(Sigp)) then 800 ng = MIN(SIZE(Gsph_c%gvec,DIM=2),SIZE(gvec_kss,DIM=2)) 801 do ig1=1,ng 802 if (ANY(Gsph_c%gvec(:,ig1)/=gvec_kss(:,ig1))) then 803 ierr=ierr+1 804 write(std_out,*)" Gsph_c, gvec_kss ",ig1,"/",ng,Gsph_c%gvec(:,ig1),gvec_kss(:,ig1) 805 end if 806 end do 807 ABI_CHECK(ierr==0,"Mismatch between Gsph_c and gvec_kss") 808 end if 809 ierr=0 810 ng = MIN(SIZE(Gsph_x%gvec,DIM=2),SIZE(gvec_kss,DIM=2)) 811 do ig1=1,ng 812 if (ANY(Gsph_x%gvec(:,ig1)/=gvec_kss(:,ig1))) then 813 ierr=ierr+1 814 write(std_out,*)" Gsph_x, gvec_kss ",ig1,"/",ng,Gsph_x%gvec(:,ig1),gvec_kss(:,ig1) 815 end if 816 end do 817 ABI_CHECK(ierr==0,"Mismatch between Gsph_x and gvec_kss") 818 !ENDDEBUG 819 820 ABI_FREE(gvec_kss) 821 ! 822 ! === Get Fourier components of the Coulombian for all q-points in the IBZ === 823 ! * If required, use a cutoff in the interaction 824 ! * Pcv%vc_sqrt contains Vc^{-1/2} 825 ! * Setup also the analytical calculation of the q->0 component 826 ! FIXME recheck ngfftf since I got different charge outside the cutoff region 827 828 if (Dtset%gw_nqlwl==0) then 829 nqlwl=1 830 ABI_MALLOC(qlwl,(3,nqlwl)) 831 qlwl(:,1)= GW_Q0_DEFAULT 832 else 833 nqlwl=Dtset%gw_nqlwl 834 ABI_MALLOC(qlwl,(3,nqlwl)) 835 qlwl(:,:)=Dtset%gw_qlwl(:,1:nqlwl) 836 end if 837 838 !The Coulomb interaction used here might have two terms : 839 !the first term generates the usual sigma self-energy, but possibly, one should subtract 840 !from it the Coulomb interaction already present in the Kohn-Sham basis, 841 !if the usefock associated to ixc is one. 842 !The latter excludes (in the present implementation) mod(Dtset%gwcalctyp,10)==5 843 nvcoul_init=1 844 call get_xclevel(Dtset%ixc,xclevel_ixc,usefock_ixc) 845 if(usefock_ixc==1)then 846 nvcoul_init=2 847 if(mod(Dtset%gwcalctyp,10)==5)then 848 write(msg,'(4a,i3,a,i3,4a,i5)')ch10,& 849 & ' The starting wavefunctions were obtained from self-consistent calculations in the planewave basis set',ch10,& 850 & ' with ixc = ',Dtset%ixc,' associated with usefock =',usefock_ixc,ch10,& 851 & ' In this case, the present implementation does not allow that the self-energy for sigma corresponds to',ch10,& 852 & ' mod(gwcalctyp,10)==5, while your gwcalctyp= ',Dtset%gwcalctyp 853 MSG_ERROR(msg) 854 endif 855 endif 856 857 do ivcoul_init=1,nvcoul_init 858 rcut = Dtset%rcut 859 icutcoul_eff=Dtset%icutcoul 860 Sigp%sigma_mixing=one 861 if( mod(Dtset%gwcalctyp,10)==5 .or. ivcoul_init==2)then 862 if(abs(Dtset%hyb_mixing)>tol8)then 863 ! Warning : the absolute value is needed, because of the singular way used to define the default for this input variable 864 Sigp%sigma_mixing=abs(Dtset%hyb_mixing) 865 else if(abs(Dtset%hyb_mixing_sr)>tol8)then 866 Sigp%sigma_mixing=abs(Dtset%hyb_mixing_sr) 867 icutcoul_eff=5 868 endif 869 if(abs(rcut)<tol6 .and. abs(Dtset%hyb_range_fock)>tol8)rcut=one/Dtset%hyb_range_fock 870 endif 871 872 !#if 1 873 if(ivcoul_init==1)then 874 875 if (Gsph_x%ng > Gsph_c%ng) then 876 call vcoul_init(Vcp,Gsph_x,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,& 877 & Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,ngfftf,comm) 878 else 879 call vcoul_init(Vcp,Gsph_c,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,& 880 & Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,ngfftf,comm) 881 end if 882 883 else 884 885 ! Use a temporary Vcp_ks to compute the Coulomb interaction already present in the Fock part of the Kohn-Sham Hamiltonian 886 if (Gsph_x%ng > Gsph_c%ng) then 887 call vcoul_init(Vcp_ks,Gsph_x,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,& 888 & Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,ngfftf,comm) 889 else 890 call vcoul_init(Vcp_ks,Gsph_c,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,& 891 & Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,ngfftf,comm) 892 end if 893 894 ! Now compute the residual Coulomb interaction 895 Vcp%vc_sqrt_resid=sqrt(Vcp%vc_sqrt**2-Sigp%sigma_mixing*Vcp_ks%vc_sqrt**2) 896 Vcp%i_sz_resid=Vcp%i_sz-Sigp%sigma_mixing*Vcp_ks%i_sz 897 ! The mixing factor has already been accounted for, so set it back to one 898 Sigp%sigma_mixing=one 899 call vcoul_free(Vcp_ks) 900 901 !DEBUG 902 write(std_out,'(a)')' setup_sigma : the residual Coulomb interaction has been computed' 903 !ENDDEBUG 904 905 endif 906 !#else 907 ! call vcoul_init(Vcp,Gsph_Max,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,ivcoul_init,Dtset%vcutgeo,& 908 !& Dtset%ecutsigx,Sigp%npwx,nqlwl,qlwl,ngfftf,comm) 909 !#endif 910 911 enddo 912 913 #if 0 914 ! Using the random q for the optical limit is one of the reasons 915 ! why sigma breaks the initial energy degeneracies. 916 Vcp%i_sz=zero 917 Vcp%vc_sqrt(1,:)=czero 918 Vcp%vcqlwl_sqrt(1,:)=czero 919 #endif 920 921 ABI_FREE(qlwl) 922 923 Sigp%ecuteps = Dtset%ecuteps 924 Sigp%ecutwfn = Dtset%ecutwfn 925 Sigp%ecutsigx = Dtset%ecutsigx 926 927 ! === Setup of the FFT mesh for the oscilator strengths === 928 ! * gwc_ngfft(7:18)==Dtset%ngfft(7:18) which is initialized before entering screening. 929 ! * Here we redefine gwc_ngfft(1:6) according to the following options : 930 ! 931 ! method==0 --> FFT grid read from fft.in (debugging purpose) 932 ! method==1 --> Normal FFT mesh 933 ! method==2 --> Slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90) 934 ! method==3 --> Doubled FFT grid, same as the the FFT for the density, 935 ! 936 ! enforce_sym==1 ==> Enforce a FFT mesh compatible with all the symmetry operation and FFT library 937 ! enforce_sym==0 ==> Find the smallest FFT grid compatbile with the library, do not care about symmetries 938 ! 939 gwc_ngfft(1:18)=Dtset%ngfft(1:18) 940 gwx_ngfft(1:18)=Dtset%ngfft(1:18) 941 942 method=2 943 if (Dtset%fftgw==00 .or. Dtset%fftgw==01) method=0 944 if (Dtset%fftgw==10 .or. Dtset%fftgw==11) method=1 945 if (Dtset%fftgw==20 .or. Dtset%fftgw==21) method=2 946 if (Dtset%fftgw==30 .or. Dtset%fftgw==31) method=3 947 enforce_sym=MOD(Dtset%fftgw,10) 948 949 ! FFT mesh for sigma_x. 950 call setmesh(gmet,Gsph_Max%gvec,gwx_ngfft,Sigp%npwvec,Sigp%npwx,Sigp%npwwfn,& 951 & gwx_nfftot,method,Sigp%mG0,Cryst,enforce_sym) 952 953 ! FFT mesh for sigma_c. 954 call setmesh(gmet,Gsph_Max%gvec,gwc_ngfft,Sigp%npwvec,Er%npwe,Sigp%npwwfn,& 955 & gwc_nfftot,method,Sigp%mG0,Cryst,enforce_sym,unit=dev_null) 956 957 !call new_setmesh(Cryst,ecut_osc,ecutwfn,nkpt,kpoints,method,Sigp%mG0,enforce_sym,gwx_ngfft,gwx_nfftot) 958 !call new_setmesh(Cryst,ecut_osc,ecutwfn,nkpt,kpoints,method,Sigp%mG0,enforce_sym,gwc_ngfft,gwc_nfftot) 959 960 ! ====================================================================== 961 ! ==== Check for presence of files with core orbitals, for PAW only ==== 962 ! ====================================================================== 963 Sigp%use_sigxcore=0 964 if (Dtset%usepaw==1.and.Dtset%gw_sigxcore==1) then 965 ii = 0 966 do itypat=1,Cryst%ntypat 967 string = Psps%filpsp(itypat) 968 fcore = "CORE_"//TRIM(basename(string)) 969 ic = INDEX (TRIM(string), "/" , back=.TRUE.) 970 if (ic>0 .and. ic<LEN_TRIM(string)) then 971 ! string defines a path, prepend path to fcore 972 fcore = Psps%filpsp(itypat)(1:ic)//TRIM(fcore) 973 end if 974 if (file_exists(fcore)) then 975 ii = ii+1 976 else 977 MSG_WARNING(sjoin("HF decoupling is required but cannot find file:", fcore)) 978 end if 979 end do 980 981 Sigp%use_sigxcore=1 982 if (ii/=Cryst%ntypat) then 983 MSG_ERROR("Files with core orbitals not found") 984 end if 985 end if ! PAW+HF decoupling 986 ! 987 ! ============================== 988 ! ==== Extrapolar technique ==== 989 ! ============================== 990 Sigp%gwcomp = Dtset%gwcomp 991 Sigp%gwencomp = Dtset%gwencomp 992 993 if (Sigp%gwcomp==1) then 994 write(msg,'(6a,e11.4,a)')ch10,& 995 & 'Using the extrapolar approximation to accelerate convergence',ch10,& 996 & 'with respect to the number of bands included',ch10,& 997 & 'with gwencomp: ',Sigp%gwencomp*Ha_eV,' [eV]' 998 call wrtout(std_out,msg,'COLL') 999 end if 1000 ! 1001 ! =================================== 1002 ! ==== Final compatibility tests ==== 1003 ! =================================== 1004 if (ANY(KS_BSt%istwfk/=1)) then 1005 MSG_WARNING('istwfk/=1 is still under development') 1006 end if 1007 1008 ltest=(KS_BSt%mband==Sigp%nbnds.and.ALL(KS_BSt%nband==Sigp%nbnds)) 1009 ABI_CHECK(ltest,'BUG in definition of KS_BSt%nband') 1010 1011 ! FIXME 1012 if (Dtset%symsigma/=0 .and. Sigp%nomegasr/=0) then 1013 if (idx_spatial_inversion(Cryst) == 0) then 1014 write(msg,'(5a)')' setup_sigma : BUG :',ch10,& 1015 & 'It is not possible to use symsigma/=0 to calculate the spectral function ',ch10,& 1016 & 'when the system does not have the spatial inversion. Please use symsigma=0 ' 1017 MSG_WARNING(msg) 1018 end if 1019 end if 1020 1021 if (mod10==SIG_GW_AC) then 1022 if (Sigp%gwcalctyp/=1) MSG_ERROR("Self-consistency with AC not implemented") 1023 if (Sigp%gwcomp==1) MSG_ERROR("AC with extrapolar technique not implemented") 1024 end if 1025 1026 call gaps_free(gaps) 1027 1028 ABI_FREE(val_indeces) 1029 1030 DBG_EXIT('COLL') 1031 1032 end subroutine setup_sigma
ABINIT/sigma_bksmask [ Functions ]
NAME
sigma_bksmask
FUNCTION
Compute tables for the distribution and the storage of the wavefunctions in the SIGMA code.
INPUTS
Dtset<type(dataset_type)>=all input variables for this dataset Sigp<sigparams_t>=Parameters governing the self-energy calculation. Kmesh <kmesh_t>=Structure describing the k-point sampling. nprocs=Total number of MPI processors my_rank=Rank of this this processor.
OUTPUT
my_spins(:)=Pointer to NULL in input. In output: list of spins treated by this node. bks_mask(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will treat this state. keep_ur(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will store this state in real space. ierr=Exit status.
PARENTS
sigma
CHILDREN
SOURCE
1262 subroutine sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,nprocs,my_spins,bks_mask,keep_ur,ierr) 1263 1264 use defs_basis 1265 use defs_datatypes 1266 use defs_abitypes 1267 use m_xmpi 1268 use m_gwdefs 1269 use m_errors 1270 use m_profiling_abi 1271 1272 use m_bz_mesh, only : kmesh_t 1273 1274 !This section has been created automatically by the script Abilint (TD). 1275 !Do not modify the following lines by hand. 1276 #undef ABI_FUNC 1277 #define ABI_FUNC 'sigma_bksmask' 1278 !End of the abilint section 1279 1280 implicit none 1281 1282 !Arguments ------------------------------------ 1283 !scalars 1284 integer,intent(in) :: my_rank,nprocs 1285 integer,intent(out) :: ierr 1286 type(Dataset_type),intent(in) :: Dtset 1287 type(sigparams_t),intent(in) :: Sigp 1288 type(kmesh_t),intent(in) :: Kmesh 1289 !arrays 1290 integer,allocatable,intent(out) :: my_spins(:) 1291 logical,intent(out) :: bks_mask(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol) 1292 logical,intent(out) :: keep_ur(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol) 1293 1294 !Local variables------------------------------- 1295 !scalars 1296 integer :: my_nspins,my_maxb,my_minb,isp,spin,band,nsppol,rank_spin 1297 logical :: store_ur 1298 !arrays 1299 integer :: tmp_spins(Sigp%nsppol),nprocs_spin(Sigp%nsppol) 1300 1301 ! ************************************************************************* 1302 1303 ierr=0; nsppol=Sigp%nsppol 1304 1305 ! List of spins for each node, number of processors per each spin 1306 ! and the MPI rank in the "spin" communicator. 1307 my_nspins=nsppol 1308 ABI_MALLOC(my_spins, (nsppol)) 1309 my_spins= [(isp, isp=1,nsppol)] 1310 nprocs_spin = nprocs; rank_spin = my_rank 1311 1312 if (nsppol==2 .and. nprocs>1) then 1313 ! Distribute spins (optimal distribution if nprocs is even) 1314 nprocs_spin(1) = nprocs/2 1315 nprocs_spin(2) = nprocs - nprocs/2 1316 my_nspins=1 1317 my_spins(1)=1 1318 if (my_rank+1>nprocs/2) then 1319 ! I will treate spin=2, compute shifted rank. 1320 my_spins(1)=2 1321 rank_spin = my_rank - nprocs/2 1322 end if 1323 end if 1324 1325 store_ur = (MODULO(Dtset%gwmem,10)==1) 1326 keep_ur=.FALSE.; bks_mask=.FALSE. 1327 1328 select case (Dtset%gwpara) 1329 case (1) 1330 ! Parallelization over transitions **without** memory distributions (Except for the spin). 1331 my_minb=1; my_maxb=Sigp%nbnds 1332 if (dtset%ucrpa>0) then 1333 bks_mask(my_minb:my_maxb,:,:)=.TRUE. 1334 if (store_ur) keep_ur(my_minb:my_maxb,:,:)=.TRUE. 1335 else 1336 do isp=1,my_nspins 1337 spin = my_spins(isp) 1338 bks_mask(my_minb:my_maxb,:,spin)=.TRUE. 1339 if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE. 1340 end do 1341 end if 1342 1343 case (2) 1344 ! Distribute bands and spin (alternating planes of bands) 1345 do isp=1,my_nspins 1346 spin = my_spins(isp) 1347 do band=1,Sigp%nbnds 1348 if (xmpi_distrib_with_replicas(band,Sigp%nbnds,rank_spin,nprocs_spin(spin))) then 1349 !if (MODULO(band-1,nprocs_spin(spin))==rank_spin) then 1350 bks_mask(band,:,spin)=.TRUE. 1351 if (store_ur) keep_ur(band,:,spin)=.TRUE. 1352 end if 1353 end do 1354 end do 1355 1356 #if 0 1357 ! Each node store the full set of occupied states to speed up Sigma_x. 1358 do isp=1,my_nspins 1359 spin = my_spins(isp) 1360 do ik_ibz=1,Kmesh%nibz 1361 ks_iv=ks_vbik(ik_ibz,spin) ! Valence index for this (k,s) 1362 bks_mask(1:ks_iv,:,spin)=.TRUE. 1363 if (store_ur) keep_ur(1:ks_iv,:,spin)=.TRUE. 1364 end do 1365 end do 1366 #endif 1367 1368 case default 1369 MSG_WARNING("Wrong value for gwpara") 1370 ierr = 1 1371 end select 1372 1373 ! Return my_spins with correct size. 1374 tmp_spins(1:my_nspins) = my_spins(1:my_nspins) 1375 1376 ABI_FREE(my_spins) 1377 ABI_MALLOC(my_spins, (my_nspins)) 1378 my_spins = tmp_spins(1:my_nspins) 1379 1380 end subroutine sigma_bksmask
ABINIT/sigma_tables [ Functions ]
NAME
sigma_tables
FUNCTION
Build symmetry tables used to speedup self-consistent GW calculations.
INPUTS
Kmesh <kmesh_t>=Structure describing the k-point sampling. [Bnd_sym(Kmesh%nibz,Sigp%nsppol)] <type(Bands_Symmetries)> SiDE EFFECTS Sigp<sigparams_t>=This routine initializes the tables: %Sigcij_tab %Sigxij_tab that are used to select the matrix elements of the self-energy that have to be calculated.
PARENTS
setup_sigma,sigma
CHILDREN
SOURCE
1062 subroutine sigma_tables(Sigp,Kmesh,Bnd_sym) 1063 1064 use defs_basis 1065 use defs_datatypes 1066 use m_gwdefs 1067 use m_profiling_abi 1068 use m_errors 1069 1070 use m_bz_mesh, only : kmesh_t 1071 use m_esymm, only : esymm_t, esymm_failed 1072 1073 !This section has been created automatically by the script Abilint (TD). 1074 !Do not modify the following lines by hand. 1075 #undef ABI_FUNC 1076 #define ABI_FUNC 'sigma_tables' 1077 !End of the abilint section 1078 1079 implicit none 1080 1081 !Arguments ------------------------------------ 1082 !scalars 1083 type(sigparams_t),intent(inout) :: Sigp 1084 type(kmesh_t),intent(in) :: Kmesh 1085 !arrays 1086 type(esymm_t),optional,intent(in) :: Bnd_sym(Kmesh%nibz,Sigp%nsppol) 1087 1088 !Local variables------------------------------- 1089 !scalars 1090 integer :: gwcalctyp,spin,ikcalc,ik_ibz,bmin,bmax,bcol,brow 1091 integer :: ii,idx_x,idx_c,irr_idx1,irr_idx2 1092 !arrays 1093 integer,allocatable :: sigc_bidx(:),sigx_bidx(:) 1094 logical :: use_sym_at(Kmesh%nibz,Sigp%nsppol) 1095 1096 ! ************************************************************************* 1097 1098 gwcalctyp=Sigp%gwcalctyp 1099 1100 ! Recreate the Sig_ij tables taking advantage of the classification of the bands. 1101 if (allocated(Sigp%Sigxij_tab)) then 1102 call sigijtab_free(Sigp%Sigxij_tab) 1103 ABI_DT_FREE(Sigp%Sigxij_tab) 1104 end if 1105 if (allocated(Sigp%Sigcij_tab)) then 1106 call sigijtab_free(Sigp%Sigcij_tab) 1107 ABI_DT_FREE(Sigp%Sigcij_tab) 1108 end if 1109 1110 ABI_DT_MALLOC(Sigp%Sigcij_tab,(Sigp%nkptgw,Sigp%nsppol)) 1111 ABI_DT_MALLOC(Sigp%Sigxij_tab,(Sigp%nkptgw,Sigp%nsppol)) 1112 1113 use_sym_at=.FALSE. 1114 if (PRESENT(Bnd_sym)) then 1115 do spin=1,Sigp%nsppol 1116 do ikcalc=1,Sigp%nkptgw 1117 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 1118 use_sym_at(ik_ibz,spin) = ( .not.esymm_failed(Bnd_sym(ik_ibz,spin)) ) 1119 end do 1120 end do 1121 end if 1122 1123 do spin=1,Sigp%nsppol 1124 do ikcalc=1,Sigp%nkptgw 1125 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 1126 1127 if (use_sym_at(ik_ibz,spin)) then 1128 if (gwcalctyp<20) then 1129 MSG_ERROR("You should not be here!") 1130 end if 1131 1132 bmin=Sigp%minbnd(ikcalc,spin); bmax=Sigp%maxbnd(ikcalc,spin) 1133 ABI_DT_MALLOC(Sigp%Sigxij_tab(ikcalc,spin)%col,(bmin:bmax)) 1134 ABI_DT_MALLOC(Sigp%Sigcij_tab(ikcalc,spin)%col,(bmin:bmax)) 1135 1136 do bcol=bmin,bmax 1137 ABI_MALLOC(sigc_bidx,(bmax-bmin+1)) 1138 ABI_MALLOC(sigx_bidx,(bmax-bmin+1)) 1139 1140 if (Bnd_sym(ik_ibz,spin)%err_status/=0) then ! Band classification failed. 1141 sigc_bidx = (/(ii,ii=bmin,bmax)/) 1142 idx_c = bmax-bmin+1 1143 sigx_bidx = (/(ii,ii=bmin,bcol)/) ! Hermitian 1144 idx_x = bcol-bmin+1 1145 else 1146 irr_idx2 = Bnd_sym(ik_ibz,spin)%b2irrep(bcol) 1147 idx_c = 0 1148 do brow=bmin,bmax 1149 irr_idx1 = Bnd_sym(ik_ibz,spin)%b2irrep(brow) 1150 if (sigma_is_herm(Sigp).and.bcol<brow) CYCLE ! Only the upper triangle for HF, SEX, or COHSEX. 1151 if (irr_idx1 == irr_idx2) then ! same character, add this row to the list. 1152 idx_c = idx_c +1 1153 sigc_bidx(idx_c) = brow 1154 end if 1155 end do 1156 idx_x = 0 1157 do brow=bmin,bcol 1158 irr_idx1 = Bnd_sym(ik_ibz,spin)%b2irrep(brow) 1159 if (bcol<brow) CYCLE ! Sig_x is always Hermitian. 1160 if (irr_idx1 == irr_idx2) then ! same character, add this row to the list. 1161 idx_x = idx_x +1 1162 sigx_bidx(idx_x) = brow 1163 end if 1164 end do 1165 end if 1166 ! 1167 ! Table for Sigma_x matrix elements taking into account symmetries of the bands. 1168 ABI_MALLOC(Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx,(idx_x)) 1169 1170 Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%size1= idx_x 1171 Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx(:) = sigx_bidx(1:idx_x) 1172 !write(std_out,*)" Sigxij_tab: ikcalc, spin, bcol ",ikcalc,spin,bcol 1173 !write(std_out,*)" size: ",idx_x,(Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx(ii),ii=1,idx_x) 1174 ! 1175 ! Table for Sigma_c matrix elements taking into account symmetries of the bands. 1176 ABI_MALLOC(Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx,(idx_c)) 1177 1178 Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%size1= idx_c 1179 Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx(:) = sigc_bidx(1:idx_c) 1180 !write(std_out,*)" Sigcij_tab: ikcalc, spin, bcol ",ikcalc,spin,bcol 1181 !write(std_out,*)" size: ",idx_c,(Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx(ii), ii=1,idx_c) 1182 1183 ABI_FREE(sigx_bidx) 1184 ABI_FREE(sigc_bidx) 1185 end do ! bcol 1186 1187 else ! Symmetries cannot be used for this (k,s). 1188 1189 bmin=Sigp%minbnd(ikcalc,spin); bmax=Sigp%maxbnd(ikcalc,spin) 1190 ABI_DT_MALLOC(Sigp%Sigcij_tab(ikcalc,spin)%col,(bmin:bmax)) 1191 ABI_DT_MALLOC(Sigp%Sigxij_tab(ikcalc,spin)%col,(bmin:bmax)) 1192 1193 if (gwcalctyp<20) then ! QP wavefunctions == KS, therefore only diagonal elements are calculated. 1194 do bcol=bmin,bmax 1195 ABI_MALLOC(Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx,(1:1)) 1196 Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%size1= 1 1197 Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx(1) = bcol 1198 ABI_MALLOC(Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx,(1:1)) 1199 Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%size1= 1 1200 Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx(1) = bcol 1201 end do 1202 else 1203 ! Use QP wavefunctions, Sigma_ij matrix is sparse but we have to classify the states in sigma. 1204 ! The only thing we can do here is filling the entire matrix taking advantage of Hermiticity (if any). 1205 do bcol=bmin,bmax 1206 ABI_MALLOC(Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx,(bcol-bmin+1)) 1207 Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%size1= bcol-bmin+1 1208 Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx(:) = (/(ii,ii=bmin,bcol)/) ! Sigma_x is Hermitian. 1209 !write(std_out,*)"Sigxij_tab: ikcalc, spin, bcol ",ikcalc,spin,bcol,Sigp%Sigxij_tab(ikcalc,spin)%col(bcol)%bidx(:) 1210 1211 ABI_MALLOC(sigc_bidx,(bmax-bmin+1)) 1212 idx_c = 0 1213 do brow=bmin,bmax 1214 if (sigma_is_herm(Sigp).and.bcol<brow) CYCLE ! Only the upper triangle of Sigc_ij is needed (SEX, COHSEX). 1215 idx_c = idx_c +1 1216 sigc_bidx(idx_c) = brow 1217 end do 1218 ABI_MALLOC(Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx,(idx_c)) 1219 Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%size1= idx_c 1220 Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx(:) = sigc_bidx(1:idx_c) 1221 ABI_FREE(sigc_bidx) 1222 !write(std_out,*)"Sigcij_tab: ikcalc, spin, bcol ",ikcalc,spin,bcol,Sigp%Sigcij_tab(ikcalc,spin)%col(bcol)%bidx(:) 1223 end do 1224 end if 1225 end if 1226 1227 end do !ikcalc 1228 end do !spin 1229 1230 end subroutine sigma_tables