Line data Source code
1 : !!****m* ABINIT/m_screening_driver
2 : !! NAME
3 : !! m_screening_driver
4 : !!
5 : !! FUNCTION
6 : !! Calculate screening and dielectric functions
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG, GMR, VO, LR, RWG, MT, RShaltaf, AS, FB)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_screening_driver
23 :
24 : use defs_basis
25 : use defs_wvltypes
26 : use m_abicore
27 : use m_dtset
28 : use m_xmpi
29 : use m_xomp
30 : use m_errors
31 : use m_abi_mixing
32 : use m_kxc
33 : use m_nctk
34 : use netcdf
35 : use libxc_functionals
36 : use m_hdr
37 : use m_dtfil
38 : use m_crystal
39 :
40 : use defs_datatypes, only : pseudopotential_type
41 : use defs_abitypes, only : MPI_type
42 : use m_time, only : timab
43 : use m_io_tools, only : open_file, file_exists, iomode_from_fname
44 : use m_fstrings, only : int2char10, sjoin, strcat, itoa, ltoa, itoa
45 : use m_energies, only : energies_type
46 : use m_numeric_tools, only : print_arr, coeffs_gausslegint, c2r
47 : use m_geometry, only : normv, vdotw, mkrdim, metric
48 : use m_gwdefs, only : GW_TOLQ0, GW_TOLQ, em1params_t, GW_Q0_DEFAULT
49 : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
50 : use m_ebands, only : ebands_t, gaps_t
51 : use m_bz_mesh, only : kmesh_t, littlegroup_t, littlegroup_free, get_ng0sh
52 : use m_kg, only : getph
53 : use m_gsphere, only : gsphere_t, setshells
54 : use m_vcoul, only : vcoul_t
55 : use m_qparticles, only : rdqps, rdgw, show_QP
56 : use m_screening, only : make_epsm1_driver, lwl_write, chi_t, chi_free, chi_new
57 : use m_io_screening, only : hscr_new, write_screening, hscr_t
58 : use m_spectra, only : spectra_t, W_EM_LF, W_EM_NLF, W_EELF
59 : use m_fftcore, only : print_ngfft
60 : use m_fft_mesh, only : rotate_FFT_mesh, cigfft, get_gfft, setmesh
61 : use m_fft, only : fourdp
62 : use m_wfd, only : wfdgw_t, wfdgw_copy, test_charge
63 : use m_wfk, only : wfk_read_eigenvalues
64 : use m_io_kss, only : make_gvec_kss
65 : use m_chi0tk, only : output_chi0sumrule
66 : use m_pawang, only : pawang_type
67 : use m_pawrad, only : pawrad_type
68 : use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize
69 : use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
70 : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
71 : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
72 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy,&
73 : & pawrhoij_free, pawrhoij_symrhoij, pawrhoij_inquire_dim
74 : use m_pawdij, only : pawdij, symdij_all
75 : use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
76 : use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free
77 : use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy
78 : use m_paw_sphharm, only : setsym_ylm
79 : use m_paw_onsite, only : pawnabla_init
80 : use m_paw_nhat, only : nhatgrid, pawmknhat
81 : use m_paw_denpot, only : pawdenpot
82 : use m_paw_init, only : pawinit, paw_gencond
83 : use m_paw_tools, only : chkpawovlp,pawprt
84 : use m_chi0, only : cchi0, cchi0q0, chi0q0_intraband
85 : use m_setvtr, only : setvtr
86 : use m_mkrho, only : prtrhomxmn
87 : use m_pspini, only : pspini
88 : use m_paw_correlations, only : pawpuxinit
89 : use m_plowannier, only : plowannier_type,init_plowannier,get_plowannier, fullbz_plowannier,destroy_plowannier
90 : use minimax_grids, only : gx_minimax_grid !, gx_get_error_message
91 : use m_pstat, only : pstat_proc
92 :
93 : implicit none
94 :
95 : private
96 : !!***
97 :
98 : public :: screening
99 : !!***
100 :
101 : contains
102 : !!***
103 :
104 : !!****f* m_screening_driver/screening
105 : !! NAME
106 : !! screening
107 : !!
108 : !! FUNCTION
109 : !! Calculate screening and dielectric functions
110 : !!
111 : !! INPUTS
112 : !! acell(3)=length scales of primitive translations (bohr)
113 : !! codvsn=code version
114 : !! Dtfil<datafiles_type)>=variables related to file names and unit numbers.
115 : !! Pawang<pawang_type)>=paw angular mesh and related data
116 : !! Pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
117 : !! Pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
118 : !! Psps <type(pseudopotential_type)>=variables related to pseudopotentials
119 : !! Before entering the first time in screening, a significant part of Psps has been initialized:
120 : !! the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid, ntypat,n1xccc,usepaw,useylm,
121 : !! and the arrays dimensioned to npsp. All the remaining components of Psps are to be initialized in
122 : !! the call to pspini. The next time the code enters screening, Psps might be identical to the
123 : !! one of the previous Dtset, in which case, no reinitialisation is scheduled in pspini.F90.
124 : !! rprim(3,3)=dimensionless real space primitive translations
125 : !!
126 : !! OUTPUT
127 : !! Output is written on the main output file.
128 : !! The symmetrical inverse dielectric matrix is stored in the _SCR file
129 : !!
130 : !! SIDE EFFECTS
131 : !! Dtset<type(dataset_type)>=all input variables for this dataset
132 : !!
133 : !! NOTES
134 : !! USE OF FFT GRIDS:
135 : !! =================
136 : !! In case of PAW:
137 : !! ---------------
138 : !! Two FFT grids are used:
139 : !! - A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
140 : !! It is defined by nfft, ngfft, mgfft, ...
141 : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
142 : !! - A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
143 : !! It is defined by nfftf, ngfftf, mgfftf, ...Total density, potentials, ... are expressed on this grid.
144 : !! In case of norm-conserving:
145 : !! ---------------------------
146 : !! - Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
147 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
148 : !!
149 : !! SOURCE
150 :
151 148 : subroutine screening(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim)
152 :
153 : !Arguments ------------------------------------
154 : !scalars
155 : character(len=8),intent(in) :: codvsn
156 : type(Datafiles_type),intent(in) :: Dtfil
157 : type(Dataset_type),intent(inout) :: Dtset
158 : type(Pawang_type),intent(inout) :: Pawang
159 : type(Pseudopotential_type),intent(inout) :: Psps
160 : !arrays
161 : real(dp),intent(in) :: acell(3),rprim(3,3)
162 : type(Pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Dtset%usepaw)
163 : type(Pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Dtset%usepaw)
164 :
165 : !Local variables ------------------------------
166 : character(len=4) :: ctype='RPA '
167 : !scalars
168 : integer,parameter :: tim_fourdp4=4,NOMEGA_PRINTED=15,master=0
169 : integer :: spin,ik_ibz,my_nbks
170 : integer :: choice,cplex,cplex_rhoij,dim_kxcg,dim_wing,ount,omp_ncpus
171 : integer :: fform_chi0,fform_em1,gnt_option,iat,ider,idir,ierr,band
172 : integer :: ifft,ii,ikbz,ikxc,initialized,iomega,ios,ipert
173 : integer :: iqibz,iqcalc,is_qeq0,isym,izero,ifirst,ilast
174 : integer :: label,mgfftf,mgfftgw
175 : integer :: nt_per_proc,work_size
176 : integer :: moved_atm_inside,moved_rhor
177 : integer :: nbcw,nbsc,nbvw,nkxc,nkxc1,n3xccc,optene,istep
178 : integer :: nfftf,nfftf_tot,nfftgw,nfftgw_tot,ngrvdw,nhatgrdim,nprocs,nspden_rhoij
179 : integer :: nscf,nzlmopt,mband
180 : integer :: optcut,optgr0,optgr1,optgr2,option,approx_type,option_test,optgrad
181 : integer :: optrad,optrhoij,psp_gencond,my_rank, ig
182 : integer :: rhoxsp_method,comm,test_type,tordering,unt_em1,unt_susc,usexcnhat, ncerr, gap_err
183 : real(dp) :: compch_fft,compch_sph,domegareal,e0,ecore,ecut_eff,ecutdg_eff
184 : real(dp) :: gsqcutc_eff,gsqcutf_eff,gsqcut_shp,omegaplasma,ucvol,vxcavg,gw_gsq,r_s
185 : real(dp) :: alpha,rhoav,factor,ec_gm,el_temp
186 : real(dp):: eff,mempercpu_mb,max_wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem
187 : integer, parameter :: epsm1_win = xmpi_undefined
188 : logical :: found,iscompatibleFFT,is_dfpt=.false.,use_tr,is_first_qcalc
189 : logical :: add_chi0_intraband,update_energies,call_pawinit
190 : character(len=10) :: string
191 : character(len=500) :: msg
192 : character(len=80) :: bar
193 296 : type(ebands_t) :: ks_ebands, qp_ebands
194 3848 : type(kmesh_t) :: Kmesh,Qmesh
195 3256 : type(vcoul_t) :: Vcp
196 7696 : type(crystal_t) :: Cryst
197 148 : type(em1params_t) :: Ep
198 : type(Energies_type) :: KS_energies
199 148 : type(gsphere_t) :: Gsph_epsG0,Gsph_wfn
200 148 : type(Hdr_type) :: Hdr_wfk,Hdr_local
201 148 : type(MPI_type) :: MPI_enreg_seq
202 : type(Pawfgr_type) :: Pawfgr
203 1036 : type(hscr_t) :: Hem1,Hchi0
204 148 : type(wfdgw_t) :: Wfd,Wfdf
205 148 : type(spectra_t) :: spectra
206 : !type(chi_t) :: chihw
207 148 : type(wvl_data) :: wvl_dummy
208 : character(len=nctk_slen) :: wing_shape
209 : !arrays
210 296 : integer :: ibocc(Dtset%nsppol),ngfft_gw(18),ngfftc(18),ngfftf(18), units(2)
211 296 : integer,allocatable :: irottb(:,:),ktabr(:,:),ktabrf(:,:),l_size_atm(:)
212 296 : integer,allocatable :: ks_vbik(:,:),ks_occ_idx(:,:),qp_vbik(:,:),nband(:,:)
213 296 : integer,allocatable :: nq_spl(:),nlmn_atm(:),gw_gfft(:,:)
214 : real(dp) :: gmet(3,3),gprimd(3,3),k0(3),qtmp(3),rmet(3,3),rprimd(3,3),tsec(2),strsxc(6)
215 592 : real(dp),allocatable :: igwene(:,:,:),chi0_sumrule(:),ec_rpa(:),rspower(:)
216 148 : real(dp),allocatable :: nhat(:,:),nhatgr(:,:,:),ph1d(:,:),ph1df(:,:)
217 148 : real(dp),allocatable :: rhog(:,:),rhor(:,:),rhor_p(:,:),rhor_kernel(:,:),taur(:,:)
218 444 : real(dp),allocatable :: z(:),zw(:),grchempottn(:,:),grewtn(:,:),grvdw(:,:),kxc(:,:),qmax(:)
219 296 : real(dp),allocatable :: ks_vhartr(:),vpsp(:),ks_vtrial(:,:),ks_vxc(:,:),xccc3d(:)
220 444 : complex(gwp),allocatable :: arr_99(:,:),kxcg(:,:),fxc_ADA(:,:,:)
221 148 : complex(dp),allocatable :: m_ks_to_qp(:,:,:,:)
222 148 : complex(dp),allocatable :: chi0_head(:,:,:), chi0_lwing(:,:,:), chi0_uwing(:,:,:)
223 148 : real(dp),allocatable :: rwork_wing(:,:,:,:), tmp_omega_wgs(:)
224 296 : complex(dp),allocatable :: chi0intra_lwing(:,:,:),chi0intra_uwing(:,:,:),chi0intra_head(:,:,:), tmp_omega(:)
225 148 : complex(gwp),allocatable,target :: chi0(:,:,:),chi0intra(:,:,:)
226 148 : complex(gwp),contiguous, pointer :: epsm1(:,:,:)
227 148 : logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:)
228 : character(len=80) :: title(2)
229 : character(len=fnlen) :: gw_fname,wfk_fname,lwl_fname
230 148 : type(littlegroup_t),pointer :: Ltg_q(:)
231 148 : type(Paw_an_type),allocatable :: Paw_an(:)
232 148 : type(Paw_ij_type),allocatable :: Paw_ij(:)
233 148 : type(Pawfgrtab_type),allocatable :: Pawfgrtab(:)
234 148 : type(Pawrhoij_type),allocatable :: Pawrhoij(:),prev_Pawrhoij(:)
235 148 : type(pawpwff_t),allocatable :: Paw_pwff(:)
236 148 : type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
237 148 : type(plowannier_type) :: wanbz,wanibz,wanibz_in
238 : real(dp) :: te_min, te_max
239 148 : type(gaps_t) :: gaps
240 148 : real(dp),allocatable :: tau_mesh(:), tau_wgs(:), iw_mesh(:), iw_wgs(:)
241 148 : real(dp),allocatable :: t2w_cos_wgs(:,:), w2t_cos_wgs(:,:), t2w_sin_wgs(:,:)
242 : real(dp) :: ft_max_error(3), cosft_duality_error
243 : !************************************************************************
244 :
245 148 : call timab(301,1,tsec) ! overall time
246 148 : call timab(302,1,tsec) ! screening(init
247 444 : units = [std_out, ab_out]
248 :
249 : write(msg,'(6a)')&
250 148 : ' SCREENING: Calculation of the susceptibility and dielectric matrices ',ch10,ch10,&
251 148 : ' Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining.',ch10,&
252 296 : ' Incorporated in ABINIT by V. Olevano, G.-M. Rignanese, and M. Torrent.'
253 148 : call wrtout(units, msg)
254 :
255 148 : if(dtset%ucrpa>0) then
256 0 : write(msg,'(6a)')ch10,&
257 0 : ' cRPA Calculation: The calculation of the polarisability is constrained (ucrpa/=0)',ch10
258 0 : call wrtout(units, msg)
259 : end if
260 : #if defined HAVE_GW_DPC
261 : if (gwp/=8) then
262 : write(msg,'(6a)')ch10,&
263 : ' Number of bytes for double precision complex /=8 ',ch10,&
264 : ' Cannot continue due to kind mismatch in BLAS library ',ch10,&
265 : ' Some BLAS interfaces are not generated by abilint '
266 : ABI_ERROR(msg)
267 : end if
268 148 : write(msg,'(a,i2,a)')'.Using double precision arithmetic ; gwpc = ',gwp,ch10
269 : #else
270 : write(msg,'(a,i2,a)')'.Using single precision arithmetic ; gwpc = ',gwp,ch10
271 : #endif
272 148 : call wrtout(units, msg)
273 :
274 : ! === Initialize MPI variables, and parallelization level ===
275 : ! gwpara: 0--> sequential run, 1--> parallelism over k-points, 2--> parallelism over bands.
276 : ! gwpara==2, each node has both fully and partially occupied states while conduction bands are divided
277 148 : comm = xmpi_world; my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
278 :
279 148 : if (my_rank == master) then
280 120 : wfk_fname = dtfil%fnamewffk
281 120 : if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then
282 0 : ABI_ERROR(msg)
283 : end if
284 : end if
285 148 : call xmpi_bcast(wfk_fname, master, comm, ierr)
286 :
287 : ! Some variables need to be initialized/nullify at start
288 148 : call KS_energies%init()
289 148 : usexcnhat=0
290 :
291 148 : call mkrdim(acell,rprim,rprimd)
292 148 : call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol)
293 :
294 : !=== Define FFT grid(s) sizes ===
295 : ! Be careful! This mesh is only used for densities and potentials. It is NOT the (usually coarser)
296 : ! GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90.
297 : ! See also NOTES in the comments at the beginning of this file.
298 : ! NOTE: The mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn.
299 :
300 148 : k0(:)=zero
301 : call pawfgr_init(Pawfgr,Dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,&
302 148 : gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=gmet,k0=k0)
303 :
304 296 : call print_ngfft([std_out], ngfftf,'Dense FFT mesh used for densities and potentials')
305 592 : nfftf_tot=PRODUCT(ngfftf(1:3))
306 :
307 : ! We can intialize MPI_enreg and fft distrib here, now ngfft are known
308 148 : call initmpi_seq(MPI_enreg_seq) ! Fake MPI_type for the sequential part.
309 148 : call MPI_enreg_seq%distribfft%init_seq('c',ngfftc(2),ngfftc(3),'all')
310 148 : call MPI_enreg_seq%distribfft%init_seq('f',ngfftf(2),ngfftf(3),'all')
311 :
312 : !=============================================
313 : !==== Open and read pseudopotential files ====
314 : !=============================================
315 148 : call pspini(Dtset,Dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,Pawrad,Pawtab,Psps,rprimd,comm_mpi=comm)
316 :
317 : ! === Initialize dimensions and basic objects ===
318 : call setup_screening(codvsn,acell,rprim,wfk_fname,Dtset,Psps,Pawtab,&
319 148 : ngfft_gw,Hdr_wfk,Hdr_local,Cryst,Kmesh,Qmesh,ks_ebands,Ltg_q,Gsph_epsG0,Gsph_wfn,Vcp,Ep,comm)
320 :
321 148 : call pstat_proc%print(_PSTAT_ARGS_)
322 148 : call timab(302,2,tsec) ! screening(init)
323 296 : call print_ngfft([std_out], ngfft_gw, header='FFT mesh used for oscillator strengths')
324 :
325 592 : nfftgw_tot=PRODUCT(ngfft_gw(1:3))
326 : mgfftgw =MAXVAL (ngfft_gw(1:3))
327 148 : nfftgw =nfftgw_tot ! no FFT //
328 :
329 : !TRYING TO RECREATE AN "ABINIT ENVIRONMENT"
330 148 : KS_energies%e_corepsp=ecore/Cryst%ucvol
331 :
332 : !==========================
333 : !=== PAW initialization ===
334 : !==========================
335 148 : if (Dtset%usepaw==1) then
336 5 : call timab(315,1,tsec) ! screening(pawin
337 :
338 5 : call chkpawovlp(Cryst%natom,Cryst%ntypat,Dtset%pawovlp,Pawtab,Cryst%rmet,Cryst%typat,Cryst%xred)
339 :
340 46 : ABI_MALLOC(Pawrhoij,(Cryst%natom))
341 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
342 5 : & nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc)
343 : call pawrhoij_alloc(Pawrhoij,cplex_rhoij,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,&
344 5 : & Cryst%typat,pawtab=Pawtab)
345 :
346 : ! Initialize values for several basic arrays stored in Pawinit
347 5 : gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2
348 :
349 : ! Test if we have to call pawinit
350 5 : call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
351 :
352 5 : if (psp_gencond==1.or.call_pawinit) then
353 1 : gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2
354 : call pawinit(dtset%effmass_free,gnt_option,gsqcut_shp,zero,Dtset%pawlcutd,Dtset%pawlmix,&
355 : & Psps%mpsang,Dtset%pawnphi,Cryst%nsym,Dtset%pawntheta,Pawang,Pawrad,&
356 1 : & Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%ixc,Dtset%usepotzero)
357 :
358 : ! Update internal values
359 1 : call paw_gencond(Dtset,gnt_option,"save",call_pawinit)
360 : else
361 4 : if (Pawtab(1)%has_kij ==1) Pawtab(1:Cryst%ntypat)%has_kij =2
362 4 : if (Pawtab(1)%has_nabla==1) Pawtab(1:Cryst%ntypat)%has_nabla=2
363 : end if
364 14 : Psps%n1xccc=MAXVAL(Pawtab(1:Cryst%ntypat)%usetcore)
365 :
366 : ! Initialize optional flags in Pawtab to zero
367 : ! (Cannot be done in Pawinit since the routine is called only if some pars. are changed)
368 14 : Pawtab(:)%has_nabla = 0
369 14 : Pawtab(:)%usepawu = 0
370 14 : Pawtab(:)%useexexch = 0
371 14 : Pawtab(:)%exchmix =zero
372 14 : Pawtab(:)%lamb_shielding = 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 5 : call pawnabla_init(Psps%mpsang,Cryst%ntypat,Pawrad,Pawtab)
377 :
378 5 : call setsym_ylm(gprimd,Pawang%l_max-1,Cryst%nsym,Dtset%pawprtvol,rprimd,Cryst%symrec,Pawang%zarot)
379 :
380 : ! Initialize and compute data for DFT+U.
381 : ! paw_dmft%use_dmft=dtset%usedmft
382 : call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,&
383 : & is_dfpt,Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,dtset%nspinor,Cryst%ntypat,dtset%optdcmagpawu,Pawang,Dtset%pawprtvol,&
384 5 : & Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa)
385 :
386 5 : if (my_rank == master) call pawtab_print(Pawtab)
387 :
388 : ! Get Pawrhoij from the header of the WFK file.
389 5 : call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij)
390 :
391 : ! Re-symmetrize rhoij.
392 : ! this call leads to a SIGFAULT, likely some pointer is not initialized correctly
393 5 : choice=1; optrhoij=1; ipert=0; idir=0
394 : !call pawrhoij_symrhoij(Pawrhoij,Pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert,Cryst%natom,&
395 : !& Cryst%nsym,Cryst%ntypat,optrhoij,Pawang,Dtset%pawprtvol,Pawtab,&
396 : !& Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat)
397 : !
398 : ! Evaluate form factors for the radial part of phi.phj-tphi.tphj ===
399 : ! rhoxsp_method=1 ! Arnaud-Alouani
400 : ! rhoxsp_method=2 ! Shiskin-Kresse
401 5 : rhoxsp_method=2
402 :
403 : ! At least for ucrpa, the Arnaud Alouani is always a better choice but needs a larger cutoff
404 5 : if(dtset%ucrpa>0) rhoxsp_method=1
405 5 : if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc
406 :
407 15 : ABI_MALLOC(gw_gfft,(3,nfftgw_tot))
408 5 : call get_gfft(ngfft_gw,(/zero,zero,zero/),gmet,gw_gsq,gw_gfft)
409 5 : ABI_FREE(gw_gfft)
410 :
411 : ! Set up q grids, make qmax 20% larger than largest expected:
412 15 : ABI_MALLOC(nq_spl,(Psps%ntypat))
413 15 : ABI_MALLOC(qmax,(Psps%ntypat))
414 14 : nq_spl = Psps%mqgrid_ff
415 14 : qmax = SQRT(gw_gsq)*1.2d0 !qmax = Psps%qgrid_ff(Psps%mqgrid_ff)
416 84 : ABI_MALLOC(Paw_pwff,(Psps%ntypat))
417 :
418 5 : call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps)
419 5 : ABI_FREE(nq_spl)
420 5 : ABI_FREE(qmax)
421 :
422 : ! Variables/arrays related to the fine FFT grid
423 20 : ABI_MALLOC(nhat,(nfftf,Dtset%nspden))
424 132482 : nhat=zero; cplex=1
425 46 : ABI_MALLOC(Pawfgrtab,(Cryst%natom))
426 5 : call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat)
427 5 : call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Dtset%nspden,Dtset%typat)
428 5 : ABI_FREE(l_size_atm)
429 5 : compch_fft=greatest_real
430 14 : usexcnhat=MAXVAL(Pawtab(:)%usexcnhat)
431 : ! * 0 --> Vloc in atomic data is Vbare (Blochl s formulation)
432 : ! * 1 --> Vloc in atomic data is VH(tnzc) (Kresse s formulation)
433 5 : write(msg,'(a,i3)')' screening : using usexcnhat = ',usexcnhat
434 5 : call wrtout(std_out, msg)
435 :
436 : ! Identify parts of the rectangular grid where the density has to be calculated.
437 5 : optcut=0;optgr0=Dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-Dtset%pawstgylm
438 5 : if (Dtset%pawcross==1) optrad=1
439 5 : if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm
440 :
441 : call nhatgrid(Cryst%atindx1,gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,&
442 5 : optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
443 :
444 15 : call timab(315,2,tsec) ! screening(pawin
445 : else
446 : ! allocate empty structure for the sake of -fcheck=all...
447 143 : ABI_MALLOC(Paw_pwff,(0))
448 143 : ABI_MALLOC(Pawrhoij,(0))
449 143 : ABI_MALLOC(Pawfgrtab,(0))
450 : end if ! End of PAW initialization.
451 :
452 : ! Consistency check and additional stuff done only for GW with PAW.
453 795 : ABI_MALLOC(Paw_onsite,(Cryst%natom))
454 148 : if (Dtset%usepaw==1) then
455 5 : if (Dtset%ecutwfn < Dtset%ecut) then
456 : write(msg,"(5a)")&
457 0 : "WARNING - ",ch10,&
458 0 : " It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,&
459 0 : " an excessive truncation of the planewave basis set can lead to unphysical results."
460 0 : call wrtout(ab_out, msg)
461 : end if
462 :
463 5 : ABI_CHECK(Dtset%useexexch==0,"LEXX not yet implemented in GW")
464 5 : ABI_CHECK(Dtset%usedmft==0,"DMFT + GW not available")
465 :
466 5 : if (Dtset%pawcross==1) then
467 0 : optgrad=1
468 : call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,Cryst%rprimd,&
469 0 : Cryst%xcart,Pawtab,Pawrad,Pawfgrtab,optgrad)
470 : end if
471 : end if
472 :
473 : ! Allocate these arrays anyway, since they are passed to subroutines.
474 148 : if (.not.allocated(nhat)) then
475 286 : ABI_MALLOC(nhat,(nfftf,0))
476 : end if
477 :
478 148 : call timab(316,1,tsec) ! screening(wfs
479 :
480 : !=====================================================
481 : !=== Prepare the distribution of the wavefunctions ===
482 : !=====================================================
483 : ! valence and partially occupied are replicate on each node while conduction bands are MPI distributed.
484 : ! This method is mandatory if gwpara==2 and/or we are using awtr==1 or the spectral method.
485 : ! If awtr==1, we evaluate chi0 taking advantage of time-reversal (speed-up~2)
486 : ! Useful indices:
487 : ! nbvw = Max. number of fully/partially occupied states over spin
488 : ! nbcw = Max. number of unoccupied states considering the spin
489 : !TODO:
490 : ! Here for semiconducting systems we have to be sure that each processor has all the
491 : ! states considered in the SCGW, moreover nbsc<nbvw
492 : ! in case of SCGW vale and conduction has to be recalculated to avoid errors
493 : ! if a metal becomes semiconductor or viceversa.
494 : ! Ideally nbvw should include only the states v such that the transition
495 : ! c-->v is taken into account in cchi0 (see GW_TOLDOCC). In the present implementation
496 :
497 592 : ABI_MALLOC(ks_occ_idx,(ks_ebands%nkpt, ks_ebands%nsppol))
498 444 : ABI_MALLOC(ks_vbik ,(ks_ebands%nkpt, ks_ebands%nsppol))
499 444 : ABI_MALLOC(qp_vbik ,(ks_ebands%nkpt, ks_ebands%nsppol))
500 :
501 148 : call ks_ebands%update_occ(Dtset%spinmagntarget, prtvol=0)
502 148 : ks_occ_idx = ks_ebands%get_occupied(tol8) ! tol8 to be consistent when the density
503 148 : ks_vbik = ks_ebands%get_valence_idx()
504 :
505 148 : ibocc(:)=MAXVAL(ks_occ_idx(:,:),DIM=1) ! Max occupied band index for each spin.
506 148 : ABI_FREE(ks_occ_idx)
507 :
508 148 : use_tr =.FALSE.; nbvw=0
509 148 : if (Dtset%gwpara==2.or.Ep%awtr==1.or.Dtset%spmeth>0) then
510 137 : use_tr = .TRUE.
511 276 : nbvw=MAXVAL(ibocc)
512 137 : nbcw=Ep%nbnds-nbvw
513 137 : write(msg,'(4a,i0,2a,i0,2a,i0,a)')ch10,&
514 137 : '- screening: taking advantage of time-reversal symmetry ',ch10,&
515 137 : '- Maximum band index for partially occupied states nbvw = ',nbvw,ch10,&
516 137 : '- Remaining bands to be divided among processors nbcw = ',nbcw,ch10,&
517 274 : '- Number of bands treated by each node ~',nbcw/nprocs,ch10
518 137 : call wrtout(ab_out, msg)
519 137 : if (Cryst%timrev/=2) then
520 0 : ABI_ERROR('Time-reversal cannot be used since cryst%timrev/=2')
521 : end if
522 : end if
523 :
524 148 : mband=Ep%nbnds
525 592 : ABI_MALLOC(nband,(Kmesh%nibz,Dtset%nsppol))
526 1159 : nband=mband
527 740 : ABI_MALLOC(bks_mask,(mband,Kmesh%nibz,Dtset%nsppol))
528 592 : ABI_MALLOC(keep_ur,(mband,Kmesh%nibz,Dtset%nsppol))
529 46746 : bks_mask=.FALSE.; keep_ur=.FALSE.
530 :
531 : ! autoparal section
532 148 : if (dtset%max_ncpus /=0) then
533 0 : ount = ab_out
534 : ! Temporary table needed to estimate memory
535 0 : ABI_MALLOC(nlmn_atm,(Cryst%natom))
536 0 : if (Dtset%usepaw==1) then
537 0 : do iat=1,Cryst%natom
538 0 : nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size
539 : end do
540 : end if
541 :
542 0 : write(ount,'(a)')"--- !Autoparal"
543 0 : write(ount,"(a)")"# Autoparal section for Screening runs"
544 0 : write(ount,"(a)") "info:"
545 0 : write(ount,"(a,i0)")" autoparal: ",dtset%autoparal
546 0 : write(ount,"(a,i0)")" max_ncpus: ",dtset%max_ncpus
547 0 : write(ount,"(a,i0)")" gwpara: ",dtset%gwpara
548 0 : write(ount,"(a,i0)")" nkpt: ",dtset%nkpt
549 0 : write(ount,"(a,i0)")" nsppol: ",dtset%nsppol
550 0 : write(ount,"(a,i0)")" nspinor: ",dtset%nspinor
551 0 : write(ount,"(a,i0)")" nbnds: ",Ep%nbnds
552 :
553 0 : work_size = nbvw * nbcw * Kmesh%nibz**2 * Dtset%nsppol
554 :
555 : ! Non-scalable memory in Mb i.e. memory that is not distribute with MPI.
556 0 : nonscal_mem = (two*gwp*Ep%npwe**2*(Ep%nomega*b2Mb)) * 1.1_dp
557 :
558 : ! List of configurations.
559 : ! Assuming an OpenMP implementation with perfect speedup!
560 0 : write(ount,"(a)")"configurations:"
561 :
562 0 : do ii=1,dtset%max_ncpus
563 0 : nt_per_proc = 0
564 0 : eff = HUGE(one)
565 0 : max_wfsmem_mb = zero
566 :
567 0 : do my_rank=0,ii-1
568 0 : call chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,ii,bks_mask,keep_ur,ierr)
569 0 : if (ierr /= 0) exit
570 0 : nt_per_proc = MAX(nt_per_proc, COUNT(bks_mask(1:nbvw,:,:)) * COUNT(bks_mask(nbvw+1:,:,:)))
571 0 : eff = MIN(eff, (one * work_size) / (ii * nt_per_proc))
572 :
573 : ! Memory needed for Fourier components ug.
574 0 : my_nbks = COUNT(bks_mask)
575 0 : ug_mem = two*gwp*Dtset%nspinor*Ep%npwwfn*my_nbks*b2Mb
576 :
577 : ! Memory needed for real space ur.
578 0 : ur_mem = two*gwp*Dtset%nspinor*nfftgw*COUNT(keep_ur)*b2Mb
579 :
580 : ! Memory needed for PAW projections Cprj
581 0 : cprj_mem = zero
582 0 : if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb
583 :
584 0 : max_wfsmem_mb = MAX(max_wfsmem_mb, ug_mem + ur_mem + cprj_mem)
585 : end do
586 0 : if (ierr /= 0) cycle
587 :
588 : ! Add the non-scalable part and increase by 10% to account for other datastructures.
589 0 : mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp
590 :
591 0 : do omp_ncpus=1,xomp_get_max_threads()
592 0 : write(ount,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus
593 0 : write(ount,"(a,i0)")" mpi_ncpus: ",ii
594 0 : write(ount,"(a,i0)")" omp_ncpus: ",omp_ncpus
595 0 : write(ount,"(a,f12.9)")" efficiency: ",eff
596 0 : write(ount,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb
597 : end do
598 : end do
599 0 : write(ount,'(a)')"..."
600 :
601 0 : ABI_FREE(nlmn_atm)
602 0 : ABI_ERROR_NODUMP("aborting now")
603 : else
604 148 : call chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,nprocs,bks_mask,keep_ur,ierr)
605 : end if
606 :
607 : ! Initialize the wf descriptor (allocate %ug and %ur if required).
608 :
609 : call wfd%init(Cryst,Pawtab,Psps,keep_ur,mband,nband,Ep%nkibz,Dtset%nsppol,bks_mask,&
610 : Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,ngfft_gw,&
611 148 : Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm)
612 :
613 148 : if (Dtset%pawcross==1) then
614 : call Wfdf%init(Cryst,Pawtab,Psps,keep_ur,mband,nband,Ep%nkibz,Dtset%nsppol,bks_mask,&
615 : Dtset%nspden,Dtset%nspinor,dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,ngfft_gw,&
616 0 : Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm)
617 : end if
618 :
619 148 : ABI_FREE(bks_mask)
620 148 : ABI_FREE(nband)
621 148 : ABI_FREE(keep_ur)
622 :
623 296 : call wfd%print([std_out])
624 : !FIXME: Rewrite the treatment of use_tr branches in cchi0 ...
625 : !Use a different nbvw for each spin.
626 : !Now use_tr means that one can use time-reversal symmetry.
627 :
628 : !==================================================
629 : !==== Read KS band structure from the KSS file ====
630 : !==================================================
631 148 : call wfd%read_wfk(wfk_fname,iomode_from_fname(wfk_fname))
632 :
633 148 : if (Dtset%pawcross==1) then
634 0 : call wfdgw_copy(Wfd, Wfdf)
635 0 : call wfdf%change_ngfft(Cryst,Psps,ngfftf)
636 : end if
637 :
638 : ! This test has been disabled (too expensive!)
639 : if (.False.) call wfd%test_ortho(Cryst,Pawtab,unit=ab_out,mode_paral="COLL")
640 :
641 148 : call pstat_proc%print(_PSTAT_ARGS_)
642 148 : call timab(316,2,tsec) ! screening(wfs
643 148 : call timab(319,1,tsec) ! screening(1)
644 :
645 148 : if (Cryst%nsym/=Dtset%nsym .and. Dtset%usepaw==1) then
646 0 : ABI_ERROR('Cryst%nsym/=Dtset%nsym, check pawinit and pawrhoij_symrhoij')
647 : end if
648 :
649 : ! Get the FFT index of $ (R^{-1}(r-\tau)) $
650 : ! S= $\transpose R^{-1}$ and k_BZ = S k_IBZ
651 : ! irottb is the FFT index of $ R^{-1} (r-\tau) $ used to symmetrize u_Sk.
652 592 : ABI_MALLOC(irottb, (nfftgw, Cryst%nsym))
653 148 : call rotate_FFT_mesh(Cryst%nsym, Cryst%symrel, Cryst%tnons, ngfft_gw, irottb, iscompatibleFFT)
654 :
655 592 : ABI_MALLOC(ktabr,(nfftgw,Kmesh%nbz))
656 7084 : do ikbz=1,Kmesh%nbz
657 6936 : isym=Kmesh%tabo(ikbz)
658 40401262 : do ifft=1,nfftgw
659 40401114 : ktabr(ifft,ikbz)=irottb(ifft,isym)
660 : end do
661 : end do
662 148 : ABI_FREE(irottb)
663 :
664 148 : if (Dtset%usepaw==1 .and. Dtset%pawcross==1) then
665 0 : ABI_MALLOC(irottb,(nfftf,Cryst%nsym))
666 0 : call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfftf,irottb,iscompatibleFFT)
667 :
668 0 : ABI_MALLOC(ktabrf,(nfftf,Kmesh%nbz))
669 0 : do ikbz=1,Kmesh%nbz
670 0 : isym=Kmesh%tabo(ikbz)
671 0 : do ifft=1,nfftf
672 0 : ktabrf(ifft,ikbz)=irottb(ifft,isym)
673 : end do
674 : end do
675 0 : ABI_FREE(irottb)
676 : else
677 148 : ABI_MALLOC(ktabrf,(0,0))
678 : end if
679 :
680 : ! Compute structure factor phases and large sphere cut-off
681 444 : ABI_MALLOC(ph1d,(2,3*(2*Dtset%mgfft+1)*Cryst%natom))
682 444 : ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*Cryst%natom))
683 148 : call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred)
684 :
685 148 : if (Psps%usepaw==1.and.Pawfgr%usefinegrid==1) then
686 2 : call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred)
687 : else
688 140015 : ph1df(:,:)=ph1d(:,:)
689 : end if
690 :
691 : ! Initialize qp_ebands using KS bands
692 : ! In case of SCGW, update qp_ebands using the QPS file.
693 148 : call ks_ebands%copy(qp_ebands)
694 148 : call timab(319,2,tsec) ! screening(1)
695 :
696 : !============================
697 : !==== Self-consistent GW ====
698 : !============================
699 148 : if (Ep%gwcalctyp>=10) then
700 49 : call timab(304,1,tsec) ! KS => QP; [wfrg]
701 :
702 : ! Initialize with KS eigenvalues and eigenfunctions.
703 294 : ABI_MALLOC(m_ks_to_qp,(Wfd%mband,Wfd%mband,Wfd%nkibz,Wfd%nsppol))
704 39764 : m_ks_to_qp = czero
705 98 : do spin=1,Wfd%nsppol
706 348 : do ik_ibz=1,Wfd%nkibz
707 3167 : do band=1,Wfd%nband(ik_ibz,spin)
708 24738 : m_ks_to_qp(band,band,:,:) = cone
709 : end do
710 : end do
711 : end do
712 :
713 : ! Read unitary transformation and QP energies.
714 : ! TODO switch on the renormalization of n in screening, QPS should report bdgw
715 196 : ABI_MALLOC(rhor_p,(nfftf,Dtset%nspden))
716 147 : ABI_MALLOC(prev_Pawrhoij,(Cryst%natom*Psps%usepaw))
717 :
718 : call rdqps(qp_ebands,Dtfil%fnameabi_qps,Dtset%usepaw,Dtset%nspden,1,nscf,&
719 49 : nfftf,ngfftf,Cryst%ucvol,Cryst,Pawtab,MPI_enreg_seq,nbsc,m_ks_to_qp,rhor_p,prev_Pawrhoij)
720 :
721 49 : ABI_FREE(rhor_p)
722 49 : ABI_FREE(prev_Pawrhoij)
723 :
724 : ! FIXME this is to preserve the old implementation for the head and the wings in ccchi0q0
725 : ! But has to be rationalized
726 49 : if (dtset%use_oldchi == 1) then
727 2981 : ks_ebands%eig = qp_ebands%eig
728 : end if
729 :
730 : ! Calculate new occ. factors and fermi level.
731 49 : call qp_ebands%update_occ(Dtset%spinmagntarget)
732 49 : qp_vbik(:,:) = qp_ebands%get_valence_idx()
733 :
734 : ! === Update only the wfg treated with GW ===
735 : ! For PAW update and re-symmetrize cprj in the full BZ, TODO add rotation in spinor space
736 49 : if (nscf/=0) call wfd%rotate(Cryst,m_ks_to_qp)
737 :
738 49 : ABI_FREE(m_ks_to_qp)
739 49 : call timab(304,2,tsec)
740 : end if ! gwcalctyp>=10
741 :
742 148 : call timab(305,1,tsec) ! screening(densit
743 : !
744 : !=== In case update the eigenvalues ===
745 : !* Either use a scissor operator or an external GW file.
746 148 : gw_fname = "__in.gw__"
747 148 : update_energies = file_exists(gw_fname)
748 :
749 148 : if (ABS(Ep%mbpt_sciss)>tol6) then
750 : write(msg,'(5a,f7.3,a)')&
751 1 : ' screening : performing a first self-consistency',ch10,&
752 1 : ' update of the energies in W by a scissor operator',ch10,&
753 2 : ' applying a scissor operator of [eV] : ',Ep%mbpt_sciss*Ha_eV,ch10
754 1 : call wrtout(units, msg)
755 1 : call qp_ebands%apply_scissors(Ep%mbpt_sciss)
756 147 : else if (update_energies) then
757 : write(msg,'(4a)')&
758 0 : ' screening : performing a first self-consistency',ch10,&
759 0 : ' update of the energies in W by a previous GW calculation via GW file: ',TRIM(gw_fname)
760 0 : call wrtout(units, msg)
761 0 : ABI_MALLOC(igwene,(qp_ebands%mband, qp_ebands%nkpt, qp_ebands%nsppol))
762 0 : call rdgw(qp_ebands,gw_fname,igwene,extrapolate=.FALSE.)
763 : !call rdgw(qp_ebands,gw_fname,igwene,extrapolate=.TRUE.)
764 0 : ABI_FREE(igwene)
765 0 : call qp_ebands%update_occ(Dtset%spinmagntarget)
766 : end if
767 :
768 : !========================
769 : !=== COMPUTE DENSITY ====
770 : !========================
771 : !* Evaluate PW part (complete charge in case of NC pseudos)
772 : !TODO this part has to be rewritten. If I decrease the tol on the occupations
773 : !I have to code some MPI stuff also if use_tr==.TRUE.
774 :
775 : !Get electronic temperature from dtset
776 148 : el_temp = merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
777 :
778 592 : ABI_MALLOC(rhor,(nfftf,Dtset%nspden))
779 592 : ABI_MALLOC(taur,(nfftf,Dtset%nspden*Dtset%usekden))
780 :
781 148 : call wfd%mkrho(cryst, psps, qp_ebands, ngfftf, nfftf, rhor)
782 148 : if (Dtset%usekden==1) call wfd%mkrho(cryst, psps, qp_ebands, ngfftf, nfftf, taur, optcalc=1)
783 :
784 148 : call timab(305,2,tsec) ! screening(densit
785 :
786 148 : nhatgrdim = 0
787 148 : if (Dtset%usepaw==1) then ! Additional computation for PAW.
788 5 : call timab(320,1,tsec) ! screening(paw
789 :
790 : ! Add the compensation charge to the PW density.
791 5 : nhatgrdim=0; if (Dtset%xclevel==2) nhatgrdim=usexcnhat*Dtset%pawnhatxc
792 5 : cplex=1; ider=2*nhatgrdim; izero=0
793 5 : if (nhatgrdim>0) then
794 5 : ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,3))
795 : else
796 12 : ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,0))
797 : end if
798 : call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,Cryst%gprimd,&
799 : Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,&
800 5 : Pawfgrtab,nhatgr,nhat,Pawrhoij,Pawrhoij,Pawtab,k0,Cryst%rprimd,Cryst%ucvol,dtset%usewvl,Cryst%xred)
801 :
802 : ! === Evaluate onsite energies, potentials, densities ===
803 : ! * Initialize variables/arrays related to the PAW spheres.
804 : ! * Initialize also lmselect (index of non-zero LM-moments of densities).
805 5 : cplex=1
806 46 : ABI_MALLOC(Paw_ij,(Cryst%natom))
807 5 : call paw_ij_nullify(Paw_ij)
808 : call paw_ij_init(Paw_ij,cplex,Dtset%nspinor,Wfd%nsppol,&
809 : Wfd%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,&
810 5 : has_dij=1,has_dijhartree=1,has_exexch_pot=1,has_pawu_occ=1)
811 :
812 5 : nkxc1=0
813 46 : ABI_MALLOC(Paw_an,(Cryst%natom))
814 5 : call paw_an_nullify(Paw_an)
815 : call paw_an_init(Paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,&
816 5 : cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=0)
817 :
818 5 : nzlmopt=-1; option=0; compch_sph=greatest_real
819 : call pawdenpot(compch_sph,el_temp,Cryst%gprimd,ipert,Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,&
820 : & Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,Paw_an,Paw_an,KS_energies%paw,Paw_ij,Pawang,&
821 : & Dtset%pawprtvol,Pawrad,Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,&
822 5 : & Dtset%xclevel,Dtset%xc_denpos,Dtset%xc_taupos,Cryst%xred,Cryst%ucvol,Psps%znuclpsp,Dtset%spinaxis)
823 5 : call timab(320,2,tsec) ! screening(paw
824 : else
825 143 : ABI_MALLOC(Paw_ij,(0))
826 143 : ABI_MALLOC(Paw_an,(0))
827 : end if ! usepaw
828 :
829 148 : call timab(321,1,tsec) ! screening(2)
830 :
831 : !JB : Should be remove : cf. l 839
832 148 : if (.not.allocated(nhatgr)) then
833 429 : ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,0))
834 : end if
835 :
836 : call test_charge(nfftf,ks_ebands%nelect,Dtset%nspden,rhor,ucvol,&
837 148 : Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,omegaplasma)
838 :
839 : !For PAW, add the compensation charge the FFT mesh, then get rho(G).
840 132625 : if (Dtset%usepaw==1) rhor(:,:)=rhor(:,:)+nhat(:,:)
841 :
842 148 : call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,rhor,ucvol=ucvol)
843 148 : if(Dtset%usekden==1)then
844 0 : call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,taur,ucvol=ucvol,optrhor=1)
845 : end if
846 :
847 148 : if (dtset%gwgamma>0 .or. dtset%gwgamma==-11) then
848 8 : ABI_MALLOC(rhor_kernel,(nfftf,Dtset%nspden))
849 : end if
850 :
851 444 : ABI_MALLOC(rhog,(2,nfftf))
852 148 : call fourdp(1,rhog,rhor(:,1),-1,MPI_enreg_seq,nfftf,1,ngfftf,tim_fourdp4)
853 :
854 : !The following steps have been gathered in the setvtr routine:
855 : !- get Ewald energy and Ewald forces
856 : !- compute local ionic pseudopotential vpsp
857 : !- eventually compute 3D core electron density xccc3d
858 : !- eventually compute vxc and vhartr
859 : !- set up ks_vtrial
860 : !**************************************************************
861 : !**** NOTE THAT Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION ****
862 : !**************************************************************
863 :
864 148 : ngrvdw=0
865 148 : ABI_MALLOC(grvdw,(3,ngrvdw))
866 444 : ABI_MALLOC(grchempottn,(3,Cryst%natom))
867 296 : ABI_MALLOC(grewtn,(3,Cryst%natom))
868 148 : nkxc=0
869 148 : if (Dtset%nspden==1) nkxc=2
870 148 : if (Dtset%nspden>=2) nkxc=3 ! check GGA and spinor that is messy !!!
871 : ! If MGGA, fxc and kxc are not available and we dont need them for the screening part (for now ...)
872 148 : if (Dtset%ixc<0 .and. libxc_functionals_ismgga()) nkxc=0
873 148 : if (nkxc/=0) then
874 592 : ABI_MALLOC(kxc,(nfftf,nkxc))
875 : end if
876 :
877 148 : n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf
878 444 : ABI_MALLOC(xccc3d,(n3xccc))
879 444 : ABI_MALLOC(ks_vhartr,(nfftf))
880 592 : ABI_MALLOC(ks_vtrial,(nfftf,Dtset%nspden))
881 296 : ABI_MALLOC(vpsp,(nfftf))
882 444 : ABI_MALLOC(ks_vxc,(nfftf,Dtset%nspden))
883 :
884 148 : optene=4; moved_atm_inside=0; moved_rhor=0; initialized=1; istep=1
885 : call setvtr(Cryst%atindx1,Dtset,KS_energies,Cryst%gmet,Cryst%gprimd,grchempottn, &
886 : grewtn,grvdw,gsqcutf_eff,istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq, &
887 : Cryst%nattyp,nfftf,ngfftf,ngrvdw,nhat,nhatgr,nhatgrdim,nkxc,Cryst%ntypat,&
888 : Psps%n1xccc,n3xccc,optene,Pawang,Pawrad,Pawrhoij,Pawtab,ph1df,Psps,rhog,rhor, &
889 : Cryst%rmet,Cryst%rprimd,strsxc,Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc, &
890 148 : vxcavg,wvl_dummy,xccc3d,Cryst%xred,taur=taur)
891 :
892 148 : if (nkxc/=0) then
893 148 : ABI_FREE(kxc)
894 : end if
895 148 : ABI_FREE(grchempottn)
896 148 : ABI_FREE(grewtn)
897 148 : ABI_FREE(grvdw)
898 148 : ABI_FREE(xccc3d)
899 :
900 : !============================
901 : !==== Compute KS PAW Dij ====
902 : !============================
903 148 : if (Dtset%usepaw==1) then
904 5 : call timab(561,1,tsec)
905 :
906 : ! Calculate unsymmetrized Dij.
907 5 : cplex=1; ipert=0; idir=0
908 : call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert,&
909 : & Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),&
910 : & Dtset%nspden,Cryst%ntypat,Paw_an,Paw_ij,Pawang,Pawfgrtab,Dtset%pawprtvol,&
911 : & Pawrad,Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,k0,Dtset%spnorbscl,&
912 : & Cryst%ucvol,dtset%cellcharge(1),ks_vtrial,ks_vxc,Cryst%xred,Dtset%znucl,&
913 5 : & nucdipmom=Dtset%nucdipmom,spinaxis=Dtset%spinaxis)
914 :
915 : ! Symmetrize KS Dij
916 : #if 0
917 : call symdij(Cryst%gprimd,Cryst%indsym,ipert,Cryst%natom,Cryst%natom,&
918 : & Cryst%nsym,Cryst%ntypat,0,Paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,&
919 : & Cryst%symrec)
920 : #else
921 : call symdij_all(Cryst%gprimd,Cryst%indsym,ipert,Cryst%natom,Cryst%natom,&
922 : & Cryst%nsym,Cryst%ntypat,Paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,&
923 5 : & Cryst%symrec)
924 : #endif
925 : ! Output of the pseudopotential strengths Dij and the augmentation occupancies Rhoij.
926 5 : call pawprt(Dtset,Cryst%natom,Paw_ij,Pawrhoij,Pawtab)
927 5 : call timab(561,2,tsec)
928 : end if
929 :
930 : ! Calculate frequency mesh.
931 : ! First omega is always zero without broadening.
932 : ! FIXME what about metals? I think we should add eta,
933 : ! this means we need to know if the system is metallic, for example using occopt
934 : ! MS Modified to account for non-zero starting frequency (19-11-2010)
935 : ! MS Modified for tangent grid (07-01-2011)
936 444 : ABI_MALLOC(Ep%omega, (Ep%nomega))
937 148 : Ep%omega(1) = CMPLX(Ep%omegaermin, zero, kind=dp)
938 :
939 148 : ep%iw_mesh_type = "None"; ep%rw_mesh_type = "None"; ep%cw_mesh_type = "None"
940 1399 : ABI_CALLOC(Ep%omega_wgs, (Ep%nomega))
941 :
942 : ! Deal with real frequencies first.
943 148 : if (Ep%nomegaer > 1) then
944 :
945 26 : if (Dtset%gw_frqre_tangrid == 0 .and. Dtset%gw_frqre_inzgrid == 0) then
946 24 : ep%rw_mesh_type = "linear"
947 24 : domegareal = (Ep%omegaermax -Ep%omegaermin) / (Ep%nomegaer -1)
948 290 : do iomega=2,Ep%nomegaer
949 290 : Ep%omega(iomega)=CMPLX(Ep%omegaermin+(iomega-1)*domegareal,zero,kind=dp)
950 : end do
951 :
952 2 : else if (Dtset%gw_frqre_tangrid == 1.and. Dtset%gw_frqre_inzgrid == 0) then
953 : ! We have tangent transformed grid
954 1 : ep%rw_mesh_type = "tangent_transform"
955 1 : ABI_WARNING('EXPERIMENTAL - Using tangent transform grid for contour deformation.')
956 1 : Ep%omegaermax = Dtset%cd_max_freq
957 1 : Ep%omegaermin = zero
958 1 : ifirst=1; ilast=Ep%nomegaer
959 1 : if (Dtset%cd_subset_freq(1)/=0) then ! Only a subset of frequencies is being calculated
960 0 : ifirst=Dtset%cd_subset_freq(1); ilast=Dtset%cd_subset_freq(2)
961 : end if
962 1 : factor = Dtset%cd_halfway_freq/TAN(pi*quarter)
963 : ! Important: here nfreqre is used because the step is set by the original grid
964 1 : domegareal=(ATAN(Ep%omegaermax/factor)*two*piinv)/(Dtset%nfreqre-1) ! Stepsize in transformed variable
965 6 : do iomega=1,Ep%nomegaer
966 6 : Ep%omega(iomega)=CMPLX(factor*TAN((iomega+ifirst-2)*domegareal*pi*half),zero,kind=dp)
967 : end do
968 1 : Ep%omegaermin = REAL(Ep%omega(1))
969 1 : Ep%omegaermax = REAL(Ep%omega(Ep%nomegaer))
970 :
971 1 : else if (Dtset%gw_frqre_tangrid==0.and.Dtset%gw_frqre_inzgrid==1) then
972 1 : ep%rw_mesh_type = "equidistant_z"
973 1 : e0=Dtset%ppmfrq; if (e0<0.1d-4) e0=omegaplasma
974 1 : domegareal=one/(Ep%nomegaer)
975 9 : do iomega=1,Ep%nomegaer
976 8 : factor = (iomega-1)*domegareal
977 9 : Ep%omega(iomega)=CMPLX(e0*factor/(one-factor),zero,kind=dp)
978 : end do
979 1 : Ep%omegaermin = REAL(Ep%omega(1))
980 1 : Ep%omegaermax = REAL(Ep%omega(Ep%nomegaer))
981 :
982 : else
983 0 : ABI_ERROR('Error in specification of real frequency grid')
984 : end if
985 : end if
986 :
987 148 : if (Ep%plasmon_pole_model .and. Ep%nomega == 2) then
988 104 : e0= Dtset%ppmfrq; if (e0 < 0.1d-4) e0 = omegaplasma
989 104 : Ep%omega(2)=CMPLX(zero,e0, kind=dp)
990 : end if
991 :
992 148 : if (Ep%analytic_continuation) then
993 : ! Negative value of dtset%nfreqim activates the minimax mesh
994 11 : ep%iw_mesh_type = "gauss_legendre"
995 11 : if (dtset%nfreqim < 0) ep%iw_mesh_type = "minimax"
996 :
997 11 : select case (ep%iw_mesh_type)
998 : case ("gauss_legendre")
999 : ! Use Gauss-Legendre quadrature method.
1000 : ! Replace $ \int_0^\infty dx f(x) $ with $ \int_0^1 dz f(1/z - 1)/z^2 $.
1001 : ! Note that the grid is not log as required by CD thus we cannot use the same SCR file.
1002 :
1003 33 : ABI_MALLOC(z, (Ep%nomegaei))
1004 22 : ABI_MALLOC(zw, (Ep%nomegaei))
1005 11 : call coeffs_gausslegint(zero, one, z, zw, Ep%nomegaei)
1006 315 : do iomega=1,Ep%nomegaei
1007 304 : Ep%omega(Ep%nomegaer + iomega) = CMPLX(zero, one/z(iomega) - one, kind=dp)
1008 315 : Ep%omega_wgs(Ep%nomegaer + iomega) = zw(iomega)
1009 : end do
1010 11 : ABI_FREE(z)
1011 11 : ABI_FREE(zw)
1012 :
1013 : case ("minimax")
1014 0 : call wrtout(std_out, sjoin("Using minimax mesh with ntau:", itoa(dtset%nfreqim)))
1015 0 : gaps = ks_ebands%get_gaps(gap_err)
1016 0 : ABI_CHECK(gap_err == 0, "gap_err")
1017 : ! ================================
1018 : ! Setup tau/omega mesh and weights
1019 : ! ================================
1020 : ! Compute min/max transition energy taking into account nsppol if any.
1021 0 : te_min = minval(gaps%cb_min - gaps%vb_max)
1022 0 : te_max = maxval(ks_ebands%eig(mband,:,:) - ks_ebands%eig(1,:,:))
1023 0 : if (te_min <= tol6) then
1024 0 : te_min = tol6
1025 0 : ABI_ERROR("System is metallic or with a very small fundamental gap!")
1026 : end if
1027 :
1028 : call gx_minimax_grid(abs(dtset%nfreqim), te_min, te_max, & ! in
1029 : tau_mesh, tau_wgs, & ! all these args are out and allocated by the routine.
1030 : iw_mesh, iw_wgs, &
1031 : t2w_cos_wgs, w2t_cos_wgs, t2w_sin_wgs, &
1032 0 : ft_max_error, cosft_duality_error, ierr)
1033 0 : ABI_CHECK(ierr == 0, "Error in gx_minimax_grid")
1034 :
1035 0 : do iomega=1,Ep%nomegaei
1036 0 : Ep%omega(Ep%nomegaer + iomega) = CMPLX(zero, iw_mesh(iomega), kind=dp)
1037 0 : Ep%omega_wgs(Ep%nomegaer + iomega) = iw_wgs(iomega)
1038 : end do
1039 :
1040 0 : ABI_FREE_NOCOUNT(tau_mesh)
1041 0 : ABI_FREE_NOCOUNT(tau_wgs)
1042 0 : ABI_FREE_NOCOUNT(iw_mesh)
1043 0 : ABI_FREE_NOCOUNT(iw_wgs)
1044 0 : ABI_FREE_NOCOUNT(t2w_cos_wgs)
1045 0 : ABI_FREE_NOCOUNT(w2t_cos_wgs)
1046 0 : ABI_FREE_NOCOUNT(t2w_sin_wgs)
1047 0 : call gaps%free()
1048 :
1049 : case default
1050 11 : ABI_ERROR(sjoin("Invalid iw_mesh_type:", ep%iw_mesh_type))
1051 : end select
1052 :
1053 137 : else if (Ep%contour_deformation .and. Dtset%cd_customnimfrqs /= 0) then
1054 1 : ep%iw_mesh_type = "custom"
1055 1 : Ep%omega(Ep%nomegaer+1)=CMPLX(zero,Dtset%cd_imfrqs(1))
1056 4 : do iomega=2,Ep%nomegaei
1057 3 : if (Dtset%cd_imfrqs(iomega) <= Dtset%cd_imfrqs(iomega-1)) then
1058 0 : ABI_ERROR(' Specified imaginary frequencies need to be strictly increasing!')
1059 : end if
1060 4 : Ep%omega(Ep%nomegaer+iomega) = CMPLX(zero, Dtset%cd_imfrqs(iomega))
1061 : end do
1062 :
1063 136 : else if (Ep%contour_deformation .and. Dtset%gw_frqim_inzgrid /= 0) then
1064 1 : ep%iw_mesh_type = "equidistant_z"
1065 1 : e0 = Dtset%ppmfrq; if (e0 < 0.1d-4) e0 = omegaplasma
1066 1 : domegareal=one/(Ep%nomegaei+1)
1067 9 : do iomega=1,Ep%nomegaei
1068 8 : factor = iomega*domegareal
1069 9 : Ep%omega(Ep%nomegaer+iomega)=CMPLX(zero,e0*factor/(one-factor),kind=dp)
1070 : end do
1071 :
1072 135 : else if (Ep%contour_deformation.and. Ep%nomegaei /= 0) then
1073 : ! This is the default for CD. see calc_sigc_cd for the integration routine.
1074 21 : ep%iw_mesh_type = "logarithmic"
1075 21 : e0=Dtset%ppmfrq; if (e0<0.1d-4) e0=omegaplasma
1076 115 : do iomega=1,Ep%nomegaei
1077 : Ep%omega(Ep%nomegaer+iomega)=CMPLX(zero,e0/(Dtset%freqim_alpha-two)&
1078 115 : * (EXP(two/(Ep%nomegaei+1)*LOG(Dtset%freqim_alpha-one)*iomega)-one),kind=dp)
1079 : end do
1080 : end if
1081 :
1082 148 : if (Dtset%cd_full_grid/=0) then
1083 : ! Full grid will be calculated. Grid values are added after the last imaginary freq.
1084 5 : do ios=1,Ep%nomegaei
1085 21 : do iomega=2,Ep%nomegaer
1086 : Ep%omega(Ep%nomegaer+Ep%nomegaei+(ios-1)*(Ep%nomegaer-1)+(iomega-1)) = &
1087 20 : CMPLX(REAL(Ep%omega(iomega)),AIMAG(Ep%omega(Ep%nomegaer+ios)))
1088 : end do
1089 : end do
1090 : end if
1091 :
1092 : ! Here we change nomega and select a slice of frequencies.
1093 436 : if (any(dtset%scr_wrange /= 0)) then
1094 : ! Consistency check.
1095 4 : call wrtout(units, sjoin(" Selecting frequency range: ", trim(ltoa(dtset%scr_wrange))))
1096 4 : ifirst = dtset%scr_wrange(1); ilast = dtset%scr_wrange(2); ii = ilast - ifirst + 1
1097 4 : ABI_CHECK_IRANGE(ifirst, 1, ep%nomega, "Invalid scr_wrange(1)")
1098 4 : ABI_CHECK_IRANGE(ilast, 1, ep%nomega, "Invalid scr_wrange(2)")
1099 4 : ABI_CHECK_IGEQ(ilast, ifirst, "scr_wrange(2) should be >= scr_wrange(1)")
1100 :
1101 : ! New number of freqs.
1102 4 : ep%nomega = ii
1103 :
1104 : ! Trasfer freqs.
1105 12 : ABI_MALLOC(tmp_omega, (ii))
1106 45 : tmp_omega(:) = ep%omega(ifirst:ilast)
1107 4 : ABI_FREE(ep%omega)
1108 4 : ABI_MOVE_ALLOC(tmp_omega, ep%omega)
1109 :
1110 : ! Trasfer weights.
1111 12 : ABI_MALLOC(tmp_omega_wgs, (ii))
1112 45 : tmp_omega_wgs(:) = ep%omega_wgs(ifirst:ilast)
1113 4 : ABI_FREE(ep%omega_wgs)
1114 4 : ABI_MOVE_ALLOC(tmp_omega_wgs, ep%omega_wgs)
1115 : end if
1116 :
1117 : ! Report frequency mesh for chi0.
1118 148 : write(msg,'(2a)')ch10,' calculating chi0 at frequencies [eV] :'
1119 148 : call wrtout(units, msg)
1120 1062 : do iomega=1,Ep%nomega
1121 914 : write(msg,'(i3,2es16.6)')iomega,Ep%omega(iomega)*Ha_eV
1122 1062 : call wrtout(units, msg)
1123 : end do
1124 :
1125 : ! Allocate chi0, wings and array for chi0_sumrule check.
1126 444 : ABI_MALLOC(chi0_sumrule, (Ep%npwe))
1127 :
1128 148 : write(msg,'(a,f12.1,a)')' Memory required for chi0 matrix= ',two*gwp*Ep%npwe**2*Ep%nI*Ep%nJ*Ep%nomega*b2Mb," [Mb]."
1129 148 : call wrtout(std_out, msg)
1130 740 : ABI_MALLOC_OR_DIE(chi0, (Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega), ierr)
1131 : !
1132 : !============================== END OF THE INITIALIZATION PART ===========================
1133 : !
1134 : !======================================================================
1135 : !==== Loop over q-points. Calculate \epsilon^{-1} and save on disc ====
1136 : !======================================================================
1137 148 : call pstat_proc%print(_PSTAT_ARGS_)
1138 148 : call timab(321,2,tsec) ! screening(2)
1139 :
1140 148 : iqcalc = 0
1141 148 : if(Dtset%plowan_compute >= 10) then
1142 : call init_plowannier(Dtset%plowan_bandf,Dtset%plowan_bandi,Dtset%plowan_compute,&
1143 : Dtset%plowan_iatom,Dtset%plowan_it,Dtset%plowan_lcalc,Dtset%plowan_natom,&
1144 : Dtset%plowan_nbl,Dtset%plowan_nt,Dtset%plowan_projcalc,Dtset%acell_orig,&
1145 : Dtset%kptns,sum(Dtset%plowan_nbl),Dtset%nimage,Dtset%nkpt,Dtset%nspinor,&
1146 0 : Dtset%nsppol,Dtset%wtk,Dtset%dmft_t2g,wanibz_in)
1147 0 : call get_plowannier(wanibz_in,wanibz,Dtset)
1148 0 : call fullbz_plowannier(Dtset,Kmesh,Cryst,Pawang,wanibz,wanbz)
1149 : end if
1150 :
1151 1082 : do iqibz=1,Qmesh%nibz
1152 934 : call timab(306,1,tsec)
1153 934 : is_first_qcalc=(iqibz==1)
1154 :
1155 : ! Selective q-point calculation.
1156 934 : found=.FALSE.; label=iqibz
1157 934 : if (Ep%nqcalc /= Ep%nqibz) then
1158 44 : do ii=1,Ep%nqcalc
1159 128 : qtmp(:)=Qmesh%ibz(:,iqibz)-Ep%qcalc(:,ii)
1160 32 : found=(normv(qtmp,gmet,'G')<GW_TOLQ)
1161 44 : if (found) then
1162 6 : label=ii; EXIT !ii
1163 : end if
1164 : end do
1165 18 : if (.not.found) CYCLE !iqibz
1166 24 : qtmp(:)=Ep%qcalc(:,1)-Qmesh%ibz(:,iqibz)
1167 6 : is_first_qcalc=(normv(qtmp,gmet,'G')<GW_TOLQ)
1168 : end if
1169 922 : iqcalc = iqcalc + 1
1170 :
1171 922 : bar=REPEAT('-',80)
1172 922 : write(msg,'(4a,1x,a,i2,a,f9.6,2(",",f9.6),3a)')ch10,ch10,bar,ch10,&
1173 4610 : ' q-point number ',label,' q = (',(Qmesh%ibz(ii,iqibz),ii=1,3),') [r.l.u.]',ch10,bar
1174 922 : call wrtout(units, msg)
1175 922 : is_qeq0 = 0; if (normv(Qmesh%ibz(:,iqibz),gmet,'G')<GW_TOLQ0) is_qeq0=1
1176 :
1177 922 : call pstat_proc%print(_PSTAT_ARGS_)
1178 922 : call timab(306,2,tsec)
1179 :
1180 922 : if (is_qeq0 == 1) then
1181 : ! Special treatment of the long wavelength limit.
1182 146 : call timab(307,1,tsec)
1183 :
1184 438 : ABI_MALLOC(chi0_head, (3,3,Ep%nomega))
1185 730 : ABI_MALLOC(chi0_lwing, (Ep%npwe*Ep%nI, Ep%nomega,3))
1186 730 : ABI_MALLOC(chi0_uwing, (Ep%npwe*Ep%nJ, Ep%nomega,3))
1187 :
1188 : call cchi0q0(use_tr,Dtset,Cryst,Ep,Psps,Kmesh,qp_ebands,ks_ebands,Gsph_epsG0,&
1189 : Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,Pawfgrtab,Paw_onsite,ktabr,ktabrf,nbvw,ngfft_gw,nfftgw,&
1190 146 : ngfftf,nfftf_tot,chi0,chi0_head,chi0_lwing,chi0_uwing,Ltg_q(iqibz),chi0_sumrule,Wfd,Wfdf,wanbz) !,mpi_enreg_seq
1191 :
1192 : !chihw = chi_new(ep%npwe, ep%nomega)
1193 : !chihw%head = chi0_head
1194 : !chihw%lwing = chi0_lwing
1195 : !chihw%uwing = chi0_uwing
1196 :
1197 : ! Add the intraband term if required and metallic occupation scheme is used.
1198 146 : add_chi0_intraband=.FALSE. !add_chi0_intraband=.TRUE.
1199 146 : if (add_chi0_intraband .and. qp_ebands%has_metal_scheme()) then
1200 :
1201 : ABI_MALLOC_OR_DIE(chi0intra,(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega), ierr)
1202 :
1203 : ABI_MALLOC(chi0intra_lwing,(Ep%npwe*Ep%nI,Ep%nomega,3))
1204 : ABI_MALLOC(chi0intra_uwing,(Ep%npwe*Ep%nJ,Ep%nomega,3))
1205 : ABI_MALLOC(chi0intra_head,(3,3,Ep%nomega))
1206 :
1207 : call chi0q0_intraband(Wfd,Cryst,Ep,Psps,qp_ebands,Gsph_epsG0,Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,use_tr,Dtset%usepawu,&
1208 : ngfft_gw,chi0intra,chi0intra_head,chi0intra_lwing,chi0intra_uwing)
1209 :
1210 : call wrtout(std_out,"Head of chi0 and chi0_intra")
1211 : do iomega=1,Ep%nomega
1212 : write(std_out,*)Ep%omega(iomega)*Ha_eV,REAL(chi0(1,1,iomega)),REAL(chi0intra(1,1,iomega))
1213 : write(std_out,*)Ep%omega(iomega)*Ha_eV,AIMAG(chi0(1,1,iomega)),AIMAG(chi0intra(1,1,iomega))
1214 : end do
1215 :
1216 : chi0 = chi0 + chi0intra
1217 : chi0_head = chi0_head + chi0intra_head
1218 : chi0_lwing = chi0_lwing + chi0intra_lwing
1219 : chi0_uwing = chi0_uwing + chi0intra_uwing
1220 :
1221 : ABI_FREE(chi0intra)
1222 : ABI_FREE(chi0intra_lwing)
1223 : ABI_FREE(chi0intra_uwing)
1224 : ABI_FREE(chi0intra_head)
1225 : end if
1226 :
1227 : if (.False.) then
1228 : lwl_fname = strcat(dtfil%filnam_ds(4), "_LWL")
1229 : call lwl_write(lwl_fname,cryst,vcp,ep%npwe,ep%nomega,gsph_epsg0%gvec,chi0,chi0_head,chi0_lwing,chi0_uwing,comm)
1230 : end if
1231 :
1232 146 : call timab(307,2,tsec)
1233 :
1234 : else
1235 : ! Calculate cchi0 for q/=0.
1236 776 : call timab(308,1,tsec)
1237 : call cchi0(use_tr,Dtset,Cryst,Qmesh%ibz(:,iqibz),Ep,Psps,Kmesh,qp_ebands,Gsph_epsG0,&
1238 : Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,nbvw,ngfft_gw,nfftgw,ngfftf,nfftf_tot,chi0,ktabr,ktabrf,&
1239 776 : Ltg_q(iqibz),chi0_sumrule,Wfd,Wfdf,wanbz)
1240 776 : call timab(308,2,tsec)
1241 : end if
1242 :
1243 : ! Print chi0(q,G,Gp,omega), then calculate epsilon and epsilon^-1 for this q-point.
1244 : ! Only master works but this part could be parallelized over frequencies.
1245 922 : call timab(309,1,tsec)
1246 :
1247 5084 : do iomega=1,MIN(Ep%nomega, NOMEGA_PRINTED)
1248 4162 : write(msg,'(1x,a,i4,a,2f9.4,a)')' chi0(G,G'') at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]'
1249 4162 : call wrtout(units, msg)
1250 4162 : write(msg,'(1x,a,i3,a,i4,a)')' chi0(q =',iqibz, ', omega =',iomega,', G,G'')'
1251 4162 : if (Ep%nqcalc /= Ep%nqibz) write(msg,'(a,i3,a,i4,a)')' chi0(q=',iqcalc,', omega=',iomega,', G,G'')'
1252 4162 : call wrtout(std_out, msg)
1253 : ! arr99 is needed to avoid the update of all the tests. Now chi0 is divided by ucvol inside (cchi0|cchi0q0).
1254 : ! TODO should be removed but GW tests have to be updated.
1255 4162 : ii = MIN(9, Ep%npwe)
1256 16648 : ABI_MALLOC(arr_99,(ii, ii))
1257 382904 : arr_99 = chi0(1:ii,1:ii,iomega) * ucvol
1258 8324 : call print_arr([ab_out], arr_99, max_r=2)
1259 8324 : call print_arr([std_out], arr_99)
1260 5084 : ABI_FREE(arr_99)
1261 : end do
1262 :
1263 922 : if (Ep%nomega > NOMEGA_PRINTED) then
1264 84 : write(msg,'(a,i3,a)')' No. of calculated frequencies > ',NOMEGA_PRINTED,', stop printing '
1265 84 : call wrtout(units, msg)
1266 : end if
1267 :
1268 : ! Write chi0 to _SUSC file
1269 : ! Master creates and write the header if this is the first q-point calculated.
1270 922 : if (Dtset%prtsuscep > 0 .and. my_rank == master) then
1271 39 : title(1)='CHI0 file: chi0'
1272 : title(2)=' '
1273 39 : if (is_qeq0 == 1) then
1274 6 : string='0'; if (Dtset%usepaw==0.and.Ep%inclvkb/=0) call int2char10(Ep%inclvkb,string)
1275 6 : title(1)=title(1)(1:21)//', calculated using inclvkb = '//string
1276 : end if
1277 :
1278 : ! Open file and write header for polarizability files.
1279 39 : if (is_first_qcalc) then
1280 8 : ikxc=0; test_type=0; tordering=1
1281 8 : hchi0 = hscr_new("polarizability",dtset,ep,hdr_local,ikxc,test_type,tordering,title,Ep%npwe,Gsph_epsG0%gvec)
1282 :
1283 8 : if (dtset%iomode == IO_MODE_ETSF) then
1284 8 : NCF_CHECK(nctk_open_create(unt_susc, nctk_ncify(dtfil%fnameabo_sus), xmpi_comm_self))
1285 8 : NCF_CHECK(cryst%ncwrite(unt_susc))
1286 8 : NCF_CHECK(qp_ebands%ncwrite(unt_susc))
1287 : else
1288 0 : unt_susc=Dtfil%unchi0
1289 0 : if (open_file(dtfil%fnameabo_sus,msg,unit=unt_susc,status='unknown',form='unformatted') /= 0) then
1290 0 : ABI_ERROR(msg)
1291 : end if
1292 : end if
1293 :
1294 8 : fform_chi0 = hchi0%fform
1295 8 : call hchi0%io(fform_chi0,2,unt_susc,xmpi_comm_self,0,Dtset%iomode)
1296 8 : call Hchi0%free()
1297 : end if
1298 :
1299 39 : call write_screening("polarizability", unt_susc, Dtset%iomode, Ep%npwe, Ep%nomega, iqcalc, chi0)
1300 :
1301 39 : if (dtset%iomode == IO_MODE_ETSF .and. is_qeq0 == 1 .and. Ep%nI == 1 .and. Ep%nJ == 1) then
1302 : ! Write head and wings to file. See cchi0 for the equations needed to build chi0(q) for q--> 0.
1303 6 : wing_shape = "two, three, number_of_coefficients_dielectric_function, number_of_frequencies_dielectric_function"
1304 : ncerr = nctk_def_arrays(unt_susc, [ &
1305 : nctkarr_t("sus_head", "dp", "two, three, three, number_of_frequencies_dielectric_function"), &
1306 : nctkarr_t("sus_upper_wing", "dp", wing_shape), &
1307 : nctkarr_t("sus_lower_wing", "dp", wing_shape) &
1308 24 : ], defmode=.True.)
1309 6 : NCF_CHECK(ncerr)
1310 :
1311 6 : NCF_CHECK(nctk_set_datamode(unt_susc))
1312 6 : NCF_CHECK(nf90_put_var(unt_susc, nctk_idname(unt_susc, "sus_head"), c2r(chi0_head)))
1313 :
1314 24 : ABI_MALLOC(rwork_wing, (2, 3, Ep%npwe * Ep%nI, Ep%nomega))
1315 51 : do iomega=1,Ep%nomega
1316 1498 : do ig=1, Ep%npwe * Ep%nI
1317 5788 : rwork_wing(1,:,ig,iomega) = real(chi0_lwing(ig,iomega,:))
1318 5833 : rwork_wing(2,:,ig,iomega) = aimag(chi0_lwing(ig,iomega,:))
1319 : end do
1320 : end do
1321 6 : NCF_CHECK(nf90_put_var(unt_susc, nctk_idname(unt_susc, "sus_lower_wing"), rwork_wing))
1322 :
1323 51 : do iomega=1,Ep%nomega
1324 1498 : do ig=1, Ep%npwe * Ep%nI
1325 5788 : rwork_wing(1,:,ig,iomega) = real(chi0_uwing(ig,iomega,:))
1326 5833 : rwork_wing(2,:,ig,iomega) = aimag(chi0_uwing(ig,iomega,:))
1327 : end do
1328 : end do
1329 6 : NCF_CHECK(nf90_put_var(unt_susc, nctk_idname(unt_susc, "sus_upper_wing"), rwork_wing))
1330 6 : ABI_FREE(rwork_wing)
1331 : end if
1332 :
1333 : end if ! is_first_qcalc
1334 :
1335 : ! Calculate the Galitskii-Migdal and RPA functionals for the correlation energy if the polarizability on a
1336 : ! Gauss-Legendre mesh along imaginary axis is available
1337 922 : if (Ep%analytic_continuation .and. Dtset%gwrpacorr>0 ) then
1338 18 : if (is_first_qcalc) then
1339 15 : ABI_MALLOC(ec_rpa,(Dtset%gwrpacorr))
1340 17 : ec_rpa(:)=zero
1341 5 : ec_gm=zero
1342 : end if
1343 18 : call calc_rpa_functional(Dtset%gwrpacorr,Dtset%gwgmcorr,label,iqibz,Ep,Vcp,Qmesh,Dtfil,gmet,chi0,comm,ec_rpa,ec_gm)
1344 18 : if (label==Ep%nqcalc) then
1345 5 : ABI_FREE(ec_rpa)
1346 : end if
1347 : end if
1348 :
1349 : ! ==========================================================
1350 : ! === Calculate RPA \tilde\epsilon^{-1} overwriting chi0 ===
1351 : ! ==========================================================
1352 922 : approx_type=0 ! RPA
1353 922 : option_test=0 ! TESTPARTICLE
1354 922 : dim_wing=0; if (is_qeq0==1) dim_wing=3
1355 :
1356 922 : if (dim_wing==0) then
1357 776 : dim_wing=1
1358 776 : if (.not.allocated(chi0_lwing)) then
1359 3104 : ABI_MALLOC(chi0_lwing,(Ep%npwe*Ep%nI,Ep%nomega,dim_wing))
1360 : end if
1361 776 : if (.not.allocated(chi0_uwing)) then
1362 3104 : ABI_MALLOC(chi0_uwing,(Ep%npwe*Ep%nJ,Ep%nomega,dim_wing))
1363 : end if
1364 776 : if (.not.allocated(chi0_head )) then
1365 2328 : ABI_MALLOC(chi0_head,(dim_wing,dim_wing,Ep%nomega))
1366 : end if
1367 776 : dim_wing=0
1368 : end if
1369 :
1370 : #if 0
1371 : ! Using the random q for the optical limit is one of the reasons
1372 : ! why sigma breaks the initial energy degeneracies.
1373 : chi0_lwing=czero
1374 : chi0_uwing=czero
1375 : chi0_head=czero
1376 : #endif
1377 :
1378 : ! Setup flags for the computation of em1
1379 : ! If the vertex is being included for the spectrum, calculate the kernel now and pass it on
1380 21184 : if (dtset%gwgamma>0) rhor_kernel = rhor
1381 :
1382 900 : select case (dtset%gwgamma)
1383 : case (0)
1384 900 : approx_type=0; option_test=0; dim_kxcg=0
1385 1800 : ABI_MALLOC(kxcg, (nfftf_tot, dim_kxcg))
1386 :
1387 : case (1, 2)
1388 : ! ALDA TDDFT kernel vertex
1389 6 : ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available")
1390 6 : ABI_WARNING('EXPERIMENTAL: Kernel is being added to screening, the SCR file will be non-standard!!')
1391 6 : ikxc=7; approx_type=1; dim_kxcg=1
1392 6 : if (Dtset%gwgamma==1) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma
1393 6 : if (Dtset%gwgamma==2) option_test=0 ! TESTPARTICLE, vertex in chi0 only
1394 18 : ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1395 : call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,rhor_kernel,&
1396 6 : Ep%npwe,dim_kxcg,kxcg,Gsph_epsG0%gvec,xmpi_comm_self)
1397 :
1398 : case (3, 4)
1399 : ! ADA non-local kernel vertex
1400 0 : ABI_CHECK(Wfd%usepaw==0,"ADA vertex + PAW not available")
1401 0 : ABI_CHECK(Wfd%nsppol==1,"ADA vertex for GWGamma not available yet for spin-polarised cases")
1402 0 : ABI_WARNING('EXPERIMENTAL: Kernel is being added to screening, the SCR file will be non-standard!!')
1403 0 : ikxc=7; approx_type=2
1404 0 : if (Dtset%gwgamma==3) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma
1405 0 : if (Dtset%gwgamma==4) option_test=0 ! TESTPARTICLE, vertex in chi0 only
1406 0 : ABI_MALLOC(fxc_ADA,(Ep%npwe,Ep%npwe,Ep%nqibz))
1407 : ! Use userrd to set kappa
1408 0 : if (Dtset%userrd==zero) Dtset%userrd = 2.1_dp
1409 : ! Set correct value of kappa (should be scaled with alpha*r_s where)
1410 : ! r_s is Wigner-Seitz radius and alpha=(4/(9*Pi))^(1/3)
1411 0 : rhoav = (omegaplasma*omegaplasma)/four_pi
1412 0 : r_s = (three/(four_pi*rhoav))**third
1413 0 : alpha = (four*ninth*piinv)**third
1414 0 : Dtset%userrd = Dtset%userrd*alpha*r_s
1415 :
1416 : call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf,Wfd%nspden,rhor_kernel,Ep%npwe,Ep%nqibz,Ep%qibz,&
1417 0 : fxc_ADA,Gsph_epsG0%gvec,xmpi_comm_self,kappa_init=Dtset%userrd)
1418 :
1419 0 : dim_kxcg=0
1420 0 : ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1421 :
1422 : case (-3, -4, -5, -6, -7, -8)
1423 : ! Bootstrap kernel and variants
1424 12 : ABI_CHECK(Dtset%usepaw==0,"GWGamma + PAW not available")
1425 12 : if (Dtset%gwgamma>-5) then
1426 8 : ABI_WARNING('EXPERIMENTAL: Bootstrap kernel is being added to screening')
1427 8 : approx_type=4
1428 4 : else if (Dtset%gwgamma>-7) then
1429 4 : ABI_WARNING('EXPERIMENTAL: Bootstrap kernel (head-only) is being added to screening')
1430 4 : approx_type=5
1431 : else
1432 0 : ABI_WARNING('EXPERIMENTAL: RPA Bootstrap kernel is being added to screening')
1433 0 : approx_type=6
1434 : end if
1435 12 : dim_kxcg=0
1436 12 : option_test=MOD(Dtset%gwgamma,2)
1437 : ! 1 -> TESTELECTRON, vertex in chi0 *and* sigma
1438 : ! 0 -> TESTPARTICLE, vertex in chi0 only
1439 24 : ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1440 :
1441 : case (-11)
1442 : ! LR+ALDA hybrid vertex kernel
1443 4 : ABI_CHECK(Dtset%usepaw==0,"GWGamma + PAW not available")
1444 4 : ikxc=7; dim_kxcg=1
1445 4 : ABI_WARNING('EXPERIMENTAL: LR+ALDA hybrid kernel is being added to screening')
1446 4 : approx_type=7
1447 4 : option_test=1 ! TESTELECTRON
1448 16 : ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1449 55308 : rhor_kernel = rhor
1450 : call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,rhor_kernel,&
1451 4 : Ep%npwe,dim_kxcg,kxcg,Gsph_epsG0%gvec,xmpi_comm_self)
1452 4 : rhoav = (omegaplasma*omegaplasma)/four_pi
1453 :
1454 : case default
1455 922 : ABI_ERROR(sjoin("Wrong gwgamma:", itoa(dtset%gwgamma)))
1456 : end select
1457 :
1458 922 : if (approx_type<2) then
1459 : ! ALDA
1460 : call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,&
1461 : approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,&
1462 906 : chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm, epsm1_win)
1463 :
1464 16 : else if (approx_type<3) then
1465 : ! ADA
1466 : call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,&
1467 : approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,&
1468 : chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm, epsm1_win, &
1469 0 : fxc_ADA=fxc_ADA(:,:,iqibz))
1470 :
1471 16 : else if (approx_type<7) then
1472 : ! Bootstrap
1473 : call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,&
1474 : approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,&
1475 12 : chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm, epsm1_win)
1476 :
1477 4 : else if (approx_type<8) then
1478 : ! LR + ALDA
1479 : call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,&
1480 : approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,&
1481 : chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm, epsm1_win, &
1482 4 : rhor=rhoav)
1483 :
1484 : else
1485 0 : ABI_ERROR(sjoin("Wrong approx_type:", itoa(approx_type)))
1486 : end if
1487 :
1488 922 : ABI_FREE(chi0_lwing)
1489 922 : ABI_FREE(chi0_uwing)
1490 922 : ABI_FREE(chi0_head)
1491 :
1492 922 : if (my_rank == master .and. is_qeq0==1) then
1493 118 : call spectra%repr(msg)
1494 118 : call wrtout(units, msg)
1495 118 : if (Ep%nomegaer > 2) then
1496 22 : call spectra%write(W_EELF ,Dtfil%fnameabo_eelf)
1497 22 : call spectra%write(W_EM_LF ,Dtfil%fnameabo_em1_lf)
1498 22 : call spectra%write(W_EM_NLF,Dtfil%fnameabo_em1_nlf)
1499 : end if
1500 : end if ! master and is_qeq0==1
1501 :
1502 : !if (is_qeq0==1) call chi_free(chihw)
1503 :
1504 922 : call spectra%free()
1505 922 : ABI_SFREE(kxcg)
1506 922 : ABI_SFREE(fxc_ADA)
1507 : !
1508 : ! Output the sum rule evaluation.
1509 : ! Vcp%vc_sqrt(:,iqibz) Contains vc^{1/2}(q,G), complex-valued due to a possible cutoff
1510 922 : epsm1 => chi0
1511 922 : call output_chi0sumrule((is_qeq0==1),iqibz,Ep%npwe,omegaplasma,chi0_sumrule,epsm1(:,:,1),Vcp%vc_sqrt(:,iqibz))
1512 :
1513 : ! If input variable npvel is larger than 0, trigger the Random Stopping Power calculation
1514 : ! Only the masternode is used
1515 922 : if (my_rank==master .and. Dtset%npvel>0) then
1516 12 : if (is_first_qcalc) then
1517 6 : ABI_MALLOC(rspower,(Dtset%npvel))
1518 82 : rspower(:)=zero
1519 : end if
1520 12 : call random_stopping_power(iqibz,Dtset%npvel,Dtset%pvelmax,Ep,Gsph_epsG0,Qmesh,Vcp,Cryst,Dtfil,epsm1,rspower)
1521 12 : if (label==Ep%nqcalc) then
1522 2 : ABI_FREE(rspower)
1523 : end if
1524 : end if
1525 :
1526 : ! Write heads and wings to main output file.
1527 922 : if (is_qeq0 == 1) then
1528 146 : write(msg,'(1x,2a)')' Heads and wings of the symmetrical epsilon^-1(G,G'') ',ch10
1529 146 : call wrtout(ab_out,msg)
1530 1056 : do iomega=1,Ep%nomega
1531 910 : write(msg,'(2x,a,i4,a,2f9.4,a)')' Upper and lower wings at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]'
1532 910 : call wrtout(ab_out, msg)
1533 1820 : call print_arr([ab_out], epsm1(1,:,iomega),max_r=9)
1534 1820 : call print_arr([ab_out], epsm1(:,1,iomega),max_r=9)
1535 1056 : call wrtout(ab_out, ch10)
1536 : end do
1537 : end if
1538 :
1539 922 : call pstat_proc%print(_PSTAT_ARGS_)
1540 922 : call timab(309,2,tsec)
1541 922 : call timab(310,1,tsec) ! wrscr
1542 :
1543 922 : if (my_rank==master) then
1544 : ! === Write the symmetrical epsilon^-1 on file ===
1545 718 : title(1)='SCR file: epsilon^-1'
1546 718 : if (is_qeq0==1) then
1547 118 : string='0'; if (Dtset%usepaw==0.and.Ep%inclvkb/=0) call int2char10(Ep%inclvkb,string)
1548 118 : title(1)=title(1)(1:21)//', calculated using inclvkb = '//string
1549 : end if
1550 718 : title(2)='TESTPARTICLE'
1551 718 : ctype='RPA'
1552 718 : title(2)(14:17)=ctype !this has to be modified
1553 :
1554 718 : if (is_first_qcalc) then
1555 : ! === Open file and write the header for the SCR file ===
1556 : ! * Here we write the RPA approximation for \tilde\epsilon^{-1}
1557 120 : ikxc=0; test_type=0; tordering=1
1558 : hem1 = hscr_new("inverse_dielectric_function",dtset,ep,hdr_local,ikxc,test_type,tordering,title,&
1559 120 : & Ep%npwe,Gsph_epsG0%gvec)
1560 120 : fform_em1 = hem1%fform
1561 120 : if (dtset%iomode == IO_MODE_ETSF) then
1562 120 : NCF_CHECK(nctk_open_create(unt_em1, nctk_ncify(dtfil%fnameabo_scr), xmpi_comm_self))
1563 120 : NCF_CHECK(cryst%ncwrite(unt_em1))
1564 120 : NCF_CHECK(qp_ebands%ncwrite(unt_em1))
1565 : else
1566 0 : unt_em1=Dtfil%unscr
1567 0 : if (open_file(dtfil%fnameabo_scr,msg,unit=unt_em1,status='unknown',form='unformatted') /= 0) then
1568 0 : ABI_ERROR(msg)
1569 : end if
1570 : end if
1571 120 : call hem1%io(fform_em1,2,unt_em1,xmpi_comm_self,0,Dtset%iomode)
1572 120 : call Hem1%free()
1573 : end if
1574 :
1575 718 : call write_screening("inverse_dielectric_function",unt_em1,Dtset%iomode,Ep%npwe,Ep%nomega,iqcalc,epsm1)
1576 : end if ! my_rank==master
1577 :
1578 1082 : call timab(310,2,tsec)
1579 : end do ! Loop over q-points
1580 :
1581 148 : if (Dtset%plowan_compute >= 10) call destroy_plowannier(wanbz)
1582 :
1583 : ! Close Files.
1584 148 : if (my_rank == master) then
1585 120 : if (dtset%iomode == IO_MODE_ETSF) then
1586 120 : NCF_CHECK(nf90_close(unt_em1))
1587 120 : if (dtset%prtsuscep > 0) then
1588 8 : NCF_CHECK(nf90_close(unt_susc))
1589 : end if
1590 : else
1591 0 : close(unt_em1)
1592 0 : if (dtset%prtsuscep > 0) close(unt_susc)
1593 : end if
1594 : end if
1595 : !
1596 : !=====================
1597 : !==== Free memory ====
1598 : !=====================
1599 148 : ABI_FREE(chi0_sumrule)
1600 148 : ABI_FREE(chi0)
1601 148 : ABI_SFREE(rhor_kernel)
1602 :
1603 148 : ABI_FREE(rhor)
1604 148 : ABI_FREE(rhog)
1605 148 : ABI_FREE(ks_vbik)
1606 148 : ABI_FREE(qp_vbik)
1607 148 : ABI_FREE(ktabr)
1608 148 : ABI_FREE(taur)
1609 148 : ABI_FREE(ks_vhartr)
1610 148 : ABI_FREE(ks_vtrial)
1611 148 : ABI_FREE(vpsp)
1612 148 : ABI_FREE(ks_vxc)
1613 148 : ABI_FREE(ph1d)
1614 148 : ABI_FREE(ph1df)
1615 148 : ABI_FREE(nhatgr)
1616 148 : ABI_FREE(nhat)
1617 148 : call pawfgr_destroy(Pawfgr)
1618 :
1619 148 : if (Dtset%usepaw==1) then ! Optional deallocation for PAW.
1620 5 : call pawrhoij_free(Pawrhoij)
1621 5 : call pawfgrtab_free(Pawfgrtab)
1622 5 : call paw_ij_free(Paw_ij)
1623 5 : call paw_an_free(Paw_an)
1624 5 : call pawpwff_free(Paw_pwff)
1625 5 : if (Dtset%pawcross==1) then
1626 0 : call paw_pwaves_lmn_free(Paw_onsite)
1627 0 : call wfdf%free()
1628 : end if
1629 : end if
1630 :
1631 179 : ABI_FREE(Pawfgrtab)
1632 157 : ABI_FREE(Paw_pwff)
1633 179 : ABI_FREE(Pawrhoij)
1634 179 : ABI_FREE(Paw_ij)
1635 179 : ABI_FREE(Paw_an)
1636 148 : ABI_FREE(ktabrf)
1637 499 : ABI_FREE(Paw_onsite)
1638 :
1639 148 : call wfd%free(); call Kmesh%free(); call Qmesh%free(); call cryst%free()
1640 148 : call Gsph_epsG0%free(); call Gsph_wfn%free(); call Vcp%free(); call Ep%free()
1641 148 : call Hdr_wfk%free(); call Hdr_local%free()
1642 148 : call ks_ebands%free(); call qp_ebands%free()
1643 148 : call destroy_mpi_enreg(MPI_enreg_seq)
1644 148 : call littlegroup_free(ltg_q)
1645 148 : ABI_FREE(Ltg_q)
1646 :
1647 148 : call timab(301,2,tsec)
1648 :
1649 888 : end subroutine screening
1650 : !!***
1651 :
1652 : !!****f* m_screening_driver/setup_screening
1653 : !! NAME
1654 : !! setup_screening
1655 : !!
1656 : !! FUNCTION
1657 : !! Initialize the Ep% data type containing the parameters used during the screening calculation.
1658 : !! as well as basic objects describing the BZ sampling .... TODO list to be completed
1659 : !!
1660 : !! INPUTS
1661 : !! wfk_fname=Name of the input WFK file.
1662 : !! acell(3)=length scales of primitive translations (Bohr).
1663 : !! rprim(3,3)=dimensionless real space primitive translations.
1664 : !! dtfil <type(datafiles_type)>=variables related to files
1665 : !!
1666 : !! OUTPUT
1667 : !! ngfft_gw(18)=Contain all needed information about the 3D FFT for the oscillator strengths.
1668 : !! See ~abinit/doc/variables/vargs.htm#ngfft
1669 : !! Ltg_q(:)<littlegroup_t>,=
1670 : !! Ep<em1params_t>=Parameters for the screening calculation.
1671 : !! Most part of it is Initialized and checked.
1672 : !! Hdr_wfk type(Hdr_type)=Header of the KSS file.
1673 : !! Cryst<crystal_t>=Definition of the unit cell and its symmetries.
1674 : !! Kmesh<kmesh_t>=Structure defining the k-point sampling (wavefunctions).
1675 : !! Qmesh<kmesh_t>=Structure defining the q-point sampling (screening)
1676 : !! Gsph_wfn<gsphere_t>=Structure defining the G-sphere for the wavefunctions (not k-dependent).
1677 : !! Gsph_epsG0<gsphere_t>=The G-sphere for the screening, enlarged to take into account for umklapps.
1678 : !! Psps <Pseudopotential_type)>=Info on pseudopotential, only for consistency check of the KSS file
1679 : !! Vcp <type vcoul_t> datatype gathering information on the coulombian cutoff technique
1680 : !! comm=MPI communicator.
1681 : !!
1682 : !! SIDE EFFECTS
1683 : !! Dtset<Dataset_type>=All input variables for this dataset.
1684 : !! %ecutwfn, %npwwfn,
1685 : !! %ecuteps, %npweps
1686 : !! might be redefined in setshells in order to close the shell.
1687 : !!
1688 : !! SOURCE
1689 :
1690 14948 : subroutine setup_screening(codvsn,acell,rprim,wfk_fname,Dtset,Psps,Pawtab,&
1691 : ngfft_gw,Hdr_wfk,Hdr_out,Cryst,Kmesh,Qmesh,ks_ebands,Ltg_q,Gsph_epsG0,Gsph_wfn,Vcp,Ep,comm)
1692 :
1693 : !Arguments ------------------------------------
1694 : !scalars
1695 : integer,intent(in) :: comm
1696 : character(len=8),intent(in) :: codvsn
1697 : character(len=fnlen),intent(in) :: wfk_fname
1698 : type(Dataset_type),intent(inout) :: Dtset !INOUT is due to setshells
1699 : type(Pseudopotential_type),intent(in) :: Psps
1700 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw)
1701 : type(em1params_t),intent(out) :: Ep
1702 : type(Hdr_type),intent(out) :: Hdr_wfk,Hdr_out
1703 : type(ebands_t),intent(out) :: ks_ebands
1704 : type(kmesh_t),intent(out) :: Kmesh,Qmesh
1705 : type(crystal_t),intent(out) :: Cryst
1706 : type(gsphere_t),intent(out) :: Gsph_epsG0,Gsph_wfn
1707 : type(vcoul_t),intent(out) :: Vcp
1708 : !arrays
1709 : integer,intent(out) :: ngfft_gw(18)
1710 : real(dp),intent(in) :: acell(3),rprim(3,3)
1711 : type(littlegroup_t),pointer :: Ltg_q(:)
1712 :
1713 : !Local variables-------------------------------
1714 : !scalars
1715 : integer,parameter :: NOMEGAGAUSS=30,NOMEGAREAL=201,pertcase0=0,master=0
1716 : integer :: bantot,ib,ibtot,ikibz,iq,iqp,isppol,ig,ng,ierr
1717 : integer :: jj,mod10,mband,ng_kss,iqbz,isym,iq_ibz,itim
1718 : integer :: timrev,use_umklp !,ncerr
1719 : integer :: npwepG0,nshepspG0,method,enforce_sym,nfftgw_tot !,spin,band,ik_ibz,
1720 : integer :: istart,iend,test_npwkss,my_rank,nprocs !ii
1721 : real(dp),parameter :: OMEGAERMAX=100.0/Ha_eV
1722 : real(dp) :: ecutepspG0,ucvol,domegareal
1723 : logical :: remove_inv,ltest,found,is_static,has_q0
1724 : character(len=500) :: msg
1725 148 : type(wvl_internal_type) :: wvl
1726 : !arrays
1727 : integer :: ng0sh_opt(3), units(2)
1728 148 : integer,allocatable :: npwarr(:)
1729 148 : integer,pointer :: gvec_kss(:,:)
1730 148 : integer,pointer :: test_gvec_kss(:,:)
1731 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),qtmp(3),sq(3),qbz(3)
1732 148 : real(dp),pointer :: energies_p(:,:,:)
1733 148 : real(dp),allocatable :: doccde(:),eigen(:),occfact(:)
1734 148 : type(Pawrhoij_type),allocatable :: Pawrhoij(:)
1735 :
1736 : ! *************************************************************************
1737 :
1738 444 : units = [std_out, ab_out]
1739 :
1740 : ! Check for calculations that are not implemented
1741 1009 : ltest = ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol) == Dtset%nband(1))
1742 148 : ABI_CHECK(ltest, 'dtset%nband(:) must be constant in the GW code.')
1743 :
1744 148 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1745 :
1746 148 : call mkrdim(acell,rprim,rprimd)
1747 148 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1748 :
1749 : ! Set up basic parameters of the calculation
1750 148 : Ep%gwcalctyp =Dtset%gwcalctyp
1751 148 : Ep%plasmon_pole_model =.TRUE.
1752 148 : Ep%analytic_continuation=.FALSE.
1753 148 : Ep%contour_deformation =.FALSE.
1754 :
1755 148 : mod10=MOD(Ep%gwcalctyp,10)
1756 148 : if (mod10/=0.and.mod10/=8) Ep%plasmon_pole_model =.FALSE.
1757 148 : if (mod10==1) Ep%analytic_continuation=.TRUE.
1758 148 : if (mod10==2.or.mod10==9) Ep%contour_deformation =.TRUE.
1759 148 : is_static=(mod10==5.or.mod10==6.or.mod10==7)
1760 :
1761 148 : Ep%nbnds =Dtset%nband(1)
1762 148 : Ep%symchi =Dtset%symchi
1763 148 : Ep%inclvkb=Dtset%inclvkb; if (Dtset%usepaw/=0) Ep%inclvkb=0
1764 148 : Ep%zcut =Dtset%zcut
1765 :
1766 148 : write(msg,'(2a,i4,2a,f10.6,a)')ch10,&
1767 148 : ' GW calculation type = ',Ep%gwcalctyp,ch10,&
1768 296 : ' zcut to avoid poles in chi0 [eV] = ',Ep%zcut*Ha_eV,ch10
1769 148 : call wrtout(std_out, msg)
1770 :
1771 148 : Ep%awtr =Dtset%awtr
1772 148 : Ep%npwe =Dtset%npweps
1773 148 : Ep%npwwfn=Dtset%npwwfn
1774 148 : Ep%npwvec=MAX(Ep%npwe,Ep%npwwfn)
1775 :
1776 148 : timrev = 2 ! This information is not reported in the header
1777 : ! 1 --> do not use time-reversal symmetry
1778 : ! 2 --> take advantage of time-reversal symmetry
1779 148 : if (any(dtset%kptopt == [3, 4])) timrev = 1
1780 :
1781 148 : if (timrev==1.and.Dtset%awtr/=0) then
1782 0 : ABI_ERROR("awtr/=0 cannot be used when time-reversal symmetry doesn't hold")
1783 : end if
1784 :
1785 : ! Read parameters from WFK and verifify them.
1786 148 : call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm)
1787 1009 : mband = MAXVAL(Hdr_wfk%nband)
1788 148 : call hdr_wfk%vs_dtset(dtset)
1789 148 : remove_inv=.FALSE.
1790 :
1791 148 : test_npwkss = 0
1792 : call make_gvec_kss(Dtset%nkpt,Dtset%kptns,Hdr_wfk%ecut_eff,Dtset%symmorphi,Dtset%nsym,Dtset%symrel,Dtset%tnons,&
1793 148 : gprimd,Dtset%prtvol,test_npwkss,test_gvec_kss,ierr)
1794 148 : ABI_CHECK(ierr==0,"Fatal error in make_gvec_kss")
1795 :
1796 444 : ABI_MALLOC(gvec_kss,(3,test_npwkss))
1797 959604 : gvec_kss = test_gvec_kss
1798 148 : ng_kss = test_npwkss
1799 :
1800 148 : if (Ep%npwvec>ng_kss) then
1801 1 : Ep%npwvec=ng_kss
1802 1 : if (Ep%npwwfn> ng_kss) Ep%npwwfn=ng_kss
1803 1 : if (Ep%npwe > ng_kss) Ep%npwe =ng_kss
1804 1 : write(msg,'(3a,3(a,i6,a))')ch10,&
1805 1 : ' Number of G-vectors found less then required. Calculation will proceed with ',ch10,&
1806 1 : ' npwvec = ',Ep%npwvec,ch10,&
1807 1 : ' npweps = ',Ep%npwe ,ch10,&
1808 2 : ' npwwfn = ',Ep%npwwfn,ch10
1809 1 : ABI_WARNING(msg)
1810 : end if
1811 :
1812 148 : ng = MIN(SIZE(gvec_kss,DIM=2),SIZE(test_gvec_kss,DIM=2))
1813 148 : ierr = 0
1814 120080 : do ig=1,ng
1815 479876 : if (ANY(gvec_kss(:,ig)/=test_gvec_kss(:,ig))) then
1816 0 : ierr=ierr+1
1817 0 : write(std_out,*)" gvec_kss ",ig,"/",ng,gvec_kss(:,ig),test_gvec_kss(:,ig)
1818 : end if
1819 : end do
1820 148 : ABI_CHECK(ierr == 0, "Mismatch between gvec_kss and test_gvec_kss")
1821 148 : ABI_FREE(test_gvec_kss)
1822 :
1823 : ! Get important dimension from Hdr_wfk
1824 : ! Check also the consistency btw Hdr_wfk and Dtset.
1825 148 : Ep%nsppol=Hdr_wfk%nsppol
1826 148 : Ep%nkibz =Hdr_wfk%nkpt
1827 :
1828 148 : if (Ep%nbnds>mband) then
1829 : write(msg,'(2a,2(a,i0))') &
1830 0 : 'Number of bands stored in WFK file is less than required. ',ch10,&
1831 0 : "WFK mband:", mband, ", chi0 nband: ", Ep%nbnds
1832 0 : ABI_ERROR(msg)
1833 : end if
1834 :
1835 148 : cryst = Hdr_wfk%get_crystal(gw_timrev=timrev, remove_inv=remove_inv)
1836 148 : call cryst%print(mode_paral='COLL')
1837 :
1838 : ! === Create basic data types for the calculation ===
1839 : ! Kmesh defines the k-point sampling for the wavefunctions.
1840 : ! Qmesh defines the q-point sampling for chi0, all possible differences k1-k2 reduced to the IBZ.
1841 : ! TODO Kmesh%bz should be in [-half, half[ but this modification is painful!
1842 :
1843 148 : call Kmesh%init(cryst, ep%nkibz, hdr_wfk%kptns, dtset%kptopt, wrap_1zone=.FALSE.)
1844 :
1845 : ! Some required information are not filled up inside kmesh_init. So doing it here, even though it is not clean
1846 1924 : Kmesh%kptrlatt(:,:) = Dtset%kptrlatt(:,:)
1847 148 : Kmesh%nshift = Dtset%nshiftk
1848 444 : ABI_MALLOC(Kmesh%shift, (3,Kmesh%nshift))
1849 740 : Kmesh%shift(:,:) = Dtset%shiftk(:,1:Dtset%nshiftk)
1850 148 : call Kmesh%print(units, header="K-mesh for the wavefunctions", prtvol=Dtset%prtvol)
1851 :
1852 : ! === Find Q-mesh ===
1853 : ! Stop if a nonzero umklapp is needed to reconstruct the BZ.
1854 : ! epsilon^-1(Sq) indeed should be symmetrized in csigme using a different expression (G-G_o is needed)
1855 148 : call qmesh%find_qmesh(Cryst, Kmesh)
1856 148 : call qmesh%print(units, "Q-mesh for the screening function", prtvol=dtset%prtvol)
1857 :
1858 7084 : do iqbz=1,Qmesh%nbz
1859 6936 : call qmesh%get_BZ_item(iqbz, qbz, iq_ibz, isym, itim)
1860 194208 : sq = (3-2*itim) * MATMUL(cryst%symrec(:,:,isym), qmesh%ibz(:,iq_ibz))
1861 27892 : if (ANY(ABS(qbz-sq) > 1.0d-4)) then
1862 : write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')&
1863 0 : ' qpoint ',qbz,' is the symmetric of ',qmesh%ibz(:,iq_ibz),ch10,&
1864 0 : ' through operation ',cryst%symrec(:,:,isym),' and itim ',itim,ch10,&
1865 0 : ' however a non zero umklapp G_o vector is required and this is not yet allowed'
1866 0 : ABI_ERROR(msg)
1867 : end if
1868 : end do
1869 :
1870 : ! Setup for long wavelength limit
1871 148 : if (Dtset%gw_nqlwl==0) then
1872 147 : Ep%nqlwl=1
1873 147 : ABI_MALLOC(Ep%qlwl,(3,Ep%nqlwl))
1874 588 : Ep%qlwl(:,1)=GW_Q0_DEFAULT ! Use default shift 0.000010, 0.000020, 0.000030
1875 : else
1876 1 : Ep%nqlwl=Dtset%gw_nqlwl
1877 3 : ABI_MALLOC(Ep%qlwl,(3,Ep%nqlwl))
1878 5 : Ep%qlwl(:,:)=Dtset%gw_qlwl(:,1:Ep%nqlwl)
1879 1 : ABI_CHECK(Ep%nqlwl==1,"nqlwl/=1 not coded")
1880 : end if
1881 : !write(std_out,*)" Using qlwl = ",Ep%qlwl
1882 :
1883 : ! Find optimal value for G-sphere enlargment due to oscillator matrix elements
1884 148 : call get_ng0sh(Kmesh%nbz,Kmesh%bz,Qmesh%nibz,Qmesh%ibz,Kmesh%nbz,Kmesh%bz,GW_TOLQ0,ng0sh_opt)
1885 148 : call wrtout(std_out,sjoin(' Optimal value for ng0sh:',ltoa(ng0sh_opt)),"COLL")
1886 :
1887 592 : Ep%mG0(:)=ng0sh_opt(:) !Ep%mG0(:) = [3, 3, 3]
1888 :
1889 : ! === In case of symmetrization, find the little group of the q"s ===
1890 : ! * For the long-wavelength limit we consider a small but finite q. However the oscillators are
1891 : ! evaluated setting q==0. Thus it is possible to take advantage of symmetries also when q --> 0.
1892 : ! * Here we calculate the enlargement of the G-sphere, npwepG0, needed to account for umklapps.
1893 : ! TODO Switch on use_umklp, write all this stuff to ab_out
1894 :
1895 148 : Ep%npwepG0 = Ep%npwe
1896 1378 : ABI_MALLOC(Ltg_q, (Qmesh%nibz))
1897 :
1898 1082 : do iq=1,Qmesh%nibz
1899 3736 : qtmp = Qmesh%ibz(:,iq); if (normv(qtmp,gmet,'G') < GW_TOLQ0) qtmp(:) = zero; use_umklp = 0
1900 1082 : call Ltg_q(iq)%init(qtmp, Kmesh%nbz, Kmesh%bz, Cryst, use_umklp, Ep%npwe, gvec=gvec_kss)
1901 : end do
1902 :
1903 148 : ecutepspG0 = Dtset%ecuteps
1904 148 : ABI_CHECK(ecutepspG0 > zero, "ecuteps must be > 0")
1905 148 : if (Ep%symchi/=0) then
1906 1117 : ecutepspG0=MAXVAL(Ltg_q(:)%max_kin_gmG0)+tol6; npwepG0=0; nshepspG0=0
1907 134 : if (my_rank == master) write(std_out,*)" Due to umklapp processes : ecutepspg0= ",ecutepspG0
1908 134 : call setshells(ecutepspG0,npwepG0,nshepspG0,Cryst%nsym,gmet,gprimd,Cryst%symrel,'eps_pG0',Cryst%ucvol)
1909 134 : Ep%npwepG0=npwepG0
1910 : end if
1911 :
1912 148 : if (Ep%npwepG0>Ep%npwvec) then
1913 : write(msg,'(3a,i5,a,i5)')&
1914 0 : ' npwepG0 > npwvec, decrease npweps or increase npwwfn. ',ch10,&
1915 0 : ' npwepG0 = ',Ep%npwepG0,' npwvec = ',Ep%npwvec
1916 0 : ABI_ERROR(msg)
1917 : end if
1918 :
1919 : ! === Create structure describing the G-sphere used for chi0/espilon and Wfns ===
1920 : ! * The cutoff is >= ecuteps to allow for umklapp
1921 148 : call Gsph_wfn%init(Cryst, Ep%npwvec, gvec=gvec_kss)
1922 :
1923 148 : call Gsph_epsG0%init(Cryst, Ep%npwepG0, gvec=gvec_kss)
1924 : !
1925 : ! =======================================================================
1926 : ! ==== Setup of the FFT mesh used for the oscillator matrix elements ====
1927 : ! =======================================================================
1928 : ! * ngfft_gw(7:18) is the same as Dtset%ngfft(7:18), initialized before entering setup_screening.
1929 : ! Here we just redefine ngfft_gw(1:6) according to the following options:
1930 : !
1931 : ! method==0 ==> FFT grid read from __fft.in__ (only for debugging purpose)
1932 : ! method==1 ==> normal FFT grid
1933 : ! method==2 ==> slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90)
1934 : ! method==3 ==> doubled FFT grid, to treat exactly the convolution defining the density,
1935 : ! Useful in sigma if ppmodel=[2,3,4] since rho(G-Gp) or to calculate matrix elements of v_Hxc.
1936 : !
1937 : ! enforce_sym==1 ==> enforce a direct space FFT mesh compatible with all symmetries operation
1938 : ! enforce_sym==0 ==> Find the smallest FFT grid compatibile with the library, do not care about symmetries
1939 : !
1940 2812 : ngfft_gw(1:18)=Dtset%ngfft(1:18); method=2
1941 148 : if (Dtset%fftgw==00 .or. Dtset%fftgw==01) method=0
1942 148 : if (Dtset%fftgw==10 .or. Dtset%fftgw==11) method=1
1943 148 : if (Dtset%fftgw==20 .or. Dtset%fftgw==21) method=2
1944 148 : if (Dtset%fftgw==30 .or. Dtset%fftgw==31) method=3
1945 148 : enforce_sym=MOD(Dtset%fftgw,10)
1946 :
1947 : ! Use npwepG0 to account for umklapps.
1948 148 : call setmesh(gmet,gvec_kss,ngfft_gw,Ep%npwvec,Ep%npwepG0,Ep%npwwfn,nfftgw_tot,method,Ep%mG0,Cryst,enforce_sym)
1949 : !call new_setmesh(Cryst,ecut_osc,ecutwfn,nkpt,kpoints,method,Ep%mG0,enforce_sym,ngfft_gw,nfftgw_tot)
1950 :
1951 148 : ABI_FREE(gvec_kss)
1952 :
1953 : ! FIXME this wont work if nqptdm/=0
1954 : call Vcp%init(Gsph_epsG0,Cryst,Qmesh,Kmesh,Dtset%gw_rcut,Dtset%gw_icutcoul,Dtset%vcutgeo,Dtset%ecuteps,Ep%npwe,Ep%nqlwl,&
1955 148 : Ep%qlwl,comm)
1956 :
1957 : #if 0
1958 : ! Using the random q for the optical limit is one of the reasons
1959 : ! why sigma breaks the initial energy degeneracies.
1960 : Vcp%i_sz=zero
1961 : Vcp%vc_sqrt(1,:)=czero
1962 : Vcp%vcqlwl_sqrt(1,:)=czero
1963 : #endif
1964 :
1965 : ! Value of scissor energy
1966 148 : Ep%mbpt_sciss=Dtset%mbpt_sciss
1967 :
1968 : ! Define the frequency mesh for epsilon according to the method used.
1969 148 : Ep%nomegaei=1
1970 148 : Ep%nomegaer=1; if (is_static) Ep%nomegaer=0
1971 148 : Ep%nomegaec=0
1972 148 : Ep%omegaermax=zero
1973 :
1974 : ! For ppmodels 2,3,4, only omega=0 is needed.
1975 148 : if (Ep%plasmon_pole_model.and.Dtset%nfreqre==1.and.Dtset%nfreqim==0) then
1976 5 : Ep%nomegaer=1; Ep%nomegaei=0
1977 5 : write(msg,'(7a)')ch10,&
1978 5 : ' The inverse dielectric matrix will be calculated on zero frequency only',ch10,&
1979 5 : ' please note that the calculated epsilon^-1 cannot be used ',ch10,&
1980 10 : ' to calculate QP corrections using plasmonpole model 1',ch10
1981 5 : call wrtout(units, msg)
1982 : end if
1983 :
1984 : ! Max number of omega along the imaginary axis
1985 148 : if (Ep%analytic_continuation.or.Ep%contour_deformation) then
1986 37 : Ep%nomegaei = abs(Dtset%nfreqim)
1987 37 : if (Dtset%gw_frqim_inzgrid==1) then
1988 1 : ABI_WARNING('iomega = z/1-z transfom grid will be used for imaginary frequency grid')
1989 : end if
1990 37 : if (Dtset%cd_customnimfrqs/=0) then
1991 1 : ABI_WARNING('Custom imaginary grid specified. Assuming experienced user.')
1992 1 : Ep%nomegaei=Dtset%cd_customnimfrqs
1993 : end if
1994 37 : if (Ep%nomegaei==-1) then
1995 0 : Ep%nomegaei=NOMEGAGAUSS
1996 0 : ABI_WARNING(sjoin('Number of imaginary frequencies set to default= ',itoa(NOMEGAGAUSS)))
1997 : end if
1998 37 : if (Ep%nomegaei==0) then
1999 3 : ABI_WARNING('nfreqim = 0! Assuming experienced user merging several frequency calculations.')
2000 : end if
2001 : end if
2002 :
2003 : ! Range and total number of real frequencies.
2004 148 : Ep%omegaermin = zero
2005 148 : if (Ep%contour_deformation) then
2006 26 : Ep%nomegaer=Dtset%nfreqre; Ep%omegaermin=Dtset%freqremin; Ep%omegaermax=Dtset%freqremax
2007 26 : if (Dtset%gw_frqre_tangrid==1) then
2008 1 : Ep%omegaermax=Dtset%cd_max_freq
2009 1 : ABI_WARNING('Tangent transfom grid will be used for real frequency grid')
2010 : end if
2011 26 : if (Dtset%gw_frqre_tangrid==1) then
2012 1 : ABI_WARNING('Tangent transfom grid will be used for real frequency grid')
2013 : end if
2014 26 : if (Ep%nomegaer==-1) then
2015 0 : Ep%nomegaer=NOMEGAREAL
2016 0 : ABI_WARNING(sjoin('Number of real frequencies set to default= ',itoa(NOMEGAREAL)))
2017 : end if
2018 26 : if (Ep%nomegaer==0) then
2019 0 : ABI_WARNING('nfreqre = 0 ! Assuming experienced user.')
2020 : end if
2021 26 : if (ABS(Ep%omegaermin)<TOL16) then
2022 25 : Ep%omegaermin=zero
2023 25 : write(msg,'(a,f8.4)')' Min real frequency set to default [Ha] = ',Ep%omegaermin
2024 25 : ABI_WARNING(msg)
2025 : end if
2026 26 : if (Ep%omegaermin>Ep%omegaermax) then
2027 0 : ABI_ERROR('freqremin > freqremax !')
2028 : end if
2029 26 : if (Ep%omegaermax<TOL16) then
2030 4 : Ep%omegaermax=OMEGAERMAX
2031 4 : write(msg,'(a,f8.4)')' Max real frequency set to default [Ha] = ',OMEGAERMAX
2032 4 : ABI_WARNING(msg)
2033 : end if
2034 : ! Check if a subset of the frequencies is to be used
2035 26 : if (Dtset%cd_subset_freq(1)/=0) then
2036 0 : istart = Dtset%cd_subset_freq(1)
2037 0 : iend = Dtset%cd_subset_freq(2)
2038 0 : if (istart>iend.or.istart<0.or.iend<0) then
2039 0 : ABI_ERROR(' check indices of cd_subset_freq!')
2040 : end if
2041 0 : write(msg,'(2(a,i0))')' Using cd_subset_freq to only do freq. from ',istart,' to ',iend
2042 0 : ABI_WARNING(msg)
2043 : ! Reset the numbers
2044 0 : if (Dtset%gw_frqre_tangrid/=1) then ! Normal equidistant grid
2045 0 : Ep%nomegaer = iend-istart+1
2046 0 : domegareal=(Ep%omegaermax-Ep%omegaermin)/(Ep%nomegaer-1)
2047 0 : Ep%omegaermin = Ep%omegaermin+(istart-1)*domegareal
2048 0 : Ep%omegaermax = Ep%omegaermin+(iend-1)*domegareal
2049 : else
2050 0 : Ep%nomegaer = iend-istart+1
2051 : end if
2052 : end if
2053 : end if
2054 :
2055 : ! Check full grid calculations
2056 148 : if (Dtset%cd_full_grid/=0) then
2057 1 : ABI_WARNING("FULL GRID IN COMPLEX PLANE CALCULATED. YOU MIGHT NOT BE ABLE TO USE SCREENING FILES!")
2058 1 : if (Dtset%cd_subset_freq(1)/=0) then
2059 0 : ABI_ERROR('cd_subset_freq cannot be used with cd_full_grid!')
2060 : end if
2061 1 : Ep%nomegaec = Ep%nomegaei*(Ep%nomegaer-1)
2062 : end if
2063 :
2064 148 : Ep%nomega=Ep%nomegaer+Ep%nomegaei+Ep%nomegaec ! Total number of frequencies.
2065 :
2066 : ! ==== Setup of the spectral method ====
2067 148 : Ep%spmeth =Dtset%spmeth; Ep%nomegasf=Dtset%nomegasf; Ep%spsmear =Dtset%spbroad
2068 :
2069 148 : if (Ep%spmeth/=0) then
2070 6 : write(msg,'(2a,i3,2a,i8)')ch10,&
2071 6 : ' setup_screening: using spectral method: ',Ep%spmeth,ch10,&
2072 12 : ' Number of frequencies for imaginary part: ',Ep%nomegasf
2073 6 : call wrtout(std_out, msg)
2074 6 : if (Ep%spmeth==2) then
2075 1 : write(msg,'(a,f8.5,a)')' Gaussian broadening = ',Ep%spsmear*Ha_eV,' [eV]'
2076 1 : call wrtout(std_out, msg)
2077 : end if
2078 : end if
2079 :
2080 148 : Ep%nI=1; Ep%nJ=1
2081 148 : if (Dtset%nspinor==2) then
2082 : !if (Dtset%usepaw==1.and.Dtset%pawspnorb>0) then
2083 : ! Ep%nI=1; Ep%nJ=4
2084 : !end if
2085 : ! For spin-spin interaction
2086 : ! Ep%nI=4; Ep%nJ=4
2087 2 : ABI_CHECK(Ep%npwepG0 == Ep%npwe, "npwepG0 must be == npwe if nspinor==2")
2088 : !ABI_CHECK(Ep%symchi == 0, "symchi/=0 and nspinor=2 not available")
2089 : end if
2090 :
2091 : ! === Enable the calculations of chi0 on user-specified q-points ===
2092 148 : Ep%nqibz=Qmesh%nibz
2093 444 : ABI_MALLOC(Ep%qibz,(3,Ep%nqibz))
2094 3884 : Ep%qibz(:,:)=Qmesh%ibz(:,:)
2095 :
2096 148 : Ep%nqcalc=Ep%nqibz
2097 148 : if (Dtset%nqptdm>0) Ep%nqcalc=Dtset%nqptdm
2098 :
2099 444 : ABI_MALLOC(Ep%qcalc,(3,Ep%nqcalc))
2100 148 : if (Ep%nqcalc/=Ep%nqibz) then
2101 3 : write(msg,'(6a)')ch10,&
2102 3 : ' Dielectric matrix will be calculated only for some ',ch10,&
2103 3 : ' selected q points provided by the user through the input variables ',ch10,&
2104 6 : ' nqptdm and qptdm'
2105 3 : call wrtout(units, msg)
2106 3 : ltest= Ep%nqcalc <= Qmesh%nibz
2107 3 : ABI_CHECK(ltest, 'nqptdm should not exceed the number of q points in the IBZ')
2108 27 : Ep%qcalc(:,:)=Dtset%qptdm(:,1:Ep%nqcalc)
2109 : ! Check whether the q-points provided are correct.
2110 9 : do iq=1,Ep%nqcalc
2111 6 : found=.FALSE.
2112 21 : do iqp=1,Qmesh%nibz
2113 84 : qtmp(:)=Ep%qcalc(:,iq)-Qmesh%ibz(:,iqp)
2114 21 : found=(normv(qtmp,gmet,'G')<GW_TOLQ)
2115 21 : if (found) EXIT
2116 : end do
2117 9 : ABI_CHECK(found, 'One or more points specified by Dtset%qptdm do not satisfy q=k1-k2')
2118 : end do
2119 : else
2120 3809 : Ep%qcalc(:,:)=Ep%qibz(:,:)
2121 : end if
2122 :
2123 : ! To write the SCR header correctly, with heads and wings, we have
2124 : ! to make sure that q==0, if present, is the first q-point in the list.
2125 148 : has_q0=.FALSE.
2126 153 : do iq=1,Ep%nqcalc
2127 153 : if (normv(Ep%qcalc(:,iq),gmet,'G')<GW_TOLQ0) then
2128 148 : has_q0=.TRUE.; EXIT
2129 : end if
2130 : end do
2131 :
2132 148 : if (has_q0.and.normv(Ep%qcalc(:,1),gmet,'G')>=GW_TOLQ0) then
2133 : write(msg,'(5a)')&
2134 0 : 'The list of q-points to be calculated contains the Gamma point, ',ch10,&
2135 0 : 'however Gamma is not the first point in the list. ' ,ch10,&
2136 0 : 'Please, change your input file accordingly. '
2137 0 : ABI_ERROR(msg)
2138 : end if
2139 :
2140 : ! === Initialize the band structure datatype ===
2141 : ! * Copy KSS energies and occupations up to Ep%nbnds==Dtset%nband(:)
2142 : ! TODO Recheck symmorphy and inversion
2143 1009 : bantot = SUM(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol))
2144 :
2145 22732 : ABI_CALLOC(doccde, (bantot))
2146 22584 : ABI_CALLOC(eigen, (bantot))
2147 22584 : ABI_CALLOC(occfact, (bantot))
2148 :
2149 : jj=0; ibtot=0
2150 298 : do isppol=1,Dtset%nsppol
2151 1159 : do ikibz=1,Dtset%nkpt
2152 37761 : do ib=1,Hdr_wfk%nband(ikibz+Dtset%nkpt*(isppol-1))
2153 36750 : ibtot=ibtot+1
2154 37611 : if (ib<=Ep%nbnds) then
2155 22288 : jj=jj+1
2156 22288 : occfact(jj)=Hdr_wfk%occ(ibtot)
2157 22288 : eigen (jj)=energies_p(ib,ikibz,isppol)
2158 : end if
2159 : end do
2160 : end do
2161 : end do
2162 148 : ABI_FREE(energies_p)
2163 :
2164 : ! Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of
2165 : ! the symmetry operations in the old GW code (symmorphy and inversion)
2166 1001 : ltest = (ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz)) < tol6))
2167 148 : if (.not. ltest) then
2168 0 : do jj=1,Kmesh%nibz
2169 0 : write(std_out, *)"wtk dtset vs kmesh:", dtset%wtk(jj), kmesh%wt(jj)
2170 : end do
2171 : end if
2172 0 : ABI_CHECK(ltest, 'Mismatch between Dtset%wtk and Kmesh%wt')
2173 :
2174 444 : ABI_MALLOC(npwarr,(Hdr_wfk%nkpt))
2175 1001 : npwarr(:)=Ep%npwwfn
2176 :
2177 : call ks_ebands%init(bantot, Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,&
2178 : doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,&
2179 : Kmesh%nibz,npwarr,Dtset%nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,&
2180 : dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
2181 148 : dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
2182 :
2183 : ! TODO modify outkss in order to calculate the eigenvalues also if NSCF calculation.
2184 : ! this fails simply because in case of NSCF occ are zero
2185 : !ltest=(ALL(ABS(occfact-ks_ebands%occ)<1.d-2))
2186 : !call assert(ltest,'difference in occfact')
2187 : !write(std_out,*)MAXVAL(ABS(occfact(:)-ks_ebands%occ(:)))
2188 :
2189 : !TODO call ebands_update_occ here
2190 : !call ks_ebands%update_occ(spinmagntarget,Dtset%prtvol)
2191 :
2192 148 : ABI_FREE(doccde)
2193 148 : ABI_FREE(eigen)
2194 148 : ABI_FREE(npwarr)
2195 :
2196 : ! Initialize abinit header for the screening part
2197 148 : call Hdr_out%init(ks_ebands,codvsn,Dtset,Pawtab,pertcase0,Psps,wvl)
2198 :
2199 : ! Get Pawrhoij from the header.
2200 475 : ABI_MALLOC(Pawrhoij,(Cryst%natom*Dtset%usepaw))
2201 148 : if (Dtset%usepaw==1) then
2202 5 : call pawrhoij_alloc(Pawrhoij,1,Dtset%nspden,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab)
2203 5 : call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij)
2204 : end if
2205 148 : call Hdr_out%update(bantot,1.0d20,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1))
2206 :
2207 148 : ABI_FREE(occfact)
2208 148 : call pawrhoij_free(Pawrhoij)
2209 179 : ABI_FREE(Pawrhoij)
2210 :
2211 : ! ==== Setup of extrapolar technique ====
2212 148 : Ep%gwcomp = Dtset%gwcomp; Ep%gwencomp = Dtset%gwencomp
2213 148 : if (Ep%gwcomp == 1) then
2214 19 : write(msg,'(a,f8.2,a)')' Using the completeness correction with gwencomp ',Ep%gwencomp*Ha_eV,' [eV] '
2215 19 : call wrtout(std_out, msg)
2216 : end if
2217 :
2218 : ! Final compatibility tests
2219 1009 : ltest = (ks_ebands%mband == Ep%nbnds .and. ALL(ks_ebands%nband == Ep%nbnds))
2220 0 : ABI_CHECK(ltest, 'BUG in definition of ks_ebands%nband')
2221 :
2222 148 : if (Ep%gwcomp==1 .and. Ep%spmeth>0) then
2223 0 : ABI_ERROR("Hilbert transform and extrapolar method are not compatible")
2224 : end if
2225 :
2226 444 : end subroutine setup_screening
2227 : !!***
2228 :
2229 : !----------------------------------------------------------------------
2230 :
2231 : !!****f* m_screening_driver/chi0_bksmask
2232 : !! NAME
2233 : !! chi0_bksmask
2234 : !!
2235 : !! FUNCTION
2236 : !! Compute tables for the distribution and the storage of the wavefunctions in the SCREENING code.
2237 : !!
2238 : !! INPUTS
2239 : !! Dtset<type(dataset_type)>=all input variables for this dataset
2240 : !! Ep<em1params_t>=Parameters for the screening calculation.
2241 : !! Kmesh <kmesh_t>=Structure describing the k-point sampling.
2242 : !! nbvw = Max. number of fully/partially occupied states over spin
2243 : !! nbcw = Max. number of unoccupied states considering the spin
2244 : !! nprocs=Total number of MPI processors
2245 : !! my_rank=Rank of this this processor.
2246 : !!
2247 : !! OUTPUT
2248 : !! bks_mask(Ep%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will treat this state.
2249 : !! keep_ur(Ep%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will store this state in real space.
2250 : !! ierr=Exit status.
2251 : !!
2252 : !! SOURCE
2253 :
2254 148 : subroutine chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,nprocs,bks_mask,keep_ur,ierr)
2255 :
2256 : !Arguments ------------------------------------
2257 : !scalars
2258 : integer,intent(in) :: my_rank,nprocs,nbvw,nbcw
2259 : integer,intent(out) :: ierr
2260 : type(Dataset_type),intent(in) :: Dtset
2261 : type(em1params_t),intent(in) :: Ep
2262 : type(kmesh_t),intent(in) :: Kmesh
2263 : !arrays
2264 : logical,intent(out) :: bks_mask(Ep%nbnds,Kmesh%nibz,Dtset%nsppol)
2265 : logical,intent(out) :: keep_ur(Ep%nbnds,Kmesh%nibz,Dtset%nsppol)
2266 :
2267 : !Local variables-------------------------------
2268 : !scalars
2269 : integer :: my_nspins,my_maxb,my_minb,isp,spin,nsppol,band,rank_spin,ib
2270 : character(len=500) :: msg
2271 : logical :: store_ur
2272 : !arrays
2273 296 : integer :: my_spins(Dtset%nsppol),nprocs_spin(Dtset%nsppol)
2274 148 : integer,allocatable :: istart(:),istop(:)
2275 :
2276 : ! *************************************************************************
2277 :
2278 148 : ierr=0; nsppol=Dtset%nsppol
2279 :
2280 448 : my_nspins=Dtset%nsppol; my_spins= [(isp,isp=1,nsppol)]
2281 :
2282 : ! List of spins for each node, number of processors per each spin
2283 : ! and the MPI rank in the "spin" communicator.
2284 298 : nprocs_spin = nprocs; rank_spin = my_rank
2285 :
2286 148 : if (nsppol==2.and.nprocs>1) then
2287 : ! Distribute spins (optimal distribution if nprocs is even)
2288 0 : nprocs_spin(1) = nprocs/2
2289 0 : nprocs_spin(2) = nprocs - nprocs/2
2290 0 : my_nspins=1; my_spins(1)=1
2291 0 : if (my_rank+1>nprocs/2) then
2292 0 : my_spins(1)=2
2293 0 : rank_spin = my_rank - nprocs/2
2294 : end if
2295 : end if
2296 :
2297 148 : store_ur = (MODULO(Dtset%gwmem,10)==1)
2298 46746 : bks_mask=.FALSE.; keep_ur=.FALSE.
2299 :
2300 180 : select case (Dtset%gwpara)
2301 : case (1)
2302 : ! Parallelization over transitions **without** memory distributions (Except for the spin).
2303 32 : my_minb=1; my_maxb=Ep%nbnds
2304 64 : do isp=1,my_nspins
2305 32 : spin = my_spins(isp)
2306 3888 : bks_mask(my_minb:my_maxb,:,spin) = .TRUE.
2307 3920 : if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE.
2308 : end do
2309 :
2310 : case (2)
2311 : ! Distribute bands and spin.
2312 234 : do isp=1,my_nspins
2313 118 : spin = my_spins(isp)
2314 :
2315 118 : if (nprocs_spin(spin) <= nbcw) then
2316 : ! Distribute nbcw empty bands among nprocs_spin (block of bands without replicas).
2317 : ! Bands are distributed in contiguous blocks because
2318 : ! this distribution is well suited for the Hilber transform
2319 : ! since each node will allocate only a smaller frequency interval
2320 : ! for the spectral function whose size scales with the number of MPI nodes.
2321 : ! Note it is now meaningless to distinguish gwcomp=0 or 1 since the workload is well balanced later on
2322 354 : ABI_MALLOC(istart,(nprocs_spin(spin)))
2323 236 : ABI_MALLOC(istop,(nprocs_spin(spin)))
2324 :
2325 118 : call xmpi_split_work2_i4b(nbcw,nprocs_spin(spin),istart,istop)
2326 :
2327 118 : my_minb = nbvw + istart(rank_spin+1)
2328 118 : my_maxb = nbvw + istop (rank_spin+1)
2329 :
2330 118 : ABI_FREE(istart)
2331 118 : ABI_FREE(istop)
2332 :
2333 118 : if (my_maxb - my_minb + 1 <= 0) then
2334 : write(msg,'(3a,2(i0,a),2a)')&
2335 0 : 'One or more processors has zero number of bands ',ch10,&
2336 0 : 'my_minb= ',my_minb,' my_maxb= ',my_maxb,ch10,&
2337 0 : 'This is a waste, decrease the number of processors.'
2338 0 : ABI_ERROR(msg)
2339 : end if
2340 :
2341 13989 : bks_mask(my_minb:my_maxb,:,spin)=.TRUE.
2342 13637 : if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE.
2343 :
2344 : else
2345 : ! New version (alternate bands with replicas if nprocs > nbcw)
2346 : ! FIXME: Fix segmentation fault with Hilbert transform.
2347 0 : do ib=1,nbcw
2348 0 : if (xmpi_distrib_with_replicas(ib,nbcw,rank_spin,nprocs_spin(spin))) then
2349 0 : band = ib + nbvw
2350 0 : bks_mask(band,:,spin)=.TRUE.
2351 0 : if (store_ur) keep_ur(band,:,spin)=.TRUE.
2352 : end if
2353 : end do
2354 : end if
2355 :
2356 : ! This is needed to have all the occupied states on each node.
2357 3511 : bks_mask(1:nbvw,:,spin) = .TRUE.
2358 3519 : if (store_ur) keep_ur(1:nbvw,:,spin)=.TRUE.
2359 : end do ! isp
2360 :
2361 : case default
2362 0 : ierr = 1
2363 148 : ABI_WARNING("Wrong value for gwpara")
2364 : end select
2365 :
2366 148 : end subroutine chi0_bksmask
2367 : !!***
2368 :
2369 : !!****f* m_screening_driver/random_stopping_power
2370 : !! NAME
2371 : !! random_stopping_power
2372 : !!
2373 : !! FUNCTION
2374 : !! Calculate the electronic random stopping power
2375 : !!
2376 : !! INPUTS
2377 : !!
2378 : !! OUTPUT
2379 : !!
2380 : !! SIDE EFFECTS
2381 : !!
2382 : !! SOURCE
2383 :
2384 12 : subroutine random_stopping_power(iqibz,npvel,pvelmax,Ep,Gsph_epsG0,Qmesh,Vcp,Cryst,Dtfil,epsm1,rspower)
2385 :
2386 : use m_splines
2387 :
2388 : !Arguments ------------------------------------
2389 : !scalars
2390 : integer,intent(in) :: iqibz,npvel
2391 : real(dp),intent(in) :: pvelmax(3)
2392 : type(em1params_t),intent(in) :: Ep
2393 : type(gsphere_t),intent(in) :: Gsph_epsG0
2394 : type(kmesh_t),intent(in) :: Qmesh
2395 : type(vcoul_t),intent(in) :: Vcp
2396 : type(crystal_t),intent(in) :: Cryst
2397 : type(Datafiles_type),intent(in) :: Dtfil
2398 : complex(gwp),intent(in) :: epsm1(Ep%npwe,Ep%npwe,Ep%nomega)
2399 : real(dp),intent(inout) :: rspower(npvel)
2400 :
2401 : !Local variables ------------------------------
2402 : integer :: ipvel,ig, iq_bz,iq_ibz,isym_q,itim_q
2403 : integer :: iomega,iomegap,nomega_re,unt_rsp, units(2)
2404 12 : integer,allocatable :: iomega_re(:)
2405 :
2406 : real(dp),parameter :: zp=1.0_dp ! Hard-coded charge of the impinging particle
2407 : real(dp) :: omega_p
2408 : real(dp) :: im_epsm1_int(1), qbz(3),qpgcart(3),qpg_red(3)
2409 24 : real(dp) :: pvel(3,npvel),pvel_norm(npvel),ypp_i(Ep%nomega),vcoul(Ep%npwe)
2410 12 : real(dp),allocatable :: im_epsm1_diag_qbz(:,:),tmp_data(:),omega_re(:)
2411 : character(len=500) :: msg
2412 : character(len=fnlen+4) :: fname
2413 : !************************************************************************
2414 :
2415 36 : units = [std_out, ab_out]
2416 :
2417 : ! First set up the velocities array from the input variables npvel and pvelmax(3)
2418 : ! Remember pvelmax is in Cartesian coordinates and so is pvel
2419 492 : do ipvel=1,npvel
2420 1920 : pvel(:,ipvel) = REAL(ipvel,dp) / REAL(npvel,dp) * pvelmax(:)
2421 1932 : pvel_norm(ipvel) = SQRT( SUM( pvel(:,ipvel)**2 ) )
2422 : enddo
2423 : !
2424 : ! Select the purely real frequency in Ep%omega
2425 12 : nomega_re=0
2426 492 : do iomega=1,Ep%nomega
2427 492 : if( AIMAG(Ep%omega(iomega)) < 1.0e-4_dp ) then
2428 480 : nomega_re=nomega_re+1
2429 : endif
2430 : enddo
2431 36 : ABI_MALLOC(omega_re,(nomega_re))
2432 36 : ABI_MALLOC(iomega_re,(nomega_re))
2433 48 : ABI_MALLOC(im_epsm1_diag_qbz,(Ep%npwe,Ep%nomega))
2434 36 : ABI_MALLOC(tmp_data,(Ep%nomega))
2435 :
2436 492 : iomegap=0
2437 492 : do iomega=1,Ep%nomega
2438 492 : if( AIMAG(Ep%omega(iomega)) < 1.0e-4_dp ) then
2439 480 : iomegap=iomegap+1
2440 480 : iomega_re(iomegap)=iomega
2441 480 : omega_re(iomegap)=REAL(Ep%omega(iomega),dp)
2442 : endif
2443 : enddo
2444 :
2445 : ! Loop over all the q-points in the full Brillouin zone and select only the
2446 : ! ones that corresponds to the correct q-point in the irreducible wedge we are
2447 : ! currently treating (index iqibz)
2448 396 : do iq_bz=1,Qmesh%nbz
2449 :
2450 : ! Perform the check and obtain the symmetry information
2451 384 : call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q)
2452 384 : if( iqibz /= iq_ibz ) cycle
2453 :
2454 : ! Apply the symmetry operation to the diagonal of epsm1
2455 2624 : do iomega=1,nomega_re
2456 41024 : do ig=1,Ep%npwe
2457 40960 : im_epsm1_diag_qbz(Gsph_epsG0%rottb(ig,itim_q,isym_q),iomega)= AIMAG( epsm1(ig,ig,iomega_re(iomega)) )
2458 : enddo
2459 : enddo
2460 : ! Apply the symmetry operation to the Coulomb interaction
2461 1024 : do ig=1,Ep%npwe
2462 1024 : vcoul(Gsph_epsG0%rottb(ig,itim_q,isym_q))=Vcp%vc_sqrt(ig,iqibz)**2
2463 : enddo
2464 :
2465 : ! Sum over G vectors
2466 1036 : do ig=1,Ep%npwe
2467 : !
2468 : ! Loop over velocities
2469 39744 : do ipvel=1,npvel
2470 :
2471 153600 : qpg_red(:) = qbz(:) + Gsph_epsG0%gvec(:,ig)
2472 : ! Transform q + G from reduced to cartesian with the symmetry operation
2473 : qpgcart(:) = two_pi * Cryst%gprimd(:,1) * qpg_red(1) &
2474 : + two_pi * Cryst%gprimd(:,2) * qpg_red(2) &
2475 153600 : + two_pi * Cryst%gprimd(:,3) * qpg_red(3)
2476 :
2477 : ! omega_p = ( q + G ) . v
2478 153600 : omega_p = DOT_PRODUCT( qpgcart(:) , pvel(:,ipvel) )
2479 :
2480 : ! Check that the calculated frequency omega_p is within the omega
2481 : ! range of epsm1 and thus that the interpolation will go fine
2482 38400 : if ( ABS(omega_p) > omega_re(nomega_re) ) then
2483 0 : write(msg,'(a,e16.4,2a,e16.4)') ' freqremax is currently ',omega_re(nomega_re),ch10,&
2484 0 : ' increase it to at least ',omega_p
2485 0 : ABI_WARNING(msg)
2486 : endif
2487 :
2488 : ! Perform the spline interpolation to obtain epsm1 at the desired omega = omega_p
2489 1612800 : tmp_data = im_epsm1_diag_qbz(ig,:)
2490 :
2491 38400 : call spline( omega_re, tmp_data, nomega_re, 1.0e+32_dp, 1.0e+32_dp, ypp_i)
2492 76800 : call splint( nomega_re, omega_re, tmp_data, ypp_i, 1, (/ ABS(omega_p) /), im_epsm1_int )
2493 :
2494 : ! Apply the odd parity of Im epsm1 in omega to recover the causal response function
2495 38400 : if (omega_p<zero) im_epsm1_int(1)=-im_epsm1_int(1)
2496 :
2497 : ! Calculate 4 * pi / |q+G|**2 * omega_p * Im{ epsm1_GG(q,omega_p) }
2498 38400 : im_epsm1_int(1) = omega_p * vcoul(ig) * im_epsm1_int(1)
2499 :
2500 : ! Accumulate the final result without the prefactor
2501 : ! (It will be included at the very end)
2502 39360 : rspower(ipvel) = rspower(ipvel) + im_epsm1_int(1)
2503 :
2504 : end do ! end of velocity loop
2505 : end do ! end G-loop
2506 :
2507 : enddo ! end of q loop in the full BZ
2508 :
2509 : ! If it is the last q, write down the result in the main output file and in a
2510 : ! separate file _RSP (for Random Stopping Power)
2511 12 : if (iqibz == Qmesh%nibz ) then
2512 :
2513 : ! Multiply by the prefactors
2514 : ! Note that this expression differs from Eq. (3.11) in Campillo PRB 58, 10307 (1998) [[cite:Campillo1998]].
2515 : ! A factor one half is missing in the paper.
2516 82 : rspower(:) = - zp**2 / ( Cryst%ucvol * Qmesh%nbz * pvel_norm(:) ) * rspower(:)
2517 :
2518 2 : write(msg,'(2a)') ch10,' ==== Random stopping power along Cartesian direction === '
2519 2 : call wrtout(units, msg)
2520 2 : write(msg,'(a,3(f12.4,2x),a)') ' ==== ',pvelmax(:),'===='
2521 2 : call wrtout(units, msg)
2522 2 : write(msg,'(a)') '# |v| (a.u.) , RSP (a.u.) '
2523 2 : call wrtout(units, msg)
2524 82 : do ipvel=1,npvel
2525 80 : write(msg,'(f16.8,4x,f16.8)') pvel_norm(ipvel),rspower(ipvel)
2526 82 : call wrtout(units, msg)
2527 : enddo
2528 2 : write(msg,'(2a)') ' ========================================================= ',ch10
2529 2 : call wrtout(units, msg)
2530 :
2531 2 : fname=TRIM(Dtfil%filnam_ds(4))//'_RSP'
2532 2 : if (open_file(fname,msg,newunit=unt_rsp,status='unknown',form='formatted') /= 0) then
2533 0 : ABI_ERROR(msg)
2534 : end if
2535 :
2536 2 : write(msg,'(a)') '# ==== Random stopping power along Cartesian direction === '
2537 2 : call wrtout(unt_rsp, msg)
2538 2 : write(msg,'(a,3(f12.4,2x))') '# ==== ',pvelmax(:)
2539 2 : call wrtout(unt_rsp, msg)
2540 2 : write(msg,'(a)') '# |v| (a.u.) , RSP (a.u.) '
2541 2 : call wrtout(unt_rsp, msg)
2542 82 : do ipvel=1,npvel
2543 80 : write(msg,'(f16.8,4x,f16.8)') pvel_norm(ipvel),rspower(ipvel)
2544 82 : call wrtout(unt_rsp,msg)
2545 : enddo
2546 2 : close(unt_rsp)
2547 : end if
2548 :
2549 12 : ABI_FREE(omega_re)
2550 12 : ABI_FREE(iomega_re)
2551 12 : ABI_FREE(im_epsm1_diag_qbz)
2552 12 : ABI_FREE(tmp_data)
2553 :
2554 12 : end subroutine random_stopping_power
2555 : !!***
2556 :
2557 : !!****f* m_screening_driver/calc_rpa_functional
2558 : !! NAME
2559 : !! calc_rpa_functional
2560 : !!
2561 : !! FUNCTION
2562 : !! Routine used to calculate the Galitskii-Migdal and RPA approximations to the correlation energy
2563 : !! from the irreducible polarizability.
2564 : !!
2565 : !! INPUTS
2566 : !! iq=index of the q-point in the array Qmesh%ibz where epsilon^-1 has to be calculated
2567 : !! Ep<em1params_t>=Structure with parameters and dimensions related to the inverse dielectric matrix.
2568 : !! Pvc<vcoul_t>=Structure gathering data on the Coulombian interaction
2569 : !! Qmesh<kmesh_t>=Data type with information on the q-sampling
2570 : !! Dtfil<Datafiles_type)>=variables related to files
2571 : !! comm=MPI communicator.
2572 : !!
2573 : !! OUTPUT
2574 : !!
2575 : !! SOURCE
2576 :
2577 18 : subroutine calc_rpa_functional(gwrpacorr,gwgmcorr,iqcalc,iq,Ep,Pvc,Qmesh,Dtfil,gmet,chi0,comm,ec_rpa,ec_gm)
2578 :
2579 : use m_hide_lapack, only : xginv, xheev
2580 :
2581 : !Arguments ------------------------------------
2582 : !scalars
2583 : integer,intent(in) :: iqcalc,iq,gwrpacorr,gwgmcorr,comm
2584 : real(dp),intent(inout) :: ec_gm
2585 : type(kmesh_t),intent(in) :: Qmesh
2586 : type(vcoul_t),intent(in) :: Pvc
2587 : type(Datafiles_type),intent(in) :: Dtfil
2588 : type(em1params_t),intent(in) :: Ep
2589 : !arrays
2590 : real(dp),intent(in) :: gmet(3,3)
2591 : real(dp),intent(inout) :: ec_rpa(gwrpacorr)
2592 : complex(gwp),intent(inout) :: chi0(Ep%npwe,Ep%npwe,Ep%nomega)
2593 :
2594 : !Local variables-------------------------------
2595 : !scalars
2596 : integer,parameter :: master = 0
2597 : integer :: ig1,ig2,ilambda,io,rank,nprocs,unt,ierr, units(2)
2598 : real(dp) :: ecorr,ecorr_gm,lambda
2599 : logical :: q_is_gamma
2600 : character(len=500) :: msg
2601 : !arrays
2602 18 : real(dp),allocatable :: z(:),zl(:),zlw(:),zw(:)
2603 18 : complex(gwp),allocatable :: chi0_diag(:),chitmp(:,:),chi0_diag_gm(:),chitmp_gm(:,:)
2604 18 : real(gwp),allocatable :: eig(:)
2605 : ! *************************************************************************
2606 :
2607 54 : units = [std_out, ab_out]
2608 :
2609 : ! initialize MPI data
2610 18 : rank = xmpi_comm_rank(comm)
2611 18 : nprocs = xmpi_comm_size(comm)
2612 :
2613 18 : ABI_CHECK(ep%iw_mesh_type == "gauss_legendre", "only gauss legendre mesh is supported")
2614 :
2615 : !if (rank==master) then ! presently only master has chi0 in screening
2616 :
2617 : ! vc_sqrt contains vc^{1/2}(q,G), complex-valued to allow for a possible cutoff
2618 18 : q_is_gamma = normv(Qmesh%ibz(:,iq),gmet,'G')<GW_TOLQ0
2619 :
2620 : ! Calculate Gauss-Legendre quadrature knots and weights for the omega integration
2621 54 : ABI_MALLOC(zw, (Ep%nomegaei))
2622 36 : ABI_MALLOC(z, (Ep%nomegaei))
2623 18 : call coeffs_gausslegint(zero, one, z, zw, Ep%nomegaei)
2624 :
2625 : ! Calculate Gauss-Legendre quadrature knots and weights for the lambda integration
2626 54 : ABI_MALLOC(zlw, (gwrpacorr))
2627 36 : ABI_MALLOC(zl, (gwrpacorr))
2628 18 : call coeffs_gausslegint(zero,one,zl,zlw,gwrpacorr)
2629 :
2630 54 : ABI_MALLOC(chi0_diag,(Ep%npwe))
2631 72 : ABI_MALLOC_OR_DIE(chitmp,(Ep%npwe,Ep%npwe), ierr)
2632 18 : if(gwgmcorr==1) then
2633 12 : ABI_MALLOC(chi0_diag_gm,(Ep%npwe))
2634 18 : ABI_MALLOC_OR_DIE(chitmp_gm,(Ep%npwe,Ep%npwe), ierr)
2635 : end if
2636 :
2637 810 : do io=2,Ep%nomega
2638 : !if (q_is_gamma) then
2639 : ! call wrtout(units, "RPA: Ignoring q==0"); cycle
2640 : !end if
2641 :
2642 810 : if(gwrpacorr==1) then ! exact integration over the coupling constant
2643 :
2644 774 : if(modulo(io-2,nprocs)/=rank) cycle ! distributing the workload
2645 :
2646 18972 : do ig2=1,Ep%npwe
2647 1194498 : do ig1=1,Ep%npwe
2648 1193724 : chitmp(ig1,ig2) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig2,iq) * chi0(ig1,ig2,io)
2649 : end do !ig1
2650 : end do !ig2
2651 :
2652 2322 : ABI_MALLOC(eig,(Ep%npwe))
2653 774 : call xheev('N','U',Ep%npwe,chitmp,eig)
2654 :
2655 18972 : do ig1=1,Ep%npwe
2656 : ec_rpa(:) = ec_rpa(:) &
2657 : & - zw(io-1) / ( z(io-1) * z(io-1) ) &
2658 36396 : & * Qmesh%wt(iq) * (-log( 1.0_dp-eig(ig1) ) - eig(ig1) ) / (2.0_dp * pi )
2659 : ec_gm = ec_gm &
2660 : & - zw(io-1) / ( z(io-1) * z(io-1) ) &
2661 18972 : & * Qmesh%wt(iq) * ( eig(ig1) / ( 1.0_dp-eig(ig1) ) - eig(ig1) ) / (2.0_dp * pi )
2662 : end do
2663 774 : ABI_FREE(eig)
2664 :
2665 : else ! numerical integration over the coupling constant
2666 :
2667 : !if(modulo( (ilambda-1)+gwrpacorr*(io-2),nprocs)/=rank) cycle ! distributing the workload
2668 :
2669 162 : do ilambda=1,gwrpacorr
2670 144 : if(modulo( (ilambda-1)+gwrpacorr*(io-2),nprocs)/=rank) cycle ! distributing the workload
2671 144 : lambda=zl(ilambda)
2672 19872 : do ig1=1,Ep%npwe
2673 19872 : chi0_diag(ig1) = Pvc%vc_sqrt(ig1,iq)**2 * chi0(ig1,ig1,io)
2674 : end do
2675 :
2676 144 : if(ilambda==1 .and. gwgmcorr==1) then ! Copy v^1/2*Chi0*v^1/2 for Galitskii-Migdal
2677 0 : chi0_diag_gm(:) = chi0_diag(:)
2678 : end if
2679 :
2680 19872 : do ig2=1,Ep%npwe
2681 2722464 : do ig1=1,Ep%npwe
2682 2702736 : chitmp(ig1,ig2) = - lambda * Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chi0(ig1,ig2,io)
2683 :
2684 2722464 : if(ilambda==1 .and. gwgmcorr==1) then ! Use lambda=1 for Galitskii-Migdal
2685 0 : chitmp_gm(ig1,ig2) = - Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chi0(ig1,ig2,io)
2686 : end if
2687 :
2688 : end do !ig1
2689 19728 : chitmp(ig2,ig2) = chitmp(ig2,ig2) + 1.0_dp
2690 :
2691 19872 : if(ilambda==1 .and. gwgmcorr==1) then ! Prepare (1-v^1/2*Chi0*v^1/2) for Galitskii-Migdal
2692 0 : chitmp_gm(ig2,ig2) = chitmp_gm(ig2,ig2) + 1.0_dp
2693 : end if
2694 :
2695 : end do !ig2
2696 144 : call xginv(chitmp(:,:),Ep%npwe)
2697 2722752 : chitmp(:,:) = matmul( chi0(:,:,io) , chitmp(:,:) )
2698 :
2699 144 : if(ilambda==1 .and. gwgmcorr==1) then ! Prepare Chi = [(1-v^1/2*Chi0*v^1/2)]^-1 * Chi0 for Galitskii-Migdal
2700 0 : call xginv(chitmp_gm(:,:),Ep%npwe)
2701 0 : chitmp_gm(:,:) = matmul( chi0(:,:,io) , chitmp_gm(:,:) )
2702 : end if
2703 :
2704 19872 : do ig1=1,Ep%npwe
2705 19728 : chi0_diag(ig1) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chitmp(ig1,ig1) - chi0_diag(ig1)
2706 :
2707 19872 : if(ilambda==1 .and. gwgmcorr==1) then ! Prepare v^1/2*Chi*v^1/2 - v^1/2*Chi0*v^1/2 for Galitskii-Migdal
2708 0 : chi0_diag_gm(ig1) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chitmp_gm(ig1,ig1) - chi0_diag_gm(ig1)
2709 : end if
2710 :
2711 : end do
2712 :
2713 19890 : do ig1=1,Ep%npwe
2714 : ec_rpa(ilambda) = ec_rpa(ilambda) &
2715 19728 : & - zw(io-1) / ( z(io-1) * z(io-1) ) * Qmesh%wt(iq) * real( chi0_diag(ig1) ) / (2.0_dp * pi )
2716 :
2717 19872 : if(ilambda==1 .and. gwgmcorr==1) then ! Integrate [v*Chi-v*Chi0](iw) dw for Galitskii-Migdal
2718 : ec_gm = ec_gm &
2719 0 : & - zw(io-1) / ( z(io-1) * z(io-1) ) * Qmesh%wt(iq) * real( chi0_diag_gm(ig1) ) / (2.0_dp * pi )
2720 : end if
2721 :
2722 : end do
2723 :
2724 : end do ! ilambda
2725 :
2726 : end if ! exact or numerical integration over the coupling constant
2727 :
2728 : end do ! io
2729 :
2730 :
2731 : ! Output the correlation energy when the last q-point to be calculated is reached
2732 : ! This would allow for a manual parallelization over q-points
2733 18 : if(iqcalc==Ep%nqcalc) then
2734 :
2735 5 : call xmpi_sum_master(ec_rpa,master,comm,ierr)
2736 5 : call xmpi_sum_master(ec_gm,master,comm,ierr)
2737 :
2738 5 : if(rank==master) then
2739 17 : ecorr = sum( zlw(:)*ec_rpa(:) )
2740 5 : ecorr_gm = ec_gm
2741 5 : if (open_file(dtfil%fnameabo_rpa, msg, newunit=unt) /=0) then
2742 0 : ABI_ERROR(msg)
2743 : end if
2744 5 : write(unt,'(a,(2x,f14.8))') '#RPA',ecorr
2745 5 : write(msg,'(2a,(2x,f14.8))') ch10,' RPA energy [Ha] :',ecorr
2746 5 : call wrtout(units, msg)
2747 5 : if(gwrpacorr>1) then
2748 9 : do ilambda=1,gwrpacorr
2749 8 : write(unt,'(i6,2x,f10.6,2x,e13.6)') ilambda,zl(ilambda),ec_rpa(ilambda)
2750 8 : write(msg,'(i6,2x,f10.6,2x,e13.6)') ilambda,zl(ilambda),ec_rpa(ilambda)
2751 9 : call wrtout(units, msg)
2752 : end do
2753 : end if
2754 5 : if(gwgmcorr==1) then ! Only exact integration over the coupling constant
2755 1 : write(unt,'(a,(2x,f14.8))') '#GM',ecorr_gm
2756 1 : write(msg,'(2a,(2x,f14.8))') ch10,' Galitskii-Migdal energy [Ha] :',ecorr_gm
2757 1 : call wrtout(units, msg)
2758 1 : write(unt,'(a1)') ' '
2759 1 : write(msg,'(a1)') ' '
2760 1 : call wrtout(units, msg)
2761 : end if
2762 5 : close(unt)
2763 : end if
2764 :
2765 : end if
2766 :
2767 18 : if(gwgmcorr==1) then
2768 6 : ABI_FREE(chitmp_gm)
2769 6 : ABI_FREE(chi0_diag_gm)
2770 : end if
2771 18 : ABI_FREE(chi0_diag)
2772 18 : ABI_FREE(chitmp)
2773 18 : ABI_FREE(zl)
2774 18 : ABI_FREE(zlw)
2775 18 : ABI_FREE(z)
2776 18 : ABI_FREE(zw)
2777 :
2778 18 : end subroutine calc_rpa_functional
2779 : !!***
2780 :
2781 6936 : end module m_screening_driver
2782 : !!***
|
