TABLE OF CONTENTS
ABINIT/sigma [ Functions ]
NAME
sigma
FUNCTION
Calculate the matrix elements of the self-energy operator.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (GMR, VO, LR, RWG, MT, MG, RShaltaf) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
acell(3)=length scales of primitive translations (bohr) codvsn=code version Dtfil<type(datafiles_type)>=variables related to files Dtset<type(dataset_type)>=all input variables for this dataset Pawang<type(pawang_type)>=paw angular mesh and related data Pawrad(ntypat*usepaw)<type(pawrad_type)>=paw radial mesh and related data Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data Psps<type(pseudopotential_type)>=variables related to pseudopotentials Before entering the first time in sigma, a significant part of Psps has been initialized : the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp. All the remaining components of Psps are to be initialized in the call to pspini. The next time the code enters screening, Psps might be identical to the one of the previous Dtset, in which case, no reinitialisation is scheduled in pspini.F90. rprim(3,3)=dimensionless real space primitive translations
OUTPUT
Converged=.TRUE. if degw are within the user-specified tolerance. Output is written on the main abinit output file. Some results are stored in external files
PARENTS
driver
NOTES
ON THE USE OF FFT GRIDS:
=================
In case of PAW:
---------------
Two FFT grids are used:
- A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
It is defined by nfft, ngfft, mgfft, ...
Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
- A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid.
In case of norm-conserving:
---------------------------
- Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
CHILDREN
calc_sigc_me,calc_sigx_me,calc_ucrpa,calc_vhxc_me,chkpawovlp classify_bands,cohsex_me,crystal_free,denfgr,destroy_mpi_enreg ebands_copy,ebands_free,ebands_interpolate_kpath,ebands_report_gap ebands_update_occ,em1results_free,energies_init,esymm_free fftdatar_write,fourdp,get_gftt,getph,gsph_free,hdr_free init_distribfft_seq,initmpi_seq,kmesh_free,kxc_ada,kxc_driver littlegroup_free,littlegroup_init,melements_free,melements_print melements_zero,melflags_reset,metric,mkdump_er,mkrdim,nhatgrid paw_an_free,paw_an_init,paw_an_nullify,paw_check_symcprj,paw_dijhf paw_gencond,paw_ij_free,paw_ij_init,paw_ij_nullify,paw_ij_print paw_mkdijexc_core,paw_pwaves_lmn_free,paw_pwaves_lmn_init,paw_qpscgw pawcprj_alloc,pawcprj_free,pawdenpot,pawdij,pawfgr_destroy,pawfgr_init pawfgrtab_free,pawfgrtab_init,pawfgrtab_print,pawinit,pawmknhat,pawprt pawpuxinit,pawpwff_free,pawpwff_init,pawrhoij_alloc,pawrhoij_copy pawrhoij_free,pawtab_get_lsize,pawtab_print,ppm_free,ppm_init prep_calc_ucrpa,print_ngfft,prtrhomxmn,pspini,rdgw,rdqps,read_rhor setsymrhoij,setup_ppmodel,setup_sigma,setvtr,show_qp,sigma_bksmask sigma_free,sigma_init,sigma_tables,sigparams_free,solve_dyson,symdij symdij_all,test_charge,timab,updt_m_lda_to_qp,vcoul_free wfd_change_ngfft,wfd_copy,wfd_distribute_bands,wfd_free,wfd_get_cprj wfd_init,wfd_mkrho,wfd_print,wfd_read_wfk,wfd_reset_ur_cprj,wfd_rotate wfd_test_ortho,write_sigma_header,write_sigma_results,wrqps,wrtout xmpi_barrier,xmpi_bcast,xmpi_sum
SOURCE
83 #if defined HAVE_CONFIG_H 84 #include "config.h" 85 #endif 86 87 #include "abi_common.h" 88 89 subroutine sigma(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim,converged) 90 91 use defs_basis 92 use m_gwdefs 93 use defs_datatypes 94 use defs_abitypes 95 use defs_wvltypes 96 use m_xmpi 97 use m_xomp 98 use m_errors 99 use m_profiling_abi 100 use m_ab7_mixing 101 use m_nctk 102 use m_kxc 103 #ifdef HAVE_NETCDF 104 use netcdf 105 #endif 106 use m_hdr 107 use libxc_functionals 108 use m_wfd 109 110 use m_numeric_tools, only : imax_loc 111 use m_fstrings, only : strcat, sjoin, itoa 112 use m_blas, only : xdotc 113 use m_io_tools, only : open_file, file_exists, iomode_from_fname 114 use m_mpinfo, only : destroy_mpi_enreg 115 use m_geometry, only : normv, mkrdim, metric 116 use m_fftcore, only : print_ngfft 117 use m_fft_mesh, only : get_gftt 118 use m_ioarr, only : fftdatar_write, read_rhor 119 use m_crystal, only : crystal_free, crystal_t 120 use m_crystal_io, only : crystal_ncwrite 121 use m_ebands, only : ebands_update_occ, ebands_copy, ebands_report_gap, get_valence_idx, get_bandenergy, & 122 & ebands_free, ebands_init, ebands_ncwrite, ebands_interpolate_kpath 123 use m_energies, only : energies_type, energies_init 124 use m_bz_mesh, only : kmesh_t, kmesh_free, littlegroup_t, littlegroup_init, littlegroup_free 125 use m_gsphere, only : gsphere_t, gsph_free 126 use m_kg, only : getph 127 use m_vcoul, only : vcoul_t, vcoul_free 128 use m_qparticles, only : wrqps, rdqps, rdgw, show_QP, updt_m_lda_to_qp 129 use m_screening, only : mkdump_er, em1results_free, epsilonm1_results 130 use m_ppmodel, only : ppm_init, ppm_free, setup_ppmodel, getem1_from_PPm, ppmodel_t 131 use m_sigma, only : sigma_init, sigma_free, sigma_ncwrite, sigma_t, sigma_get_exene, & 132 & write_sigma_header, write_sigma_results 133 use m_dyson_solver, only : solve_dyson 134 use m_esymm, only : esymm_t, esymm_free, esymm_failed 135 use m_melemts, only : melflags_reset, melements_print, melements_free, melflags_t, melements_t, melements_zero 136 use m_pawang, only : pawang_type 137 use m_pawrad, only : pawrad_type 138 use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize 139 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 140 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_print 141 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print 142 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, pawrhoij_get_nspden, symrhoij 143 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, paw_overlap 144 use m_pawdij, only : pawdij, symdij_all 145 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 146 use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free 147 use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free 148 use m_paw_slater, only : paw_mkdijexc_core, paw_dijhf 149 use m_paw_dmft, only : paw_dmft_type 150 151 !This section has been created automatically by the script Abilint (TD). 152 !Do not modify the following lines by hand. 153 #undef ABI_FUNC 154 #define ABI_FUNC 'sigma' 155 use interfaces_14_hidewrite 156 use interfaces_18_timing 157 use interfaces_51_manage_mpi 158 use interfaces_53_ffts 159 use interfaces_64_psp 160 use interfaces_65_paw 161 use interfaces_67_common 162 use interfaces_69_wfdesc 163 use interfaces_70_gw 164 !End of the abilint section 165 166 implicit none 167 168 !Arguments ------------------------------------ 169 !scalars 170 logical,intent(out) :: converged 171 character(len=6),intent(in) :: codvsn 172 type(Datafiles_type),intent(in) :: Dtfil 173 type(Dataset_type),intent(inout) :: Dtset 174 type(Pawang_type),intent(inout) :: Pawang 175 type(Pseudopotential_type),intent(inout) :: Psps 176 !arrays 177 real(dp),intent(in) :: acell(3),rprim(3,3) 178 type(Pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Psps%usepaw) 179 type(Pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Psps%usepaw) 180 181 !Local variables------------------------------- 182 !scalars 183 integer,parameter :: level40=40,tim_fourdp5=5,master=0,cplex1=1 184 integer :: approx_type,b1gw,b2gw,choice,cplex,cplex_dij,band 185 integer :: dim_kxcg,gwcalctyp,gnt_option,has_dijU,has_dijso,iab,bmin,bmax,irr_idx1,irr_idx2 186 integer :: iat,ib,ib1,ib2,ic,id_required,ider,idir,ii,ik,ierr,ount 187 integer :: ik_bz,ikcalc,ik_ibz,ikxc,ipert,npw_k,omp_ncpus 188 integer :: isp,is_idx,istep,itypat,itypatcor,izero,jj,first_band,last_band 189 integer :: ks_iv,lcor,lmn2_size_max,mband,my_nband 190 integer :: mgfftf,mod10,moved_atm_inside,moved_rhor,n3xccc !,mgfft 191 integer :: nbsc,ndij,ndim,nfftf,nfftf_tot,nkcalc,gwc_nfft,gwc_nfftot,gwx_nfft,gwx_nfftot 192 integer :: ngrvdw,nhatgrdim,nkxc,nkxc1,nprocs,nscf,nspden_rhoij,nzlmopt,optene 193 integer :: optcut,optgr0,optgr1,optgr2,option,option_test,option_dij,optrad,optrhoij,psp_gencond 194 integer :: my_rank,rhoxsp_method,comm,use_aerhor,use_umklp,usexcnhat 195 integer :: ioe0j,spin,io,jb,nomega_sigc 196 integer :: temp_unt,ncid 197 integer :: work_size,nstates_per_proc,my_nbks 198 !integer :: jb_qp,ib_ks,ks_irr 199 real(dp) :: compch_fft,compch_sph,r_s,rhoav,alpha,opt_ecut 200 real(dp) :: drude_plsmf,my_plsmf,ecore,ecut_eff,ecutdg_eff,ehartree 201 real(dp) :: ex_energy,gsqcutc_eff,gsqcutf_eff,gsqcut_shp,norm,oldefermi 202 real(dp) :: ucvol,vxcavg,vxcavg_qp 203 real(dp) :: gwc_gsq,gwx_gsq,gw_gsq 204 real(dp):: eff,mempercpu_mb,max_wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem 205 complex(dpc) :: max_degw,cdummy 206 logical :: use_paw_aeur,dbg_mode,pole_screening,call_pawinit 207 character(len=500) :: msg 208 character(len=fnlen) :: wfk_fname,pawden_fname 209 type(kmesh_t) :: Kmesh,Qmesh 210 type(ebands_t) :: KS_BSt,QP_BSt 211 type(vcoul_t) :: Vcp 212 type(crystal_t) :: Cryst 213 type(Energies_type) :: KS_energies,QP_energies 214 type(Epsilonm1_results) :: Er 215 type(gsphere_t) :: Gsph_Max,Gsph_x,Gsph_c 216 type(hdr_type) :: Hdr_wfk,Hdr_sigma,hdr_rhor 217 type(melflags_t) :: KS_mflags,QP_mflags 218 type(melements_t) :: KS_me,QP_me 219 type(MPI_type) :: MPI_enreg_seq 220 type(paw_dmft_type) :: Paw_dmft 221 type(pawfgr_type) :: Pawfgr 222 type(ppmodel_t) :: PPm 223 type(sigparams_t) :: Sigp 224 type(sigma_t) :: Sr 225 type(wfd_t),target :: Wfd,Wfdf 226 type(wvl_data) :: Wvl 227 !arrays 228 integer :: gwc_ngfft(18),ngfftc(18),ngfftf(18),gwx_ngfft(18) 229 integer,allocatable :: nq_spl(:),nlmn_atm(:),my_spins(:) 230 integer,allocatable :: tmp_gfft(:,:),ks_vbik(:,:),nband(:,:),l_size_atm(:),qp_vbik(:,:) 231 integer,allocatable :: tmp_kstab(:,:,:),ks_irreptab(:,:,:),qp_irreptab(:,:,:),my_band_list(:) 232 real(dp),parameter :: k0(3)=zero 233 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),strsxc(6),tsec(2) 234 real(dp),allocatable :: grchempottn(:,:),grewtn(:,:),grvdw(:,:),qmax(:) 235 real(dp),allocatable :: ks_nhat(:,:),ks_nhatgr(:,:,:),ks_rhog(:,:) 236 real(dp),allocatable :: ks_rhor(:,:),ks_vhartr(:),ks_vtrial(:,:),ks_vxc(:,:) 237 real(dp),allocatable :: ks_taug(:,:),ks_taur(:,:) 238 real(dp),allocatable :: kxc(:,:),qp_kxc(:,:),ph1d(:,:),ph1df(:,:) 239 real(dp),allocatable :: prev_rhor(:,:),prev_taur(:,:),qp_nhat(:,:) 240 real(dp),allocatable :: qp_nhatgr(:,:,:),qp_rhog(:,:),qp_rhor_paw(:,:) 241 real(dp),allocatable :: qp_rhor_n_one(:,:),qp_rhor_nt_one(:,:) 242 real(dp),allocatable :: qp_rhor(:,:),qp_vhartr(:),qp_vtrial(:,:),qp_vxc(:,:) 243 real(dp),allocatable :: qp_taur(:,:),qp_taug(:,:),igwene(:,:,:) 244 real(dp),allocatable :: vpsp(:),xccc3d(:),dijexc_core(:,:,:),dij_hf(:,:,:) 245 !real(dp),allocatable :: osoc_bks(:, :, :) 246 real(dp),allocatable :: ks_aepaw_rhor(:,:) !,ks_n_one_rhor(:,:),ks_nt_one_rhor(:,:) 247 complex(dpc) :: ovlp(2) 248 complex(dpc),allocatable :: ctmp(:,:),hbare(:,:,:,:) 249 complex(dpc),target,allocatable :: sigcme(:,:,:,:,:) 250 complex(dpc),allocatable :: hlda(:,:,:,:),htmp(:,:,:,:),uks2qp(:,:) 251 complex(gwpc),allocatable :: kxcg(:,:),fxc_ADA(:,:,:) 252 complex(gwpc),ABI_CONTIGUOUS pointer :: ug1(:) 253 !complex(dpc),pointer :: sigcme_p(:,:,:,:) 254 complex(dpc),allocatable :: sigcme_k(:,:,:,:) 255 complex(dpc), allocatable :: rhot1_q_m(:,:,:,:,:,:,:) 256 complex(dpc), allocatable :: M1_q_m(:,:,:,:,:,:,:) 257 logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:),bmask(:) 258 type(esymm_t),target,allocatable :: KS_sym(:,:) 259 type(esymm_t),pointer :: QP_sym(:,:) 260 type(pawcprj_type),allocatable :: Cp1(:,:),Cp2(:,:) 261 type(littlegroup_t),allocatable :: Ltg_k(:) 262 type(Paw_an_type),allocatable :: KS_paw_an(:),QP_paw_an(:) 263 type(Paw_ij_type),allocatable :: KS_paw_ij(:),QP_paw_ij(:) 264 type(Pawfgrtab_type),allocatable :: Pawfgrtab(:) 265 type(Pawrhoij_type),allocatable :: KS_Pawrhoij(:),QP_pawrhoij(:),prev_Pawrhoij(:),tmp_pawrhoij(:) 266 type(pawpwff_t),allocatable :: Paw_pwff(:) 267 type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:) 268 269 !************************************************************************ 270 271 DBG_ENTER('COLL') 272 273 call timab(401,1,tsec) ! sigma(Total) 274 call timab(402,1,tsec) ! sigma(Init1) 275 276 write(msg,'(7a)')& 277 & ' SIGMA: Calculation of the GW corrections ',ch10,ch10,& 278 & ' Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining.',ch10,& 279 & ' Incorporated in ABINIT by V. Olevano, G.-M. Rignanese, and M. Torrent.' 280 call wrtout(std_out,msg,'COLL') 281 call wrtout(ab_out,msg,'COLL') 282 283 if(dtset%ucrpa>0) then 284 write(msg,'(6a)')ch10,& 285 & ' cRPA Calculation: Calculation of the screened Coulomb interaction (ucrpa/=0) ',ch10 286 call wrtout(ab_out, msg,'COLL') 287 call wrtout(std_out,msg,'COLL') 288 end if 289 290 #if defined HAVE_GW_DPC 291 if (gwpc/=8) then 292 write(msg,'(6a)')ch10,& 293 & 'Number of bytes for double precision complex /=8 ',ch10,& 294 & 'Cannot continue due to kind mismatch in BLAS library ',ch10,& 295 & 'Some BLAS interfaces are not generated by abilint ' 296 MSG_ERROR(msg) 297 end if 298 write(msg,'(a,i2,a)')'.Using double precision arithmetic ; gwpc = ',gwpc,ch10 299 #else 300 write(msg,'(a,i2,a)')'.Using single precision arithmetic ; gwpc = ',gwpc,ch10 301 #endif 302 call wrtout(std_out,msg,'COLL') 303 call wrtout(ab_out,msg,'COLL') 304 305 gwcalctyp=Dtset%gwcalctyp 306 mod10 =MOD(Dtset%gwcalctyp,10) 307 308 ! Perform some additional checks for hybrid functional calculations 309 if(mod10==5) then 310 if (Dtset%ixc_sigma<0 .and. .not.libxc_functionals_check()) then 311 msg='Hybrid functional calculations require the compilation with LIBXC library' 312 MSG_ERROR(msg) 313 end if 314 ! XG 20171116 : I do not agree with this condition, as one might like to do a one-shot hybrid functional calculation 315 ! on top of a LDA/GGA calculation ... give the power (and risks) to the user ! 316 ! if(gwcalctyp<10) then 317 ! msg='gwcalctyp should enforce update of the energies and/or wavefunctions when performing hybrid functional calculation' 318 ! MSG_ERROR(msg) 319 ! end if 320 if(Dtset%usepaw==1) then 321 msg='PAW version of hybrid functional calculations is not implemented' 322 MSG_ERROR(msg) 323 end if 324 end if 325 326 !=== Initialize MPI variables, and parallelization level === 327 !* gwpara 0--> sequential run, 1--> parallelism over k-points, 2--> parallelism over bands. 328 !* In case of gwpara==1 memory is not parallelized. 329 !* If gwpara==2, bands are divided among processors but each proc has all the states where GW corrections are required. 330 comm = xmpi_world; my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 331 converged = .FALSE. 332 333 if (my_rank == master) then 334 wfk_fname = dtfil%fnamewffk 335 if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then 336 MSG_ERROR(msg) 337 end if 338 end if 339 call xmpi_bcast(wfk_fname, master, comm, ierr) 340 341 ! === Some variables need to be initialized/nullify at start === 342 call energies_init(KS_energies) 343 usexcnhat=0 344 call mkrdim(acell,rprim,rprimd) 345 call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol) 346 ! 347 ! === Define FFT grid(s) sizes === 348 ! * Be careful! This mesh is only used for densities, potentials and the matrix elements of v_Hxc. It is NOT the 349 ! (usually coarser) GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90. 350 ! See also NOTES in the comments at the beginning of this file. 351 ! NOTE: This mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn. 352 353 call pawfgr_init(Pawfgr,Dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,& 354 & gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=gmet,k0=k0) 355 356 ! Fake MPI_type for the sequential part. 357 call initmpi_seq(MPI_enreg_seq) 358 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all') 359 call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfftf(2),ngfftf(3),'all') 360 361 call print_ngfft(ngfftf,header='Dense FFT mesh used for densities and potentials') 362 nfftf_tot=PRODUCT(ngfftf(1:3)) 363 364 ! Open and read pseudopotential files === 365 call pspini(Dtset,Dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,Pawrad,Pawtab,Psps,rprimd,comm_mpi=comm) 366 367 call timab(402,2,tsec) ! Init1 368 ! 369 ! =============================================== 370 ! ==== Initialize Sigp, Er and basic objects ==== 371 ! =============================================== 372 ! * Sigp is completetly initialized here. 373 ! * Er is only initialized with dimensions, (SCR|SUSC) file is read in mkdump_Er 374 call timab(403,1,tsec) ! setup_sigma 375 376 call setup_sigma(codvsn,wfk_fname,acell,rprim,ngfftf,Dtset,Dtfil,Psps,Pawtab,& 377 & gwx_ngfft,gwc_ngfft,Hdr_wfk,Hdr_sigma,Cryst,Kmesh,Qmesh,KS_BSt,Gsph_Max,Gsph_x,Gsph_c,Vcp,Er,Sigp,comm) 378 379 call timab(403,2,tsec) ! setup_sigma 380 call timab(402,1,tsec) ! Init1 381 382 !XG090617 Please, do not remove this write, unless you have checked 383 !that the code executes correctly on max+g95 (especially, Tv5#70). 384 !It is one more a silly write, perhaps needed because the compiler does not treat correctly non-nullified pointers. 385 if (sigma_needs_w(Sigp) .and. my_rank==master) then 386 write(std_out,*)' screening after setup_sigma : Er%Hscr%headform=',Er%Hscr%headform 387 end if 388 !END XG090617 389 390 pole_screening = .FALSE. 391 if (Er%fform==2002) then 392 pole_screening = .TRUE. 393 MSG_WARNING(' EXPERIMENTAL - Using a pole-fit screening!') 394 end if 395 396 call print_ngfft(gwc_ngfft,header='FFT mesh for oscillator strengths used for Sigma_c') 397 call print_ngfft(gwx_ngfft,header='FFT mesh for oscillator strengths used for Sigma_x') 398 399 b1gw=Sigp%minbdgw 400 b2gw=Sigp%maxbdgw 401 402 gwc_nfftot=PRODUCT(gwc_ngfft(1:3)) 403 gwc_nfft =gwc_nfftot !no FFT // 404 405 gwx_nfftot=PRODUCT(gwx_ngfft(1:3)) 406 gwx_nfft =gwx_nfftot !no FFT // 407 ! 408 !TRYING TO RECREATE AN "ABINIT ENVIRONMENT" 409 KS_energies%e_corepsp=ecore/Cryst%ucvol 410 411 !=== Calculate KS occupation numbers and ks_vbk(nkibz,nsppol) ==== 412 !* ks_vbk gives the (valence|last Fermi band) index for each k and spin. 413 !* spinmagntarget is passed to fermi.F90 to fix the problem with newocc in case of magnetic metals 414 ABI_MALLOC(ks_vbik,(KS_BSt%nkpt,KS_BSt%nsppol)) 415 ABI_MALLOC(qp_vbik,(KS_BSt%nkpt,KS_BSt%nsppol)) 416 417 !call ebands_update_occ(KS_BSt,Dtset%spinmagntarget,prtvol=0) 418 ks_vbik(:,:) = get_valence_idx(KS_BSt) 419 420 ! ============================ 421 ! ==== PAW initialization ==== 422 ! ============================ 423 if (Dtset%usepaw==1) then 424 call chkpawovlp(Cryst%natom,Cryst%ntypat,Dtset%pawovlp,Pawtab,Cryst%rmet,Cryst%typat,Cryst%xred) 425 426 cplex_dij=Dtset%nspinor; cplex=1; ndij=1 427 428 ABI_DT_MALLOC(KS_Pawrhoij,(Cryst%natom)) 429 nspden_rhoij=pawrhoij_get_nspden(Dtset%nspden,Dtset%nspinor,Dtset%pawspnorb) 430 call pawrhoij_alloc(KS_Pawrhoij,Dtset%pawcpxocc,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab) 431 432 ! Initialize values for several basic arrays 433 gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2 434 435 ! Test if we have to call pawinit 436 call paw_gencond(Dtset,gnt_option,"test",call_pawinit) 437 438 if (psp_gencond==1.or. call_pawinit) then 439 call timab(553,1,tsec) 440 gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2 441 call pawinit(gnt_option,gsqcut_shp,zero,Dtset%pawlcutd,Dtset%pawlmix,& 442 & Psps%mpsang,Dtset%pawnphi,Cryst%nsym,Dtset%pawntheta,Pawang,Pawrad,& 443 & Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%xclevel,Dtset%usepotzero) 444 call timab(553,2,tsec) 445 446 ! Update internal values 447 call paw_gencond(Dtset,gnt_option,"save",call_pawinit) 448 449 else 450 if (Pawtab(1)%has_kij ==1) Pawtab(1:Cryst%ntypat)%has_kij =2 451 if (Pawtab(1)%has_nabla==1) Pawtab(1:Cryst%ntypat)%has_nabla=2 452 end if 453 Psps%n1xccc=MAXVAL(Pawtab(1:Cryst%ntypat)%usetcore) 454 455 ! Initialize optional flags in Pawtab to zero 456 ! (Cannot be done in Pawinit since the routine is called only if some pars. are changed) 457 Pawtab(:)%has_nabla = 0 458 Pawtab(:)%usepawu = 0 459 Pawtab(:)%useexexch = 0 460 Pawtab(:)%exchmix =zero 461 462 call setsymrhoij(gprimd,Pawang%l_max-1,Cryst%nsym,Dtset%pawprtvol,Cryst%rprimd,Cryst%symrec,Pawang%zarot) 463 464 ! Initialize and compute data for LDA+U 465 Paw_dmft%use_dmft=Dtset%usedmft 466 if (Dtset%usepawu>0.or.Dtset%useexexch>0) then 467 call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,& 468 & Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,Cryst%ntypat,Pawang,Dtset%pawprtvol,& 469 & Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa) 470 end if 471 472 if (my_rank == master) call pawtab_print(Pawtab) 473 474 ! Get Pawrhoij from the header of the WFK file. 475 call pawrhoij_copy(Hdr_wfk%pawrhoij,KS_Pawrhoij) 476 477 ! Re-symmetrize symrhoij. 478 ! FIXME this call leads to a SIGFAULT, likely some pointer is not initialized correctly 479 choice=1; optrhoij=1; ipert=0; idir=0 480 ! call symrhoij(KS_Pawrhoij,KS_Pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert,& 481 ! & Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij,Pawang,Dtset%pawprtvol,& 482 ! & Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat) 483 484 ! Evaluate form factor of radial part of phi.phj-tphi.tphj. 485 rhoxsp_method=1 ! Arnaud-Alouani 486 ! rhoxsp_method=2 ! Shiskin-Kresse 487 if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc 488 489 ! The q-grid must contain the FFT mesh used for sigma_c and the G-sphere for the exchange part. 490 ! We use the FFT mesh for sigma_c since COHSEX and the extrapolar method require oscillator 491 ! strengths on the FFT mesh. 492 ABI_MALLOC(tmp_gfft,(3,gwc_nfftot)) 493 call get_gftt(gwc_ngfft,k0,gmet,gwc_gsq,tmp_gfft) 494 ABI_FREE(tmp_gfft) 495 496 gwx_gsq = Dtset%ecutsigx/(two*pi**2) 497 ! allocate(tmp_gfft(3,gwx_nfftot)); q0=zero 498 ! call get_gftt(gwx_ngfft,q0,gmet,gwx_gsq,tmp_gfft) 499 ! deallocate(tmp_gfft) 500 gw_gsq = MAX(gwx_gsq,gwc_gsq) 501 502 ! * Set up q-grid, make qmax 20% larger than largest expected. 503 ABI_MALLOC(nq_spl,(Psps%ntypat)) 504 ABI_MALLOC(qmax,(Psps%ntypat)) 505 qmax = SQRT(gw_gsq)*1.2d0 ! qmax=Psps%qgrid_ff(Psps%mqgrid_ff) 506 nq_spl = Psps%mqgrid_ff 507 ! write(std_out,*)"using nq_spl",nq_spl,"qmax=",qmax 508 ABI_DT_MALLOC(Paw_pwff,(Psps%ntypat)) 509 510 call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps) 511 512 ABI_FREE(nq_spl) 513 ABI_FREE(qmax) 514 ! 515 ! Variables/arrays related to the fine FFT grid === 516 ABI_MALLOC(ks_nhat,(nfftf,Dtset%nspden)) 517 ks_nhat=zero 518 ABI_DT_MALLOC(Pawfgrtab,(Cryst%natom)) 519 call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat) 520 cplex=1 521 call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Dtset%nspden,Dtset%typat) 522 ABI_FREE(l_size_atm) 523 compch_fft=greatest_real 524 usexcnhat=MAXVAL(Pawtab(:)%usexcnhat) 525 ! * 0 if Vloc in atomic data is Vbare (Blochl s formulation) 526 ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation) 527 write(msg,'(a,i2)')' sigma : using usexcnhat = ',usexcnhat 528 call wrtout(std_out,msg,'COLL') 529 ! 530 ! Identify parts of the rectangular grid where the density has to be calculated === 531 optcut=0; optgr0=Dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-Dtset%pawstgylm 532 if (Dtset%pawcross==1) optrad=1 533 if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm 534 535 call nhatgrid(Cryst%atindx1,gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,& 536 & optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 537 538 call pawfgrtab_print(Pawfgrtab,Cryst%natom,unit=std_out,prtvol=Dtset%pawprtvol) 539 else 540 ABI_DT_MALLOC(Paw_pwff,(0)) 541 ABI_DT_MALLOC(Pawfgrtab,(0)) 542 end if !End of PAW Initialization 543 544 ! Consistency check and additional stuff done only for GW with PAW. 545 ABI_DT_MALLOC(Paw_onsite,(0)) 546 547 if (Dtset%usepaw==1) then 548 if (Dtset%ecutwfn < Dtset%ecut) then 549 write(msg,"(5a)")& 550 & "WARNING - ",ch10,& 551 & " It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,& 552 & " an excessive truncation of the planewave basis set can lead to unphysical results." 553 MSG_WARNING(msg) 554 end if 555 556 ABI_CHECK(Dtset%useexexch==0,"LEXX not yet implemented in GW") 557 ABI_CHECK(Paw_dmft%use_dmft==0,"DMFT + GW not available") 558 559 ! Optionally read core orbitals from file and calculate $ \<\phi_i|Sigma_x^\core|\phi_j\> $ for the HF decoupling. 560 if (Sigp%use_sigxcore==1) then 561 lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size) 562 ABI_MALLOC(dijexc_core,(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat)) 563 564 call paw_mkdijexc_core(ndij,cplex_dij,lmn2_size_max,Cryst,Pawtab,Pawrad,dijexc_core,Dtset%prtvol,Psps%filpsp) 565 end if ! HF decoupling 566 567 if (Dtset%pawcross==1) then 568 if (allocated(Paw_onsite) ) then 569 ABI_DT_FREE(Paw_onsite) 570 end if 571 ABI_DT_MALLOC(Paw_onsite,(Cryst%natom)) 572 call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,& 573 & Cryst%rprimd,Cryst%xcart,Pawtab,Pawrad,Pawfgrtab) 574 end if 575 end if 576 577 ! Allocate these arrays anyway, since they are passed to subroutines. 578 if (.not.allocated(ks_nhat)) then 579 ABI_MALLOC(ks_nhat,(nfftf,0)) 580 end if 581 if (.not.allocated(dijexc_core)) then 582 ABI_MALLOC(dijexc_core,(1,1,0)) 583 end if 584 585 ! ================================================== 586 ! ==== Read KS band structure from the KSS file ==== 587 ! ================================================== 588 ! 589 ! Initialize Wfd, allocate wavefunctions and precalculate tables to do the FFT using the coarse gwc_ngfft. 590 mband=Sigp%nbnds 591 ABI_MALLOC(bks_mask,(mband,Kmesh%nibz,Sigp%nsppol)) 592 ABI_MALLOC(keep_ur ,(mband,Kmesh%nibz,Sigp%nsppol)) 593 keep_ur=.FALSE.; bks_mask=.FALSE. 594 595 ABI_MALLOC(nband,(Kmesh%nibz,Sigp%nsppol)) 596 nband=mband 597 598 ! autoparal section 599 if (dtset%max_ncpus/=0) then 600 ount =ab_out 601 ! Temporary table needed to estimate memory 602 ABI_MALLOC(nlmn_atm,(Cryst%natom)) 603 if (Dtset%usepaw==1) then 604 do iat=1,Cryst%natom 605 nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size 606 end do 607 end if 608 609 write(ount,'(a)')"--- !Autoparal" 610 write(ount,"(a)")'#Autoparal section for Sigma runs.' 611 write(ount,"(a)") "info:" 612 write(ount,"(a,i0)")" autoparal: ",dtset%autoparal 613 write(ount,"(a,i0)")" max_ncpus: ",dtset%max_ncpus 614 write(ount,"(a,i0)")" gwpara: ",dtset%gwpara 615 write(ount,"(a,i0)")" nkpt: ",dtset%nkpt 616 write(ount,"(a,i0)")" nsppol: ",dtset%nsppol 617 write(ount,"(a,i0)")" nspinor: ",dtset%nspinor 618 write(ount,"(a,i0)")" nbnds: ",Sigp%nbnds 619 620 work_size = mband * Kmesh%nibz * Sigp%nsppol 621 622 ! Non-scalable memory in Mb i.e. memory that is not distribute with MPI. 623 nonscal_mem = (two*gwpc*Er%npwe**2*Er%nomega*(Er%mqmem+1)*b2Mb) * 1.1_dp 624 625 ! List of configurations. 626 ! Assuming an OpenMP implementation with perfect speedup! 627 write(ount,"(a)")"configurations:" 628 do ii=1,dtset%max_ncpus 629 nstates_per_proc = 0 630 eff = HUGE(one) 631 max_wfsmem_mb = zero 632 633 do my_rank=0,ii-1 634 call sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,ii,my_spins,bks_mask,keep_ur,ierr) 635 ABI_FREE(my_spins) 636 if (ierr /= 0) exit 637 my_nbks = COUNT(bks_mask) 638 nstates_per_proc = MAX(nstates_per_proc, my_nbks) 639 eff = MIN(eff, (one * work_size) / (ii * nstates_per_proc)) 640 641 ! Memory needed for Fourier components ug. 642 ug_mem = two*gwpc*Dtset%nspinor*Sigp%npwwfn*my_nbks*b2Mb 643 ! Memory needed for real space ur (use gwc_nfft, instead of gwx_nfft) 644 ur_mem = two*gwpc*Dtset%nspinor*gwc_nfft*COUNT(keep_ur)*b2Mb 645 ! Memory needed for PAW projections Cprj 646 cprj_mem = zero 647 if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb 648 max_wfsmem_mb = MAX(max_wfsmem_mb, ug_mem + ur_mem + cprj_mem) 649 end do 650 if (ierr /= 0) cycle 651 652 ! Add the non-scalable part and increase by 10% to account for other datastructures. 653 mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp 654 do omp_ncpus=1,xomp_get_max_threads() 655 write(ount,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus 656 write(ount,"(a,i0)")" mpi_ncpus: ",ii 657 write(ount,"(a,i0)")" omp_ncpus: ",omp_ncpus 658 write(ount,"(a,f12.9)")" efficiency: ",eff 659 write(ount,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb 660 end do 661 end do 662 write(ount,'(a)')"..." 663 664 ABI_FREE(nlmn_atm) 665 MSG_ERROR_NODUMP("aborting now") 666 667 else 668 call sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,nprocs,my_spins,bks_mask,keep_ur,ierr) 669 ABI_CHECK(ierr==0, "Error in sigma_bksmask") 670 end if 671 672 ! Then each node owns the wavefunctions where GW corrections are required. 673 do isp=1,SIZE(my_spins) 674 spin = my_spins(isp) 675 do ikcalc=1,Sigp%nkptgw 676 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW 677 ii=Sigp%minbnd(ikcalc,spin); jj=Sigp%maxbnd(ikcalc,spin) 678 bks_mask(ii:jj,ik_ibz,spin) = .TRUE. 679 if (MODULO(Dtset%gwmem,10)==1) keep_ur(ii:jj,ik_ibz,spin)=.TRUE. 680 end do 681 end do 682 683 ABI_FREE(my_spins) 684 opt_ecut=Dtset%ecutwfn 685 !opt_ecut=zero 686 687 call wfd_init(Wfd,Cryst,Pawtab,Psps,keep_ur,Dtset%paral_kgb,Sigp%npwwfn,mband,nband,Kmesh%nibz,Sigp%nsppol,bks_mask,& 688 & Dtset%nspden,Dtset%nspinor,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,& 689 & Gsph_Max%gvec,Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,opt_ecut=opt_ecut) 690 691 if (Dtset%pawcross==1) then 692 call wfd_init(Wfdf,Cryst,Pawtab,Psps,keep_ur,Dtset%paral_kgb,Sigp%npwwfn,mband,nband,Kmesh%nibz,Sigp%nsppol,bks_mask,& 693 & Dtset%nspden,Dtset%nspinor,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,& 694 & Gsph_Max%gvec,Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,opt_ecut=opt_ecut) 695 end if 696 697 ABI_FREE(bks_mask) 698 ABI_FREE(nband) 699 ABI_FREE(keep_ur) 700 701 call timab(402,2,tsec) ! sigma(Init1) 702 call timab(404,1,tsec) ! rdkss 703 704 call wfd_read_wfk(Wfd,wfk_fname,iomode_from_fname(wfk_fname)) 705 706 if (Dtset%pawcross==1) then 707 call wfd_copy(Wfd,Wfdf) 708 call wfd_change_ngfft(Wfdf,Cryst,Psps,ngfftf) 709 end if 710 711 ! This test has been disabled (too expensive!) 712 if (.False.) call wfd_test_ortho(Wfd,Cryst,Pawtab,unit=ab_out,mode_paral="COLL") 713 714 call timab(404,2,tsec) ! rdkss 715 call timab(405,1,tsec) ! Init2 716 717 !Debugging section. 718 !if (.TRUE.) then 719 if (.FALSE.) then 720 ! 721 if (.FALSE..and.Wfd%usepaw==1) then 722 ABI_DT_MALLOC(Cp1,(Wfd%natom,Wfd%nspinor)) 723 call pawcprj_alloc(Cp1,0,Wfd%nlmn_atm) 724 ABI_DT_MALLOC(Cp2,(Wfd%natom,Wfd%nspinor)) 725 call pawcprj_alloc(Cp2,0,Wfd%nlmn_atm) 726 727 call wfd_change_ngfft(Wfd,Cryst,Psps,ngfftf) 728 729 do spin=1,Wfd%nsppol 730 do ik_bz=1,Kmesh%nbz 731 ik_ibz=Kmesh%tab(ik_bz) 732 do band=1,Wfd%nband(ik_ibz,spin) 733 call paw_check_symcprj(Wfd,ik_bz,band,spin,1,Cryst,Kmesh,Psps,Pawtab,Pawang,Cp1) 734 call paw_check_symcprj(Wfd,ik_bz,band,spin,2,Cryst,Kmesh,Psps,Pawtab,Pawang,Cp2) 735 736 do iat=1,Cryst%natom 737 do isp=1,Wfd%nspinor 738 write(789,'(3i2,/,(f8.4))') band,ik_bz,spin,Cp1(iat,isp)%cp 739 write(790,'(3i2,/,(f8.4))') band,ik_bz,spin,Cp2(iat,isp)%cp 740 write(791,'(3i2,/,(f8.4))') band,ik_bz,spin,Cp1(iat,isp)%cp(1,:)**2 + Cp1(iat,isp)%cp(2,:)**2 741 write(792,'(3i2,/,(f8.4))') band,ik_bz,spin,Cp2(iat,isp)%cp(1,:)**2 + Cp2(iat,isp)%cp(2,:)**2 742 end do 743 end do 744 end do 745 end do 746 end do 747 748 call pawcprj_free(Cp1) 749 ABI_DT_FREE(Cp1) 750 call pawcprj_free(Cp2) 751 ABI_DT_FREE(Cp2) 752 end if 753 754 end if 755 756 ! ============================================================== 757 ! ==== Find little group of the k-points for GW corrections ==== 758 ! ============================================================== 759 ! * The little group is calculated only if sys_sigma. 760 ! * If use_umklp==1 then also symmetries requiring an umklapp to preserve k_gw are included. 761 ! 762 ABI_DT_MALLOC(Ltg_k,(Sigp%nkptgw)) 763 use_umklp=1 764 do ikcalc=1,Sigp%nkptgw 765 if (Sigp%symsigma/=0) then 766 call littlegroup_init(Sigp%kptgw(:,ikcalc),Qmesh,Cryst,use_umklp,Ltg_k(ikcalc),0) 767 end if 768 end do 769 ! 770 !=== Compute structure factor phases and large sphere cut-off === 771 !WARNING cannot use Dtset%mgfft, this has to be checked better 772 !mgfft=MAXVAL(ngfftc(:)) 773 !allocate(ph1d(2,3*(2*mgfft+1)*Cryst%natom),ph1df(2,3*(2*mgfftf+1)*Cryst%natom)) 774 !write(std_out,*)' CHECK ',Dtset%mgfftdg,mgfftf 775 !if (Dtset%mgfftdg/=mgfftf) then 776 !write(std_out,*)"WARNING Dtset%mgfftf /= mgfftf" 777 !write(std_out,*)'HACKING Dtset%mgfftf' 778 !Dtset%mgfftdg=mgfftf 779 !end if 780 ABI_MALLOC(ph1d,(2,3*(2*Dtset%mgfft+1)*Cryst%natom)) 781 ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*Cryst%natom)) 782 783 call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred) 784 785 if (Psps%usepaw==1.and.Pawfgr%usefinegrid==1) then 786 call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred) 787 else 788 ph1df(:,:)=ph1d(:,:) 789 end if 790 791 !=================================================================================== 792 !==== Classify the GW wavefunctions according to the irreducible representation ==== 793 !=================================================================================== 794 !* Warning still under development. 795 !* Only for SCGW. 796 !bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw 797 798 ABI_DT_MALLOC(KS_sym,(Wfd%nkibz,Wfd%nsppol)) 799 800 if (Sigp%symsigma==1.and.gwcalctyp>=20) then 801 ! call check_zarot(Gsph_c%ng,Cryst,gwc_ngfft,Gsph_c%gvec,Psps,Pawang,Gsph_c%rottb,Gsph_c%rottbm1) 802 use_paw_aeur=.FALSE. ! should pass ngfftf but the dense mesh is not forced to be symmetric 803 do spin=1,Wfd%nsppol 804 do ikcalc=1,Sigp%nkptgw 805 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 806 first_band = Sigp%minbnd(ikcalc,spin) 807 last_band = Sigp%maxbnd(ikcalc,spin) 808 ! call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,Cryst,KS_BSt,Pawtab,Pawrad,Pawang,Psps,& 809 call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,Wfd%ngfft,Cryst,KS_BSt,Pawtab,Pawrad,Pawang,Psps,& 810 & Dtset%tolsym,KS_sym(ik_ibz,spin)) 811 end do 812 end do 813 ! Recreate the Sig_ij tables taking advantage of the classification of the bands. 814 call sigma_tables(Sigp,Kmesh,KS_sym) 815 end if 816 817 call timab(405,2,tsec) ! Init2 818 call timab(406,1,tsec) ! make_vhxc 819 820 !=========================== 821 !=== COMPUTE THE DENSITY === 822 !=========================== 823 ! * Evaluate the planewave part (complete charge in case of NC pseudos). 824 ABI_MALLOC(ks_rhor,(nfftf,Dtset%nspden)) 825 ABI_MALLOC(ks_taur,(nfftf,Dtset%nspden*Dtset%usekden)) 826 827 call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,KS_BSt,ngfftf,nfftf,ks_rhor) 828 829 if (Dtset%usekden==1) then 830 call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,KS_BSt,ngfftf,nfftf,ks_taur,optcalc=1) 831 end if 832 833 !======================================== 834 !==== Additional computation for PAW ==== 835 !======================================== 836 nhatgrdim=0 837 if (Dtset%usepaw==1) then 838 ! Calculate the compensation charge nhat. 839 if (Dtset%xclevel==2) nhatgrdim=usexcnhat*Dtset%pawnhatxc 840 cplex=1; ider=2*nhatgrdim; izero=0 841 if (nhatgrdim>0) then 842 ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim)) 843 end if 844 if (nhatgrdim==0) then 845 ABI_MALLOC(ks_nhatgr,(0,0,0)) 846 end if 847 848 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,Cryst%gprimd,& 849 & Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,& 850 & Pawfgrtab,ks_nhatgr,ks_nhat,KS_Pawrhoij,KS_Pawrhoij,Pawtab,k0,Cryst%rprimd,& 851 & Cryst%ucvol,dtset%usewvl,Cryst%xred) 852 853 !if (nhatgrdim==0) then 854 ! ABI_FREE(ks_nhatgr) 855 !end if 856 857 ! === Evaluate onsite energies, potentials, densities === 858 ! * Initialize variables/arrays related to the PAW spheres. 859 ! * Initialize also lmselect (index of non-zero LM-moments of densities). 860 ABI_DT_MALLOC(KS_paw_ij,(Cryst%natom)) 861 ! cplex=1;cplex_dij=Dtset%nspinor 862 has_dijso=Dtset%pawspnorb; has_dijU=Dtset%usepawu 863 call paw_ij_nullify(KS_paw_ij) 864 call paw_ij_init(KS_paw_ij,cplex,Dtset%nspinor,Dtset%nsppol,& 865 & Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,& 866 & has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=1,has_dijxc_hat=1,has_dijxc_val=1,& 867 & has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1) 868 869 nkxc1=0 870 ABI_DT_MALLOC(KS_paw_an,(Cryst%natom)) 871 call paw_an_nullify(KS_paw_an) 872 call paw_an_init(KS_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,Dtset%nspden,& 873 & cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=1) 874 ! 875 ! Calculate onsite vxc with and without core charge. 876 nzlmopt=-1; option=0; compch_sph=greatest_real 877 call pawdenpot(compch_sph,KS_energies%e_paw,KS_energies%e_pawdc,ipert,& 878 & Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,& 879 & Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,KS_Paw_an,KS_Paw_an,KS_paw_ij,& 880 & Pawang,Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,& 881 & Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Cryst%ucvol,Psps%znuclpsp) 882 883 ! TO BE REMOVED ASAP !!! 884 _IBM6("Silly write for IBM6") 885 else 886 ABI_MALLOC(ks_nhatgr,(0,0,0)) 887 ABI_DT_MALLOC(KS_paw_ij,(0)) 888 ABI_DT_MALLOC(KS_paw_an,(0)) 889 end if !PAW 890 891 !if (.not.allocated(ks_nhatgr)) then 892 ! ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,0)) 893 !end if 894 895 call test_charge(nfftf,KS_BSt%nelect,Dtset%nspden,ks_rhor,Cryst%ucvol,& 896 & Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf) 897 ! 898 !For PAW, add the compensation charge on the FFT mesh, then get rho(G). 899 if (Dtset%usepaw==1) ks_rhor=ks_rhor+ks_nhat 900 901 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_rhor,ucvol=ucvol) 902 if(Dtset%usekden==1)then 903 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_taur,optrhor=1,ucvol=ucvol) 904 end if 905 906 ABI_MALLOC(ks_rhog,(2,nfftf)) 907 ABI_MALLOC(ks_taug,(2,nfftf*Dtset%usekden)) 908 call fourdp(1,ks_rhog,ks_rhor(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp5) 909 if(Dtset%usekden==1)then 910 call fourdp(1,ks_taug,ks_taur(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp5) 911 end if 912 913 !The following steps have been gathered in the setvtr routine: 914 !- get Ewald energy and Ewald forces 915 !- compute local ionic pseudopotential vpsp 916 !- eventually compute 3D core electron density xccc3d 917 !- eventually compute vxc and vhartr 918 !- set up ks_vtrial 919 ! 920 !******************************************************************* 921 !**** NOTE THAT HERE Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION **** 922 !******************************************************************* 923 924 ngrvdw=0 925 ABI_MALLOC(grvdw,(3,ngrvdw)) 926 ABI_MALLOC(grchempottn,(3,Cryst%natom)) 927 ABI_MALLOC(grewtn,(3,Cryst%natom)) 928 nkxc=0 929 if (Dtset%nspden==1) nkxc=2 930 if (Dtset%nspden>=2) nkxc=3 ! check GGA and spinor, quite a messy part!!! 931 !In case of MGGA, fxc and kxc are not available and we dont need them 932 !for the screening part (for now ...) 933 if (Dtset%ixc<0.and.libxc_functionals_ismgga()) nkxc=0 934 if (nkxc/=0) then 935 ABI_MALLOC(kxc,(nfftf,nkxc)) 936 end if 937 938 n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf 939 ABI_MALLOC(xccc3d,(n3xccc)) 940 ABI_MALLOC(ks_vhartr,(nfftf)) 941 ABI_MALLOC(ks_vtrial,(nfftf,Dtset%nspden)) 942 ABI_MALLOC(vpsp,(nfftf)) 943 ABI_MALLOC(ks_vxc,(nfftf,Dtset%nspden)) 944 945 optene=4; moved_atm_inside=0; moved_rhor=0; istep=1 946 947 call setvtr(Cryst%atindx1,Dtset,KS_energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,& 948 & istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,& 949 & Cryst%nattyp,nfftf,ngfftf,ngrvdw,ks_nhat,ks_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,& 950 & optene,pawrad,Pawtab,ph1df,Psps,ks_rhog,ks_rhor,Cryst%rmet,Cryst%rprimd,strsxc,& 951 & Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,Wvl,xccc3d,Cryst%xred,taug=ks_taug,taur=ks_taur) 952 953 !============================ 954 !==== Compute KS PAW Dij ==== 955 !============================ 956 if (Dtset%usepaw==1) then 957 !TO BE REMOVED 958 _IBM6("Another silly write for IBM6") 959 call timab(561,1,tsec) 960 961 ! Calculate the unsymmetrized Dij. 962 ipert=0; idir=0 963 call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert,& 964 & Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),& 965 & Dtset%nspden,Cryst%ntypat,KS_paw_an,KS_paw_ij,Pawang,Pawfgrtab,& 966 & Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,& 967 & k0,Dtset%spnorbscl,Cryst%ucvol,dtset%charge,ks_vtrial,ks_vxc,Cryst%xred,& 968 & nucdipmom=Dtset%nucdipmom) 969 970 ! Symmetrize KS Dij 971 #if 0 972 call symdij(Cryst%gprimd,Cryst%indsym,ipert,& 973 & Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,0,KS_paw_ij,Pawang,& 974 & Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec) 975 #else 976 call symdij_all(Cryst%gprimd,Cryst%indsym,ipert,& 977 & Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,KS_paw_ij,Pawang,& 978 & Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec) 979 #endif 980 ! 981 ! Output the pseudopotential strengths Dij and the augmentation occupancies Rhoij. 982 call pawprt(Dtset,Cryst%natom,KS_paw_ij,KS_Pawrhoij,Pawtab) 983 call timab(561,2,tsec) 984 end if 985 986 call timab(406,2,tsec) ! make_vhxc 987 988 !=== Calculate Vxc(b1,b2,k,s)=<b1,k,s|v_{xc}|b2,k,s> for all the states included in GW === 989 ! * ks_vxcvalme is calculated without NLCC, ks_vxcme contains NLCC (if any) 990 ! * This part is parallelized within MPI_COMM_WORD since each node has all GW wavefunctions. 991 ! * ks_vUme is zero unless we are using LDA+U as starting point, see calc_vHxc_braket 992 ! * Note that vH matrix elements are calculated using the true uncutted interaction. 993 call timab(407,1,tsec) ! vHxc_me 994 995 call melflags_reset(KS_mflags) 996 KS_mflags%has_vhartree=1 997 KS_mflags%has_vxc =1 998 KS_mflags%has_vxcval =1 999 if (Dtset%usepawu>0 ) KS_mflags%has_vu =1 1000 if (Dtset%useexexch>0 ) KS_mflags%has_lexexch=1 1001 if (Sigp%use_sigxcore==1) KS_mflags%has_sxcore =1 1002 if (gwcalctyp<10 ) KS_mflags%only_diago =1 ! off-diagonal elements only for SC on wavefunctions. 1003 1004 if (.FALSE.) then ! quick and dirty hack to test HF contribution. 1005 MSG_WARNING("testing on-site HF") 1006 lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size) 1007 ABI_MALLOC(dij_hf,(cplex_dij*lmn2_size_max,ndij,Cryst%natom)) 1008 call paw_dijhf(ndij,cplex_dij,lmn2_size_max,Cryst%natom,Cryst%ntypat,Pawtab,Pawrad,Pawang,& 1009 & KS_Pawrhoij,dij_hf,Dtset%prtvol) 1010 1011 do iat=1,Cryst%natom 1012 itypat = Cryst%typat(iat) 1013 ii = Pawtab(itypat)%lmn2_size 1014 KS_Paw_ij(iat)%dijxc(:,:) = dij_hf(1:cplex_dij*ii,:,iat) 1015 KS_Paw_ij(iat)%dijxc_hat(:,:) = zero 1016 end do 1017 ABI_FREE(dij_hf) 1018 1019 ! option_dij=3 1020 ! call symdij(Cryst%gprimd,Cryst%indsym,ipert,& 1021 ! & Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,& 1022 ! & KS_paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,& 1023 ! & Cryst%symafm,Cryst%symrec) 1024 end if 1025 1026 ABI_MALLOC(tmp_kstab,(2,Wfd%nkibz,Wfd%nsppol)) 1027 tmp_kstab=0 1028 do spin=1,Sigp%nsppol 1029 do ikcalc=1,Sigp%nkptgw ! No spin dependent! 1030 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 1031 tmp_kstab(1,ik_ibz,spin)=Sigp%minbnd(ikcalc,spin) 1032 tmp_kstab(2,ik_ibz,spin)=Sigp%maxbnd(ikcalc,spin) 1033 end do 1034 end do 1035 1036 call calc_vhxc_me(Wfd,KS_mflags,KS_me,Cryst,Dtset,nfftf,ngfftf,& 1037 & ks_vtrial,ks_vhartr,ks_vxc,Psps,Pawtab,KS_paw_an,Pawang,Pawfgrtab,KS_paw_ij,dijexc_core,& 1038 & ks_rhor,usexcnhat,ks_nhat,ks_nhatgr,nhatgrdim,tmp_kstab,taug=ks_taug,taur=ks_taur) 1039 ABI_FREE(tmp_kstab) 1040 1041 !#ifdef DEV_HAVE_SCGW_SYM 1042 !Set KS matrix elements connecting different irreps to zero. Do not touch unknown bands!. 1043 if (gwcalctyp>=20 .and. Sigp%symsigma > 0) then 1044 bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw 1045 ABI_MALLOC(ks_irreptab,(bmin:bmax,Kmesh%nibz,Sigp%nsppol)) 1046 ks_irreptab=0 1047 do spin=1,Sigp%nsppol 1048 do ikcalc=1,Sigp%nkptgw 1049 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 1050 first_band = Sigp%minbnd(ikcalc,spin) 1051 last_band = Sigp%maxbnd(ikcalc,spin) 1052 if (.not.esymm_failed(KS_sym(ik_ibz,spin))) then 1053 ks_irreptab(first_band:last_band,ik_ibz,spin) = KS_sym(ik_ibz,spin)%b2irrep(first_band:last_band) 1054 ! ks_irreptab(bmin:bmax,ik_ibz,spin) = KS_sym(ik_ibz,spin)%b2irrep(bmin:bmax) 1055 end if 1056 end do 1057 end do 1058 call melements_zero(KS_me,ks_irreptab) 1059 ABI_FREE(ks_irreptab) 1060 end if 1061 !#endif 1062 1063 call melements_print(KS_me,header="Matrix elements in the KS basis set",prtvol=Dtset%prtvol) 1064 ! 1065 ! If possible, calculate the EXX energy from the between the frozen core 1066 ! and the valence electrons using KS wavefunctions 1067 ! MG: BE careful here, since ex_energy is meaningful only is all occupied states are calculated. 1068 if( KS_mflags%has_sxcore ==1 ) then 1069 ! TODO 1070 !ex_energy = mels_get_exene_core(KS_me,kmesh,KS_BSt) 1071 ex_energy=zero 1072 do spin=1,Sigp%nsppol 1073 do ik=1,Kmesh%nibz 1074 do ib=b1gw,b2gw 1075 if (Sigp%nsig_ab==1) then 1076 ex_energy = ex_energy + half*KS_BSt%occ(ib,ik,spin)*Kmesh%wt(ik)*KS_me%sxcore(ib,ib,ik,spin) 1077 else 1078 ex_energy = ex_energy + half*KS_BSt%occ(ib,ik,spin)*Kmesh%wt(ik)*SUM(KS_me%sxcore(ib,ib,ik,:)) 1079 end if 1080 end do 1081 end do 1082 end do 1083 write(msg,'(a,2(es16.6,a))')' CORE Exchange energy with KS wavefunctions: ',ex_energy,' Ha ,',ex_energy*Ha_eV,' eV' 1084 call wrtout(std_out,msg,'COLL') 1085 end if 1086 1087 call timab(407,2,tsec) ! vHxc_me 1088 1089 call timab(408,1,tsec) ! hqp_init 1090 1091 ! Do not break this coding! When gwcalctyp>10, the order of the bands can be interexchanged after 1092 ! the diagonalization. Therefore, we have to correctly assign the matrix elements to the corresponding 1093 ! bands and we cannot skip the following even though it looks unuseful. 1094 if (gwcalctyp>=10) then 1095 call wrtout(std_out,ch10//' *************** KS Energies *******************','COLL') 1096 end if 1097 1098 !=== QP_BSt stores energies and occ. used for the calculation === 1099 ! * Initialize QP_BSt with KS values. 1100 ! * In case of SC update QP_BSt using the QPS file. 1101 call ebands_copy(KS_BSt,QP_BSt) 1102 1103 ABI_MALLOC(qp_rhor,(nfftf,Dtset%nspden)) 1104 ABI_MALLOC(qp_taur,(nfftf,Dtset%nspden*Dtset%usekden)) 1105 QP_sym => KS_sym 1106 1107 if (gwcalctyp<10) then 1108 ! one-shot GW, just do a copy of the KS density. 1109 qp_rhor=ks_rhor 1110 if(Dtset%usekden==1)qp_taur=ks_taur 1111 QP_sym => KS_sym 1112 else 1113 ! Self-consistent GW. 1114 ! Read the unitary matrix and the QP energies of the previous step from the QPS file. 1115 call energies_init(QP_energies) 1116 QP_energies%e_corepsp=ecore/Cryst%ucvol 1117 1118 ! m_lda_to_qp(ib,jb,k,s) := <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}> 1119 ! Initialize the QP amplitudes with KS wavefunctions. 1120 ABI_MALLOC(Sr%m_lda_to_qp,(Sigp%nbnds,Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)) 1121 Sr%m_lda_to_qp=czero 1122 do ib=1,Sigp%nbnds 1123 Sr%m_lda_to_qp(ib,ib,:,:)=cone 1124 end do 1125 1126 ! Now read m_lda_to_qp and update the energies in QP_BSt. 1127 ! TODO switch on the renormalization of n in sigma. 1128 ABI_MALLOC(prev_rhor,(nfftf,Dtset%nspden)) 1129 ABI_MALLOC(prev_taur,(nfftf,Dtset%nspden*Dtset%usekden)) 1130 ABI_DT_MALLOC(prev_Pawrhoij,(Cryst%natom*Psps%usepaw)) 1131 1132 call rdqps(QP_BSt,Dtfil%fnameabi_qps,Dtset%usepaw,Dtset%nspden,1,nscf,& 1133 & nfftf,ngfftf,Cryst%ucvol,Dtset%paral_kgb,Cryst,Pawtab,MPI_enreg_seq,nbsc,Sr%m_lda_to_qp,prev_rhor,prev_Pawrhoij) 1134 1135 ! Find the irreps associated to the QP amplitudes starting from the analogous table for the KS states. 1136 ! bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw 1137 ! allocate(qp_irreptab(bmin:bmax,Kmesh%nibz,Sigp%nsppol)) 1138 ! qp_irreptab=0 1139 ! !qp_irreptab=ks_irreptab 1140 1141 ! do jb_qp=bmin,bmax 1142 ! do ib_ks=bmin,bmax 1143 ! if (ABS(Sr%m_lda_to_qp(ib_ks,jb_qp,ik_ibz,spin)) > tol12) then ! jb_qp has same the same character as ib_ks. 1144 ! ks_irr = ks_irreptab(ib_ks,ib_ks,ik_ibz,spin) 1145 ! qp_irreptab(jb_qp,jb_qp,ik_ibz,spin) = ks_irr 1146 ! do ii=bmin,bmax 1147 ! if (ks_irr == ks_irreptab(ii,ib_ks,ik_ibz,spin)) then 1148 ! qp_irreptab(jb_qp,ii,ik_ibz,spin) = ks_irr 1149 ! end if 1150 ! end do 1151 ! end if 1152 ! end do 1153 ! end do 1154 1155 if (nscf==0) prev_rhor=ks_rhor 1156 if (nscf==0 .and. Dtset%usekden==1) prev_taur=ks_taur 1157 1158 if (nscf>0.and.gwcalctyp>=20.and.wfd_iam_master(Wfd)) then 1159 ! Print the unitary transformation on std_out. 1160 call show_QP(QP_BSt,Sr%m_lda_to_qp,fromb=Sigp%minbdgw,tob=Sigp%maxbdgw,unit=std_out,tolmat=0.001_dp) 1161 end if 1162 1163 ! Compute QP wfg as linear combination of KS states === 1164 ! * Wfd%ug is modified inside calc_wf_qp 1165 ! * For PAW, update also the on-site projections. 1166 ! * WARNING the first dimension of MPI_enreg MUST be Kmesh%nibz 1167 ! TODO here we should use nbsc instead of nbnds 1168 1169 call wfd_rotate(Wfd,Cryst,Sr%m_lda_to_qp) 1170 1171 ! Reinit the storage mode of Wfd as ug have been changed === 1172 ! Update also the wavefunctions for GW corrections on each processor 1173 call wfd_reset_ur_cprj(Wfd) 1174 1175 ! This test has been disabled (too expensive!) 1176 if (.False.) call wfd_test_ortho(Wfd, Cryst, Pawtab, unit=std_out) 1177 1178 ! Compute QP occupation numbers. 1179 call wrtout(std_out,'sigma: calculating QP occupation numbers:','COLL') 1180 call ebands_update_occ(QP_BSt,Dtset%spinmagntarget,prtvol=0) 1181 qp_vbik(:,:) = get_valence_idx(QP_BSt) 1182 1183 ! #ifdef DEV_HAVE_SCGW_SYM 1184 ! Calculate the irreducible representations of the new QP amplitdues. 1185 if (Sigp%symsigma==1.and.gwcalctyp>=20) then 1186 ABI_DT_MALLOC(QP_sym,(Wfd%nkibz,Wfd%nsppol)) 1187 use_paw_aeur=.FALSE. ! should pass ngfftf but the dense mesh is not forced to be symmetric 1188 do spin=1,Wfd%nsppol 1189 do ikcalc=1,Sigp%nkptgw 1190 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 1191 ! Quick fix for SCGW+symm TODO fix properly! 1192 first_band = Sigp%minbnd(ikcalc,spin) 1193 last_band = Sigp%maxbnd(ikcalc,spin) 1194 ! first_band = MINVAL(Sigp%minbnd(:,spin)) 1195 ! last_band = MAXVAL(Sigp%maxbnd(:,spin)) 1196 ! call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,Cryst,QP_BSt,Pawtab,Pawrad,Pawang,Psps,& 1197 call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,Wfd%ngfft,Cryst,QP_BSt,Pawtab,Pawrad,Pawang,Psps,& 1198 & Dtset%tolsym,QP_sym(ik_ibz,spin)) 1199 end do 1200 end do 1201 1202 ! Recreate the Sig_ij tables taking advantage of the classification of the bands. 1203 call sigma_tables(Sigp,Kmesh,QP_sym) 1204 end if 1205 ! #endif 1206 1207 ! Compute QP density using the updated wfg. 1208 call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,QP_BSt,ngfftf,nfftf,qp_rhor) 1209 if (Dtset%usekden==1) then 1210 call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,QP_BSt,ngfftf,nfftf,qp_taur,optcalc=1) 1211 end if 1212 1213 ! ======================================== 1214 ! ==== QP self-consistent GW with PAW ==== 1215 ! ======================================== 1216 if (Dtset%usepaw==1) then 1217 ABI_MALLOC(qp_nhat,(nfftf,Dtset%nspden)) 1218 nhatgrdim=0; if (Dtset%xclevel==2) nhatgrdim=usexcnhat 1219 ABI_MALLOC(qp_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim)) 1220 1221 ABI_DT_MALLOC(QP_pawrhoij,(Cryst%natom)) 1222 ABI_DT_MALLOC(QP_paw_ij,(Cryst%natom)) 1223 ABI_DT_MALLOC(QP_paw_an,(Cryst%natom)) 1224 1225 ! Calculate new QP quantities: nhat, nhatgr, rho_ij, paw_ij, and paw_an. 1226 call paw_qpscgw(Wfd,nscf,nfftf,ngfftf,Dtset,Cryst,Kmesh,Psps,QP_BSt,& 1227 & Pawang,Pawrad,Pawtab,Pawfgrtab,prev_Pawrhoij,& 1228 & QP_pawrhoij,QP_paw_ij,QP_paw_an,QP_energies,qp_nhat,nhatgrdim,qp_nhatgr,compch_sph,compch_fft) 1229 else 1230 ABI_MALLOC(qp_nhat,(0,0)) 1231 ABI_MALLOC(qp_nhatgr,(0,0,0)) 1232 ABI_DT_MALLOC(QP_pawrhoij,(0)) 1233 ABI_DT_MALLOC(QP_paw_ij,(0)) 1234 ABI_DT_MALLOC(QP_paw_an,(0)) 1235 end if 1236 1237 ! here I should renormalize the density 1238 call test_charge(nfftf,KS_BSt%nelect,Dtset%nspden,qp_rhor,Cryst%ucvol,& 1239 & Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf) 1240 1241 if (Dtset%usepaw==1) qp_rhor(:,:)=qp_rhor(:,:)+qp_nhat(:,:) ! Add the "hat" term. 1242 1243 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,qp_rhor,ucvol=ucvol) 1244 if(Dtset%usekden==1) then 1245 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,qp_taur,optrhor=1,ucvol=ucvol) 1246 end if 1247 1248 ! Simple mixing of the PW density to damp oscillations in the Hartree potential. 1249 if (nscf>0 .and. (ABS(Dtset%rhoqpmix-one)>tol12) ) then 1250 write(msg,'(2a,f6.3)')ch10,' sigma: mixing QP densities using rhoqpmix= ',Dtset%rhoqpmix 1251 call wrtout(std_out,msg,'COLL') 1252 qp_rhor = prev_rhor + Dtset%rhoqpmix*(qp_rhor-prev_rhor) 1253 if(Dtset%usekden==1)qp_taur = prev_taur + Dtset%rhoqpmix*(qp_taur-prev_taur) ! mix taur. 1254 end if 1255 1256 ABI_FREE(prev_rhor) 1257 ABI_FREE(prev_taur) 1258 if (Psps%usepaw==1.and.nscf>0) then 1259 call pawrhoij_free(prev_pawrhoij) 1260 end if 1261 ABI_DT_FREE(prev_pawrhoij) 1262 1263 ABI_MALLOC(qp_rhog,(2,nfftf)) 1264 ABI_MALLOC(qp_taug,(2,nfftf*Dtset%usekden)) 1265 call fourdp(1,qp_rhog,qp_rhor(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp5) 1266 if(Dtset%usekden==1)call fourdp(1,qp_taug,qp_taur(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp5) 1267 1268 ! =========================================== 1269 ! ==== Optional output of the QP density ==== 1270 ! =========================================== 1271 if (Dtset%prtden/=0.and.wfd_iam_master(Wfd)) then 1272 call fftdatar_write("qp_rhor",dtfil%fnameabo_qp_den,dtset%iomode,hdr_sigma,& 1273 cryst,ngfftf,cplex1,nfftf,dtset%nspden,qp_rhor,mpi_enreg_seq,ebands=qp_bst) 1274 end if 1275 1276 ! =========================================== 1277 ! === Optional output of the full QP density 1278 ! =========================================== 1279 if (Wfd%usepaw==1.and.Dtset%prtden==2) then 1280 ABI_MALLOC(qp_rhor_paw ,(nfftf,Wfd%nspden)) 1281 ABI_MALLOC(qp_rhor_n_one ,(nfftf,Wfd%nspden)) 1282 ABI_MALLOC(qp_rhor_nt_one,(nfftf,Wfd%nspden)) 1283 1284 call denfgr(Cryst%atindx1,Cryst%gmet,comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,qp_nhat,& 1285 & Wfd%nspinor,Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,QP_pawrhoij,Pawtab,Dtset%prtvol,& 1286 & qp_rhor,qp_rhor_paw,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 1287 1288 ABI_FREE(qp_rhor_n_one) 1289 ABI_FREE(qp_rhor_nt_one) 1290 if (Dtset%prtvol>9) then ! Print a normalisation check 1291 norm = SUM(qp_rhor_paw(:,1))*Cryst%ucvol/PRODUCT(Pawfgr%ngfft(1:3)) 1292 write(msg,'(a,F8.4)') ' QUASIPARTICLE DENSITY CALCULATED - NORM OF DENSITY: ',norm 1293 call wrtout(std_out,msg,'PERS') 1294 end if 1295 1296 ! Write the density to file 1297 if (my_rank==master) then 1298 call fftdatar_write("qp_pawrhor",dtfil%fnameabo_qp_pawden,dtset%iomode,hdr_sigma,& 1299 cryst,ngfftf,cplex1,nfftf,dtset%nspden,qp_rhor_paw,mpi_enreg_seq,ebands=qp_bst) 1300 end if 1301 ABI_FREE(qp_rhor_paw) 1302 end if 1303 1304 nkxc=0 1305 if (Dtset%nspden==1) nkxc=2 1306 if (Dtset%nspden>=2) nkxc=3 !check GGA and spinor that is messy !!! 1307 !In case of MGGA, fxc and kxc are not available and we dont need them 1308 !for the screening part (for now ...) 1309 if (Dtset%ixc<0.and.libxc_functionals_ismgga()) nkxc=0 1310 if (nkxc/=0) then 1311 ABI_MALLOC(qp_kxc,(nfftf,nkxc)) 1312 end if 1313 1314 ! **** NOTE THAT Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION **** 1315 n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf 1316 ABI_MALLOC(qp_vhartr,(nfftf)) 1317 ABI_MALLOC(qp_vtrial,(nfftf,Dtset%nspden)) 1318 ABI_MALLOC(qp_vxc,(nfftf,Dtset%nspden)) 1319 1320 optene=4; moved_atm_inside=0; moved_rhor=0; istep=1 1321 1322 call setvtr(Cryst%atindx1,Dtset,QP_energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,& 1323 & istep,qp_kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,& 1324 & Cryst%nattyp,nfftf,ngfftf,ngrvdw,qp_nhat,qp_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,& 1325 & optene,pawrad,Pawtab,ph1df,Psps,qp_rhog,qp_rhor,Cryst%rmet,Cryst%rprimd,strsxc,& 1326 & Cryst%ucvol,usexcnhat,qp_vhartr,vpsp,qp_vtrial,qp_vxc,vxcavg_qp,Wvl,& 1327 & xccc3d,Cryst%xred,taug=qp_taug,taur=qp_taur) 1328 1329 if (allocated(qp_kxc)) then 1330 ABI_FREE(qp_kxc) 1331 end if 1332 1333 if (Dtset%usepaw==1) then 1334 call timab(561,1,tsec) 1335 1336 ! Compute QP Dij 1337 ipert=0; idir=0 1338 call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert,& 1339 & Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),& 1340 & Dtset%nspden,Cryst%ntypat,QP_paw_an,QP_paw_ij,Pawang,Pawfgrtab,& 1341 & Dtset%pawprtvol,Pawrad,QP_pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,& 1342 & k0,Dtset%spnorbscl,Cryst%ucvol,dtset%charge,qp_vtrial,qp_vxc,Cryst%xred,& 1343 & nucdipmom=Dtset%nucdipmom) 1344 1345 ! Symmetrize total Dij 1346 option_dij=0 1347 1348 #if 0 1349 call symdij(Cryst%gprimd,Cryst%indsym,ipert,& 1350 & Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,& 1351 & QP_paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,& 1352 & Cryst%symafm,Cryst%symrec) 1353 #else 1354 call symdij_all(Cryst%gprimd,Cryst%indsym,ipert,& 1355 & Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,& 1356 & QP_paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,& 1357 & Cryst%symafm,Cryst%symrec) 1358 #endif 1359 1360 ! Output the QP pseudopotential strengths Dij and the augmentation occupancies Rhoij. 1361 call pawprt(Dtset,Cryst%natom,QP_paw_ij,QP_Pawrhoij,Pawtab) 1362 call timab(561,2,tsec) 1363 end if 1364 1365 ehartree=half*SUM(qp_rhor(:,1)*qp_vhartr(:))/DBLE(nfftf)*Cryst%ucvol 1366 1367 write(msg,'(a,80a)')ch10,('-',ii=1,80) 1368 call wrtout(ab_out,msg,'COLL') 1369 write(msg,'(5a,f9.4,3a,es21.14,2a,es21.14)')ch10,& 1370 & ' QP results after the unitary transformation in the KS subspace: ',ch10,ch10,& 1371 & ' Number of electrons = ',qp_rhog(1,1)*Cryst%ucvol,ch10,ch10,& 1372 & ' QP Band energy [Ha] = ',get_bandenergy(QP_BSt),ch10,& 1373 & ' QP Hartree energy [Ha] = ',ehartree 1374 call wrtout(ab_out,msg,'COLL') 1375 write(msg,'(a,80a)')ch10,('-',ii=1,80) 1376 call wrtout(ab_out,msg,'COLL') 1377 ! 1378 ! TODO Since plasmonpole model 2-3-4 depend on the Fourier components of the density 1379 ! in case of self-consistency we might calculate here the ppm coefficients using qp_rhor 1380 end if ! gwcalctyp>=10 1381 ! 1382 !=== KS hamiltonian hlda(b1,b1,k,s)= <b1,k,s|H_s|b1,k,s> === 1383 ABI_MALLOC(hlda,(b1gw:b2gw,b1gw:b2gw,Kmesh%nibz,Sigp%nsppol*Sigp%nsig_ab)) 1384 hlda = czero 1385 1386 if (Dtset%nspinor==1) then 1387 do spin=1,Sigp%nsppol 1388 do ik=1,Kmesh%nibz 1389 do ib=b1gw,b2gw 1390 hlda(ib,ib,ik,spin) = KS_BSt%eig(ib,ik,spin) 1391 end do 1392 end do 1393 end do 1394 else 1395 ! Spinorial case 1396 ! * Note that here vxc contains the contribution of the core. 1397 ! * Scale ovlp if orthonormalization is not satisfied as npwwfn might be < npwvec. 1398 ! TODO add spin-orbit case and gwpara 2 1399 ABI_MALLOC(my_band_list, (wfd%mband)) 1400 ABI_MALLOC(bmask, (wfd%mband)) 1401 bmask = .False.; bmask(b1gw:b2gw) = .True. 1402 1403 if (Wfd%usepaw==1) then 1404 ABI_DT_MALLOC(Cp1,(Wfd%natom,Wfd%nspinor)) 1405 call pawcprj_alloc(Cp1,0,Wfd%nlmn_atm) 1406 end if 1407 1408 do spin=1,Sigp%nsppol 1409 do ik_ibz=1,Kmesh%nibz 1410 ! Distribute bands in [b1gw, b2gw] range 1411 call wfd_distribute_bands(wfd, ik_ibz, spin, my_nband, my_band_list, bmask=bmask) 1412 if (my_nband == 0) cycle 1413 npw_k = Wfd%npwarr(ik_ibz) 1414 do ii=1,my_nband ! ib=b1gw,b2gw in sequential 1415 ib = my_band_list(ii) 1416 ug1 => Wfd%Wave(ib,ik_ibz,spin)%ug 1417 cdummy = xdotc(npw_k*Wfd%nspinor,ug1,1,ug1,1) 1418 ovlp(1) = REAL(cdummy); ovlp(2) = AIMAG(cdummy) 1419 1420 if (Psps%usepaw==1) then 1421 call wfd_get_cprj(Wfd,ib,ik_ibz,spin,Cryst,Cp1,sorted=.FALSE.) 1422 ovlp = ovlp + paw_overlap(Cp1,Cp1,Cryst%typat,Pawtab) 1423 end if 1424 ! write(std_out,*)ovlp(1),ovlp(2) 1425 norm=DBLE(ovlp(1)+ovlp(2)) 1426 ovlp(1)=DBLE(ovlp(1)/norm); ovlp(2)=DBLE(ovlp(2)/norm) ! ovlp(2)=cone-ovlp(1) 1427 hlda(ib,ib,ik_ibz,1) = KS_BSt%eig(ib,ik_ibz,1)*ovlp(1)-KS_me%vxc(ib,ib,ik_ibz,3) 1428 hlda(ib,ib,ik_ibz,2) = KS_BSt%eig(ib,ik_ibz,1)*ovlp(2)-KS_me%vxc(ib,ib,ik_ibz,4) 1429 hlda(ib,ib,ik_ibz,3) = KS_me%vxc(ib,ib,ik_ibz,3) 1430 hlda(ib,ib,ik_ibz,4) = KS_me%vxc(ib,ib,ik_ibz,4) 1431 end do 1432 end do 1433 end do 1434 1435 call xmpi_sum(hlda, wfd%comm, ierr) 1436 1437 ABI_FREE(my_band_list) 1438 ABI_FREE(bmask) 1439 if (Wfd%usepaw==1) then 1440 call pawcprj_free(Cp1) 1441 ABI_DT_FREE(Cp1) 1442 end if 1443 end if 1444 1445 !=== Initialize Sigma results === 1446 !TODO it is better if we use ragged arrays indexed by the k-point 1447 call sigma_init(Sigp,Kmesh%nibz,Dtset%usepawu,Sr) 1448 1449 ! Setup of the bare Hamiltonian := T + v_{loc} + v_{nl} + v_H. 1450 ! * The representation depends wheter we are updating the wfs or not. 1451 ! * ks_vUme is zero unless we are using LDA+U as starting point, see calc_vHxc_braket 1452 ! * Note that vH matrix elements are calculated using the true uncutted interaction. 1453 1454 if (gwcalctyp<10) then 1455 ! For one-shot GW use the KS representation. 1456 Sr%hhartree=hlda-KS_me%vxcval 1457 ! Additional goodies for PAW 1458 ! * LDA +U Hamiltonian 1459 ! * LEXX. 1460 ! * Core contribution estimated using Fock exchange. 1461 if (Dtset%usepaw==1) then 1462 if (Sigp%use_sigxcore==1) Sr%hhartree=hlda - (KS_me%vxc - KS_me%sxcore) 1463 if (Dtset%usepawu>0) Sr%hhartree=Sr%hhartree-KS_me%vu 1464 if (Dtset%useexexch>0) then 1465 MSG_ERROR("useexexch > 0 not implemented") 1466 Sr%hhartree = Sr%hhartree - KS_me%vlexx 1467 end if 1468 end if 1469 else 1470 ! Self-consistent on energies and|or wavefunctions. 1471 ! * For NC get the bare Hamiltonian $H_{bare}= T+v_{loc}+ v_{nl}$ in the KS representation 1472 ! * For PAW, calculate the matrix elements of h0, store also the new Dij in QP_Paw_ij. 1473 ! * h0 is defined as T+vH[tn+nhat+tnZc] + vxc[tnc] + dij_eff and 1474 ! dij_eff = dij^0 + dij^hartree + dij^xc-dij^xc_val + dijhat - dijhat_val. 1475 ! In the above expression tn, tnhat are QP quantities. 1476 if (Dtset%usepaw==0) then 1477 ABI_MALLOC(hbare,(b1gw:b2gw,b1gw:b2gw,Kmesh%nibz,Sigp%nsppol*Sigp%nsig_ab)) 1478 hbare=hlda-KS_me%vhartree-KS_me%vxcval 1479 1480 ! Change basis from KS to QP, hbare is overwritten: A_{QP} = U^\dagger A_{KS} U 1481 ABI_MALLOC(htmp,(b1gw:b2gw,b1gw:b2gw,Kmesh%nibz,Sigp%nsppol*Sigp%nsig_ab)) 1482 ABI_MALLOC(ctmp,(b1gw:b2gw,b1gw:b2gw)) 1483 ABI_MALLOC(uks2qp,(b1gw:b2gw,b1gw:b2gw)) 1484 1485 htmp=hbare; hbare=czero 1486 1487 do spin=1,Sigp%nsppol 1488 do ik=1,Kmesh%nibz 1489 uks2qp(:,:) = Sr%m_lda_to_qp(b1gw:b2gw,b1gw:b2gw,ik,spin) 1490 do iab=1,Sigp%nsig_ab 1491 is_idx=spin; if (Sigp%nsig_ab>1) is_idx=iab 1492 ctmp(:,:)=MATMUL(htmp(:,:,ik,is_idx),uks2qp(:,:)) 1493 hbare(:,:,ik,is_idx)=MATMUL(TRANSPOSE(CONJG(uks2qp)),ctmp) 1494 end do 1495 end do 1496 end do 1497 ABI_FREE(htmp) 1498 ABI_FREE(ctmp) 1499 ABI_FREE(uks2qp) 1500 end if ! usepaw==0 1501 1502 ! Calculate the QP matrix elements 1503 ! This part is parallelized within MPI_COMM_WORD since each node has all GW wavefunctions. 1504 ! For PAW, we have to construct the new bare Hamiltonian. 1505 call wrtout(std_out,ch10//' *************** QP Energies *******************','COLL') 1506 1507 call melflags_reset(QP_mflags) 1508 ! if (gwcalctyp<20) QP_mflags%only_diago=1 ! For e-only, no need of off-diagonal elements. 1509 QP_mflags%has_vhartree=1 1510 if (Dtset%usepaw==1) QP_mflags%has_hbare =1 1511 ! QP_mflags%has_vxc =1 1512 ! QP_mflags%has_vxcval =1 1513 ! if (Sigp%gwcalctyp >100) QP_mflags%has_vxcval_hybrid=1 1514 if (mod10==5 .and. & 1515 & (Dtset%ixc_sigma==-402 .or. Dtset%ixc_sigma==-406 .or. Dtset%ixc_sigma==-427 .or. Dtset%ixc_sigma==-428 .or.& 1516 & Dtset%ixc_sigma==-456 .or. Dtset%ixc_sigma==41 .or. Dtset%ixc_sigma==42)) then 1517 QP_mflags%has_vxcval_hybrid=1 1518 end if 1519 ! if (Sigp%use_sigxcore==1) QP_mflags%has_sxcore =1 1520 ! if (Dtset%usepawu>0) QP_mflags%has_vu =1 1521 ! if (Dtset%useexexch>0) QP_mflags%has_lexexch=1 1522 1523 ABI_MALLOC(tmp_kstab,(2,Wfd%nkibz,Wfd%nsppol)) 1524 tmp_kstab=0 1525 do spin=1,Sigp%nsppol 1526 do ikcalc=1,Sigp%nkptgw ! No spin dependent! 1527 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 1528 tmp_kstab(1,ik_ibz,spin)=Sigp%minbnd(ikcalc,spin) 1529 tmp_kstab(2,ik_ibz,spin)=Sigp%maxbnd(ikcalc,spin) 1530 end do 1531 end do 1532 1533 call calc_vhxc_me(Wfd,QP_mflags,QP_me,Cryst,Dtset,nfftf,ngfftf,& 1534 & qp_vtrial,qp_vhartr,qp_vxc,Psps,Pawtab,QP_paw_an,Pawang,Pawfgrtab,QP_paw_ij,dijexc_core,& 1535 & qp_rhor,usexcnhat,qp_nhat,qp_nhatgr,nhatgrdim,tmp_kstab,taug=qp_taug,taur=qp_taur) 1536 ABI_FREE(tmp_kstab) 1537 1538 ! #ifdef DEV_HAVE_SCGW_SYM 1539 if (gwcalctyp>=20 .and. Sigp%symsigma>0) then 1540 bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw 1541 ABI_MALLOC(qp_irreptab,(bmin:bmax,Kmesh%nibz,Sigp%nsppol)) 1542 qp_irreptab=0 1543 do spin=1,Sigp%nsppol 1544 do ikcalc=1,Sigp%nkptgw 1545 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 1546 first_band = Sigp%minbnd(ikcalc,spin) 1547 last_band = Sigp%maxbnd(ikcalc,spin) 1548 if (.not.esymm_failed(QP_sym(ik_ibz,spin))) then 1549 qp_irreptab(first_band:last_band,ik_ibz,spin) = QP_sym(ik_ibz,spin)%b2irrep(first_band:last_band) 1550 ! qp_irreptab(bmin:bmax,ik_ibz,spin) = QP_sym(ik_ibz,spin)%b2irrep(bmin:bmax) 1551 end if 1552 end do 1553 end do 1554 call melements_zero(QP_me,qp_irreptab) 1555 ABI_FREE(qp_irreptab) 1556 end if 1557 ! #endif 1558 1559 call melements_print(QP_me,header="Matrix elements in the QP basis set",prtvol=Dtset%prtvol) 1560 1561 ! * Output the QP pseudopotential strengths Dij and the augmentation occupancies Rhoij. 1562 if (Dtset%usepaw==1) then 1563 call wrtout(std_out," *** After calc_vHxc_braket *** ",'COLL') 1564 ! TODO terminate the implementation of this routine. 1565 call paw_ij_print(QP_Paw_ij,unit=std_out,pawprtvol=Dtset%pawprtvol,pawspnorb=Dtset%pawspnorb,mode_paral="COLL") 1566 call pawprt(Dtset,Cryst%natom,QP_paw_ij,QP_Pawrhoij,Pawtab) 1567 end if 1568 1569 if (Dtset%usepaw==0) then 1570 ! GA : We have an odd bug here. I have to unroll this loop, otherwise it 1571 ! might cause segfault when running on several nodes. 1572 ! 1573 ! Sr%hhartree = hbare + QP_me%vhartree 1574 if(QP_mflags%has_vxcval_hybrid==0) then 1575 do spin=1,Sigp%nsppol*Sr%nsig_ab 1576 do ikcalc=1,Sr%nkibz 1577 do ib1=b1gw,b2gw 1578 do ib2=b1gw,b2gw 1579 Sr%hhartree(ib2,ib1,ikcalc,spin) = hbare(ib2,ib1,ikcalc,spin) + QP_me%vhartree(ib2,ib1,ikcalc,spin) 1580 end do 1581 end do 1582 end do 1583 end do 1584 else 1585 do spin=1,Sigp%nsppol*Sr%nsig_ab 1586 do ikcalc=1,Sr%nkibz 1587 do ib1=b1gw,b2gw 1588 do ib2=b1gw,b2gw 1589 Sr%hhartree(ib2,ib1,ikcalc,spin) = hbare(ib2,ib1,ikcalc,spin) + & 1590 & QP_me%vhartree(ib2,ib1,ikcalc,spin) + & 1591 & QP_me%vxcval_hybrid(ib2,ib1,ikcalc,spin) 1592 end do 1593 end do 1594 end do 1595 end do 1596 end if 1597 else 1598 Sr%hhartree=QP_me%hbare 1599 end if 1600 1601 ! #ifdef DEV_HAVE_SCGW_SYM 1602 if (gwcalctyp>=20 .and. Sigp%symsigma > 0) then 1603 ! bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw 1604 do spin=1,Sigp%nsppol 1605 do ik_ibz=1,Kmesh%nibz 1606 if (.not.esymm_failed(QP_sym(ik_ibz,spin))) then 1607 bmin=Sigp%minbnd(ik_ibz,spin); bmax=Sigp%minbnd(ik_ibz,spin) 1608 do ib2=bmin,bmax 1609 irr_idx2 = QP_sym(ik_ibz,spin)%b2irrep(ib2) 1610 do ib1=bmin,bmax 1611 irr_idx1 = QP_sym(ik_ibz,spin)%b2irrep(ib1) 1612 if (irr_idx1/=irr_idx2 .and. ALL((/irr_idx1,irr_idx2/)/=0) ) Sr%hhartree(ib1,ib2,ik_ibz,spin) = czero 1613 end do 1614 end do 1615 end if 1616 end do 1617 end do 1618 end if 1619 ! #endif 1620 1621 ABI_FREE(qp_rhog) 1622 ABI_FREE(qp_vhartr) 1623 ABI_FREE(qp_vxc) 1624 ABI_FREE(qp_taug) 1625 ABI_FREE(qp_nhat) 1626 ABI_FREE(qp_nhatgr) 1627 call melements_free(QP_me) 1628 end if ! gwcalctyp<10 1629 1630 ! Free some memory 1631 if (allocated(hbare)) then 1632 ABI_FREE(hbare) 1633 end if 1634 if (allocated(hlda )) then 1635 ABI_FREE(hlda) 1636 end if 1637 1638 ! Prepare the storage of QP amplitudes and energies === 1639 ! Initialize with KS wavefunctions and energies. 1640 Sr%eigvec_qp=czero; Sr%en_qp_diago=zero 1641 do ib=1,Sigp%nbnds 1642 Sr%en_qp_diago(ib,:,:)=KS_BSt%eig(ib,:,:) 1643 Sr%eigvec_qp(ib,ib,:,:)=cone 1644 end do 1645 1646 ! Store <n,k,s|V_xc[n_val]|n,k,s> and <n,k,s|V_U|n,k,s> === 1647 ! Note that we store the matrix elements of V_xc in the KS basis set, not in the QP basis set 1648 ! Matrix elements of V_U are zero unless we are using LDA+U as starting point 1649 do ib=b1gw,b2gw 1650 Sr%vxcme(ib,:,:)=KS_me%vxcval(ib,ib,:,:) 1651 if (Dtset%usepawu>0) Sr%vUme (ib,:,:)=KS_me%vu(ib,ib,:,:) 1652 end do 1653 1654 ! Initial guess for the GW energies 1655 ! Save the energies of the previous iteration. 1656 do spin=1,Sigp%nsppol 1657 do ik=1,Kmesh%nibz 1658 do ib=1,Sigp%nbnds 1659 Sr%e0 (ib,ik,spin)=QP_BSt%eig(ib,ik,spin) 1660 Sr%egw(ib,ik,spin)=QP_BSt%eig(ib,ik,spin) 1661 end do 1662 Sr%e0gap(ik,spin)=zero 1663 ks_iv=ks_vbik(ik,spin) 1664 if (Sigp%nbnds>=ks_iv+1) Sr%e0gap(ik,spin)=Sr%e0(ks_iv+1,ik,spin)-Sr%e0(ks_iv,ik,spin) 1665 end do 1666 end do 1667 ! 1668 !=== If required apply a scissor operator or update the energies === 1669 !TODO check if other Sr entries have to be updated 1670 !moreover this part should be done only in case of semiconductors 1671 !FIXME To me it makes more sense if we apply the scissor to KS_BS but I have to RECHECK csigme 1672 if (ABS(Sigp%mbpt_sciss)>tol6) then 1673 write(msg,'(6a,f10.5,2a)')ch10,& 1674 & ' sigma : performing a first self-consistency',ch10,& 1675 & ' update of the energies in G by a scissor operator',ch10, & 1676 & ' applying a scissor operator of ',Sigp%mbpt_sciss*Ha_eV,' [eV] ',ch10 1677 call wrtout(std_out,msg,'COLL') 1678 call wrtout(ab_out,msg,'COLL') 1679 do spin=1,Sigp%nsppol 1680 do ik=1,Kmesh%nibz 1681 ks_iv=ks_vbik(ik,spin) 1682 if (Sigp%nbnds>=ks_iv+1) then 1683 Sr%egw (ks_iv+1:Sigp%nbnds,ik,spin) = Sr%egw (ks_iv+1:Sigp%nbnds,ik,spin)+Sigp%mbpt_sciss 1684 QP_BSt%eig(ks_iv+1:Sigp%nbnds,ik,spin) = QP_BSt%eig(ks_iv+1:Sigp%nbnds,ik,spin)+Sigp%mbpt_sciss 1685 end if 1686 end do 1687 end do 1688 ! %call apply_scissor(QP_BSt,Sigp%mbpt_sciss) 1689 else if (.FALSE.) then 1690 write(msg,'(4a)')ch10,& 1691 & ' sigma : performing a first self-consistency',ch10,& 1692 & ' update of the energies in G by a previous GW calculation' 1693 call wrtout(std_out,msg,'COLL') 1694 call wrtout(ab_out,msg,'COLL') 1695 ! TODO Recheck this part, is not clear to me! 1696 ABI_MALLOC(igwene,(QP_Bst%mband,QP_Bst%nkpt,QP_Bst%nsppol)) 1697 call rdgw(QP_BSt,'__in.gw__',igwene,extrapolate=.TRUE.) 1698 ABI_FREE(igwene) 1699 Sr%egw=QP_BSt%eig 1700 ! 1701 ! * Recalculate the new fermi level. 1702 call ebands_update_occ(QP_BSt,Dtset%spinmagntarget,prtvol=0) 1703 end if 1704 ! 1705 !In case of AC refer all the energies wrt to the fermi level 1706 !Take care because results from ppmodel cannot be used for AC 1707 !FIXME check ks_energy or qp_energy (in case of SCGW?) 1708 1709 if (mod10==SIG_GW_AC) then 1710 ! All these quantities will be passed to csigme 1711 ! if I skipped the self-consistent part then here I have to use fermi 1712 QP_BSt%eig = QP_BSt%eig -QP_BSt%fermie 1713 Sr%egw = Sr%egw-QP_BSt%fermie 1714 Sr%e0 = Sr%e0 -QP_BSt%fermie 1715 oldefermi=QP_BSt%fermie 1716 ! TODO Recheck fermi 1717 ! Clean EVERYTHING in particulare the treatment of E fermi 1718 QP_BSt%fermie=zero 1719 end if 1720 ! 1721 ! === Setup frequencies around the KS\QP eigenvalues to compute Sigma derivatives (notice the spin) === 1722 ! TODO it is better using an odd Sr%nomega4sd so that the KS\QP eigenvalue is in the middle 1723 ioe0j=Sr%nomega4sd/2+1 1724 do spin=1,Sigp%nsppol 1725 do ikcalc=1,Sigp%nkptgw 1726 ib1=Sigp%minbnd(ikcalc,spin) 1727 ib2=Sigp%maxbnd(ikcalc,spin) 1728 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 1729 do jb=ib1,ib2 1730 do io=1,Sr%nomega4sd 1731 Sr%omega4sd(jb,ik_ibz,io,spin)=Sr%egw(jb,ik_ibz,spin)+Sigp%deltae*(io-ioe0j) 1732 end do 1733 end do 1734 end do 1735 end do 1736 1737 call timab(408,2,tsec) ! hqp_init 1738 call timab(409,1,tsec) ! getW 1739 1740 ! === Get epsilon^{-1} either from the _SCR or the _SUSC file and store it in Er%epsm1 === 1741 ! * If Er%mqmem==0, allocate and read a single q-slice inside csigme. 1742 ! TODO Er%nomega should be initialized so that only the frequencies really needed are stored in memory 1743 1744 ! TODO: The same piece of code is present in screening. 1745 if (sigma_needs_w(Sigp)) then 1746 select case (dtset%gwgamma) 1747 case (0) 1748 id_required=4; ikxc=0; approx_type=0; option_test=0; dim_kxcg=0 1749 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1750 1751 case (1, 2) 1752 ! ALDA TDDFT kernel vertex 1753 !ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available") 1754 ABI_CHECK(Er%ID==0,"Er%ID should be 0") 1755 1756 if (Dtset%usepaw==1) then 1757 ! If we have PAW, we need the full density on the fine grid 1758 ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden)) 1759 if (Dtset%getpawden==0 .and. Dtset%irdpawden==0) then 1760 MSG_ERROR("Must use get/irdpawden to provide a _PAWDEN file!") 1761 end if 1762 call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin),"COLL") 1763 if (.not. file_exists(dtfil%filpawdensin)) then 1764 MSG_ERROR(sjoin("Missing file:", dtfil%filpawdensin)) 1765 end if 1766 1767 ABI_DT_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw)) 1768 call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, & 1769 ks_aepaw_rhor, hdr_rhor, tmp_pawrhoij, wfd%comm) 1770 1771 call hdr_free(hdr_rhor) 1772 call pawrhoij_free(tmp_pawrhoij) 1773 ABI_DT_FREE(tmp_pawrhoij) 1774 end if ! Dtset%usepaw==1 1775 1776 id_required=4; ikxc=7; approx_type=1; dim_kxcg=1 1777 if (Dtset%gwgamma==1) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma 1778 if (Dtset%gwgamma==2) option_test=0 ! TESTPARTICLE, vertex in chi0 only 1779 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1780 1781 dbg_mode=.FALSE. 1782 if (Dtset%usepaw==1) then 1783 ! Use PAW all-electron density 1784 call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,ks_aepaw_rhor,& 1785 Er%npwe,dim_kxcg,kxcg,Gsph_c%gvec,xmpi_comm_self,dbg_mode=dbg_mode) 1786 else 1787 ! Norm-conserving 1788 call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,ks_rhor,& 1789 Er%npwe,dim_kxcg,kxcg,Gsph_c%gvec,xmpi_comm_self,dbg_mode=dbg_mode) 1790 end if 1791 1792 case (3, 4) 1793 ! ADA non-local kernel vertex 1794 !ABI_CHECK(Dtset%usepaw==0,"GWGamma + PAW not available") 1795 ABI_CHECK(Er%ID==0,"Er%ID should be 0") 1796 ABI_CHECK(Sigp%nsppol==1,"ADA vertex for GWGamma not available yet for spin-polarised cases") 1797 1798 if (Dtset%usepaw==1) then ! If we have PAW, we need the full density on the fine grid 1799 ABI_MALLOC(ks_aepaw_rhor,(nfftf,Sigp%nsppol)) 1800 if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then 1801 MSG_ERROR("Must use get/irdpawden to provide a _PAWDEN file!") 1802 end if 1803 call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin),"COLL") 1804 if (.not. file_exists(dtfil%filpawdensin)) then 1805 MSG_ERROR(sjoin("Missing file:", dtfil%filpawdensin)) 1806 end if 1807 1808 ABI_DT_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw)) 1809 1810 call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, & 1811 ks_aepaw_rhor, hdr_rhor, tmp_pawrhoij, wfd%comm) 1812 1813 call hdr_free(hdr_rhor) 1814 call pawrhoij_free(tmp_pawrhoij) 1815 ABI_DT_FREE(tmp_pawrhoij) 1816 end if ! Dtset%usepaw==1 1817 1818 id_required=4; ikxc=7; approx_type=2; 1819 if (Dtset%gwgamma==3) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma 1820 if (Dtset%gwgamma==4) option_test=0 ! TESTPARTICLE, vertex in chi0 only 1821 ABI_MALLOC(fxc_ADA,(Er%npwe,Er%npwe,Er%nqibz)) 1822 ! Use userrd to set kappa 1823 if (Dtset%userrd==zero) Dtset%userrd = 2.1_dp 1824 ! Set correct value of kappa (should be scaled with alpha*r_s where) 1825 ! r_s is Wigner-Seitz radius and alpha=(4/(9*Pi))^(1/3) 1826 rhoav = (drude_plsmf*drude_plsmf)/four_pi 1827 r_s = (three/(four_pi*rhoav))**third 1828 alpha = (four*ninth*piinv)**third 1829 Dtset%userrd = Dtset%userrd/(alpha*r_s) 1830 1831 dbg_mode=.TRUE. 1832 if (Dtset%usepaw==1) then 1833 ! Use PAW all-electron density 1834 call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,& 1835 ks_aepaw_rhor,Er%npwe,Er%nqibz,Er%qibz,& 1836 fxc_ADA,Gsph_c%gvec,xmpi_comm_self,kappa_init=Dtset%userrd,dbg_mode=dbg_mode) 1837 else 1838 ! Norm conserving 1839 call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,& 1840 ks_rhor,Er%npwe,Er%nqibz,Er%qibz,& 1841 fxc_ADA,Gsph_c%gvec,xmpi_comm_self,kappa_init=Dtset%userrd,dbg_mode=dbg_mode) 1842 end if 1843 1844 dim_kxcg = 0 1845 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1846 1847 ! @WC: bootstrap -- 1848 case (-3, -4, -5, -6, -7, -8) 1849 ! ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available") 1850 ABI_CHECK(Er%ID==0,"Er%ID should be 0") 1851 1852 if (Dtset%usepaw==1) then 1853 ! If we have PAW, we need the full density on the fine grid 1854 ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden)) 1855 if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then 1856 MSG_ERROR("Must use get/irdpawden to provide a _PAWDEN file!") 1857 end if 1858 call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin),"COLL") 1859 if (.not. file_exists(dtfil%filpawdensin)) then 1860 MSG_ERROR(sjoin("Missing file:", dtfil%filpawdensin)) 1861 end if 1862 1863 ABI_DT_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw)) 1864 1865 call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, & 1866 ks_aepaw_rhor, hdr_rhor, tmp_pawrhoij, wfd%comm) 1867 1868 call hdr_free(hdr_rhor) 1869 call pawrhoij_free(tmp_pawrhoij) 1870 ABI_DT_FREE(tmp_pawrhoij) 1871 end if ! Dtset%usepaw==1 1872 1873 id_required=4; ikxc=7; dim_kxcg=0 1874 1875 if (dtset%gwgamma>-5) then 1876 approx_type=4 ! full fxc(G,G') 1877 else if (dtset%gwgamma>-7) then 1878 approx_type=5 ! fxc(0,0) one-shot 1879 else 1880 approx_type=6 ! rpa-type bootstrap 1881 end if 1882 1883 option_test=MOD(Dtset%gwgamma,2) 1884 ! 1 -> TESTELECTRON, vertex in chi0 *and* sigma 1885 ! 0 -> TESTPARTICLE, vertex in chi0 only 1886 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1887 ! --@WC 1888 1889 case default 1890 MSG_ERROR(sjoin("Wrong gwgamma:", itoa(dtset%gwgamma))) 1891 end select 1892 1893 ! Set plasma frequency 1894 my_plsmf = drude_plsmf; if (Dtset%ppmfrq>tol6) my_plsmf = Dtset%ppmfrq 1895 Dtset%ppmfrq = my_plsmf 1896 1897 ! if(dtset%ucrpa==0) then 1898 if (Dtset%gwgamma<3) then 1899 call mkdump_Er(Er,Vcp,Er%npwe,Gsph_c%gvec,dim_kxcg,kxcg,id_required,& 1900 & approx_type,ikxc,option_test,Dtfil%fnameabo_scr,Dtset%iomode,& 1901 & nfftf_tot,ngfftf,comm) 1902 else 1903 call mkdump_Er(Er,Vcp,Er%npwe,Gsph_c%gvec,dim_kxcg,kxcg,id_required,& 1904 & approx_type,ikxc,option_test,Dtfil%fnameabo_scr,Dtset%iomode,& 1905 & nfftf_tot,ngfftf,comm,fxc_ADA) 1906 end if 1907 ! end if 1908 if (allocated(kxcg)) then 1909 ABI_FREE(kxcg) 1910 end if 1911 if (allocated(fxc_ADA)) then 1912 ABI_FREE(fxc_ADA) 1913 end if 1914 if (allocated(ks_aepaw_rhor)) then 1915 ABI_FREE(ks_aepaw_rhor) 1916 end if 1917 end if 1918 ! 1919 ! ================================================ 1920 ! ==== Calculate plasmonpole model parameters ==== 1921 ! ================================================ 1922 ! TODO Maybe its better if we use mqmem as input variable 1923 use_aerhor=0 1924 ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden*use_aerhor)) 1925 1926 if (sigma_needs_ppm(Sigp)) then 1927 my_plsmf=drude_plsmf; if (Dtset%ppmfrq>tol6) my_plsmf=Dtset%ppmfrq 1928 call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,Sigp%ppmodel,my_plsmf,Dtset%gw_invalid_freq) 1929 1930 ! PPm%force_plsmf= force_ppmfrq ! this line to change the plasme frequency in HL expression. 1931 1932 if (Wfd%usepaw==1 .and. Ppm%userho==1) then 1933 ! * For PAW and ppmodel 2-3-4 we need the AE rho(G) without compensation charge. 1934 ! * It would be possible to calculate rho(G) using Paw_pwff, though. It should be faster but 1935 ! results will depend on the expression used for the matrix elements. This approach is safer. 1936 use_aerhor=1 1937 ABI_FREE(ks_aepaw_rhor) 1938 ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden)) 1939 1940 ! Check if input density file is available, otherwise compute 1941 pawden_fname = strcat(Dtfil%filnam_ds(3), '_PAWDEN') 1942 call wrtout(std_out,sjoin('Checking for existence of:',pawden_fname)) 1943 if (file_exists(pawden_fname)) then 1944 ! Read density from file 1945 ABI_DT_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw)) 1946 1947 call read_rhor(pawden_fname, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, & 1948 ks_aepaw_rhor, hdr_rhor, tmp_pawrhoij, wfd%comm) 1949 1950 call hdr_free(hdr_rhor) 1951 call pawrhoij_free(tmp_pawrhoij) 1952 ABI_DT_FREE(tmp_pawrhoij) 1953 else 1954 ! Have to calculate PAW AW rhor from scratch 1955 ABI_MALLOC(qp_rhor_n_one,(pawfgr%nfft,Dtset%nspden)) 1956 ABI_MALLOC(qp_rhor_nt_one,(pawfgr%nfft,Dtset%nspden)) 1957 1958 ! FIXME 1959 MSG_WARNING(" denfgr in sigma seems to produce wrong results") 1960 write(std_out,*)" input tilde ks_rhor integrates: ",SUM(ks_rhor(:,1))*Cryst%ucvol/nfftf 1961 1962 call denfgr(Cryst%atindx1,Cryst%gmet,comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,ks_nhat,Wfd%nspinor,& 1963 Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,KS_Pawrhoij,Pawtab,Dtset%prtvol, & 1964 ks_rhor,ks_aepaw_rhor,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 1965 1966 ABI_FREE(qp_rhor_n_one) 1967 ABI_FREE(qp_rhor_nt_one) 1968 end if 1969 1970 write(msg,'(a,f8.4)')' sigma: PAW AE density used for PPmodel integrates to: ',SUM(ks_aepaw_rhor(:,1))*Cryst%ucvol/nfftf 1971 call wrtout(std_out,msg,'PERS') 1972 1973 if (Er%mqmem/=0) then 1974 ! Calculate ppmodel parameters for all q-points. 1975 call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,nfftf,Gsph_c%gvec,ngfftf,ks_aepaw_rhor(:,1)) 1976 end if 1977 1978 else 1979 ! NC or PAW with PPmodel 1. 1980 if (Er%mqmem/=0) then 1981 ! Calculate ppmodel parameters for all q-points 1982 call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,nfftf,Gsph_c%gvec,ngfftf,ks_rhor(:,1)) 1983 end if 1984 end if ! PAW or NC PPm and/or needs density 1985 end if ! sigma_needs_ppm 1986 1987 call timab(409,2,tsec) ! getW 1988 1989 if (wfd_iam_master(Wfd)) then 1990 ! Write info on the run on ab_out, then open files to store final results. 1991 call ebands_report_gap(KS_BSt,header='KS Band Gaps',unit=ab_out) 1992 if(dtset%ucrpa==0) then 1993 call write_sigma_header(Sigp,Er,Cryst,Kmesh,Qmesh) 1994 end if 1995 1996 if (open_file(Dtfil%fnameabo_gw,msg,unit=unt_gw,status='unknown',form='formatted') /=0) then 1997 MSG_ERROR(msg) 1998 end if 1999 write(unt_gw,*)Sigp%nkptgw,Sigp%nsppol 2000 2001 if (open_file(Dtfil%fnameabo_gwdiag,msg,unit=unt_gwdiag,status='unknown',form='formatted') /= 0) then 2002 MSG_ERROR(msg) 2003 end if 2004 write(unt_gwdiag,*)Sigp%nkptgw,Sigp%nsppol 2005 2006 if (open_file(Dtfil%fnameabo_sig,msg,unit=unt_sig,status='unknown',form='formatted') /= 0) then 2007 MSG_ERROR(msg) 2008 end if 2009 if (open_file(Dtfil%fnameabo_sgr,msg,unit=unt_sgr,status='unknown',form='formatted') /= 0) then 2010 MSG_ERROR(msg) 2011 end if 2012 2013 if (mod10==SIG_GW_AC) then 2014 ! Sigma along the imaginary axis. 2015 if (open_file(Dtfil%fnameabo_sgm,msg,unit=unt_sgm,status='unknown',form='formatted') /= 0) then 2016 MSG_ERROR(msg) 2017 end if 2018 end if 2019 end if 2020 2021 !======================================================================= 2022 !==== Calculate self-energy and output the results for each k-point ==== 2023 !======================================================================= 2024 !* Here it would be possible to calculate the QP correction for the same k-point using a PPmodel 2025 !in the first iteration just to improve the initial guess and CD or AC in the second step. Useful 2026 !if the KS level is far from the expected QP results. Previously it was possible to do such 2027 !calculation by simply listing the same k-point twice in kptgw. Now this trick is not allowed anymore. 2028 !Everything, indeed, should be done in a clean and transparent way inside csigme. 2029 2030 call wfd_print(Wfd,mode_paral='PERS') 2031 2032 call wrtout(std_out,sigma_type_from_key(mod10),'COLL') 2033 2034 if (gwcalctyp<10) then 2035 msg = " Perturbative Calculation" 2036 if (gwcalctyp==1) msg = " Newton Raphson method " 2037 else if (gwcalctyp<20) then 2038 msg = " Self-Consistent on Energies only" 2039 else 2040 msg = " Self-Consistent on Energies and Wavefunctions" 2041 end if 2042 if(Dtset%ucrpa==0) then 2043 call wrtout(std_out,msg,'COLL') 2044 end if 2045 ! 2046 !================================================= 2047 !==== Calculate the matrix elements of Sigma ===== 2048 !================================================= 2049 2050 nomega_sigc=Sr%nomega_r+Sr%nomega4sd; if (mod10==SIG_GW_AC) nomega_sigc=Sr%nomega_i 2051 2052 ! min and max band indeces for GW corrections. 2053 ib1=Sigp%minbdgw; ib2=Sigp%maxbdgw 2054 2055 !MG TODO: I don't like the fact that ib1 and ib2 are redefined here because this 2056 ! prevents me from refactoring the code. In particular I want to store the self-energy 2057 ! results inside the sigma_results datatypes hence one needs to know all the dimensions 2058 ! at the beginning of the execution (e.g. in setup_sigma) so that one can easily allocate the arrays in the type. 2059 if(Dtset%ucrpa>=1) then 2060 !Read the band 2061 if (open_file("forlb.ovlp",msg,newunit=temp_unt,form="formatted", status="unknown") /= 0) then 2062 MSG_ERROR(msg) 2063 end if 2064 rewind(temp_unt) 2065 read(temp_unt,*) 2066 read(temp_unt,*) msg, ib1, ib2 2067 close(temp_unt) 2068 end if 2069 2070 ABI_MALLOC(sigcme,(nomega_sigc,ib1:ib2,ib1:ib2,Sigp%nkptgw,Sigp%nsppol*Sigp%nsig_ab)) 2071 sigcme=czero 2072 2073 ! if (.False. .and. psps%usepaw == 0 .and. wfd%nspinor == 1 .and. any(dtset%so_psp /= 0)) then 2074 ! call wrtout(std_out, "Computing SOC contribution with first-order perturbation theory") 2075 ! ABI_MALLOC(bks_mask, (wfd%mband, wfd%nkibz, wfd%nsppol)) 2076 ! bks_mask = .False. 2077 ! do spin=1,wfd%nsppol 2078 ! do ikcalc=1,Sigp%nkptgw 2079 ! ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW 2080 ! ii=Sigp%minbnd(ikcalc, spin); jj=Sigp%maxbnd(ikcalc, spin) 2081 ! bks_mask(ii:jj, ik_ibz, spin) = .True. 2082 ! end do 2083 ! end do 2084 ! 2085 ! call wfd_get_socpert(wfd, cryst, psps, pawtab, bks_mask, osoc_bks) 2086 ! ABI_FREE(bks_mask) 2087 ! ABI_FREE(osoc_bks) 2088 ! end if 2089 2090 !========================================================== 2091 !==== Exchange part using the dense gwx_ngfft FFT mesh ==== 2092 !========================================================== 2093 !TODO : load distribution is not optimal if band parallelism is used. 2094 !but this problem was also affecting the old implementation. 2095 call timab(421,1,tsec) ! calc_sigx_me 2096 2097 ! This routine shows how the wavefunctions are distributed. 2098 !call wfd_show_bkstab(Wfd,unit=std_out) 2099 2100 if(Dtset%ucrpa>=1) then 2101 ! !Calculation of the Rho_n for calc_ucrpa 2102 call wrtout(std_out,sjoin("begin of Ucrpa calc for a nb of kpoints of: ",itoa(Sigp%nkptgw)),"COLL") 2103 ! Wannier basis: rhot1_q_m will need an atom index in the near future. 2104 ic=0 2105 do itypat=1,cryst%ntypat 2106 if(pawtab(itypat)%lpawu.ne.-1) then 2107 ndim=2*dtset%lpawu(itypat)+1 2108 ic=ic+1 2109 itypatcor=itypat 2110 lcor=dtset%lpawu(itypat) 2111 end if 2112 end do 2113 if(ic>1) then 2114 MSG_ERROR("number of correlated species is larger than one") 2115 end if 2116 ABI_ALLOCATE(rhot1_q_m,(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,ndim,ndim,sigp%npwx,Qmesh%nibz)) 2117 ABI_ALLOCATE(M1_q_m,(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,ndim,ndim,sigp%npwx,Qmesh%nibz)) 2118 2119 M1_q_m=czero 2120 rhot1_q_m=czero 2121 2122 ! do ikcalc=1,Sigp%nkptgw 2123 2124 ! if(cryst%nsym==1) nkcalc=Kmesh%nbz 2125 ! if(cryst%nsym>1) nkcalc=Kmesh%nibz 2126 nkcalc=Kmesh%nbz 2127 do ikcalc=1,nkcalc ! for the oscillator strengh, spins are identical without SOC 2128 ! if(Sigp%nkptgw/=Kmesh%nbz) then 2129 ! write(msg,'(6a)')ch10,& 2130 !& ' nkptgw and nbz differs: this is not allowed to compute U in cRPA ' 2131 ! MSG_ERROR(msg) 2132 ! endif 2133 write(std_out,*) "ikcalc",ikcalc,size(Sigp%kptgw2bz),nkcalc 2134 2135 if(mod(ikcalc-1,nprocs)==Wfd%my_rank) then 2136 if(cryst%nsym==1) ik_ibz=Kmesh%tab(ikcalc) ! Index of the irred k-point for GW 2137 if(cryst%nsym>1) ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW 2138 2139 call prep_calc_ucrpa(ik_ibz,ikcalc,itypatcor,ib1,ib2,Cryst,QP_bst,Sigp,Gsph_x,Vcp,& 2140 & Kmesh,Qmesh,lcor,M1_q_m,Pawtab,Pawang,Paw_pwff,& 2141 & Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,gwx_ngfft,& 2142 & ngfftf,Dtset%prtvol,Dtset%pawcross,rhot1_q_m) 2143 end if 2144 end do 2145 2146 call xmpi_sum(rhot1_q_m,Wfd%comm,ierr) 2147 call xmpi_sum(M1_q_m,Wfd%comm,ierr) 2148 ! if(Cryst%nsym==1) then 2149 M1_q_m=M1_q_m/Kmesh%nbz/Wfd%nsppol 2150 rhot1_q_m=rhot1_q_m/Kmesh%nbz/Wfd%nsppol 2151 ! endif 2152 ! Calculation of U in cRPA: need to treat with a different cutoff the 2153 ! bare coulomb and the screened coulomb interaction. 2154 call calc_ucrpa(itypatcor,cryst,Kmesh,lcor,M1_q_m,Qmesh,Er%npwe,sigp%npwx,& 2155 & Cryst%nsym,rhot1_q_m,Sigp%nomegasr,Sigp%minomega_r,Sigp%maxomega_r,ib1,ib2,'Gsum',Cryst%ucvol,Wfd,Er%fname) 2156 2157 ABI_DEALLOCATE(rhot1_q_m) 2158 ABI_DEALLOCATE(M1_q_m) 2159 2160 else 2161 do ikcalc=1,Sigp%nkptgw 2162 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW 2163 ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point) 2164 ib2=MAXVAL(Sigp%maxbnd(ikcalc,:)) 2165 2166 call calc_sigx_me(ik_ibz,ikcalc,ib1,ib2,Cryst,QP_bst,Sigp,Sr,Gsph_x,Vcp,Kmesh,Qmesh,Ltg_k(ikcalc),& 2167 & Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,& 2168 & gwx_ngfft,ngfftf,Dtset%prtvol,Dtset%pawcross) 2169 end do 2170 2171 ! for the time being, do not remove this barrier! 2172 call xmpi_barrier(Wfd%comm) 2173 2174 call timab(421,2,tsec) ! calc_sigx_me 2175 2176 ! ========================================================== 2177 ! ==== Correlation part using the coarse gwc_ngfft mesh ==== 2178 ! ========================================================== 2179 if (mod10/=SIG_HF) then 2180 do ikcalc=1,Sigp%nkptgw 2181 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW 2182 ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point) 2183 ib2=MAXVAL(Sigp%maxbnd(ikcalc,:)) 2184 2185 ! sigcme_p => sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:) ! Causes annoying Fortran runtime warning on abiref. 2186 ABI_ALLOCATE(sigcme_k,(nomega_sigc,ib2-ib1+1,ib2-ib1+1,Sigp%nsppol*Sigp%nsig_ab)) 2187 sigcme_k=czero 2188 if (any(mod10 == [SIG_SEX, SIG_COHSEX])) then 2189 ! Calculate static COHSEX or SEX using the coarse gwc_ngfft mesh. 2190 call cohsex_me(ik_ibz,ikcalc,nomega_sigc,ib1,ib2,Cryst,QP_BSt,Sigp,Sr,Er,Gsph_c,Vcp,Kmesh,Qmesh,& 2191 & Ltg_k(ikcalc),Pawtab,Pawang,Paw_pwff,Psps,Wfd,QP_sym,& 2192 !& gwc_ngfft,Dtset%iomode,Dtset%prtvol,sigcme_p) 2193 & gwc_ngfft,Dtset%iomode,Dtset%prtvol,sigcme_k) 2194 else 2195 ! Compute correlated part using the coarse gwc_ngfft mesh. 2196 call calc_sigc_me(ik_ibz,ikcalc,nomega_sigc,ib1,ib2,Dtset,Cryst,QP_BSt,Sigp,Sr,Er,Gsph_Max,Gsph_c,Vcp,Kmesh,Qmesh,& 2197 & Ltg_k(ikcalc),PPm,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,& 2198 !& gwc_ngfft,ngfftf,nfftf,ks_rhor,use_aerhor,ks_aepaw_rhor,sigcme_p) 2199 & gwc_ngfft,ngfftf,nfftf,ks_rhor,use_aerhor,ks_aepaw_rhor,sigcme_k) 2200 end if 2201 sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:)=sigcme_k 2202 ABI_DEALLOCATE(sigcme_k) 2203 2204 end do 2205 end if 2206 2207 call xmpi_barrier(Wfd%comm) 2208 2209 ! ===================================================== 2210 ! ==== Solve Dyson equation storing results in Sr% ==== 2211 ! ===================================================== 2212 ! * Use perturbative approach or AC to find QP corrections. 2213 ! * If qp-GW, diagonalize also H0+Sigma in the KS basis set to get the 2214 ! new QP amplitudes and energies (Sr%eigvec_qp and Sr%en_qp_diago. 2215 ! TODO AC with spinor not implemented yet. 2216 ! TODO Diagonalization of Sigma+hhartre with AC is wrong. 2217 ! 2218 call timab(425,1,tsec) ! solve_dyson 2219 do ikcalc=1,Sigp%nkptgw 2220 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW 2221 ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indeces for GW corrections (for this k-point) 2222 ib2=MAXVAL(Sigp%maxbnd(ikcalc,:)) 2223 2224 ! sigcme_p => sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:) ! Causes annoying Fortran runtime warning on abiref. 2225 ABI_ALLOCATE(sigcme_k,(nomega_sigc,ib2-ib1+1,ib2-ib1+1,Sigp%nsppol*Sigp%nsig_ab)) 2226 sigcme_k=sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:) 2227 2228 ! call solve_dyson(ikcalc,ib1,ib2,nomega_sigc,Sigp,Kmesh,sigcme_p,QP_BSt%eig,Sr,Dtset%prtvol,Dtfil,Wfd%comm) 2229 call solve_dyson(ikcalc,ib1,ib2,nomega_sigc,Sigp,Kmesh,sigcme_k,QP_BSt%eig,Sr,Dtset%prtvol,Dtfil,Wfd%comm) 2230 ABI_DEALLOCATE(sigcme_k) 2231 ! 2232 ! Calculate direct gap for each spin and print out final results. 2233 ! We use the valence index of the KS system because we still do not know 2234 ! all the QP corrections. Ideally one should use the QP valence index 2235 do spin=1,Sigp%nsppol 2236 if ( Sigp%maxbnd(ikcalc,spin) >= ks_vbik(ik_ibz,spin)+1 .and. & 2237 & Sigp%minbnd(ikcalc,spin) <= ks_vbik(ik_ibz,spin) ) then 2238 ks_iv=ks_vbik(ik_ibz,spin) 2239 Sr%egwgap (ik_ibz,spin)= Sr%egw(ks_iv+1,ik_ibz,spin) - Sr%egw(ks_iv,ik_ibz,spin) 2240 Sr%degwgap(ik_ibz,spin)= Sr%degw(ks_iv+1,ik_ibz,spin) - Sr%degw(ks_iv,ik_ibz,spin) 2241 else 2242 ! The "gap" cannot be computed 2243 Sr%e0gap(ik_ibz,spin)=zero; Sr%egwgap (ik_ibz,spin)=zero; Sr%degwgap(ik_ibz,spin)=zero 2244 end if 2245 end do 2246 2247 if (wfd_iam_master(Wfd)) call write_sigma_results(ikcalc,ik_ibz,Sigp,Sr,KS_BSt) 2248 end do !ikcalc 2249 2250 call timab(425,2,tsec) ! solve_dyson 2251 call timab(426,1,tsec) ! finalize 2252 2253 ! Update the energies in QP_BSt 2254 ! If QPSCGW, use diagonalized eigenvalues otherwise perturbative results. 2255 if (gwcalctyp>=10) then 2256 do ib=1,Sigp%nbnds 2257 QP_BSt%eig(ib,:,:)=Sr%en_qp_diago(ib,:,:) 2258 end do 2259 else 2260 QP_BSt%eig=Sr%egw 2261 end if 2262 2263 ! ================================================================================ 2264 ! ==== This part is done only if all k-points in the IBZ have been calculated ==== 2265 ! ================================================================================ 2266 if (Sigp%nkptgw==Kmesh%nibz) then 2267 2268 ! Recalculate new occupations and Fermi level. 2269 call ebands_update_occ(QP_BSt,Dtset%spinmagntarget,prtvol=Dtset%prtvol) 2270 qp_vbik(:,:) = get_valence_idx(QP_BSt) 2271 2272 write(msg,'(2a,3x,2(es16.6,a))')ch10,' New Fermi energy : ',QP_BSt%fermie,' Ha ,',QP_BSt%fermie*Ha_eV,' eV' 2273 call wrtout(std_out,msg,'COLL') 2274 call wrtout(ab_out,msg,'COLL') 2275 2276 ! === If all k-points and all occupied bands are calculated, output EXX === 2277 ! FIXME here be careful about the check on ks_vbik in case of metals 2278 ! if (my_rank==master.and.Sigp%nkptgw==Kmesh%nibz.and.ALL(Sigp%minbnd(:)==1).and.ALL(Sigp%maxbnd(:)>=MAXVAL(nbv(:)))) then 2279 ! if (ALL(Sigp%minbnd==1).and. ALL(Sigp%maxbnd>=MAXVAL(MAXVAL(ks_vbik(:,:),DIM=1))) ) then 2280 if (ALL(Sigp%minbnd==1).and. ALL(Sigp%maxbnd>=ks_vbik) ) then ! FIXME here the two arrays use a different indexing. 2281 2282 ex_energy = sigma_get_exene(Sr,Kmesh,QP_BSt) 2283 write(msg,'(a,2(es16.6,a))')' New Exchange energy : ',ex_energy,' Ha ,',ex_energy*Ha_eV,' eV' 2284 call wrtout(std_out,msg,'COLL') 2285 call wrtout(ab_out,msg,'COLL') 2286 end if 2287 2288 ! Report the QP gaps (Fundamental and Optical) 2289 call ebands_report_gap(QP_BSt,header='QP Band Gaps',unit=ab_out) 2290 2291 ! Band structure interpolation from QP energies computed on the k-mesh. 2292 if (nint(dtset%einterp(1)) /= 0 .and. all(sigp%minbdgw == sigp%minbnd) .and. all(sigp%maxbdgw == sigp%maxbnd)) then 2293 call ebands_interpolate_kpath(QP_BSt, dtset, cryst, [sigp%minbdgw, sigp%maxbdgw], dtfil%filnam_ds(4), comm) 2294 end if 2295 end if ! Sigp%nkptgw==Kmesh%nibz 2296 ! 2297 ! Write SCF data in case of self-consistent calculation === 2298 ! Save Sr%en_qp_diago, Sr%eigvec_qp and m_lda_to_qp in the _QPS file. 2299 ! Note that in the first iteration qp_rhor contains KS rhor, then the mixed rhor. 2300 if (gwcalctyp>=10) then 2301 ! Calculate the new m_lda_to_qp 2302 call updt_m_lda_to_qp(Sigp,Kmesh,nscf,Sr,Sr%m_lda_to_qp) 2303 2304 if (wfd_iam_master(Wfd)) then 2305 call wrqps(Dtfil%fnameabo_qps,Sigp,Cryst,Kmesh,Psps,Pawtab,QP_Pawrhoij,& 2306 & Dtset%nspden,nscf,nfftf,ngfftf,Sr,QP_BSt,Sr%m_lda_to_qp,qp_rhor) 2307 end if 2308 2309 ! Report the MAX variation for each kptgw and spin 2310 call wrtout(ab_out,ch10//' Convergence of QP corrections ','COLL') 2311 converged=.TRUE. 2312 do spin=1,Sigp%nsppol 2313 write(msg,'(a,i2,a)')' >>>>> For spin ',spin,' <<<<< ' 2314 call wrtout(ab_out,msg,'COLL') 2315 do ikcalc=1,Sigp%nkptgw 2316 ib1 = Sigp%minbnd(ikcalc,spin); ib2 = Sigp%maxbnd(ikcalc,spin) 2317 ik_bz = Sigp%kptgw2bz(ikcalc); ik_ibz = Kmesh%tab(ik_bz) 2318 ii = imax_loc(ABS(Sr%degw(ib1:ib2,ik_ibz,spin))) 2319 max_degw = Sr%degw(ii,ik_ibz,spin) 2320 write(msg,('(a,i3,a,2f8.3,a,i3)'))'. kptgw no:',ikcalc,'; Maximum DeltaE = (',max_degw*Ha_eV,') for band index:',ii 2321 call wrtout(ab_out,msg,'COLL') 2322 converged = (converged .and. ABS(max_degw) < Dtset%gw_toldfeig) 2323 end do 2324 end do 2325 end if 2326 2327 ! ============================================ 2328 ! ==== Save the GW results in NETCDF file ==== 2329 ! ============================================ 2330 #ifdef HAVE_NETCDF 2331 if (wfd_iam_master(Wfd)) then 2332 NCF_CHECK(nctk_open_create(ncid, strcat(dtfil%filnam_ds(4), '_SIGRES.nc'), xmpi_comm_self)) 2333 NCF_CHECK(nctk_defnwrite_ivars(ncid, ["sigres_version"], [1])) 2334 NCF_CHECK(crystal_ncwrite(Cryst, ncid)) 2335 NCF_CHECK(ebands_ncwrite(KS_Bst, ncid)) 2336 NCF_CHECK(sigma_ncwrite(Sigp, Er, Sr, ncid)) 2337 ! Add qp_rhor. Note that qp_rhor == ks_rhor if wavefunctions are not updated. 2338 !ncerr = nctk_write_datar("qp_rhor",path,ngfft,cplex,nfft,nspden,& 2339 ! comm_fft,fftn3_distrib,ffti3_local,datar,action) 2340 NCF_CHECK(nf90_close(ncid)) 2341 end if 2342 #endif 2343 2344 end if ! ucrpa 2345 2346 !----------------------------- END OF THE CALCULATION ------------------------ 2347 ! 2348 !===================== 2349 !==== Close Files ==== 2350 !===================== 2351 if (wfd_iam_master(Wfd)) then 2352 close(unt_gw ) 2353 close(unt_gwdiag) 2354 close(unt_sig) 2355 close(unt_sgr) 2356 if (mod10==SIG_GW_AC) close(unt_sgm) 2357 end if 2358 ! 2359 !=============================================== 2360 !==== Free arrays and local data structures ==== 2361 !=============================================== 2362 ABI_FREE(ks_vbik) 2363 ABI_FREE(qp_vbik) 2364 ABI_FREE(ph1d) 2365 ABI_FREE(ph1df) 2366 ABI_FREE(qp_rhor) 2367 ABI_FREE(ks_rhor) 2368 ABI_FREE(ks_rhog) 2369 ABI_FREE(qp_taur) 2370 ABI_FREE(ks_taur) 2371 ABI_FREE(ks_taug) 2372 ABI_FREE(ks_vhartr) 2373 ABI_FREE(ks_vtrial) 2374 ABI_FREE(vpsp) 2375 ABI_FREE(ks_vxc) 2376 ABI_FREE(xccc3d) 2377 ABI_FREE(grchempottn) 2378 ABI_FREE(grewtn) 2379 ABI_FREE(grvdw) 2380 ABI_FREE(sigcme) 2381 2382 if (allocated(kxc)) then 2383 ABI_FREE(kxc) 2384 end if 2385 if (allocated(qp_vtrial)) then 2386 ABI_FREE(qp_vtrial) 2387 end if 2388 2389 ABI_FREE(ks_nhat) 2390 ABI_FREE(ks_nhatgr) 2391 ABI_FREE(dijexc_core) 2392 ABI_FREE(ks_aepaw_rhor) 2393 call pawfgr_destroy(Pawfgr) 2394 2395 ! Deallocation for PAW. 2396 if (Dtset%usepaw==1) then 2397 call pawrhoij_free(KS_Pawrhoij) 2398 ABI_DT_FREE(KS_Pawrhoij) 2399 call pawfgrtab_free(Pawfgrtab) 2400 call paw_ij_free(KS_paw_ij) 2401 call paw_an_free(KS_paw_an) 2402 call pawpwff_free(Paw_pwff) 2403 if (gwcalctyp>=10) then 2404 call pawrhoij_free(QP_pawrhoij) 2405 call paw_ij_free(QP_paw_ij) 2406 call paw_an_free(QP_paw_an) 2407 end if 2408 if (Dtset%pawcross==1) then 2409 call paw_pwaves_lmn_free(Paw_onsite) 2410 Wfdf%bks_comm = xmpi_comm_null 2411 call wfd_free(Wfdf) 2412 end if 2413 end if 2414 ABI_DT_FREE(Paw_onsite) 2415 ABI_DT_FREE(Pawfgrtab) 2416 ABI_DT_FREE(Paw_pwff) 2417 ABI_DT_FREE(KS_paw_ij) 2418 ABI_DT_FREE(KS_paw_an) 2419 if (gwcalctyp>=10) then 2420 ABI_DT_FREE(QP_pawrhoij) 2421 ABI_DT_FREE(QP_paw_an) 2422 ABI_DT_FREE(QP_paw_ij) 2423 end if 2424 2425 call wfd_free(Wfd) 2426 call destroy_mpi_enreg(MPI_enreg_seq) 2427 call littlegroup_free(Ltg_k) 2428 ABI_DT_FREE(Ltg_k) 2429 call kmesh_free(Kmesh) 2430 call kmesh_free(Qmesh) 2431 call gsph_free(Gsph_Max) 2432 call gsph_free(Gsph_x) 2433 call gsph_free(Gsph_c) 2434 call vcoul_free(Vcp) 2435 call crystal_free(Cryst) 2436 call sigma_free(Sr) 2437 call em1results_free(Er) 2438 call ppm_free(PPm) 2439 call hdr_free(Hdr_sigma) 2440 call hdr_free(Hdr_wfk) 2441 call ebands_free(KS_BSt) 2442 call ebands_free(QP_BSt) 2443 call melements_free(KS_me) 2444 call esymm_free(KS_sym) 2445 ABI_DT_FREE(KS_sym) 2446 2447 if (Sigp%symsigma==1.and.gwcalctyp>=20) then 2448 call esymm_free(QP_sym) 2449 ABI_DT_FREE(QP_sym) 2450 end if 2451 2452 call sigparams_free(Sigp) 2453 2454 call timab(426,2,tsec) ! finalize 2455 call timab(401,2,tsec) 2456 2457 DBG_EXIT('COLL') 2458 2459 end subroutine sigma