TABLE OF CONTENTS
ABINIT/screening [ Functions ]
NAME
screening
FUNCTION
Calculate screening and dielectric functions
COPYRIGHT
Copyright (C) 2001-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<datafiles_type)>=variables related to file names and unit numbers. Pawang<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 screening, 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
Output is written on the main output file. The symmetrical inverse dielectric matrix is stored in the _SCR file
SIDE EFFECTS
Dtset<type(dataset_type)>=all input variables for this dataset
NOTES
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.
PARENTS
driver
CHILDREN
apply_scissor,calc_rpa_functional,cchi0,cchi0q0,chi0_bksmask chi0q0_intraband,chi_free,chkpawovlp,coeffs_gausslegint,crystal_free destroy_mpi_enreg,ebands_copy,ebands_free,ebands_update_occ em1params_free,energies_init,fourdp,get_gftt,getph,gsph_free,hdr_free hscr_free,hscr_io,init_distribfft_seq,initmpi_seq,kmesh_free,kxc_ada kxc_driver,littlegroup_free,lwl_write,make_epsm1_driver,metric,mkrdim nhatgrid,output_chi0sumrule,paw_an_free,paw_an_init,paw_an_nullify paw_gencond,paw_ij_free,paw_ij_init,paw_ij_nullify,paw_pwaves_lmn_free paw_pwaves_lmn_init,pawdenpot,pawdij,pawfgr_destroy,pawfgr_init pawfgrtab_free,pawfgrtab_init,pawinit,pawmknhat,pawnabla_init,pawprt pawpuxinit,pawpwff_free,pawpwff_init,pawrhoij_alloc,pawrhoij_copy pawrhoij_free,pawtab_get_lsize,pawtab_print,print_arr,print_ngfft prtrhomxmn,pspini,random_stopping_power,rdgw,rdqps,rotate_fft_mesh setsymrhoij,setup_screening,setvtr,spectra_free,spectra_repr spectra_write,symdij,symdij_all,test_charge,timab,vcoul_free wfd_change_ngfft,wfd_copy,wfd_free,wfd_init,wfd_mkrho,wfd_print wfd_read_wfk,wfd_rotate,wfd_test_ortho,write_screening,wrtout xmpi_bcast
SOURCE
79 #if defined HAVE_CONFIG_H 80 #include "config.h" 81 #endif 82 83 #include "abi_common.h" 84 85 subroutine screening(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim) 86 87 use defs_basis 88 use defs_datatypes 89 use defs_abitypes 90 use defs_wvltypes 91 use m_profiling_abi 92 use m_xmpi 93 use m_xomp 94 use m_errors 95 use m_ab7_mixing 96 use m_kxc 97 use m_nctk 98 #ifdef HAVE_NETCDF 99 use netcdf 100 #endif 101 use libxc_functionals 102 use m_hdr 103 104 use m_io_tools, only : open_file, file_exists, iomode_from_fname 105 use m_fstrings, only : int2char10, sjoin, strcat, itoa 106 use m_energies, only : energies_type, energies_init 107 use m_numeric_tools, only : print_arr, iseven, coeffs_gausslegint 108 use m_geometry, only : normv, vdotw, mkrdim, metric 109 use m_gwdefs, only : GW_TOLQ0, GW_TOLQ, em1params_free, em1params_t, GW_Q0_DEFAULT 110 use m_mpinfo, only : destroy_mpi_enreg 111 use m_crystal, only : crystal_free, crystal_t 112 #ifdef HAVE_NETCDF 113 use m_crystal_io, only : crystal_ncwrite 114 #endif 115 use m_ebands, only : ebands_update_occ, ebands_copy, get_valence_idx, get_occupied, apply_scissor, & 116 & ebands_free, ebands_has_metal_scheme, ebands_ncwrite 117 use m_bz_mesh, only : kmesh_t, kmesh_free, littlegroup_t, littlegroup_free 118 use m_kg, only : getph 119 use m_gsphere, only : gsph_free, gsphere_t 120 use m_vcoul, only : vcoul_t, vcoul_free 121 use m_qparticles, only : rdqps, rdgw, show_QP 122 use m_screening, only : make_epsm1_driver, lwl_write, chi_t, chi_free, chi_new 123 use m_io_screening, only : hscr_new, hscr_io, write_screening, hscr_free, hscr_t 124 use m_spectra, only : spectra_t, spectra_write, spectra_repr, spectra_free, W_EM_LF, W_EM_NLF, W_EELF 125 use m_fftcore, only : print_ngfft 126 use m_fft_mesh, only : rotate_FFT_mesh, cigfft, get_gftt 127 use m_wfd, only : wfd_init, wfd_free, wfd_nullify, wfd_print, wfd_t, wfd_rotate, wfd_test_ortho,& 128 & wfd_read_wfk, wfd_test_ortho, wfd_copy, wfd_change_ngfft 129 use m_chi0, only : output_chi0sumrule 130 use m_pawang, only : pawang_type 131 use m_pawrad, only : pawrad_type 132 use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize 133 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 134 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 135 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free 136 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy,& 137 & pawrhoij_free, symrhoij, pawrhoij_get_nspden 138 use m_pawdij, only : pawdij, symdij_all 139 use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free 140 use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free 141 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 142 143 !This section has been created automatically by the script Abilint (TD). 144 !Do not modify the following lines by hand. 145 #undef ABI_FUNC 146 #define ABI_FUNC 'screening' 147 use interfaces_14_hidewrite 148 use interfaces_18_timing 149 use interfaces_51_manage_mpi 150 use interfaces_53_ffts 151 use interfaces_64_psp 152 use interfaces_65_paw 153 use interfaces_67_common 154 use interfaces_69_wfdesc 155 use interfaces_70_gw 156 !End of the abilint section 157 158 implicit none 159 160 !Arguments ------------------------------------ 161 !scalars 162 character(len=6),intent(in) :: codvsn 163 type(Datafiles_type),intent(in) :: Dtfil 164 type(Dataset_type),intent(inout) :: Dtset 165 type(Pawang_type),intent(inout) :: Pawang 166 type(Pseudopotential_type),intent(inout) :: Psps 167 !arrays 168 real(dp),intent(in) :: acell(3),rprim(3,3) 169 type(Pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Dtset%usepaw) 170 type(Pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Dtset%usepaw) 171 172 !Local variables ------------------------------ 173 character(len=4) :: ctype='RPA ' 174 !scalars 175 integer,parameter :: level30=30,tim_fourdp4=4,NOMEGA_PRINTED=15,master=0 176 integer :: spin,ik_ibz,my_nbks 177 integer :: choice,cplex,dim_kxcg,dim_wing,ount,omp_ncpus 178 integer :: fform_chi0,fform_em1,gnt_option,iat,ider,idir,ierr,band 179 integer :: ifft,ii,ikbz,ikxc,initialized,iomega,ios,ipert 180 integer :: iqibz,iqcalc,is_qeq0,isym,izero,ifirst,ilast 181 integer :: label,mgfftf,mgfftgw 182 integer :: nt_per_proc,work_size 183 integer :: moved_atm_inside,moved_rhor 184 integer :: nbcw,nbsc,nbvw,nkxc,nkxc1,n3xccc,optene,istep 185 integer :: nfftf,nfftf_tot,nfftgw,nfftgw_tot,ngrvdw,nhatgrdim,nprocs,nspden_rhoij 186 integer :: nscf,nzlmopt,mband 187 integer :: optcut,optgr0,optgr1,optgr2,option,approx_type,option_test,optgrad 188 integer :: optrad,optrhoij,psp_gencond,my_rank 189 integer :: rhoxsp_method,comm,test_type,tordering,unt_em1,unt_susc,usexcnhat 190 real(dp) :: compch_fft,compch_sph,domegareal,e0,ecore,ecut_eff,ecutdg_eff 191 real(dp) :: gsqcutc_eff,gsqcutf_eff,gsqcut_shp,omegaplasma,ucvol,vxcavg,gw_gsq,r_s 192 real(dp) :: alpha,rhoav,opt_ecut,factor 193 real(dp):: eff,mempercpu_mb,max_wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem 194 logical :: found,iscompatibleFFT,use_tr,is_first_qcalc 195 logical :: add_chi0_intraband,update_energies,call_pawinit 196 character(len=10) :: string 197 character(len=500) :: msg 198 character(len=80) :: bar 199 type(ebands_t) :: KS_BSt,QP_BSt 200 type(kmesh_t) :: Kmesh,Qmesh 201 type(vcoul_t) :: Vcp 202 type(crystal_t) :: Cryst 203 type(em1params_t) :: Ep 204 type(Energies_type) :: KS_energies 205 type(gsphere_t) :: Gsph_epsG0,Gsph_wfn 206 type(Hdr_type) :: Hdr_wfk,Hdr_local 207 type(MPI_type) :: MPI_enreg_seq 208 type(Pawfgr_type) :: Pawfgr 209 type(hscr_t) :: Hem1,Hchi0 210 type(wfd_t) :: Wfd,Wfdf 211 type(spectra_t) :: spectra 212 type(chi_t) :: chihw 213 type(wvl_data) :: wvl_dummy 214 !arrays 215 integer :: ibocc(Dtset%nsppol),ngfft_gw(18),ngfftc(18),ngfftf(18) 216 integer,allocatable :: irottb(:,:),ktabr(:,:),ktabrf(:,:),l_size_atm(:) 217 integer,allocatable :: ks_vbik(:,:),ks_occ_idx(:,:),qp_vbik(:,:),nband(:,:) 218 integer,allocatable :: nq_spl(:),nlmn_atm(:),gw_gfft(:,:) 219 real(dp) :: gmet(3,3),gprimd(3,3),k0(3),qtmp(3),rmet(3,3),rprimd(3,3),tsec(2),strsxc(6) 220 real(dp),allocatable :: igwene(:,:,:),chi0_sumrule(:),ec_rpa(:),rspower(:) 221 real(dp),allocatable :: nhat(:,:),nhatgr(:,:,:),ph1d(:,:),ph1df(:,:) 222 real(dp),allocatable :: rhog(:,:),rhor(:,:),rhor_p(:,:),rhor_kernel(:,:),taug(:,:),taur(:,:) 223 real(dp),allocatable :: z(:),zw(:),grchempottn(:,:),grewtn(:,:),grvdw(:,:),kxc(:,:),qmax(:) 224 real(dp),allocatable :: ks_vhartr(:),vpsp(:),ks_vtrial(:,:),ks_vxc(:,:),xccc3d(:) 225 complex(gwpc),allocatable :: arr_99(:,:),kxcg(:,:),fxc_ADA(:,:,:) 226 complex(dpc),allocatable :: m_lda_to_qp(:,:,:,:) 227 complex(dpc),allocatable :: chi0_lwing(:,:,:),chi0_uwing(:,:,:),chi0_head(:,:,:) 228 complex(dpc),allocatable :: chi0intra_lwing(:,:,:),chi0intra_uwing(:,:,:),chi0intra_head(:,:,:) 229 complex(gwpc),allocatable,target :: chi0(:,:,:),chi0intra(:,:,:) 230 complex(gwpc),ABI_CONTIGUOUS pointer :: epsm1(:,:,:) 231 logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:) 232 character(len=80) :: title(2) 233 character(len=fnlen) :: gw_fname,wfk_fname,lwl_fname 234 type(littlegroup_t),pointer :: Ltg_q(:) 235 type(Paw_an_type),allocatable :: Paw_an(:) 236 type(Paw_ij_type),allocatable :: Paw_ij(:) 237 type(Pawfgrtab_type),allocatable :: Pawfgrtab(:) 238 type(Pawrhoij_type),allocatable :: Pawrhoij(:),prev_Pawrhoij(:) 239 type(pawpwff_t),allocatable :: Paw_pwff(:) 240 type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:) 241 242 !************************************************************************ 243 244 DBG_ENTER("COLL") 245 246 call timab(301,1,tsec) ! overall time 247 call timab(302,1,tsec) ! screening(init 248 249 write(msg,'(6a)')& 250 & ' SCREENING: Calculation of the susceptibility and dielectric matrices ',ch10,ch10,& 251 & ' Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining.',ch10,& 252 & ' Incorporated in ABINIT by V. Olevano, G.-M. Rignanese, and M. Torrent.' 253 call wrtout(ab_out, msg,'COLL') 254 call wrtout(std_out,msg,'COLL') 255 256 if(dtset%ucrpa>0) then 257 write(msg,'(6a)')ch10,& 258 & ' cRPA Calculation: The calculation of the polarisability is constrained (ucrpa/=0)',ch10 259 call wrtout(ab_out, msg,'COLL') 260 call wrtout(std_out,msg,'COLL') 261 end if 262 #if defined HAVE_GW_DPC 263 if (gwpc/=8) then 264 write(msg,'(6a)')ch10,& 265 & ' Number of bytes for double precision complex /=8 ',ch10,& 266 & ' Cannot continue due to kind mismatch in BLAS library ',ch10,& 267 & ' Some BLAS interfaces are not generated by abilint ' 268 MSG_ERROR(msg) 269 end if 270 write(msg,'(a,i2,a)')'.Using double precision arithmetic ; gwpc = ',gwpc,ch10 271 #else 272 write(msg,'(a,i2,a)')'.Using single precision arithmetic ; gwpc = ',gwpc,ch10 273 #endif 274 call wrtout(std_out,msg,'COLL') 275 call wrtout(ab_out, msg,'COLL') 276 277 ! === Initialize MPI variables, and parallelization level === 278 ! gwpara: 0--> sequential run, 1--> parallelism over k-points, 2--> parallelism over bands. 279 ! gwpara==2, each node has both fully and partially occupied states while conduction bands are divided 280 comm = xmpi_world; my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 281 282 if (my_rank == master) then 283 wfk_fname = dtfil%fnamewffk 284 if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then 285 MSG_ERROR(msg) 286 end if 287 end if 288 call xmpi_bcast(wfk_fname, master, comm, ierr) 289 290 ! Some variables need to be initialized/nullify at start 291 call energies_init(KS_energies) 292 usexcnhat=0 293 294 call mkrdim(acell,rprim,rprimd) 295 call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol) 296 297 !=== Define FFT grid(s) sizes === 298 ! Be careful! This mesh is only used for densities and potentials. It is NOT the (usually coarser) 299 ! GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90. 300 ! See also NOTES in the comments at the beginning of this file. 301 ! NOTE: The mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn. 302 303 k0(:)=zero 304 call pawfgr_init(Pawfgr,Dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,& 305 & gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=gmet,k0=k0) 306 307 call print_ngfft(ngfftf,'Dense FFT mesh used for densities and potentials') 308 nfftf_tot=PRODUCT(ngfftf(1:3)) 309 310 ! We can intialize MPI_enreg and fft distrib here, now ngfft are known 311 call initmpi_seq(MPI_enreg_seq) ! Fake MPI_type for the sequential part. 312 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all') 313 call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfftf(2),ngfftf(3),'all') 314 315 !============================================= 316 !==== Open and read pseudopotential files ==== 317 !============================================= 318 call pspini(Dtset,Dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,Pawrad,Pawtab,Psps,rprimd,comm_mpi=comm) 319 320 !=== Initialize dimensions and basic objects === 321 call setup_screening(codvsn,acell,rprim,ngfftf,wfk_fname,dtfil,Dtset,Psps,Pawtab,& 322 & ngfft_gw,Hdr_wfk,Hdr_local,Cryst,Kmesh,Qmesh,KS_BSt,Ltg_q,Gsph_epsG0,Gsph_wfn,Vcp,Ep,comm) 323 324 call timab(302,2,tsec) ! screening(init) 325 call print_ngfft(ngfft_gw,'FFT mesh used for oscillator strengths') 326 327 nfftgw_tot=PRODUCT(ngfft_gw(1:3)) 328 mgfftgw =MAXVAL (ngfft_gw(1:3)) 329 nfftgw =nfftgw_tot ! no FFT // 330 331 !TRYING TO RECREATE AN "ABINIT ENVIRONMENT" 332 KS_energies%e_corepsp=ecore/Cryst%ucvol 333 334 !========================== 335 !=== PAW initialization === 336 !========================== 337 if (Dtset%usepaw==1) then 338 call timab(315,1,tsec) ! screening(pawin 339 340 call chkpawovlp(Cryst%natom,Cryst%ntypat,Dtset%pawovlp,Pawtab,Cryst%rmet,Cryst%typat,Cryst%xred) 341 342 ABI_DT_MALLOC(Pawrhoij,(Cryst%natom)) 343 nspden_rhoij=pawrhoij_get_nspden(Dtset%nspden,Dtset%nspinor,Dtset%pawspnorb) 344 call pawrhoij_alloc(Pawrhoij,Dtset%pawcpxocc,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,& 345 & Cryst%typat,pawtab=Pawtab) 346 347 ! Initialize values for several basic arrays stored in Pawinit 348 gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2 349 350 ! Test if we have to call pawinit 351 call paw_gencond(Dtset,gnt_option,"test",call_pawinit) 352 353 if (psp_gencond==1.or.call_pawinit) then 354 gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2 355 call pawinit(gnt_option,gsqcut_shp,zero,Dtset%pawlcutd,Dtset%pawlmix,& 356 & Psps%mpsang,Dtset%pawnphi,Cryst%nsym,Dtset%pawntheta,Pawang,Pawrad,& 357 & Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%xclevel,Dtset%usepotzero) 358 359 ! Update internal values 360 call paw_gencond(Dtset,gnt_option,"save",call_pawinit) 361 else 362 if (Pawtab(1)%has_kij ==1) Pawtab(1:Cryst%ntypat)%has_kij =2 363 if (Pawtab(1)%has_nabla==1) Pawtab(1:Cryst%ntypat)%has_nabla=2 364 end if 365 Psps%n1xccc=MAXVAL(Pawtab(1:Cryst%ntypat)%usetcore) 366 367 ! Initialize optional flags in Pawtab to zero 368 ! (Cannot be done in Pawinit since the routine is called only if some pars. are changed) 369 Pawtab(:)%has_nabla = 0 370 Pawtab(:)%usepawu = 0 371 Pawtab(:)%useexexch = 0 372 Pawtab(:)%exchmix =zero 373 ! 374 ! * Evaluate <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for the long wavelength limit. 375 ! TODO solve problem with memory leak and clean this part as well as the associated flag 376 call pawnabla_init(Psps%mpsang,Cryst%ntypat,Pawrad,Pawtab) 377 378 call setsymrhoij(gprimd,Pawang%l_max-1,Cryst%nsym,Dtset%pawprtvol,rprimd,Cryst%symrec,Pawang%zarot) 379 380 ! * Initialize and compute data for LDA+U. 381 ! paw_dmft%use_dmft=dtset%usedmft 382 if (Dtset%usepawu>0.or.Dtset%useexexch>0) then 383 call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,& 384 & Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,Cryst%ntypat,Pawang,Dtset%pawprtvol,& 385 & Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa) 386 end if 387 388 if (my_rank == master) call pawtab_print(Pawtab) 389 390 ! Get Pawrhoij from the header of the WFK file. 391 call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij) 392 393 ! Re-symmetrize symrhoij. 394 ! this call leads to a SIGFAULT, likely some pointer is not initialized correctly 395 choice=1; optrhoij=1; ipert=0; idir=0 396 ! call symrhoij(Pawrhoij,Pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert,Cryst%natom,& 397 ! & Cryst%nsym,Cryst%ntypat,optrhoij,Pawang,Dtset%pawprtvol,Pawtab,& 398 ! & Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat) 399 ! 400 ! Evaluate form factors for the radial part of phi.phj-tphi.tphj === 401 ! rhoxsp_method=1 ! Arnaud-Alouani 402 ! rhoxsp_method=2 ! Shiskin-Kresse 403 rhoxsp_method=2 404 405 ! At least for ucrpa, the Arnaud Alouani is always a better choice but needs a larger cutoff 406 if(dtset%ucrpa>0) rhoxsp_method=1 407 if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc 408 409 ABI_MALLOC(gw_gfft,(3,nfftgw_tot)) 410 call get_gftt(ngfft_gw,(/zero,zero,zero/),gmet,gw_gsq,gw_gfft) 411 ABI_FREE(gw_gfft) 412 413 ! Set up q grids, make qmax 20% larger than largest expected: 414 ABI_MALLOC(nq_spl,(Psps%ntypat)) 415 ABI_MALLOC(qmax,(Psps%ntypat)) 416 nq_spl = Psps%mqgrid_ff 417 qmax = SQRT(gw_gsq)*1.2d0 !qmax = Psps%qgrid_ff(Psps%mqgrid_ff) 418 ABI_DT_MALLOC(Paw_pwff,(Psps%ntypat)) 419 420 call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps) 421 ABI_FREE(nq_spl) 422 ABI_FREE(qmax) 423 424 ! Variables/arrays related to the fine FFT grid 425 ABI_MALLOC(nhat,(nfftf,Dtset%nspden)) 426 nhat=zero; cplex=1 427 ABI_DT_MALLOC(Pawfgrtab,(Cryst%natom)) 428 call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat) 429 call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Dtset%nspden,Dtset%typat) 430 ABI_FREE(l_size_atm) 431 compch_fft=greatest_real 432 usexcnhat=MAXVAL(Pawtab(:)%usexcnhat) 433 ! * 0 --> Vloc in atomic data is Vbare (Blochl s formulation) 434 ! * 1 --> Vloc in atomic data is VH(tnzc) (Kresse s formulation) 435 write(msg,'(a,i3)')' screening : using usexcnhat = ',usexcnhat 436 call wrtout(std_out,msg,'COLL') 437 ! 438 ! Identify parts of the rectangular grid where the density has to be calculated. 439 optcut=0;optgr0=Dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-Dtset%pawstgylm 440 if (Dtset%pawcross==1) optrad=1 441 if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm 442 443 call nhatgrid(Cryst%atindx1,gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,& 444 & optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 445 446 call timab(315,2,tsec) ! screening(pawin 447 else 448 ! allocate empty structure for the sake of -fcheck=all... 449 ABI_DT_MALLOC(Paw_pwff,(0)) 450 ABI_DT_MALLOC(Pawrhoij,(0)) 451 ABI_DT_MALLOC(Pawfgrtab,(0)) 452 end if ! End of PAW initialization. 453 454 ! Consistency check and additional stuff done only for GW with PAW. 455 ABI_DT_MALLOC(Paw_onsite,(Cryst%natom)) 456 if (Dtset%usepaw==1) then 457 if (Dtset%ecutwfn < Dtset%ecut) then 458 write(msg,"(5a)")& 459 & "WARNING - ",ch10,& 460 & " It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,& 461 & " an excessive truncation of the planewave basis set can lead to unphysical results." 462 call wrtout(ab_out,msg,'COLL') 463 end if 464 465 ABI_CHECK(Dtset%useexexch==0,"LEXX not yet implemented in GW") 466 ABI_CHECK(Dtset%usedmft==0,"DMFT + GW not available") 467 468 if (Dtset%pawcross==1) then 469 optgrad=1 470 call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,Cryst%rprimd,& 471 & Cryst%xcart,Pawtab,Pawrad,Pawfgrtab,optgrad) 472 end if 473 end if 474 475 !Allocate these arrays anyway, since they are passed to subroutines. 476 if (.not.allocated(nhat)) then 477 ABI_MALLOC(nhat,(nfftf,0)) 478 end if 479 480 call timab(316,1,tsec) ! screening(wfs 481 482 !===================================================== 483 !=== Prepare the distribution of the wavefunctions === 484 !===================================================== 485 ! valence and partially occupied are replicate on each node while conduction bands are MPI distributed. 486 ! This method is mandatory if gwpara==2 and/or we are using awtr==1 or the spectral method. 487 ! If awtr==1, we evaluate chi0 taking advantage of time-reversal (speed-up~2) 488 ! Useful indeces: 489 ! nbvw = Max. number of fully/partially occupied states over spin 490 ! nbcw = Max. number of unoccupied states considering the spin 491 !TODO: 492 ! Here for semiconducting systems we have to be sure that each processor has all the 493 ! states considered in the SCGW, moreover nbsc<nbvw 494 ! in case of SCGW vale and conduction has to be recalculated to avoid errors 495 ! if a metal becomes semiconductor or viceversa. 496 ! Ideally nbvw should include only the states v such that the transition 497 ! c-->v is taken into account in cchi0 (see GW_TOLDOCC). In the present implementation 498 499 ABI_MALLOC(ks_occ_idx,(KS_BSt%nkpt,KS_BSt%nsppol)) 500 ABI_MALLOC(ks_vbik ,(KS_BSt%nkpt,KS_BSt%nsppol)) 501 ABI_MALLOC(qp_vbik ,(KS_BSt%nkpt,KS_BSt%nsppol)) 502 503 call ebands_update_occ(KS_BSt,Dtset%spinmagntarget,prtvol=0) 504 ks_occ_idx = get_occupied(KS_BSt,tol8) ! tol8 to be consistent when the density 505 ks_vbik = get_valence_idx(KS_BSt) 506 507 ibocc(:)=MAXVAL(ks_occ_idx(:,:),DIM=1) ! Max occupied band index for each spin. 508 ABI_FREE(ks_occ_idx) 509 510 use_tr=.FALSE.; nbvw=0 511 if (Dtset%gwpara==2.or.Ep%awtr==1.or.Dtset%spmeth>0) then 512 use_tr=.TRUE. 513 nbvw=MAXVAL(ibocc) 514 nbcw=Ep%nbnds-nbvw 515 write(msg,'(4a,i0,2a,i0,2a,i0,a)')ch10,& 516 & '- screening: taking advantage of time-reversal symmetry ',ch10,& 517 & '- Maximum band index for partially occupied states nbvw = ',nbvw,ch10,& 518 & '- Remaining bands to be divided among processors nbcw = ',nbcw,ch10,& 519 & '- Number of bands treated by each node ~',nbcw/nprocs,ch10 520 call wrtout(ab_out,msg,'COLL') 521 if (Cryst%timrev/=2) then 522 MSG_ERROR('Time-reversal cannot be used since cryst%timrev/=2') 523 end if 524 end if 525 526 mband=Ep%nbnds 527 ABI_MALLOC(nband,(Kmesh%nibz,Dtset%nsppol)) 528 nband=mband 529 ABI_MALLOC(bks_mask,(mband,Kmesh%nibz,Dtset%nsppol)) 530 ABI_MALLOC(keep_ur,(mband,Kmesh%nibz,Dtset%nsppol)) 531 bks_mask=.FALSE.; keep_ur=.FALSE. 532 533 ! autoparal section 534 if (dtset%max_ncpus /=0) then 535 ount = ab_out 536 ! Temporary table needed to estimate memory 537 ABI_MALLOC(nlmn_atm,(Cryst%natom)) 538 if (Dtset%usepaw==1) then 539 do iat=1,Cryst%natom 540 nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size 541 end do 542 end if 543 544 write(ount,'(a)')"--- !Autoparal" 545 write(ount,"(a)")"# Autoparal section for Screening runs" 546 write(ount,"(a)") "info:" 547 write(ount,"(a,i0)")" autoparal: ",dtset%autoparal 548 write(ount,"(a,i0)")" max_ncpus: ",dtset%max_ncpus 549 write(ount,"(a,i0)")" gwpara: ",dtset%gwpara 550 write(ount,"(a,i0)")" nkpt: ",dtset%nkpt 551 write(ount,"(a,i0)")" nsppol: ",dtset%nsppol 552 write(ount,"(a,i0)")" nspinor: ",dtset%nspinor 553 write(ount,"(a,i0)")" nbnds: ",Ep%nbnds 554 555 work_size = nbvw * nbcw * Kmesh%nibz**2 * Dtset%nsppol 556 557 ! Non-scalable memory in Mb i.e. memory that is not distribute with MPI. 558 nonscal_mem = (two*gwpc*Ep%npwe**2*(Ep%nomega*b2Mb)) * 1.1_dp 559 560 ! List of configurations. 561 ! Assuming an OpenMP implementation with perfect speedup! 562 write(ount,"(a)")"configurations:" 563 564 do ii=1,dtset%max_ncpus 565 nt_per_proc = 0 566 eff = HUGE(one) 567 max_wfsmem_mb = zero 568 569 do my_rank=0,ii-1 570 call chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,ii,bks_mask,keep_ur,ierr) 571 if (ierr /= 0) exit 572 nt_per_proc = MAX(nt_per_proc, COUNT(bks_mask(1:nbvw,:,:)) * COUNT(bks_mask(nbvw+1:,:,:))) 573 eff = MIN(eff, (one * work_size) / (ii * nt_per_proc)) 574 575 ! Memory needed for Fourier components ug. 576 my_nbks = COUNT(bks_mask) 577 ug_mem = two*gwpc*Dtset%nspinor*Ep%npwwfn*my_nbks*b2Mb 578 579 ! Memory needed for real space ur. 580 ur_mem = two*gwpc*Dtset%nspinor*nfftgw*COUNT(keep_ur)*b2Mb 581 582 ! Memory needed for PAW projections Cprj 583 cprj_mem = zero 584 if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb 585 586 max_wfsmem_mb = MAX(max_wfsmem_mb, ug_mem + ur_mem + cprj_mem) 587 end do 588 if (ierr /= 0) cycle 589 590 ! Add the non-scalable part and increase by 10% to account for other datastructures. 591 mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp 592 593 do omp_ncpus=1,xomp_get_max_threads() 594 write(ount,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus 595 write(ount,"(a,i0)")" mpi_ncpus: ",ii 596 write(ount,"(a,i0)")" omp_ncpus: ",omp_ncpus 597 write(ount,"(a,f12.9)")" efficiency: ",eff 598 write(ount,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb 599 end do 600 end do 601 write(ount,'(a)')"..." 602 603 ABI_FREE(nlmn_atm) 604 MSG_ERROR_NODUMP("aborting now") 605 else 606 call chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,nprocs,bks_mask,keep_ur,ierr) 607 end if 608 609 ! Initialize the Wf_info object (allocate %ug and %ur if required). 610 opt_ecut=Dtset%ecutwfn 611 !opt_ecut=zero 612 613 !call gsph_init(Gsph_wfn,Cryst,Ep%npwvec,gvec=gvec_kss) 614 615 call wfd_init(Wfd,Cryst,Pawtab,Psps,keep_ur,Dtset%paral_kgb,Ep%npwwfn,mband,nband,Ep%nkibz,Dtset%nsppol,bks_mask,& 616 & Dtset%nspden,Dtset%nspinor,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,ngfft_gw,& 617 & Gsph_wfn%gvec,Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,opt_ecut=opt_ecut) 618 619 if (Dtset%pawcross==1) then 620 call wfd_init(Wfdf,Cryst,Pawtab,Psps,keep_ur,Dtset%paral_kgb,Ep%npwwfn,mband,nband,Ep%nkibz,Dtset%nsppol,bks_mask,& 621 & Dtset%nspden,Dtset%nspinor,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,ngfft_gw,& 622 & Gsph_wfn%gvec,Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,opt_ecut=opt_ecut) 623 end if 624 625 ABI_FREE(bks_mask) 626 ABI_FREE(nband) 627 ABI_FREE(keep_ur) 628 629 call wfd_print(Wfd,mode_paral='PERS') 630 !FIXME: Rewrite the treatment of use_tr branches in cchi0 ... 631 !Use a different nbvw for each spin. 632 !Now use_tr means that one can use time-reversal symmetry. 633 634 !================================================== 635 !==== Read KS band structure from the KSS file ==== 636 !================================================== 637 call wfd_read_wfk(Wfd,wfk_fname,iomode_from_fname(wfk_fname)) 638 639 if (Dtset%pawcross==1) then 640 call wfd_copy(Wfd,Wfdf) 641 call wfd_change_ngfft(Wfdf,Cryst,Psps,ngfftf) 642 end if 643 644 ! This test has been disabled (too expensive!) 645 if (.False.) call wfd_test_ortho(Wfd,Cryst,Pawtab,unit=ab_out,mode_paral="COLL") 646 647 call timab(316,2,tsec) ! screening(wfs 648 call timab(319,1,tsec) ! screening(1) 649 650 if (Cryst%nsym/=Dtset%nsym .and. Dtset%usepaw==1) then 651 MSG_ERROR('Cryst%nsym/=Dtset%nsym, check pawinit and symrhoij') 652 end if 653 654 ! Get the FFT index of $ (R^{-1}(r-\tau)) $ 655 ! S= $\transpose R^{-1}$ and k_BZ = S k_IBZ 656 ! irottb is the FFT index of $ R^{-1} (r-\tau) $ used to symmetrize u_Sk. 657 ABI_MALLOC(irottb,(nfftgw,Cryst%nsym)) 658 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfft_gw,irottb,iscompatibleFFT) 659 660 ABI_MALLOC(ktabr,(nfftgw,Kmesh%nbz)) 661 do ikbz=1,Kmesh%nbz 662 isym=Kmesh%tabo(ikbz) 663 do ifft=1,nfftgw 664 ktabr(ifft,ikbz)=irottb(ifft,isym) 665 end do 666 end do 667 ABI_FREE(irottb) 668 669 if (Dtset%usepaw==1 .and. Dtset%pawcross==1) then 670 ABI_MALLOC(irottb,(nfftf,Cryst%nsym)) 671 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfftf,irottb,iscompatibleFFT) 672 673 ABI_MALLOC(ktabrf,(nfftf,Kmesh%nbz)) 674 do ikbz=1,Kmesh%nbz 675 isym=Kmesh%tabo(ikbz) 676 do ifft=1,nfftf 677 ktabrf(ifft,ikbz)=irottb(ifft,isym) 678 end do 679 end do 680 ABI_FREE(irottb) 681 else 682 ABI_MALLOC(ktabrf,(0,0)) 683 end if 684 685 !=== Compute structure factor phases and large sphere cut-off === 686 !WARNING cannot use Dtset%mgfft, this has to be checked better 687 !mgfft=MAXVAL(ngfftc(:)) 688 !allocate(ph1d(2,3*(2*mgfft+1)*Cryst%natom),ph1df(2,3*(2*mgfftf+1)*Cryst%natom)) 689 write(std_out,*)' CHECK ',Dtset%mgfftdg,mgfftf 690 !if (Dtset%mgfftdg/=mgfftf) write(std_out,*)"WARNING Dtset%mgfftf /= mgfftf" 691 ABI_MALLOC(ph1d,(2,3*(2*Dtset%mgfft+1)*Cryst%natom)) 692 ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*Cryst%natom)) 693 call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred) 694 695 if (Psps%usepaw==1.and.Pawfgr%usefinegrid==1) then 696 call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred) 697 else 698 ph1df(:,:)=ph1d(:,:) 699 end if 700 701 ! Initialize QP_BSt using KS bands 702 ! In case of SCGW, update QP_BSt using the QPS file. 703 call ebands_copy(KS_BSt,QP_BSt) 704 705 call timab(319,2,tsec) ! screening(1) 706 707 !============================ 708 !==== Self-consistent GW ==== 709 !============================ 710 if (Ep%gwcalctyp>=10) then 711 call timab(304,1,tsec) ! KS => QP; [wfrg] 712 713 ! Initialize with KS eigenvalues and eigenfunctions. 714 ABI_MALLOC(m_lda_to_qp,(Wfd%mband,Wfd%mband,Wfd%nkibz,Wfd%nsppol)) 715 m_lda_to_qp = czero 716 do spin=1,Wfd%nsppol 717 do ik_ibz=1,Wfd%nkibz 718 do band=1,Wfd%nband(ik_ibz,spin) 719 m_lda_to_qp(band,band,:,:) = cone 720 end do 721 end do 722 end do 723 724 ! Read unitary transformation and QP energies. 725 ! TODO switch on the renormalization of n in screening, QPS should report bdgw 726 ABI_MALLOC(rhor_p,(nfftf,Dtset%nspden)) 727 ABI_DT_MALLOC(prev_Pawrhoij,(Cryst%natom*Psps%usepaw)) 728 729 call rdqps(QP_BSt,Dtfil%fnameabi_qps,Dtset%usepaw,Dtset%nspden,1,nscf,& 730 & nfftf,ngfftf,Cryst%ucvol,Dtset%paral_kgb,Cryst,Pawtab,MPI_enreg_seq,nbsc,m_lda_to_qp,rhor_p,prev_Pawrhoij) 731 732 ABI_FREE(rhor_p) 733 ABI_DT_FREE(prev_Pawrhoij) 734 735 ! FIXME this is to preserve the old implementation for the head and the wings in ccchi0q0 736 ! But has to be rationalized 737 KS_BSt%eig=QP_BSt%eig 738 739 ! Calculate new occ. factors and fermi level. 740 call ebands_update_occ(QP_BSt,Dtset%spinmagntarget) 741 qp_vbik(:,:) = get_valence_idx(QP_BSt) 742 743 ! === Update only the wfg treated with GW === 744 ! For PAW update and re-symmetrize cprj in the full BZ, TODO add rotation in spinor space 745 if (nscf/=0) call wfd_rotate(Wfd,Cryst,m_lda_to_qp) 746 747 ABI_FREE(m_lda_to_qp) 748 call timab(304,2,tsec) 749 end if ! gwcalctyp>=10 750 751 call timab(305,1,tsec) ! screening(densit 752 ! 753 !=== In case update the eigenvalues === 754 !* Either use a scissor operator or an external GW file. 755 gw_fname = "__in.gw__" 756 update_energies = file_exists(gw_fname) 757 758 if (ABS(Ep%mbpt_sciss)>tol6) then 759 write(msg,'(5a,f7.3,a)')& 760 & ' screening : performing a first self-consistency',ch10,& 761 & ' update of the energies in W by a scissor operator',ch10,& 762 & ' applying a scissor operator of [eV] : ',Ep%mbpt_sciss*Ha_eV,ch10 763 call wrtout(std_out,msg,'COLL') 764 call wrtout(ab_out,msg,'COLL') 765 call apply_scissor(QP_BSt,Ep%mbpt_sciss) 766 else if (update_energies) then 767 write(msg,'(4a)')& 768 & ' screening : performing a first self-consistency',ch10,& 769 & ' update of the energies in W by a previous GW calculation via GW file: ',TRIM(gw_fname) 770 call wrtout(std_out,msg,'COLL') 771 call wrtout(ab_out,msg,'COLL') 772 ABI_MALLOC(igwene,(QP_Bst%mband,QP_Bst%nkpt,QP_Bst%nsppol)) 773 call rdgw(QP_BSt,gw_fname,igwene,extrapolate=.FALSE.) 774 ! call rdgw(QP_BSt,gw_fname,igwene,extrapolate=.TRUE.) 775 ABI_FREE(igwene) 776 call ebands_update_occ(QP_BSt,Dtset%spinmagntarget) 777 end if 778 779 !======================== 780 !=== COMPUTE DENSITY ==== 781 !======================== 782 !* Evaluate PW part (complete charge in case of NC pseudos) 783 !TODO this part has to be rewritten. If I decrease the tol on the occupations 784 !I have to code some MPI stuff also if use_tr==.TRUE. 785 786 ABI_MALLOC(rhor,(nfftf,Dtset%nspden)) 787 ABI_MALLOC(taur,(nfftf,Dtset%nspden*Dtset%usekden)) 788 789 call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,QP_BSt,ngfftf,nfftf,rhor) 790 if (Dtset%usekden==1) then 791 call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,QP_BSt,ngfftf,nfftf,taur,optcalc=1) 792 end if 793 794 call timab(305,2,tsec) ! screening(densit 795 796 nhatgrdim = 0 797 if (Dtset%usepaw==1) then ! Additional computation for PAW. 798 call timab(320,1,tsec) ! screening(paw 799 800 ! Add the compensation charge to the PW density. 801 nhatgrdim=0; if (Dtset%xclevel==2) nhatgrdim=usexcnhat*Dtset%pawnhatxc 802 cplex=1; ider=2*nhatgrdim; izero=0 803 if (nhatgrdim>0) then 804 ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,3)) 805 else 806 ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,0)) 807 end if 808 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,Cryst%gprimd,& 809 & Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,& 810 & Pawfgrtab,nhatgr,nhat,Pawrhoij,Pawrhoij,Pawtab,k0,Cryst%rprimd,Cryst%ucvol,dtset%usewvl,Cryst%xred) 811 812 ! === Evaluate onsite energies, potentials, densities === 813 ! * Initialize variables/arrays related to the PAW spheres. 814 ! * Initialize also lmselect (index of non-zero LM-moments of densities). 815 cplex=1 816 ABI_DT_MALLOC(Paw_ij,(Cryst%natom)) 817 call paw_ij_nullify(Paw_ij) 818 call paw_ij_init(Paw_ij,cplex,Dtset%nspinor,Wfd%nsppol,& 819 & Wfd%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,& 820 & has_dij=1,has_dijhartree=1,has_exexch_pot=1,has_pawu_occ=1) 821 822 nkxc1=0 823 ABI_DT_MALLOC(Paw_an,(Cryst%natom)) 824 call paw_an_nullify(Paw_an) 825 call paw_an_init(Paw_an,Cryst%natom,Cryst%ntypat,nkxc1,Dtset%nspden,& 826 & cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=0) 827 828 nzlmopt=-1; option=0; compch_sph=greatest_real 829 call pawdenpot(compch_sph,KS_energies%e_paw,KS_energies%e_pawdc,ipert,Dtset%ixc,& 830 & Cryst%natom,Cryst%natom,Dtset%nspden,Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,Paw_an,Paw_an,& 831 & Paw_ij,Pawang,Dtset%pawprtvol,Pawrad,Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,& 832 & Dtset%xclevel,Dtset%xc_denpos,Cryst%ucvol,Psps%znuclpsp) 833 call timab(320,2,tsec) ! screening(paw 834 else 835 ABI_DT_MALLOC(Paw_ij,(0)) 836 ABI_DT_MALLOC(Paw_an,(0)) 837 end if ! usepaw 838 839 call timab(321,1,tsec) ! screening(2) 840 841 !JB : Should be remove : cf. l 839 842 if (.not.allocated(nhatgr)) then 843 ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,0)) 844 end if 845 846 call test_charge(nfftf,KS_BSt%nelect,Dtset%nspden,rhor,ucvol,& 847 & Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,omegaplasma) 848 849 !For PAW, add the compensation charge the FFT mesh, then get rho(G). 850 if (Dtset%usepaw==1) rhor(:,:)=rhor(:,:)+nhat(:,:) 851 852 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,rhor,ucvol=ucvol) 853 if(Dtset%usekden==1)then 854 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,taur,ucvol=ucvol,optrhor=1) 855 end if 856 857 if (dtset%gwgamma>0) then 858 ABI_MALLOC(rhor_kernel,(nfftf,Dtset%nspden)) 859 end if 860 861 ABI_MALLOC(rhog,(2,nfftf)) 862 ABI_MALLOC(taug,(2,nfftf*Dtset%usekden)) 863 call fourdp(1,rhog,rhor(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp4) 864 if(Dtset%usekden==1)then 865 call fourdp(1,taug,taur(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp4) 866 end if 867 868 !The following steps have been gathered in the setvtr routine: 869 !- get Ewald energy and Ewald forces 870 !- compute local ionic pseudopotential vpsp 871 !- eventually compute 3D core electron density xccc3d 872 !- eventually compute vxc and vhartr 873 !- set up ks_vtrial 874 !************************************************************** 875 !**** NOTE THAT Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION **** 876 !************************************************************** 877 878 ngrvdw=0 879 ABI_MALLOC(grvdw,(3,ngrvdw)) 880 ABI_MALLOC(grchempottn,(3,Cryst%natom)) 881 ABI_MALLOC(grewtn,(3,Cryst%natom)) 882 nkxc=0 883 if (Dtset%nspden==1) nkxc=2 884 if (Dtset%nspden>=2) nkxc=3 ! check GGA and spinor that is messy !!! 885 ! If MGGA, fxc and kxc are not available and we dont need them for the screening part (for now ...) 886 if (Dtset%ixc<0 .and. libxc_functionals_ismgga()) nkxc=0 887 if (nkxc/=0) then 888 ABI_MALLOC(kxc,(nfftf,nkxc)) 889 end if 890 891 n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf 892 ABI_MALLOC(xccc3d,(n3xccc)) 893 ABI_MALLOC(ks_vhartr,(nfftf)) 894 ABI_MALLOC(ks_vtrial,(nfftf,Dtset%nspden)) 895 ABI_MALLOC(vpsp,(nfftf)) 896 ABI_MALLOC(ks_vxc,(nfftf,Dtset%nspden)) 897 898 optene=4; moved_atm_inside=0; moved_rhor=0; initialized=1; istep=1 899 call setvtr(Cryst%atindx1,Dtset,KS_energies,Cryst%gmet,Cryst%gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,istep,kxc,mgfftf,& 900 & moved_atm_inside,moved_rhor,MPI_enreg_seq, & 901 & Cryst%nattyp,nfftf,ngfftf,ngrvdw,nhat,nhatgr,nhatgrdim,nkxc,Cryst%ntypat,& 902 & Psps%n1xccc,n3xccc,optene,pawrad,Pawtab,ph1df,Psps,rhog,rhor,Cryst%rmet,Cryst%rprimd,strsxc,Cryst%ucvol,usexcnhat,& 903 & ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,wvl_dummy,xccc3d,Cryst%xred,taug=taug,taur=taur) 904 905 if (nkxc/=0) then 906 ABI_FREE(kxc) 907 end if 908 ABI_FREE(grchempottn) 909 ABI_FREE(grewtn) 910 ABI_FREE(grvdw) 911 ABI_FREE(xccc3d) 912 913 !============================ 914 !==== Compute KS PAW Dij ==== 915 !============================ 916 if (Dtset%usepaw==1) then 917 call timab(561,1,tsec) 918 919 ! Calculate unsymmetrized Dij. 920 cplex=1; ipert=0; idir=0 921 call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert,& 922 & Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),& 923 & Dtset%nspden,Cryst%ntypat,Paw_an,Paw_ij,Pawang,Pawfgrtab,Dtset%pawprtvol,& 924 & Pawrad,Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,k0,Dtset%spnorbscl,& 925 & Cryst%ucvol,dtset%charge,ks_vtrial,ks_vxc,Cryst%xred,& 926 & nucdipmom=Dtset%nucdipmom) 927 928 ! Symmetrize KS Dij 929 #if 0 930 call symdij(Cryst%gprimd,Cryst%indsym,ipert,Cryst%natom,Cryst%natom,& 931 & Cryst%nsym,Cryst%ntypat,0,Paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,& 932 & Cryst%symrec) 933 #else 934 call symdij_all(Cryst%gprimd,Cryst%indsym,ipert,Cryst%natom,Cryst%natom,& 935 & Cryst%nsym,Cryst%ntypat,Paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,& 936 & Cryst%symrec) 937 #endif 938 ! 939 ! Output of the pseudopotential strengths Dij and the augmentation occupancies Rhoij. 940 call pawprt(Dtset,Cryst%natom,Paw_ij,Pawrhoij,Pawtab) 941 call timab(561,2,tsec) 942 end if 943 ! 944 !=== Calculate the frequency mesh === 945 !* First omega is always zero without broadening. 946 !FIXME what about metals? I think we should add eta, this means we need to know if the system is metallic, for example using occopt 947 !MS Modified to account for non-zero starting frequency (19-11-2010) 948 !MS Modified for tangent grid (07-01-2011) 949 ABI_MALLOC(Ep%omega,(Ep%nomega)) 950 Ep%omega(1)=CMPLX(Ep%omegaermin,zero,kind=dpc) 951 952 if (Ep%nomegaer>1) then ! Avoid division by zero. 953 if (Dtset%gw_frqre_tangrid==0.and.Dtset%gw_frqre_inzgrid==0) then 954 domegareal=(Ep%omegaermax-Ep%omegaermin)/(Ep%nomegaer-1) 955 do iomega=2,Ep%nomegaer 956 Ep%omega(iomega)=CMPLX(Ep%omegaermin+(iomega-1)*domegareal,zero,kind=dpc) 957 end do 958 else if (Dtset%gw_frqre_tangrid==1.and.Dtset%gw_frqre_inzgrid==0) then ! We have tangent transformed grid 959 MSG_WARNING('EXPERIMENTAL - Using tangent transform grid for contour deformation.') 960 Ep%omegaermax = Dtset%cd_max_freq 961 Ep%omegaermin = zero 962 ifirst=1; ilast=Ep%nomegaer 963 if (Dtset%cd_subset_freq(1)/=0) then ! Only a subset of frequencies is being calculated 964 ifirst=Dtset%cd_subset_freq(1); ilast=Dtset%cd_subset_freq(2) 965 end if 966 factor = Dtset%cd_halfway_freq/TAN(pi*quarter) 967 ! *Important*-here nfreqre is used because the step is set by the original grid 968 domegareal=(ATAN(Ep%omegaermax/factor)*two*piinv)/(Dtset%nfreqre-1) ! Stepsize in transformed variable 969 do iomega=1,Ep%nomegaer 970 Ep%omega(iomega)=CMPLX(factor*TAN((iomega+ifirst-2)*domegareal*pi*half),zero,kind=dpc) 971 end do 972 Ep%omegaermin = REAL(Ep%omega(1)) 973 Ep%omegaermax = REAL(Ep%omega(Ep%nomegaer)) 974 else if (Dtset%gw_frqre_tangrid==0.and.Dtset%gw_frqre_inzgrid==1) then 975 e0=Dtset%ppmfrq; if (e0<0.1d-4) e0=omegaplasma 976 domegareal=one/(Ep%nomegaer) 977 do iomega=1,Ep%nomegaer 978 factor = (iomega-1)*domegareal 979 Ep%omega(iomega)=CMPLX(e0*factor/(one-factor),zero,kind=dpc) 980 end do 981 Ep%omegaermin = REAL(Ep%omega(1)) 982 Ep%omegaermax = REAL(Ep%omega(Ep%nomegaer)) 983 else 984 MSG_ERROR('Error in specification of real frequency grid') 985 end if 986 end if 987 988 if (Ep%plasmon_pole_model.and.Ep%nomega==2) then 989 e0=Dtset%ppmfrq; if (e0<0.1d-4) e0=omegaplasma 990 Ep%omega(2)=CMPLX(zero,e0,kind=dpc) 991 end if 992 ! 993 !=== For AC, use Gauss-Legendre quadrature method === 994 !* Replace $ \int_0^\infty dx f(x) $ with $ \int_0^1 dz f(1/z - 1)/z^2 $. 995 !* Note that the grid is not log as required by CD, so we cannot use the same SCR file. 996 if (Ep%analytic_continuation) then 997 ABI_MALLOC(z,(Ep%nomegaei)) 998 ABI_MALLOC(zw,(Ep%nomegaei)) 999 call coeffs_gausslegint(zero,one,z,zw,Ep%nomegaei) 1000 do iomega=1,Ep%nomegaei 1001 Ep%omega(Ep%nomegaer+iomega)=CMPLX(zero,one/z(iomega)-one,kind=dpc) 1002 end do 1003 ABI_FREE(z) 1004 ABI_FREE(zw) 1005 else if (Ep%contour_deformation.and.(Dtset%cd_customnimfrqs/=0)) then 1006 Ep%omega(Ep%nomegaer+1)=CMPLX(zero,Dtset%cd_imfrqs(1)) 1007 do iomega=2,Ep%nomegaei 1008 if (Dtset%cd_imfrqs(iomega)<=Dtset%cd_imfrqs(iomega-1)) then 1009 MSG_ERROR(' Specified imaginary frequencies need to be strictly increasing!') 1010 end if 1011 Ep%omega(Ep%nomegaer+iomega)=CMPLX(zero,Dtset%cd_imfrqs(iomega)) 1012 end do 1013 else if (Ep%contour_deformation.and.(Dtset%gw_frqim_inzgrid/=0)) then 1014 e0=Dtset%ppmfrq; if (e0<0.1d-4) e0=omegaplasma 1015 domegareal=one/(Ep%nomegaei+1) 1016 do iomega=1,Ep%nomegaei 1017 factor = iomega*domegareal 1018 Ep%omega(Ep%nomegaer+iomega)=CMPLX(zero,e0*factor/(one-factor),kind=dpc) 1019 end do 1020 else if (Ep%contour_deformation.and.(Ep%nomegaei/=0)) then 1021 e0=Dtset%ppmfrq; if (e0<0.1d-4) e0=omegaplasma 1022 do iomega=1,Ep%nomegaei 1023 Ep%omega(Ep%nomegaer+iomega)=CMPLX(zero,e0/(Dtset%freqim_alpha-two)& 1024 & *(EXP(two/(Ep%nomegaei+1)*LOG(Dtset%freqim_alpha-one)*iomega)-one),kind=dpc) 1025 end do 1026 end if 1027 1028 if (Dtset%cd_full_grid/=0) then ! Full grid will be calculated 1029 ! Grid values are added after the last imaginary freq. 1030 do ios=1,Ep%nomegaei 1031 do iomega=2,Ep%nomegaer 1032 Ep%omega(Ep%nomegaer+Ep%nomegaei+(ios-1)*(Ep%nomegaer-1)+(iomega-1)) = & 1033 & CMPLX(REAL(Ep%omega(iomega)),AIMAG(Ep%omega(Ep%nomegaer+ios))) 1034 end do 1035 end do 1036 end if 1037 1038 ! Report frequency mesh for chi0. 1039 write(msg,'(2a)')ch10,' calculating chi0 at frequencies [eV] :' 1040 call wrtout(std_out,msg,'COLL'); call wrtout(ab_out,msg,'COLL') 1041 do iomega=1,Ep%nomega 1042 write(msg,'(i3,2es16.6)')iomega,Ep%omega(iomega)*Ha_eV 1043 call wrtout(std_out,msg,'COLL'); call wrtout(ab_out,msg,'COLL') 1044 end do 1045 1046 ! Allocate chi0, wings and array for chi0_sumrule check. 1047 ABI_MALLOC(chi0_sumrule,(Ep%npwe)) 1048 1049 write(msg,'(a,f12.1,a)')' Memory required for chi0 matrix= ',two*gwpc*Ep%npwe**2*Ep%nI*Ep%nJ*Ep%nomega*b2Mb," [Mb]." 1050 call wrtout(std_out,msg,'COLL') 1051 ABI_STAT_MALLOC(chi0,(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega), ierr) 1052 ABI_CHECK(ierr==0, "Out of memory in chi0") 1053 ! 1054 !============================== END OF THE INITIALIZATION PART =========================== 1055 ! 1056 !====================================================================== 1057 !==== Loop over q-points. Calculate \epsilon^{-1} and save on disc ==== 1058 !====================================================================== 1059 call timab(321,2,tsec) ! screening(2) 1060 1061 iqcalc = 0 1062 do iqibz=1,Qmesh%nibz 1063 call timab(306,1,tsec) 1064 is_first_qcalc=(iqibz==1) 1065 1066 ! Selective q-point calculation. 1067 found=.FALSE.; label=iqibz 1068 if (Ep%nqcalc/=Ep%nqibz) then 1069 do ii=1,Ep%nqcalc 1070 qtmp(:)=Qmesh%ibz(:,iqibz)-Ep%qcalc(:,ii) 1071 found=(normv(qtmp,gmet,'G')<GW_TOLQ) 1072 if (found) then 1073 label=ii; EXIT !ii 1074 end if 1075 end do 1076 if (.not.found) CYCLE !iqibz 1077 qtmp(:)=Ep%qcalc(:,1)-Qmesh%ibz(:,iqibz) 1078 is_first_qcalc=(normv(qtmp,gmet,'G')<GW_TOLQ) 1079 end if 1080 iqcalc = iqcalc + 1 1081 1082 bar=REPEAT('-',80) 1083 write(msg,'(4a,1x,a,i2,a,f9.6,2(",",f9.6),3a)')ch10,ch10,bar,ch10,& 1084 & ' q-point number ',label,' q = (',(Qmesh%ibz(ii,iqibz),ii=1,3),') [r.l.u.]',ch10,bar 1085 call wrtout(std_out,msg,'COLL') 1086 call wrtout(ab_out,msg,'COLL') 1087 is_qeq0 = 0; if (normv(Qmesh%ibz(:,iqibz),gmet,'G')<GW_TOLQ0) is_qeq0=1 1088 1089 call timab(306,2,tsec) 1090 1091 if (is_qeq0==1) then 1092 ! Special treatment of the long wavelength limit. 1093 call timab(307,1,tsec) 1094 1095 ABI_MALLOC(chi0_lwing,(Ep%npwe*Ep%nI,Ep%nomega,3)) 1096 ABI_MALLOC(chi0_uwing,(Ep%npwe*Ep%nJ,Ep%nomega,3)) 1097 ABI_MALLOC(chi0_head,(3,3,Ep%nomega)) 1098 1099 call cchi0q0(use_tr,Dtset,Cryst,Ep,Psps,Kmesh,QP_BSt,KS_BSt,Gsph_epsG0,& 1100 & Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,Pawfgrtab,Paw_onsite,ktabr,ktabrf,nbvw,ngfft_gw,nfftgw,& 1101 & ngfftf,nfftf_tot,chi0,chi0_head,chi0_lwing,chi0_uwing,Ltg_q(iqibz),chi0_sumrule,Wfd,Wfdf) 1102 1103 chihw = chi_new(ep%npwe, ep%nomega) 1104 chihw%head = chi0_head 1105 chihw%lwing = chi0_lwing 1106 chihw%uwing = chi0_uwing 1107 1108 ! Add the intraband term if required and metallic occupation scheme is used. 1109 add_chi0_intraband=.FALSE. !add_chi0_intraband=.TRUE. 1110 if (add_chi0_intraband .and. ebands_has_metal_scheme(QP_BSt)) then 1111 1112 ABI_STAT_MALLOC(chi0intra,(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega), ierr) 1113 ABI_CHECK(ierr==0, "Out of memory in chi0intra") 1114 1115 ABI_MALLOC(chi0intra_lwing,(Ep%npwe*Ep%nI,Ep%nomega,3)) 1116 ABI_MALLOC(chi0intra_uwing,(Ep%npwe*Ep%nJ,Ep%nomega,3)) 1117 ABI_MALLOC(chi0intra_head,(3,3,Ep%nomega)) 1118 1119 call chi0q0_intraband(Wfd,Cryst,Ep,Psps,QP_BSt,Gsph_epsG0,Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,use_tr,Dtset%usepawu,& 1120 ngfft_gw,chi0intra,chi0intra_head,chi0intra_lwing,chi0intra_uwing) 1121 1122 call wrtout(std_out,"Head of chi0 and chi0_intra","COLL") 1123 do iomega=1,Ep%nomega 1124 write(std_out,*)Ep%omega(iomega)*Ha_eV,REAL(chi0(1,1,iomega)),REAL(chi0intra(1,1,iomega)) 1125 write(std_out,*)Ep%omega(iomega)*Ha_eV,AIMAG(chi0(1,1,iomega)),AIMAG(chi0intra(1,1,iomega)) 1126 end do 1127 1128 chi0 = chi0 + chi0intra 1129 chi0_head = chi0_head + chi0intra_head 1130 chi0_lwing = chi0_lwing + chi0intra_lwing 1131 chi0_uwing = chi0_uwing + chi0intra_uwing 1132 1133 ABI_FREE(chi0intra) 1134 ABI_FREE(chi0intra_lwing) 1135 ABI_FREE(chi0intra_uwing) 1136 ABI_FREE(chi0intra_head) 1137 end if 1138 1139 if (.False.) then 1140 lwl_fname = strcat(dtfil%filnam_ds(4), "_LWL") 1141 call lwl_write(lwl_fname,cryst,vcp,ep%npwe,ep%nomega,gsph_epsg0%gvec,chi0,chi0_head,chi0_lwing,chi0_uwing,comm) 1142 end if 1143 1144 call timab(307,2,tsec) 1145 1146 else 1147 ! Calculate cchi0 for q/=0. 1148 call timab(308,1,tsec) 1149 1150 call cchi0(use_tr,Dtset,Cryst,Qmesh%ibz(:,iqibz),Ep,Psps,Kmesh,QP_BSt,Gsph_epsG0,& 1151 & Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,nbvw,ngfft_gw,nfftgw,ngfftf,nfftf_tot,chi0,ktabr,ktabrf,& 1152 & Ltg_q(iqibz),chi0_sumrule,Wfd,Wfdf) 1153 1154 call timab(308,2,tsec) 1155 end if 1156 1157 ! Print chi0(q,G,Gp,omega), then calculate epsilon and epsilon^-1 for this q-point. 1158 ! Only master works but this part could be parallelized over frequencies. 1159 call timab(309,1,tsec) 1160 1161 do iomega=1,MIN(Ep%nomega,NOMEGA_PRINTED) 1162 write(msg,'(1x,a,i4,a,2f9.4,a)')' chi0(G,G'') at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]' 1163 call wrtout(std_out,msg,'COLL') 1164 call wrtout(ab_out,msg,'COLL') 1165 write(msg,'(1x,a,i3,a,i4,a)')' chi0(q =',iqibz, ', omega =',iomega,', G,G'')' 1166 if (Ep%nqcalc/=Ep%nqibz) write(msg,'(a,i3,a,i4,a)')' chi0(q=',iqcalc,', omega=',iomega,', G,G'')' 1167 call wrtout(std_out,msg,'COLL') 1168 ! arr99 is needed to avoid the update of all the tests. Now chi0 is divided by ucvol inside (cchi0|cchi0q0). 1169 ! TODO should be removed but GW tests have to be updated. 1170 ii = MIN(9,Ep%npwe) 1171 ABI_MALLOC(arr_99,(ii,ii)) 1172 arr_99 = chi0(1:ii,1:ii,iomega)*ucvol 1173 call print_arr(arr_99,max_r=2,unit=ab_out) 1174 call print_arr(arr_99,unit=std_out) 1175 ABI_FREE(arr_99) 1176 !call print_arr(chi0(:,:,iomega),max_r=2,unit=ab_out) 1177 !call print_arr(chi0(:,:,iomega),unit=std_out) 1178 end do 1179 1180 if (Ep%nomega>NOMEGA_PRINTED) then 1181 write(msg,'(a,i3,a)')' No. of calculated frequencies > ',NOMEGA_PRINTED,', stop printing ' 1182 call wrtout(std_out,msg,'COLL') 1183 call wrtout(ab_out,msg,'COLL') 1184 end if 1185 1186 ! Write chi0 to _SUSC file 1187 ! Master creates and write the header if this is the first q-point calculated. 1188 if (Dtset%prtsuscep>0 .and. my_rank==master) then 1189 title(1)='CHI0 file: chi0' 1190 title(2)=' ' 1191 if (is_qeq0==1) then 1192 string='0'; if (Dtset%usepaw==0.and.Ep%inclvkb/=0) call int2char10(Ep%inclvkb,string) 1193 title(1)=title(1)(1:21)//', calculated using inclvkb = '//string 1194 end if 1195 1196 ! Open file and write header for polarizability files. 1197 if (is_first_qcalc) then 1198 ikxc=0; test_type=0; tordering=1 1199 hchi0 = hscr_new("polarizability",dtset,ep,hdr_local,ikxc,test_type,tordering,title,Ep%npwe,Gsph_epsG0%gvec) 1200 if (dtset%iomode == IO_MODE_ETSF) then 1201 #ifdef HAVE_NETCDF 1202 NCF_CHECK(nctk_open_create(unt_susc, nctk_ncify(dtfil%fnameabo_sus), xmpi_comm_self)) 1203 NCF_CHECK(crystal_ncwrite(cryst, unt_susc)) 1204 NCF_CHECK(ebands_ncwrite(QP_BSt, unt_susc)) 1205 #endif 1206 else 1207 unt_susc=Dtfil%unchi0 1208 if (open_file(dtfil%fnameabo_sus,msg,unit=unt_susc,status='unknown',form='unformatted') /= 0) then 1209 MSG_ERROR(msg) 1210 end if 1211 end if 1212 fform_chi0 = hchi0%fform 1213 call hscr_io(hchi0,fform_chi0,2,unt_susc,xmpi_comm_self,0,Dtset%iomode) 1214 call hscr_free(Hchi0) 1215 end if 1216 call write_screening("polarizability",unt_susc,Dtset%iomode,Ep%npwe,Ep%nomega,iqcalc,chi0) 1217 end if 1218 1219 ! Calculate the RPA functional correlation energy if the polarizability on a 1220 ! Gauss-Legendre mesh along imaginary axis is available 1221 if (Ep%analytic_continuation .and. Dtset%gwrpacorr>0 ) then 1222 if (is_first_qcalc) then 1223 ABI_MALLOC(ec_rpa,(Dtset%gwrpacorr)) 1224 ec_rpa(:)=zero 1225 end if 1226 call calc_rpa_functional(Dtset%gwrpacorr,label,iqibz,Ep,Vcp,Qmesh,Dtfil,gmet,chi0,comm,ec_rpa) 1227 if (label==Ep%nqcalc) then 1228 ABI_FREE(ec_rpa) 1229 end if 1230 end if 1231 1232 ! ========================================================== 1233 ! === Calculate RPA \tilde\epsilon^{-1} overwriting chi0 === 1234 ! ========================================================== 1235 approx_type=0 ! RPA 1236 option_test=0 ! TESTPARTICLE 1237 dim_wing=0; if (is_qeq0==1) dim_wing=3 1238 1239 if (dim_wing==0) then 1240 dim_wing=1 1241 if (.not.allocated(chi0_lwing)) then 1242 ABI_MALLOC(chi0_lwing,(Ep%npwe*Ep%nI,Ep%nomega,dim_wing)) 1243 end if 1244 if (.not.allocated(chi0_uwing)) then 1245 ABI_MALLOC(chi0_uwing,(Ep%npwe*Ep%nJ,Ep%nomega,dim_wing)) 1246 end if 1247 if (.not.allocated(chi0_head )) then 1248 ABI_MALLOC(chi0_head,(dim_wing,dim_wing,Ep%nomega)) 1249 end if 1250 dim_wing=0 1251 end if 1252 1253 #if 0 1254 ! Using the random q for the optical limit is one of the reasons 1255 ! why sigma breaks the initial energy degeneracies. 1256 chi0_lwing=czero 1257 chi0_uwing=czero 1258 chi0_head=czero 1259 #endif 1260 1261 ! Setup flags for the computation of em1 1262 ! If the vertex is being included for the spectrum, calculate the kernel now and pass it on 1263 if (dtset%gwgamma>0) rhor_kernel = rhor 1264 1265 select case (dtset%gwgamma) 1266 case (0) 1267 approx_type=0; option_test=0; dim_kxcg=0 1268 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1269 1270 case (1, 2) 1271 ! ALDA TDDFT kernel vertex 1272 ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available") 1273 MSG_WARNING('EXPERIMENTAL: Kernel is being added to screening, the SCR file will be non-standard!!') 1274 ikxc=7; approx_type=1; dim_kxcg=1 1275 if (Dtset%gwgamma==1) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma 1276 if (Dtset%gwgamma==2) option_test=0 ! TESTPARTICLE, vertex in chi0 only 1277 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1278 call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,rhor_kernel,& 1279 Ep%npwe,dim_kxcg,kxcg,Gsph_epsG0%gvec,xmpi_comm_self) 1280 1281 case (3, 4) 1282 ! ADA non-local kernel vertex 1283 ABI_CHECK(Wfd%usepaw==0,"ADA vertex + PAW not available") 1284 ABI_CHECK(Wfd%nsppol==1,"ADA vertex for GWGamma not available yet for spin-polarised cases") 1285 MSG_WARNING('EXPERIMENTAL: Kernel is being added to screening, the SCR file will be non-standard!!') 1286 ikxc=7; approx_type=2 1287 if (Dtset%gwgamma==3) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma 1288 if (Dtset%gwgamma==4) option_test=0 ! TESTPARTICLE, vertex in chi0 only 1289 ABI_MALLOC(fxc_ADA,(Ep%npwe,Ep%npwe,Ep%nqibz)) 1290 ! Use userrd to set kappa 1291 if (Dtset%userrd==zero) Dtset%userrd = 2.1_dp 1292 ! Set correct value of kappa (should be scaled with alpha*r_s where) 1293 ! r_s is Wigner-Seitz radius and alpha=(4/(9*Pi))^(1/3) 1294 rhoav = (omegaplasma*omegaplasma)/four_pi 1295 r_s = (three/(four_pi*rhoav))**third 1296 alpha = (four*ninth*piinv)**third 1297 Dtset%userrd = Dtset%userrd*alpha*r_s 1298 1299 call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf,Wfd%nspden,rhor_kernel,Ep%npwe,Ep%nqibz,Ep%qibz,& 1300 fxc_ADA,Gsph_epsG0%gvec,xmpi_comm_self,kappa_init=Dtset%userrd) 1301 1302 dim_kxcg=0 1303 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1304 1305 ! @WC: bootstrap -- 1306 case (-3, -4, -5, -6, -7, -8) 1307 ABI_CHECK(Dtset%usepaw==0,"GWGamma + PAW not available") 1308 if (Dtset%gwgamma>-5) then 1309 MSG_WARNING('EXPERIMENTAL: Bootstrap kernel is being added to screening') 1310 approx_type=4 1311 else if (Dtset%gwgamma>-7) then 1312 MSG_WARNING('EXPERIMENTAL: Bootstrap kernel (head-only) is being added to screening') 1313 approx_type=5 1314 else 1315 MSG_WARNING('EXPERIMENTAL: Bootstrap kernel (RPA-type, head-only) is being added to screening') 1316 approx_type=6 1317 end if 1318 dim_kxcg=0 1319 option_test=MOD(Dtset%gwgamma,2) 1320 ! 1 -> TESTELECTRON, vertex in chi0 *and* sigma 1321 ! 0 -> TESTPARTICLE, vertex in chi0 only 1322 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1323 !--@WC 1324 1325 case default 1326 MSG_ERROR(sjoin("Wrong gwgamma:", itoa(dtset%gwgamma))) 1327 end select 1328 1329 if (approx_type<2) then 1330 call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,& 1331 approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,& 1332 chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm) 1333 else if (approx_type<3) then !@WC: ADA 1334 call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,& 1335 approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,& 1336 chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm,fxc_ADA=fxc_ADA(:,:,iqibz)) 1337 else if (approx_type<7) then !@WC: bootstrap 1338 call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,& 1339 approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,& 1340 chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm) 1341 else 1342 MSG_ERROR(sjoin("Wrong approx_type:", itoa(approx_type))) 1343 end if 1344 1345 ABI_FREE(chi0_lwing) 1346 ABI_FREE(chi0_uwing) 1347 ABI_FREE(chi0_head) 1348 1349 if (my_rank==master .and. is_qeq0==1) then 1350 call spectra_repr(spectra,msg) 1351 call wrtout(std_out,msg,'COLL') 1352 call wrtout(ab_out,msg,'COLL') 1353 1354 if (Ep%nomegaer>2) then 1355 call spectra_write(spectra,W_EELF ,Dtfil%fnameabo_eelf) 1356 call spectra_write(spectra,W_EM_LF ,Dtfil%fnameabo_em1_lf) 1357 call spectra_write(spectra,W_EM_NLF,Dtfil%fnameabo_em1_nlf) 1358 end if 1359 end if ! master and is_qeq0==1 1360 1361 if (is_qeq0==1) call chi_free(chihw) 1362 1363 call spectra_free(spectra) 1364 if (allocated(kxcg)) then 1365 ABI_FREE(kxcg) 1366 end if 1367 if (allocated(fxc_ADA)) then 1368 ABI_FREE(fxc_ADA) 1369 end if 1370 ! 1371 ! Output the sum rule evaluation. 1372 ! Vcp%vc_sqrt(:,iqibz) Contains vc^{1/2}(q,G), complex-valued due to a possible cutoff 1373 epsm1 => chi0 1374 call output_chi0sumrule((is_qeq0==1),iqibz,Ep%npwe,omegaplasma,chi0_sumrule,epsm1(:,:,1),Vcp%vc_sqrt(:,iqibz)) 1375 1376 ! If input variable npvel is larger than 0, trigger the Random Stopping Power calculation 1377 ! Only the masternode is used 1378 if (my_rank==master .and. Dtset%npvel>0) then 1379 if (is_first_qcalc) then 1380 ABI_MALLOC(rspower,(Dtset%npvel)) 1381 rspower(:)=zero 1382 end if 1383 call random_stopping_power(iqibz,Dtset%npvel,Dtset%pvelmax,Ep,Gsph_epsG0,Qmesh,Vcp,Cryst,Dtfil,epsm1,rspower) 1384 if (label==Ep%nqcalc) then 1385 ABI_FREE(rspower) 1386 end if 1387 end if 1388 1389 ! Write heads and wings to main output file. 1390 if (is_qeq0==1) then 1391 write(msg,'(1x,2a)')' Heads and wings of the symmetrical epsilon^-1(G,G'') ',ch10 1392 call wrtout(ab_out,msg,'COLL') 1393 do iomega=1,Ep%nomega 1394 write(msg,'(2x,a,i4,a,2f9.4,a)')' Upper and lower wings at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]' 1395 call wrtout(ab_out,msg,'COLL') 1396 call print_arr(epsm1(1,:,iomega),max_r=9,unit=ab_out) 1397 call print_arr(epsm1(:,1,iomega),max_r=9,unit=ab_out) 1398 call wrtout(ab_out,ch10,'COLL') 1399 end do 1400 end if 1401 1402 call timab(309,2,tsec) 1403 call timab(310,1,tsec) ! wrscr 1404 1405 if (my_rank==master) then 1406 ! === Write the symmetrical epsilon^-1 on file === 1407 title(1)='SCR file: epsilon^-1' 1408 if (is_qeq0==1) then 1409 string='0'; if (Dtset%usepaw==0.and.Ep%inclvkb/=0) call int2char10(Ep%inclvkb,string) 1410 title(1)=title(1)(1:21)//', calculated using inclvkb = '//string 1411 end if 1412 title(2)='TESTPARTICLE' 1413 ctype='RPA' 1414 title(2)(14:17)=ctype !this has to be modified 1415 1416 if (is_first_qcalc) then 1417 ! === Open file and write the header for the SCR file === 1418 ! * Here we write the RPA approximation for \tilde\epsilon^{-1} 1419 ikxc=0; test_type=0; tordering=1 1420 hem1 = hscr_new("inverse_dielectric_function",dtset,ep,hdr_local,ikxc,test_type,tordering,title,& 1421 & Ep%npwe,Gsph_epsG0%gvec) 1422 fform_em1 = hem1%fform 1423 if (dtset%iomode == IO_MODE_ETSF) then 1424 #ifdef HAVE_NETCDF 1425 NCF_CHECK(nctk_open_create(unt_em1, nctk_ncify(dtfil%fnameabo_scr), xmpi_comm_self)) 1426 NCF_CHECK(crystal_ncwrite(cryst, unt_em1)) 1427 NCF_CHECK(ebands_ncwrite(QP_BSt, unt_em1)) 1428 #endif 1429 else 1430 unt_em1=Dtfil%unscr 1431 if (open_file(dtfil%fnameabo_scr,msg,unit=unt_em1,status='unknown',form='unformatted') /= 0) then 1432 MSG_ERROR(msg) 1433 end if 1434 end if 1435 call hscr_io(hem1,fform_em1,2,unt_em1,xmpi_comm_self,0,Dtset%iomode) 1436 call hscr_free(Hem1) 1437 end if 1438 1439 call write_screening("inverse_dielectric_function",unt_em1,Dtset%iomode,Ep%npwe,Ep%nomega,iqcalc,epsm1) 1440 end if ! my_rank==master 1441 1442 call timab(310,2,tsec) 1443 end do ! Loop over q-points 1444 1445 !----------------------------- END OF THE CALCULATION ------------------------ 1446 1447 ! Close Files. 1448 if (my_rank==master) then 1449 if (dtset%iomode == IO_MODE_ETSF) then 1450 #ifdef HAVE_NETCDF 1451 NCF_CHECK(nf90_close(unt_em1)) 1452 if (dtset%prtsuscep>0) then 1453 NCF_CHECK(nf90_close(unt_susc)) 1454 end if 1455 #endif 1456 else 1457 close(unt_em1) 1458 if (dtset%prtsuscep>0) close(unt_susc) 1459 end if 1460 end if 1461 ! 1462 !===================== 1463 !==== Free memory ==== 1464 !===================== 1465 ABI_FREE(chi0_sumrule) 1466 ABI_FREE(chi0) 1467 if (allocated(rhor_kernel)) then 1468 ABI_FREE(rhor_kernel) 1469 end if 1470 1471 ABI_FREE(rhor) 1472 ABI_FREE(rhog) 1473 ABI_FREE(ks_vbik) 1474 ABI_FREE(qp_vbik) 1475 ABI_FREE(ktabr) 1476 ABI_FREE(taur) 1477 ABI_FREE(taug) 1478 ABI_FREE(ks_vhartr) 1479 ABI_FREE(ks_vtrial) 1480 ABI_FREE(vpsp) 1481 ABI_FREE(ks_vxc) 1482 ABI_FREE(ph1d) 1483 ABI_FREE(ph1df) 1484 ABI_FREE(nhatgr) 1485 ABI_FREE(nhat) 1486 call pawfgr_destroy(Pawfgr) 1487 1488 if (Dtset%usepaw==1) then ! Optional deallocation for PAW. 1489 call pawrhoij_free(Pawrhoij) 1490 call pawfgrtab_free(Pawfgrtab) 1491 call paw_ij_free(Paw_ij) 1492 call paw_an_free(Paw_an) 1493 call pawpwff_free(Paw_pwff) 1494 if (Dtset%pawcross==1) then 1495 call paw_pwaves_lmn_free(Paw_onsite) 1496 Wfdf%bks_comm = xmpi_comm_null 1497 call wfd_free(Wfdf) 1498 end if 1499 end if 1500 1501 ABI_DT_FREE(Pawfgrtab) 1502 ABI_DT_FREE(Paw_pwff) 1503 ABI_DT_FREE(Pawrhoij) 1504 ABI_DT_FREE(Paw_ij) 1505 ABI_DT_FREE(Paw_an) 1506 ABI_FREE(ktabrf) 1507 ABI_DT_FREE(Paw_onsite) 1508 1509 call wfd_free(Wfd) 1510 call kmesh_free(Kmesh) 1511 call kmesh_free(Qmesh) 1512 call crystal_free(Cryst) 1513 call gsph_free(Gsph_epsG0) 1514 call gsph_free(Gsph_wfn) 1515 call vcoul_free(Vcp) 1516 call em1params_free(Ep) 1517 call hdr_free(Hdr_wfk) 1518 call hdr_free(Hdr_local) 1519 call ebands_free(KS_BSt) 1520 call ebands_free(QP_BSt) 1521 call littlegroup_free(Ltg_q) 1522 call destroy_mpi_enreg(MPI_enreg_seq) 1523 ABI_DT_FREE(Ltg_q) 1524 1525 call timab(301,2,tsec) 1526 1527 DBG_EXIT("COLL") 1528 1529 end subroutine screening