Branch data 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 : 150 : 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 : 300 : type(ebands_t) :: ks_ebands, qp_ebands
194 [ + + + + : 3900 : type(kmesh_t) :: Kmesh,Qmesh
+ + + + ]
195 [ + + + + : 3300 : type(vcoul_t) :: Vcp
+ + + + +
+ ]
196 [ + + + + : 7800 : type(crystal_t) :: Cryst
+ + + + +
+ + + + +
+ + + + ]
197 : 150 : type(em1params_t) :: Ep
198 : : type(Energies_type) :: KS_energies
199 : 150 : type(gsphere_t) :: Gsph_epsG0,Gsph_wfn
200 : 150 : type(Hdr_type) :: Hdr_wfk,Hdr_local
201 : 150 : type(MPI_type) :: MPI_enreg_seq
202 : : type(Pawfgr_type) :: Pawfgr
203 [ + + + + ]: 1050 : type(hscr_t) :: Hem1,Hchi0
204 : 150 : type(wfdgw_t) :: Wfd,Wfdf
205 : 150 : type(spectra_t) :: spectra
206 : : !type(chi_t) :: chihw
207 : 150 : type(wvl_data) :: wvl_dummy
208 : : character(len=nctk_slen) :: wing_shape
209 : : !arrays
210 : 300 : integer :: ibocc(Dtset%nsppol),ngfft_gw(18),ngfftc(18),ngfftf(18), units(2)
211 : 300 : integer,allocatable :: irottb(:,:),ktabr(:,:),ktabrf(:,:),l_size_atm(:)
212 : 300 : integer,allocatable :: ks_vbik(:,:),ks_occ_idx(:,:),qp_vbik(:,:),nband(:,:)
213 : 300 : 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 : 600 : real(dp),allocatable :: igwene(:,:,:),chi0_sumrule(:),ec_rpa(:),rspower(:)
216 : 150 : real(dp),allocatable :: nhat(:,:),nhatgr(:,:,:),ph1d(:,:),ph1df(:,:)
217 : 150 : real(dp),allocatable :: rhog(:,:),rhor(:,:),rhor_p(:,:),rhor_kernel(:,:),taur(:,:)
218 : 450 : real(dp),allocatable :: z(:),zw(:),grchempottn(:,:),grewtn(:,:),grvdw(:,:),kxc(:,:),qmax(:)
219 : 300 : real(dp),allocatable :: ks_vhartr(:),vpsp(:),ks_vtrial(:,:),ks_vxc(:,:),xccc3d(:)
220 : 450 : complex(gwp),allocatable :: arr_99(:,:),kxcg(:,:),fxc_ADA(:,:,:)
221 : 150 : complex(dp),allocatable :: m_ks_to_qp(:,:,:,:)
222 : 150 : complex(dp),allocatable :: chi0_head(:,:,:), chi0_lwing(:,:,:), chi0_uwing(:,:,:)
223 : 150 : real(dp),allocatable :: rwork_wing(:,:,:,:), tmp_omega_wgs(:)
224 : 300 : complex(dp),allocatable :: chi0intra_lwing(:,:,:),chi0intra_uwing(:,:,:),chi0intra_head(:,:,:), tmp_omega(:)
225 : 150 : complex(gwp),allocatable,target :: chi0(:,:,:),chi0intra(:,:,:)
226 : 150 : complex(gwp),contiguous, pointer :: epsm1(:,:,:)
227 : 150 : logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:)
228 : : character(len=80) :: title(2)
229 : : character(len=fnlen) :: gw_fname,wfk_fname,lwl_fname
230 : 150 : type(littlegroup_t),pointer :: Ltg_q(:)
231 : 150 : type(Paw_an_type),allocatable :: Paw_an(:)
232 : 150 : type(Paw_ij_type),allocatable :: Paw_ij(:)
233 : 150 : type(Pawfgrtab_type),allocatable :: Pawfgrtab(:)
234 : 150 : type(Pawrhoij_type),allocatable :: Pawrhoij(:),prev_Pawrhoij(:)
235 : 150 : type(pawpwff_t),allocatable :: Paw_pwff(:)
236 : 150 : type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
237 : 150 : type(plowannier_type) :: wanbz,wanibz,wanibz_in
238 : : real(dp) :: te_min, te_max
239 : 150 : type(gaps_t) :: gaps
240 : 150 : real(dp),allocatable :: tau_mesh(:), tau_wgs(:), iw_mesh(:), iw_wgs(:)
241 : 150 : 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 : 150 : call timab(301,1,tsec) ! overall time
246 : 150 : call timab(302,1,tsec) ! screening(init
247 [ + + ]: 450 : units = [std_out, ab_out]
248 : :
249 : : write(msg,'(6a)')&
250 : 150 : ' SCREENING: Calculation of the susceptibility and dielectric matrices ',ch10,ch10,&
251 : 150 : ' Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining.',ch10,&
252 : 300 : ' Incorporated in ABINIT by V. Olevano, G.-M. Rignanese, and M. Torrent.'
253 : 150 : call wrtout(units, msg)
254 : :
255 [ - + ]: 150 : 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 : 150 : 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 : 150 : 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 : 150 : comm = xmpi_world; my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
278 : :
279 [ + + ]: 150 : if (my_rank == master) then
280 : 122 : wfk_fname = dtfil%fnamewffk
281 [ - + ]: 122 : if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then
282 : 0 : ABI_ERROR(msg)
283 : : end if
284 : : end if
285 : 150 : call xmpi_bcast(wfk_fname, master, comm, ierr)
286 : :
287 : : ! Some variables need to be initialized/nullify at start
288 : 150 : call KS_energies%init()
289 : 150 : usexcnhat=0
290 : :
291 : 150 : call mkrdim(acell,rprim,rprimd)
292 : 150 : 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 : 150 : k0(:)=zero
301 : : call pawfgr_init(Pawfgr,Dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,&
302 : 150 : gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=gmet,k0=k0)
303 : :
304 [ + + ]: 300 : call print_ngfft([std_out], ngfftf,'Dense FFT mesh used for densities and potentials')
305 [ + + ]: 600 : nfftf_tot=PRODUCT(ngfftf(1:3))
306 : :
307 : : ! We can intialize MPI_enreg and fft distrib here, now ngfft are known
308 : 150 : call initmpi_seq(MPI_enreg_seq) ! Fake MPI_type for the sequential part.
309 : 150 : call MPI_enreg_seq%distribfft%init_seq('c',ngfftc(2),ngfftc(3),'all')
310 : 150 : call MPI_enreg_seq%distribfft%init_seq('f',ngfftf(2),ngfftf(3),'all')
311 : :
312 : : !=============================================
313 : : !==== Open and read pseudopotential files ====
314 : : !=============================================
315 : 150 : 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 : 150 : ngfft_gw,Hdr_wfk,Hdr_local,Cryst,Kmesh,Qmesh,ks_ebands,Ltg_q,Gsph_epsG0,Gsph_wfn,Vcp,Ep,comm)
320 : :
321 : 150 : call pstat_proc%print(_PSTAT_ARGS_)
322 : 150 : call timab(302,2,tsec) ! screening(init)
323 [ + + ]: 300 : call print_ngfft([std_out], ngfft_gw, header='FFT mesh used for oscillator strengths')
324 : :
325 [ + + ]: 600 : nfftgw_tot=PRODUCT(ngfft_gw(1:3))
326 : : mgfftgw =MAXVAL (ngfft_gw(1:3))
327 : 150 : nfftgw =nfftgw_tot ! no FFT //
328 : :
329 : : !TRYING TO RECREATE AN "ABINIT ENVIRONMENT"
330 : 150 : KS_energies%e_corepsp=ecore/Cryst%ucvol
331 : :
332 : : !==========================
333 : : !=== PAW initialization ===
334 : : !==========================
335 [ + + ]: 150 : if (Dtset%usepaw==1) then
336 : 6 : call timab(315,1,tsec) ! screening(pawin
337 : :
338 : 6 : call chkpawovlp(Cryst%natom,Cryst%ntypat,Dtset%pawovlp,Pawtab,Cryst%rmet,Cryst%typat,Cryst%xred)
339 : :
340 [ - + + - : 51 : ABI_MALLOC(Pawrhoij,(Cryst%natom))
+ - + - -
+ - + - +
+ + + - -
+ - + - +
- + - + -
+ ]
341 : : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
342 : 6 : & nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc)
343 : : call pawrhoij_alloc(Pawrhoij,cplex_rhoij,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,&
344 : 6 : & Cryst%typat,pawtab=Pawtab)
345 : :
346 : : ! Initialize values for several basic arrays stored in Pawinit
347 [ + - + - : 6 : 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 : 6 : call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
351 : :
352 [ + + - + ]: 6 : 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 [ - + - - ]: 5 : if (Pawtab(1)%has_kij ==1) Pawtab(1:Cryst%ntypat)%has_kij =2
362 [ - + - - ]: 5 : if (Pawtab(1)%has_nabla==1) Pawtab(1:Cryst%ntypat)%has_nabla=2
363 : : end if
364 [ + + ]: 17 : 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 [ + + ]: 17 : Pawtab(:)%has_nabla = 0
369 [ + + ]: 17 : Pawtab(:)%usepawu = 0
370 [ + + ]: 17 : Pawtab(:)%useexexch = 0
371 [ + + ]: 17 : Pawtab(:)%exchmix =zero
372 [ + + ]: 17 : 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 : 6 : call pawnabla_init(Psps%mpsang,Cryst%ntypat,Pawrad,Pawtab)
377 : :
378 : 6 : 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 : 6 : & Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa)
385 : :
386 [ + - ]: 6 : if (my_rank == master) call pawtab_print(Pawtab)
387 : :
388 : : ! Get Pawrhoij from the header of the WFK file.
389 : 6 : 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 : 6 : 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 : 6 : rhoxsp_method=2
402 : :
403 : : ! At least for ucrpa, the Arnaud Alouani is always a better choice but needs a larger cutoff
404 [ - + ]: 6 : if(dtset%ucrpa>0) rhoxsp_method=1
405 [ + + ]: 6 : if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc
406 : :
407 [ - + + - : 18 : ABI_MALLOC(gw_gfft,(3,nfftgw_tot))
+ - + - -
+ - + ]
408 : 6 : call get_gfft(ngfft_gw,(/zero,zero,zero/),gmet,gw_gsq,gw_gfft)
409 [ - + ]: 6 : ABI_FREE(gw_gfft)
410 : :
411 : : ! Set up q grids, make qmax 20% larger than largest expected:
412 [ - + + - : 18 : ABI_MALLOC(nq_spl,(Psps%ntypat))
+ - + - -
+ - + ]
413 [ - + + - : 18 : ABI_MALLOC(qmax,(Psps%ntypat))
+ - + - -
+ - + ]
414 [ + + ]: 17 : nq_spl = Psps%mqgrid_ff
415 [ + + ]: 17 : qmax = SQRT(gw_gsq)*1.2d0 !qmax = Psps%qgrid_ff(Psps%mqgrid_ff)
416 [ - + + - : 101 : ABI_MALLOC(Paw_pwff,(Psps%ntypat))
+ - + - -
+ - + - +
+ + + + +
+ + - - +
- + ]
417 : :
418 : 6 : call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps)
419 [ - + ]: 6 : ABI_FREE(nq_spl)
420 [ - + ]: 6 : ABI_FREE(qmax)
421 : :
422 : : ! Variables/arrays related to the fine FFT grid
423 [ - + + - : 24 : ABI_MALLOC(nhat,(nfftf,Dtset%nspden))
- + + - +
- + - - +
- + ]
424 [ + + + + ]: 159484 : nhat=zero; cplex=1
425 [ - + + - : 51 : ABI_MALLOC(Pawfgrtab,(Cryst%natom))
+ - + - -
+ - + - +
+ + + - -
+ - + - +
- + - + -
+ - + -
+ ]
426 : 6 : call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat)
427 : 6 : call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Dtset%nspden,Dtset%typat)
428 [ - + ]: 6 : ABI_FREE(l_size_atm)
429 : 6 : compch_fft=greatest_real
430 [ + + ]: 17 : 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 : 6 : write(msg,'(a,i3)')' screening : using usexcnhat = ',usexcnhat
434 : 6 : call wrtout(std_out, msg)
435 : :
436 : : ! Identify parts of the rectangular grid where the density has to be calculated.
437 : 6 : optcut=0;optgr0=Dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-Dtset%pawstgylm
438 [ - + ]: 6 : if (Dtset%pawcross==1) optrad=1
439 [ + + + - ]: 6 : 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 : 6 : optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
443 : :
444 : 18 : call timab(315,2,tsec) ! screening(pawin
445 : : else
446 : : ! allocate empty structure for the sake of -fcheck=all...
447 [ - + - + ]: 144 : ABI_MALLOC(Paw_pwff,(0))
448 [ - + - + ]: 144 : ABI_MALLOC(Pawrhoij,(0))
449 [ - + - + ]: 144 : 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 [ - + + - : 805 : ABI_MALLOC(Paw_onsite,(Cryst%natom))
+ - + - -
+ - + - +
+ + + - -
+ - + - +
- + - + ]
454 [ + + ]: 150 : if (Dtset%usepaw==1) then
455 [ - + ]: 6 : 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 [ - + ]: 6 : ABI_CHECK(Dtset%useexexch==0,"LEXX not yet implemented in GW")
464 [ - + ]: 6 : ABI_CHECK(Dtset%usedmft==0,"DMFT + GW not available")
465 : :
466 [ - + ]: 6 : 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 [ + + ]: 150 : if (.not.allocated(nhat)) then
475 [ - + + - : 288 : ABI_MALLOC(nhat,(nfftf,0))
- + - + -
+ ]
476 : : end if
477 : :
478 : 150 : 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 [ - + + - : 600 : ABI_MALLOC(ks_occ_idx,(ks_ebands%nkpt, ks_ebands%nsppol))
- + + - +
- + - - +
- + - + ]
498 [ - + + - : 450 : ABI_MALLOC(ks_vbik ,(ks_ebands%nkpt, ks_ebands%nsppol))
- + + - -
+ - + -
+ ]
499 [ - + + - : 450 : ABI_MALLOC(qp_vbik ,(ks_ebands%nkpt, ks_ebands%nsppol))
- + + - -
+ - + ]
500 : :
501 : 150 : call ks_ebands%update_occ(Dtset%spinmagntarget, prtvol=0)
502 [ + - + - : 150 : ks_occ_idx = ks_ebands%get_occupied(tol8) ! tol8 to be consistent when the density
- + - - ]
503 [ + - + - : 150 : ks_vbik = ks_ebands%get_valence_idx()
- + - - ]
504 : :
505 : 150 : ibocc(:)=MAXVAL(ks_occ_idx(:,:),DIM=1) ! Max occupied band index for each spin.
506 [ - + ]: 150 : ABI_FREE(ks_occ_idx)
507 : :
508 : 150 : use_tr =.FALSE.; nbvw=0
509 [ + + + + : 150 : if (Dtset%gwpara==2.or.Ep%awtr==1.or.Dtset%spmeth>0) then
- + ]
510 : 139 : use_tr = .TRUE.
511 [ + + ]: 280 : nbvw=MAXVAL(ibocc)
512 : 139 : nbcw=Ep%nbnds-nbvw
513 : 139 : write(msg,'(4a,i0,2a,i0,2a,i0,a)')ch10,&
514 : 139 : '- screening: taking advantage of time-reversal symmetry ',ch10,&
515 : 139 : '- Maximum band index for partially occupied states nbvw = ',nbvw,ch10,&
516 : 139 : '- Remaining bands to be divided among processors nbcw = ',nbcw,ch10,&
517 : 278 : '- Number of bands treated by each node ~',nbcw/nprocs,ch10
518 : 139 : call wrtout(ab_out, msg)
519 [ - + ]: 139 : 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 : 150 : mband=Ep%nbnds
525 [ - + + - : 600 : ABI_MALLOC(nband,(Kmesh%nibz,Dtset%nsppol))
- + + - +
- + - - +
- + ]
526 [ + + + + ]: 1168 : nband=mband
527 [ - + + - : 750 : ABI_MALLOC(bks_mask,(mband,Kmesh%nibz,Dtset%nsppol))
- + + - -
+ + - + -
+ - - + -
+ ]
528 [ - + + - : 600 : ABI_MALLOC(keep_ur,(mband,Kmesh%nibz,Dtset%nsppol))
- + + - -
+ + - - +
- + ]
529 [ + + + + : 46886 : bks_mask=.FALSE.; keep_ur=.FALSE.
+ + + + +
+ + + ]
530 : :
531 : : ! autoparal section
532 [ - + ]: 150 : 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 : 150 : 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 : 150 : Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm)
612 : :
613 [ - + ]: 150 : 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 [ - + ]: 150 : ABI_FREE(bks_mask)
620 [ - + ]: 150 : ABI_FREE(nband)
621 [ - + ]: 150 : ABI_FREE(keep_ur)
622 : :
623 [ + + ]: 300 : 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 : 150 : call wfd%read_wfk(wfk_fname,iomode_from_fname(wfk_fname))
632 : :
633 [ - + ]: 150 : 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 : 150 : call pstat_proc%print(_PSTAT_ARGS_)
642 : 150 : call timab(316,2,tsec) ! screening(wfs
643 : 150 : call timab(319,1,tsec) ! screening(1)
644 : :
645 [ - + - - ]: 150 : 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 [ - + + - : 600 : ABI_MALLOC(irottb, (nfftgw, Cryst%nsym))
- + + - +
- + - - +
- + ]
653 : 150 : call rotate_FFT_mesh(Cryst%nsym, Cryst%symrel, Cryst%tnons, ngfft_gw, irottb, iscompatibleFFT)
654 : :
655 [ - + + - : 600 : ABI_MALLOC(ktabr,(nfftgw,Kmesh%nbz))
- + + - +
- + - - +
- + ]
656 [ + + ]: 7126 : do ikbz=1,Kmesh%nbz
657 : 6976 : isym=Kmesh%tabo(ikbz)
658 [ + + ]: 40451632 : do ifft=1,nfftgw
659 : 40451482 : ktabr(ifft,ikbz)=irottb(ifft,isym)
660 : : end do
661 : : end do
662 [ - + ]: 150 : ABI_FREE(irottb)
663 : :
664 [ + + - + ]: 150 : 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 [ - + ]: 150 : ABI_MALLOC(ktabrf,(0,0))
678 : : end if
679 : :
680 : : ! Compute structure factor phases and large sphere cut-off
681 [ + - - + : 450 : ABI_MALLOC(ph1d,(2,3*(2*Dtset%mgfft+1)*Cryst%natom))
+ - + - +
- - + - +
- + ]
682 [ + - - + : 450 : ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*Cryst%natom))
+ - + - +
- - + - +
- + ]
683 : 150 : call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred)
684 : :
685 [ + + + + ]: 150 : if (Psps%usepaw==1.and.Pawfgr%usefinegrid==1) then
686 : 3 : call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred)
687 : : else
688 [ + + + + ]: 140574 : 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 : 150 : call ks_ebands%copy(qp_ebands)
694 : 150 : call timab(319,2,tsec) ! screening(1)
695 : :
696 : : !============================
697 : : !==== Self-consistent GW ====
698 : : !============================
699 [ + + ]: 150 : 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 : 150 : 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 : 150 : gw_fname = "__in.gw__"
747 : 150 : update_energies = file_exists(gw_fname)
748 : :
749 [ + + ]: 150 : 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 [ - + ]: 149 : 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 [ - + - - : 150 : el_temp = merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
- - ]
777 : :
778 [ - + + - : 600 : ABI_MALLOC(rhor,(nfftf,Dtset%nspden))
- + + - +
- + - - +
- + ]
779 [ - + + - : 600 : ABI_MALLOC(taur,(nfftf,Dtset%nspden*Dtset%usekden))
+ - - - +
- - + - +
- + ]
780 : :
781 : 150 : call wfd%mkrho(cryst, psps, qp_ebands, ngfftf, nfftf, rhor)
782 [ - + ]: 150 : if (Dtset%usekden==1) call wfd%mkrho(cryst, psps, qp_ebands, ngfftf, nfftf, taur, optcalc=1)
783 : :
784 : 150 : call timab(305,2,tsec) ! screening(densit
785 : :
786 : 150 : nhatgrdim = 0
787 [ + + ]: 150 : if (Dtset%usepaw==1) then ! Additional computation for PAW.
788 : 6 : call timab(320,1,tsec) ! screening(paw
789 : :
790 : : ! Add the compensation charge to the PW density.
791 [ + + ]: 6 : nhatgrdim=0; if (Dtset%xclevel==2) nhatgrdim=usexcnhat*Dtset%pawnhatxc
792 : 6 : cplex=1; ider=2*nhatgrdim; izero=0
793 [ + + ]: 6 : if (nhatgrdim>0) then
794 [ - + + - : 5 : ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,3))
- + + - +
- + - + -
- + - + -
+ ]
795 : : else
796 [ - + + - : 15 : 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 : 6 : 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 : 6 : cplex=1
806 [ - + + - : 51 : ABI_MALLOC(Paw_ij,(Cryst%natom))
+ - + - -
+ - + - +
+ + + - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ ]
807 : 6 : 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 : 6 : has_dij=1,has_dijhartree=1,has_exexch_pot=1,has_pawu_occ=1)
811 : :
812 : 6 : nkxc1=0
813 [ - + + - : 51 : ABI_MALLOC(Paw_an,(Cryst%natom))
+ - + - -
+ - + - +
+ + + - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + ]
814 : 6 : call paw_an_nullify(Paw_an)
815 : : call paw_an_init(Paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,&
816 : 6 : cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=0)
817 : :
818 : 6 : 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 : 6 : & Dtset%xclevel,Dtset%xc_denpos,Dtset%xc_taupos,Cryst%xred,Cryst%ucvol,Psps%znuclpsp,Dtset%spinaxis)
823 : 6 : call timab(320,2,tsec) ! screening(paw
824 : : else
825 [ - + - + ]: 144 : ABI_MALLOC(Paw_ij,(0))
826 [ - + - + ]: 144 : ABI_MALLOC(Paw_an,(0))
827 : : end if ! usepaw
828 : :
829 : 150 : call timab(321,1,tsec) ! screening(2)
830 : :
831 : : !JB : Should be remove : cf. l 839
832 [ + + ]: 150 : if (.not.allocated(nhatgr)) then
833 [ - + + - : 432 : ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,0))
- + + - -
+ - + -
+ ]
834 : : end if
835 : :
836 : : call test_charge(nfftf,ks_ebands%nelect,Dtset%nspden,rhor,ucvol,&
837 : 150 : 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 [ + + + + : 159628 : if (Dtset%usepaw==1) rhor(:,:)=rhor(:,:)+nhat(:,:)
+ + ]
841 : :
842 : 150 : call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,rhor,ucvol=ucvol)
843 [ - + ]: 150 : 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 [ + + ]: 150 : if (dtset%gwgamma>0 .or. dtset%gwgamma==-11) then
848 [ - + + - : 8 : ABI_MALLOC(rhor_kernel,(nfftf,Dtset%nspden))
- + + - +
- + - - +
- + ]
849 : : end if
850 : :
851 [ - + + - : 450 : ABI_MALLOC(rhog,(2,nfftf))
+ - + - -
+ - + ]
852 : 150 : 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 : 150 : ngrvdw=0
865 [ - + ]: 150 : ABI_MALLOC(grvdw,(3,ngrvdw))
866 [ - + + - : 450 : ABI_MALLOC(grchempottn,(3,Cryst%natom))
+ - + - -
+ - + ]
867 [ - + + - : 300 : ABI_MALLOC(grewtn,(3,Cryst%natom))
- + - + ]
868 : 150 : nkxc=0
869 [ + + ]: 150 : if (Dtset%nspden==1) nkxc=2
870 [ + + ]: 150 : 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 [ + + - + ]: 150 : if (Dtset%ixc<0 .and. libxc_functionals_ismgga()) nkxc=0
873 [ + - ]: 150 : if (nkxc/=0) then
874 [ - + + - : 600 : ABI_MALLOC(kxc,(nfftf,nkxc))
- + + - +
- + - - +
- + ]
875 : : end if
876 : :
877 [ + + ]: 150 : n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf
878 [ + + + - : 450 : ABI_MALLOC(xccc3d,(n3xccc))
+ - + + -
+ - + ]
879 [ - + + - : 450 : ABI_MALLOC(ks_vhartr,(nfftf))
+ - + - -
+ - + ]
880 [ - + + - : 600 : ABI_MALLOC(ks_vtrial,(nfftf,Dtset%nspden))
- + + - +
- + - - +
- + ]
881 [ - + + - : 300 : ABI_MALLOC(vpsp,(nfftf))
- + - + ]
882 [ - + + - : 450 : ABI_MALLOC(ks_vxc,(nfftf,Dtset%nspden))
- + + - -
+ - + -
+ ]
883 : :
884 : 150 : 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 : 150 : vxcavg,wvl_dummy,xccc3d,Cryst%xred,taur=taur)
891 : :
892 [ + - ]: 150 : if (nkxc/=0) then
893 [ - + ]: 150 : ABI_FREE(kxc)
894 : : end if
895 [ - + ]: 150 : ABI_FREE(grchempottn)
896 [ - + ]: 150 : ABI_FREE(grewtn)
897 [ - + ]: 150 : ABI_FREE(grvdw)
898 [ - + ]: 150 : ABI_FREE(xccc3d)
899 : :
900 : : !============================
901 : : !==== Compute KS PAW Dij ====
902 : : !============================
903 [ + + ]: 150 : if (Dtset%usepaw==1) then
904 : 6 : call timab(561,1,tsec)
905 : :
906 : : ! Calculate unsymmetrized Dij.
907 : 6 : 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 : 6 : & 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 : 6 : & Cryst%symrec)
924 : : #endif
925 : : ! Output of the pseudopotential strengths Dij and the augmentation occupancies Rhoij.
926 : 6 : call pawprt(Dtset,Cryst%natom,Paw_ij,Pawrhoij,Pawtab)
927 : 6 : 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 [ - + + - : 450 : ABI_MALLOC(Ep%omega, (Ep%nomega))
+ - + - -
+ - + -
+ ]
937 : 150 : Ep%omega(1) = CMPLX(Ep%omegaermin, zero, kind=dp)
938 : :
939 : 150 : ep%iw_mesh_type = "None"; ep%rw_mesh_type = "None"; ep%cw_mesh_type = "None"
940 [ - + + - : 1409 : ABI_CALLOC(Ep%omega_wgs, (Ep%nomega))
+ - + - -
+ - + - +
+ + ]
941 : :
942 : : ! Deal with real frequencies first.
943 [ + + ]: 150 : 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 : : ! Ep%omegaermax (freqremax) has been set in setup_screening
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 [ + + + + ]: 150 : if (Ep%plasmon_pole_model .and. Ep%nomega == 2) then
988 [ + + ]: 106 : e0= Dtset%ppmfrq; if (e0 < 0.1d-4) e0 = omegaplasma
989 : 106 : Ep%omega(2)=CMPLX(zero,e0, kind=dp)
990 : : end if
991 : :
992 [ + + ]: 150 : 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 [ + + + + ]: 139 : 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 [ + + + + ]: 138 : 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 [ + + + + ]: 137 : 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 [ + + ]: 150 : 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 [ + + + + : 442 : 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 : 150 : write(msg,'(2a)')ch10,' calculating chi0 at frequencies [eV] :'
1119 : 150 : call wrtout(units, msg)
1120 [ + + ]: 1068 : do iomega=1,Ep%nomega
1121 : 918 : write(msg,'(i3,2es16.6)')iomega,Ep%omega(iomega)*Ha_eV
1122 : 1068 : call wrtout(units, msg)
1123 : : end do
1124 : :
1125 : : ! Allocate chi0, wings and array for chi0_sumrule check.
1126 [ - + + - : 450 : ABI_MALLOC(chi0_sumrule, (Ep%npwe))
+ - + - -
+ - + ]
1127 : :
1128 : 150 : 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 : 150 : call wrtout(std_out, msg)
1130 [ - + + - : 750 : 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 : 150 : call pstat_proc%print(_PSTAT_ARGS_)
1138 : 150 : call timab(321,2,tsec) ! screening(2)
1139 : :
1140 : 150 : iqcalc = 0
1141 [ - + ]: 150 : 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 [ + + ]: 1093 : do iqibz=1,Qmesh%nibz
1152 : 943 : call timab(306,1,tsec)
1153 : 943 : is_first_qcalc=(iqibz==1)
1154 : :
1155 : : ! Selective q-point calculation.
1156 : 943 : found=.FALSE.; label=iqibz
1157 [ + + ]: 943 : 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 : 931 : iqcalc = iqcalc + 1
1170 : :
1171 : 931 : bar=REPEAT('-',80)
1172 : 931 : write(msg,'(4a,1x,a,i2,a,f9.6,2(",",f9.6),3a)')ch10,ch10,bar,ch10,&
1173 [ + - + + ]: 4655 : ' q-point number ',label,' q = (',(Qmesh%ibz(ii,iqibz),ii=1,3),') [r.l.u.]',ch10,bar
1174 : 931 : call wrtout(units, msg)
1175 [ + + ]: 931 : is_qeq0 = 0; if (normv(Qmesh%ibz(:,iqibz),gmet,'G')<GW_TOLQ0) is_qeq0=1
1176 : :
1177 : 931 : call pstat_proc%print(_PSTAT_ARGS_)
1178 : 931 : call timab(306,2,tsec)
1179 : :
1180 [ + + ]: 931 : if (is_qeq0 == 1) then
1181 : : ! Special treatment of the long wavelength limit.
1182 : 148 : call timab(307,1,tsec)
1183 : :
1184 [ - + + - : 444 : ABI_MALLOC(chi0_head, (3,3,Ep%nomega))
+ - + - -
+ - + -
+ ]
1185 [ - + + - : 740 : ABI_MALLOC(chi0_lwing, (Ep%npwe*Ep%nI, Ep%nomega,3))
- + + - +
- + - + -
- + - + -
+ ]
1186 [ - + + - : 740 : 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 : 148 : 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 : 148 : add_chi0_intraband=.FALSE. !add_chi0_intraband=.TRUE.
1199 : 148 : 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 : 148 : call timab(307,2,tsec)
1233 : :
1234 : : else
1235 : : ! Calculate cchi0 for q/=0.
1236 : 783 : 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 : 783 : Ltg_q(iqibz),chi0_sumrule,Wfd,Wfdf,wanbz)
1240 : 783 : 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 : 931 : call timab(309,1,tsec)
1246 : :
1247 [ + + ]: 5111 : do iomega=1,MIN(Ep%nomega, NOMEGA_PRINTED)
1248 : 4180 : write(msg,'(1x,a,i4,a,2f9.4,a)')' chi0(G,G'') at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]'
1249 : 4180 : call wrtout(units, msg)
1250 : 4180 : write(msg,'(1x,a,i3,a,i4,a)')' chi0(q =',iqibz, ', omega =',iomega,', G,G'')'
1251 [ + + ]: 4180 : if (Ep%nqcalc /= Ep%nqibz) write(msg,'(a,i3,a,i4,a)')' chi0(q=',iqcalc,', omega=',iomega,', G,G'')'
1252 : 4180 : 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 : 4180 : ii = MIN(9, Ep%npwe)
1256 [ - + + - : 16720 : ABI_MALLOC(arr_99,(ii, ii))
- + + - +
- + - - +
- + - + ]
1257 [ - + + + : 384560 : arr_99 = chi0(1:ii,1:ii,iomega) * ucvol
+ + ]
1258 [ + + ]: 8360 : call print_arr([ab_out], arr_99, max_r=2)
1259 [ + + ]: 8360 : call print_arr([std_out], arr_99)
1260 [ - + ]: 5111 : ABI_FREE(arr_99)
1261 : : end do
1262 : :
1263 [ + + ]: 931 : 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 [ + + + - ]: 931 : 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 [ + + + + ]: 931 : 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 : 931 : approx_type=0 ! RPA
1353 : 931 : option_test=0 ! TESTPARTICLE
1354 [ + + ]: 931 : dim_wing=0; if (is_qeq0==1) dim_wing=3
1355 : :
1356 [ + + ]: 931 : if (dim_wing==0) then
1357 : 783 : dim_wing=1
1358 [ + - ]: 783 : if (.not.allocated(chi0_lwing)) then
1359 [ - + + - : 3132 : ABI_MALLOC(chi0_lwing,(Ep%npwe*Ep%nI,Ep%nomega,dim_wing))
- + + - +
- + - - +
- + - + ]
1360 : : end if
1361 [ + - ]: 783 : if (.not.allocated(chi0_uwing)) then
1362 [ - + + - : 3132 : ABI_MALLOC(chi0_uwing,(Ep%npwe*Ep%nJ,Ep%nomega,dim_wing))
- + + - +
- + - - +
- + - + ]
1363 : : end if
1364 [ + - ]: 783 : if (.not.allocated(chi0_head )) then
1365 [ - + + - : 2349 : ABI_MALLOC(chi0_head,(dim_wing,dim_wing,Ep%nomega))
+ - + - -
+ - + -
+ ]
1366 : : end if
1367 : 783 : 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 [ + + + - : 21193 : if (dtset%gwgamma>0) rhor_kernel = rhor
+ - - + -
- + + +
+ ]
1381 : :
1382 : 909 : select case (dtset%gwgamma)
1383 : : case (0)
1384 : 909 : approx_type=0; option_test=0; dim_kxcg=0
1385 [ - + + - : 1818 : 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 [ + + - + : 931 : ABI_ERROR(sjoin("Wrong gwgamma:", itoa(dtset%gwgamma)))
+ - ]
1456 : : end select
1457 : :
1458 [ + + ]: 931 : 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 : 915 : 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 [ - + ]: 931 : ABI_FREE(chi0_lwing)
1489 [ - + ]: 931 : ABI_FREE(chi0_uwing)
1490 [ - + ]: 931 : ABI_FREE(chi0_head)
1491 : :
1492 [ + + + + ]: 931 : if (my_rank == master .and. is_qeq0==1) then
1493 : 120 : call spectra%repr(msg)
1494 : 120 : call wrtout(units, msg)
1495 [ + + ]: 120 : 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 : 931 : call spectra%free()
1505 [ + - ]: 931 : ABI_SFREE(kxcg)
1506 [ - + ]: 931 : 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 : 931 : epsm1 => chi0
1511 : 931 : 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 [ + + + + ]: 931 : 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 [ + + ]: 931 : if (is_qeq0 == 1) then
1528 : 148 : write(msg,'(1x,2a)')' Heads and wings of the symmetrical epsilon^-1(G,G'') ',ch10
1529 : 148 : call wrtout(ab_out,msg)
1530 [ + + ]: 1062 : do iomega=1,Ep%nomega
1531 : 914 : 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 : 914 : call wrtout(ab_out, msg)
1533 [ + + ]: 1828 : call print_arr([ab_out], epsm1(1,:,iomega),max_r=9)
1534 [ + + ]: 1828 : call print_arr([ab_out], epsm1(:,1,iomega),max_r=9)
1535 : 1062 : call wrtout(ab_out, ch10)
1536 : : end do
1537 : : end if
1538 : :
1539 : 931 : call pstat_proc%print(_PSTAT_ARGS_)
1540 : 931 : call timab(309,2,tsec)
1541 : 931 : call timab(310,1,tsec) ! wrscr
1542 : :
1543 [ + + ]: 931 : if (my_rank==master) then
1544 : : ! === Write the symmetrical epsilon^-1 on file ===
1545 : 727 : title(1)='SCR file: epsilon^-1'
1546 [ + + ]: 727 : if (is_qeq0==1) then
1547 [ + + + + ]: 120 : string='0'; if (Dtset%usepaw==0.and.Ep%inclvkb/=0) call int2char10(Ep%inclvkb,string)
1548 : 120 : title(1)=title(1)(1:21)//', calculated using inclvkb = '//string
1549 : : end if
1550 : 727 : title(2)='TESTPARTICLE'
1551 : 727 : ctype='RPA'
1552 : 727 : title(2)(14:17)=ctype !this has to be modified
1553 : :
1554 [ + + ]: 727 : 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 : 122 : ikxc=0; test_type=0; tordering=1
1558 : : hem1 = hscr_new("inverse_dielectric_function",dtset,ep,hdr_local,ikxc,test_type,tordering,title,&
1559 [ - + - + : 122 : & Ep%npwe,Gsph_epsG0%gvec)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
- - - - -
- - ]
1560 : 122 : fform_em1 = hem1%fform
1561 [ + - ]: 122 : if (dtset%iomode == IO_MODE_ETSF) then
1562 [ - + ]: 122 : NCF_CHECK(nctk_open_create(unt_em1, nctk_ncify(dtfil%fnameabo_scr), xmpi_comm_self))
1563 [ - + ]: 122 : NCF_CHECK(cryst%ncwrite(unt_em1))
1564 [ - + ]: 122 : 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 : 122 : call hem1%io(fform_em1,2,unt_em1,xmpi_comm_self,0,Dtset%iomode)
1572 : 122 : call Hem1%free()
1573 : : end if
1574 : :
1575 : 727 : call write_screening("inverse_dielectric_function",unt_em1,Dtset%iomode,Ep%npwe,Ep%nomega,iqcalc,epsm1)
1576 : : end if ! my_rank==master
1577 : :
1578 : 1093 : call timab(310,2,tsec)
1579 : : end do ! Loop over q-points
1580 : :
1581 [ - + ]: 150 : if (Dtset%plowan_compute >= 10) call destroy_plowannier(wanbz)
1582 : :
1583 : : ! Close Files.
1584 [ + + ]: 150 : if (my_rank == master) then
1585 [ + - ]: 122 : if (dtset%iomode == IO_MODE_ETSF) then
1586 [ - + ]: 122 : NCF_CHECK(nf90_close(unt_em1))
1587 [ + + ]: 122 : 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 [ - + ]: 150 : ABI_FREE(chi0_sumrule)
1600 [ - + ]: 150 : ABI_FREE(chi0)
1601 [ + + ]: 150 : ABI_SFREE(rhor_kernel)
1602 : :
1603 [ - + ]: 150 : ABI_FREE(rhor)
1604 [ - + ]: 150 : ABI_FREE(rhog)
1605 [ - + ]: 150 : ABI_FREE(ks_vbik)
1606 [ - + ]: 150 : ABI_FREE(qp_vbik)
1607 [ - + ]: 150 : ABI_FREE(ktabr)
1608 [ - + ]: 150 : ABI_FREE(taur)
1609 [ - + ]: 150 : ABI_FREE(ks_vhartr)
1610 [ - + ]: 150 : ABI_FREE(ks_vtrial)
1611 [ - + ]: 150 : ABI_FREE(vpsp)
1612 [ - + ]: 150 : ABI_FREE(ks_vxc)
1613 [ - + ]: 150 : ABI_FREE(ph1d)
1614 [ - + ]: 150 : ABI_FREE(ph1df)
1615 [ - + ]: 150 : ABI_FREE(nhatgr)
1616 [ - + ]: 150 : ABI_FREE(nhat)
1617 : 150 : call pawfgr_destroy(Pawfgr)
1618 : :
1619 [ + + ]: 150 : if (Dtset%usepaw==1) then ! Optional deallocation for PAW.
1620 : 6 : call pawrhoij_free(Pawrhoij)
1621 : 6 : call pawfgrtab_free(Pawfgrtab)
1622 : 6 : call paw_ij_free(Paw_ij)
1623 : 6 : call paw_an_free(Paw_an)
1624 : 6 : call pawpwff_free(Paw_pwff)
1625 [ - + ]: 6 : 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 [ + - + + : 183 : ABI_FREE(Pawfgrtab)
- + - + -
+ - + - +
- + - + -
+ - + ]
1632 [ + - + + : 161 : ABI_FREE(Paw_pwff)
- + - + -
+ ]
1633 [ + - + + : 183 : ABI_FREE(Pawrhoij)
- + - + -
+ - + - +
- + - + ]
1634 [ + - + + : 183 : ABI_FREE(Paw_ij)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + ]
1635 [ + - + + : 183 : ABI_FREE(Paw_an)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
1636 [ - + ]: 150 : ABI_FREE(ktabrf)
1637 [ + - + + : 505 : ABI_FREE(Paw_onsite)
- + - + -
+ - + - +
- + ]
1638 : :
1639 : 150 : call wfd%free(); call Kmesh%free(); call Qmesh%free(); call cryst%free()
1640 : 150 : call Gsph_epsG0%free(); call Gsph_wfn%free(); call Vcp%free(); call Ep%free()
1641 : 150 : call Hdr_wfk%free(); call Hdr_local%free()
1642 : 150 : call ks_ebands%free(); call qp_ebands%free()
1643 : 150 : call destroy_mpi_enreg(MPI_enreg_seq)
1644 : 150 : call littlegroup_free(ltg_q)
1645 [ - + ]: 150 : ABI_FREE(Ltg_q)
1646 : :
1647 : 150 : call timab(301,2,tsec)
1648 : :
1649 [ - + - + : 900 : 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 [ - + - + : 15150 : 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),parameter :: OMEGAERMAX_TANGRID=36.74930883_dp ! 1000 eV, former default of cd_max_freq
1723 : : real(dp) :: ecutepspG0,ucvol,domegareal
1724 : : logical :: remove_inv,ltest,found,is_static,has_q0
1725 : : character(len=500) :: msg
1726 : 150 : type(wvl_internal_type) :: wvl
1727 : : !arrays
1728 : : integer :: ng0sh_opt(3), units(2)
1729 : 150 : integer,allocatable :: npwarr(:)
1730 : 150 : integer,pointer :: gvec_kss(:,:)
1731 : 150 : integer,pointer :: test_gvec_kss(:,:)
1732 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),qtmp(3),sq(3),qbz(3)
1733 : 150 : real(dp),pointer :: energies_p(:,:,:)
1734 : 150 : real(dp),allocatable :: doccde(:),eigen(:),occfact(:)
1735 : 150 : type(Pawrhoij_type),allocatable :: Pawrhoij(:)
1736 : :
1737 : : ! *************************************************************************
1738 : :
1739 [ + + ]: 450 : units = [std_out, ab_out]
1740 : :
1741 : : ! Check for calculations that are not implemented
1742 [ + + + - ]: 1016 : ltest = ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol) == Dtset%nband(1))
1743 [ - + ]: 150 : ABI_CHECK(ltest, 'dtset%nband(:) must be constant in the GW code.')
1744 : :
1745 : 150 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1746 : :
1747 : 150 : call mkrdim(acell,rprim,rprimd)
1748 : 150 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1749 : :
1750 : : ! Set up basic parameters of the calculation
1751 : 150 : Ep%gwcalctyp =Dtset%gwcalctyp
1752 : 150 : Ep%plasmon_pole_model =.TRUE.
1753 : 150 : Ep%analytic_continuation=.FALSE.
1754 : 150 : Ep%contour_deformation =.FALSE.
1755 : :
1756 : 150 : mod10=MOD(Ep%gwcalctyp,10)
1757 [ + + ]: 150 : if (mod10/=0.and.mod10/=8) Ep%plasmon_pole_model =.FALSE.
1758 [ + + ]: 150 : if (mod10==1) Ep%analytic_continuation=.TRUE.
1759 [ + + ]: 150 : if (mod10==2.or.mod10==9) Ep%contour_deformation =.TRUE.
1760 : 150 : is_static=(mod10==5.or.mod10==6.or.mod10==7)
1761 : :
1762 : 150 : Ep%nbnds =Dtset%nband(1)
1763 : 150 : Ep%symchi =Dtset%symchi
1764 [ + + ]: 150 : Ep%inclvkb=Dtset%inclvkb; if (Dtset%usepaw/=0) Ep%inclvkb=0
1765 : 150 : Ep%zcut =Dtset%zcut
1766 : :
1767 : 150 : write(msg,'(2a,i4,2a,f10.6,a)')ch10,&
1768 : 150 : ' GW calculation type = ',Ep%gwcalctyp,ch10,&
1769 : 300 : ' zcut to avoid poles in chi0 [eV] = ',Ep%zcut*Ha_eV,ch10
1770 : 150 : call wrtout(std_out, msg)
1771 : :
1772 : 150 : Ep%awtr =Dtset%awtr
1773 : 150 : Ep%npwe =Dtset%npweps
1774 : 150 : Ep%npwwfn=Dtset%npwwfn
1775 : 150 : Ep%npwvec=MAX(Ep%npwe,Ep%npwwfn)
1776 : :
1777 : 150 : timrev = 2 ! This information is not reported in the header
1778 : : ! 1 --> do not use time-reversal symmetry
1779 : : ! 2 --> take advantage of time-reversal symmetry
1780 [ + + ]: 150 : if (any(dtset%kptopt == [3, 4])) timrev = 1
1781 : :
1782 [ + + - + ]: 150 : if (timrev==1.and.Dtset%awtr/=0) then
1783 : 0 : ABI_ERROR("awtr/=0 cannot be used when time-reversal symmetry doesn't hold")
1784 : : end if
1785 : :
1786 : : ! Read parameters from WFK and verifify them.
1787 : 150 : call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm)
1788 [ + + ]: 1016 : mband = MAXVAL(Hdr_wfk%nband)
1789 : 150 : call hdr_wfk%vs_dtset(dtset)
1790 : 150 : remove_inv=.FALSE.
1791 : :
1792 : 150 : test_npwkss = 0
1793 : : call make_gvec_kss(Dtset%nkpt,Dtset%kptns,Hdr_wfk%ecut_eff,Dtset%symmorphi,Dtset%nsym,Dtset%symrel,Dtset%tnons,&
1794 : 150 : gprimd,Dtset%prtvol,test_npwkss,test_gvec_kss,ierr)
1795 [ - + ]: 150 : ABI_CHECK(ierr==0,"Fatal error in make_gvec_kss")
1796 : :
1797 [ - + + - : 450 : ABI_MALLOC(gvec_kss,(3,test_npwkss))
+ - + - -
+ - + ]
1798 [ + - + + : 964502 : gvec_kss = test_gvec_kss
+ + + + +
+ ]
1799 : 150 : ng_kss = test_npwkss
1800 : :
1801 [ + + ]: 150 : if (Ep%npwvec>ng_kss) then
1802 : 1 : Ep%npwvec=ng_kss
1803 [ + - ]: 1 : if (Ep%npwwfn> ng_kss) Ep%npwwfn=ng_kss
1804 [ - + ]: 1 : if (Ep%npwe > ng_kss) Ep%npwe =ng_kss
1805 : 1 : write(msg,'(3a,3(a,i6,a))')ch10,&
1806 : 1 : ' Number of G-vectors found less then required. Calculation will proceed with ',ch10,&
1807 : 1 : ' npwvec = ',Ep%npwvec,ch10,&
1808 : 1 : ' npweps = ',Ep%npwe ,ch10,&
1809 : 2 : ' npwwfn = ',Ep%npwwfn,ch10
1810 : 1 : ABI_WARNING(msg)
1811 : : end if
1812 : :
1813 : 150 : ng = MIN(SIZE(gvec_kss,DIM=2),SIZE(test_gvec_kss,DIM=2))
1814 : 150 : ierr = 0
1815 [ + + ]: 120694 : do ig=1,ng
1816 [ + + + - : 482326 : if (ANY(gvec_kss(:,ig)/=test_gvec_kss(:,ig))) then
- + ]
1817 : 0 : ierr=ierr+1
1818 : 0 : write(std_out,*)" gvec_kss ",ig,"/",ng,gvec_kss(:,ig),test_gvec_kss(:,ig)
1819 : : end if
1820 : : end do
1821 [ - + ]: 150 : ABI_CHECK(ierr == 0, "Mismatch between gvec_kss and test_gvec_kss")
1822 [ - + ]: 150 : ABI_FREE(test_gvec_kss)
1823 : :
1824 : : ! Get important dimension from Hdr_wfk
1825 : : ! Check also the consistency btw Hdr_wfk and Dtset.
1826 : 150 : Ep%nsppol=Hdr_wfk%nsppol
1827 : 150 : Ep%nkibz =Hdr_wfk%nkpt
1828 : :
1829 [ - + ]: 150 : if (Ep%nbnds>mband) then
1830 : : write(msg,'(2a,2(a,i0))') &
1831 : 0 : 'Number of bands stored in WFK file is less than required. ',ch10,&
1832 : 0 : "WFK mband:", mband, ", chi0 nband: ", Ep%nbnds
1833 : 0 : ABI_ERROR(msg)
1834 : : end if
1835 : :
1836 [ - + - + : 150 : cryst = Hdr_wfk%get_crystal(gw_timrev=timrev, remove_inv=remove_inv)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + ]
1837 : 150 : call cryst%print(mode_paral='COLL')
1838 : :
1839 : : ! === Create basic data types for the calculation ===
1840 : : ! Kmesh defines the k-point sampling for the wavefunctions.
1841 : : ! Qmesh defines the q-point sampling for chi0, all possible differences k1-k2 reduced to the IBZ.
1842 : : ! TODO Kmesh%bz should be in [-half, half[ but this modification is painful!
1843 : :
1844 : 150 : call Kmesh%init(cryst, ep%nkibz, hdr_wfk%kptns, dtset%kptopt, wrap_1zone=.FALSE.)
1845 : :
1846 : : ! Some required information are not filled up inside kmesh_init. So doing it here, even though it is not clean
1847 [ + + + + ]: 1950 : Kmesh%kptrlatt(:,:) = Dtset%kptrlatt(:,:)
1848 : 150 : Kmesh%nshift = Dtset%nshiftk
1849 [ - + + - : 450 : ABI_MALLOC(Kmesh%shift, (3,Kmesh%nshift))
+ - + - -
+ - + -
+ ]
1850 [ + + + + ]: 750 : Kmesh%shift(:,:) = Dtset%shiftk(:,1:Dtset%nshiftk)
1851 : 150 : call Kmesh%print(units, header="K-mesh for the wavefunctions", prtvol=Dtset%prtvol)
1852 : :
1853 : : ! === Find Q-mesh ===
1854 : : ! Stop if a nonzero umklapp is needed to reconstruct the BZ.
1855 : : ! epsilon^-1(Sq) indeed should be symmetrized in csigme using a different expression (G-G_o is needed)
1856 : 150 : call qmesh%find_qmesh(Cryst, Kmesh)
1857 : 150 : call qmesh%print(units, "Q-mesh for the screening function", prtvol=dtset%prtvol)
1858 : :
1859 [ + + ]: 7126 : do iqbz=1,Qmesh%nbz
1860 : 6976 : call qmesh%get_BZ_item(iqbz, qbz, iq_ibz, isym, itim)
1861 [ + + + + : 195328 : sq = (3-2*itim) * MATMUL(cryst%symrec(:,:,isym), qmesh%ibz(:,iq_ibz))
+ + + + +
- + - +
+ ]
1862 [ + + + - : 28054 : if (ANY(ABS(qbz-sq) > 1.0d-4)) then
- + ]
1863 : : write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')&
1864 : 0 : ' qpoint ',qbz,' is the symmetric of ',qmesh%ibz(:,iq_ibz),ch10,&
1865 : 0 : ' through operation ',cryst%symrec(:,:,isym),' and itim ',itim,ch10,&
1866 : 0 : ' however a non zero umklapp G_o vector is required and this is not yet allowed'
1867 : 0 : ABI_ERROR(msg)
1868 : : end if
1869 : : end do
1870 : :
1871 : : ! Setup for long wavelength limit
1872 [ + + ]: 150 : if (Dtset%gw_nqlwl==0) then
1873 : 149 : Ep%nqlwl=1
1874 [ - + - + ]: 149 : ABI_MALLOC(Ep%qlwl,(3,Ep%nqlwl))
1875 [ + + ]: 596 : Ep%qlwl(:,1)=GW_Q0_DEFAULT ! Use default shift 0.000010, 0.000020, 0.000030
1876 : : else
1877 : 1 : Ep%nqlwl=Dtset%gw_nqlwl
1878 [ - + + - : 3 : ABI_MALLOC(Ep%qlwl,(3,Ep%nqlwl))
+ - + - -
+ - + -
+ ]
1879 [ + + + + ]: 5 : Ep%qlwl(:,:)=Dtset%gw_qlwl(:,1:Ep%nqlwl)
1880 [ - + ]: 1 : ABI_CHECK(Ep%nqlwl==1,"nqlwl/=1 not coded")
1881 : : end if
1882 : : !write(std_out,*)" Using qlwl = ",Ep%qlwl
1883 : :
1884 : : ! Find optimal value for G-sphere enlargment due to oscillator matrix elements
1885 : 150 : call get_ng0sh(Kmesh%nbz,Kmesh%bz,Qmesh%nibz,Qmesh%ibz,Kmesh%nbz,Kmesh%bz,GW_TOLQ0,ng0sh_opt)
1886 : 150 : call wrtout(std_out,sjoin(' Optimal value for ng0sh:',ltoa(ng0sh_opt)),"COLL")
1887 : :
1888 [ + + ]: 600 : Ep%mG0(:)=ng0sh_opt(:) !Ep%mG0(:) = [3, 3, 3]
1889 : :
1890 : : ! === In case of symmetrization, find the little group of the q"s ===
1891 : : ! * For the long-wavelength limit we consider a small but finite q. However the oscillators are
1892 : : ! evaluated setting q==0. Thus it is possible to take advantage of symmetries also when q --> 0.
1893 : : ! * Here we calculate the enlargement of the G-sphere, npwepG0, needed to account for umklapps.
1894 : : ! TODO Switch on use_umklp, write all this stuff to ab_out
1895 : :
1896 : 150 : Ep%npwepG0 = Ep%npwe
1897 [ - + + - : 1393 : ABI_MALLOC(Ltg_q, (Qmesh%nibz))
+ - + - -
+ - + +
+ ]
1898 : :
1899 [ + + ]: 1093 : do iq=1,Qmesh%nibz
1900 [ + + + + ]: 3772 : qtmp = Qmesh%ibz(:,iq); if (normv(qtmp,gmet,'G') < GW_TOLQ0) qtmp(:) = zero; use_umklp = 0
1901 : 1093 : call Ltg_q(iq)%init(qtmp, Kmesh%nbz, Kmesh%bz, Cryst, use_umklp, Ep%npwe, gvec=gvec_kss)
1902 : : end do
1903 : :
1904 : 150 : ecutepspG0 = Dtset%ecuteps
1905 [ - + ]: 150 : ABI_CHECK(ecutepspG0 > zero, "ecuteps must be > 0")
1906 [ + + ]: 150 : if (Ep%symchi/=0) then
1907 [ + - - + : 1130 : ecutepspG0=MAXVAL(Ltg_q(:)%max_kin_gmG0)+tol6; npwepG0=0; nshepspG0=0
- - + + ]
1908 [ + + ]: 136 : if (my_rank == master) write(std_out,*)" Due to umklapp processes : ecutepspg0= ",ecutepspG0
1909 : 136 : call setshells(ecutepspG0,npwepG0,nshepspG0,Cryst%nsym,gmet,gprimd,Cryst%symrel,'eps_pG0',Cryst%ucvol)
1910 : 136 : Ep%npwepG0=npwepG0
1911 : : end if
1912 : :
1913 [ - + ]: 150 : if (Ep%npwepG0>Ep%npwvec) then
1914 : : write(msg,'(3a,i5,a,i5)')&
1915 : 0 : ' npwepG0 > npwvec, decrease npweps or increase npwwfn. ',ch10,&
1916 : 0 : ' npwepG0 = ',Ep%npwepG0,' npwvec = ',Ep%npwvec
1917 : 0 : ABI_ERROR(msg)
1918 : : end if
1919 : :
1920 : : ! === Create structure describing the G-sphere used for chi0/espilon and Wfns ===
1921 : : ! * The cutoff is >= ecuteps to allow for umklapp
1922 [ - + ]: 150 : call Gsph_wfn%init(Cryst, Ep%npwvec, gvec=gvec_kss)
1923 : :
1924 [ - + ]: 150 : call Gsph_epsG0%init(Cryst, Ep%npwepG0, gvec=gvec_kss)
1925 : : !
1926 : : ! =======================================================================
1927 : : ! ==== Setup of the FFT mesh used for the oscillator matrix elements ====
1928 : : ! =======================================================================
1929 : : ! * ngfft_gw(7:18) is the same as Dtset%ngfft(7:18), initialized before entering setup_screening.
1930 : : ! Here we just redefine ngfft_gw(1:6) according to the following options:
1931 : : !
1932 : : ! method==0 ==> FFT grid read from __fft.in__ (only for debugging purpose)
1933 : : ! method==1 ==> normal FFT grid
1934 : : ! method==2 ==> slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90)
1935 : : ! method==3 ==> doubled FFT grid, to treat exactly the convolution defining the density,
1936 : : ! Useful in sigma if ppmodel=[2,3,4] since rho(G-Gp) or to calculate matrix elements of v_Hxc.
1937 : : !
1938 : : ! enforce_sym==1 ==> enforce a direct space FFT mesh compatible with all symmetries operation
1939 : : ! enforce_sym==0 ==> Find the smallest FFT grid compatibile with the library, do not care about symmetries
1940 : : !
1941 [ + + ]: 2850 : ngfft_gw(1:18)=Dtset%ngfft(1:18); method=2
1942 [ - + ]: 150 : if (Dtset%fftgw==00 .or. Dtset%fftgw==01) method=0
1943 [ + + ]: 150 : if (Dtset%fftgw==10 .or. Dtset%fftgw==11) method=1
1944 [ + + ]: 150 : if (Dtset%fftgw==20 .or. Dtset%fftgw==21) method=2
1945 [ + + ]: 150 : if (Dtset%fftgw==30 .or. Dtset%fftgw==31) method=3
1946 : 150 : enforce_sym=MOD(Dtset%fftgw,10)
1947 : :
1948 : : ! Use npwepG0 to account for umklapps.
1949 [ - + ]: 150 : call setmesh(gmet,gvec_kss,ngfft_gw,Ep%npwvec,Ep%npwepG0,Ep%npwwfn,nfftgw_tot,method,Ep%mG0,Cryst,enforce_sym)
1950 : : !call new_setmesh(Cryst,ecut_osc,ecutwfn,nkpt,kpoints,method,Ep%mG0,enforce_sym,ngfft_gw,nfftgw_tot)
1951 : :
1952 [ - + ]: 150 : ABI_FREE(gvec_kss)
1953 : :
1954 : : ! FIXME this wont work if nqptdm/=0
1955 : : call Vcp%init(Gsph_epsG0,Cryst,Qmesh,Kmesh,Dtset%gw_rcut,Dtset%gw_icutcoul,Dtset%vcutgeo,Ep%npwe,Ep%nqlwl,&
1956 : 150 : Ep%qlwl,comm)
1957 : :
1958 : : #if 0
1959 : : ! Using the random q for the optical limit is one of the reasons
1960 : : ! why sigma breaks the initial energy degeneracies.
1961 : : Vcp%i_sz=zero
1962 : : Vcp%vc_sqrt(1,:)=czero
1963 : : Vcp%vcqlwl_sqrt(1,:)=czero
1964 : : #endif
1965 : :
1966 : : ! Value of scissor energy
1967 : 150 : Ep%mbpt_sciss=Dtset%mbpt_sciss
1968 : :
1969 : : ! Define the frequency mesh for epsilon according to the method used.
1970 : 150 : Ep%nomegaei=1
1971 [ + + ]: 150 : Ep%nomegaer=1; if (is_static) Ep%nomegaer=0
1972 : 150 : Ep%nomegaec=0
1973 : 150 : Ep%omegaermax=zero
1974 : :
1975 : : ! For ppmodels 2,3,4, only omega=0 is needed.
1976 [ + + + + : 150 : if (Ep%plasmon_pole_model.and.Dtset%nfreqre==1.and.Dtset%nfreqim==0) then
+ - ]
1977 : 5 : Ep%nomegaer=1; Ep%nomegaei=0
1978 : 5 : write(msg,'(7a)')ch10,&
1979 : 5 : ' The inverse dielectric matrix will be calculated on zero frequency only',ch10,&
1980 : 5 : ' please note that the calculated epsilon^-1 cannot be used ',ch10,&
1981 : 10 : ' to calculate QP corrections using plasmonpole model 1',ch10
1982 : 5 : call wrtout(units, msg)
1983 : : end if
1984 : :
1985 : : ! Max number of omega along the imaginary axis
1986 [ + + + + ]: 150 : if (Ep%analytic_continuation.or.Ep%contour_deformation) then
1987 : 37 : Ep%nomegaei = abs(Dtset%nfreqim)
1988 [ + + ]: 37 : if (Dtset%gw_frqim_inzgrid==1) then
1989 : 1 : ABI_WARNING('iomega = z/1-z transfom grid will be used for imaginary frequency grid')
1990 : : end if
1991 [ + + ]: 37 : if (Dtset%cd_customnimfrqs/=0) then
1992 : 1 : ABI_WARNING('Custom imaginary grid specified. Assuming experienced user.')
1993 : 1 : Ep%nomegaei=Dtset%cd_customnimfrqs
1994 : : end if
1995 [ - + ]: 37 : if (Ep%nomegaei==-1) then
1996 : 0 : Ep%nomegaei=NOMEGAGAUSS
1997 : 0 : ABI_WARNING(sjoin('Number of imaginary frequencies set to default= ',itoa(NOMEGAGAUSS)))
1998 : : end if
1999 [ + + ]: 37 : if (Ep%nomegaei==0) then
2000 : 3 : ABI_WARNING('nfreqim = 0! Assuming experienced user merging several frequency calculations.')
2001 : : end if
2002 : : end if
2003 : :
2004 : : ! Range and total number of real frequencies.
2005 : 150 : Ep%omegaermin = zero
2006 [ + + ]: 150 : if (Ep%contour_deformation) then
2007 : 26 : Ep%nomegaer=Dtset%nfreqre; Ep%omegaermin=Dtset%freqremin; Ep%omegaermax=Dtset%freqremax
2008 [ + + ]: 26 : if (Dtset%gw_frqre_tangrid==1) then
2009 : 1 : ABI_WARNING('Tangent transform grid will be used for real frequency grid')
2010 [ - + ]: 1 : if (Ep%omegaermin>TOL16) then
2011 : 0 : ABI_WARNING('freqremin is ignored with gw_frqre_tangrid=1, the tangent grid always starts at zero. Use cd_subset_freq to split the grid.')
2012 : 0 : Ep%omegaermin=zero
2013 : : end if
2014 [ - + ]: 1 : if (Ep%omegaermax<TOL16) then
2015 : 0 : Ep%omegaermax=OMEGAERMAX_TANGRID
2016 : 0 : write(msg,'(a,f8.4)')' Max real frequency of tangent grid set to default [Ha] = ',OMEGAERMAX_TANGRID
2017 : 0 : ABI_WARNING(msg)
2018 : : end if
2019 : : end if
2020 [ - + ]: 26 : if (Ep%nomegaer==-1) then
2021 : 0 : Ep%nomegaer=NOMEGAREAL
2022 : 0 : ABI_WARNING(sjoin('Number of real frequencies set to default= ',itoa(NOMEGAREAL)))
2023 : : end if
2024 [ - + ]: 26 : if (Ep%nomegaer==0) then
2025 : 0 : ABI_WARNING('nfreqre = 0 ! Assuming experienced user.')
2026 : : end if
2027 [ + + ]: 26 : if (ABS(Ep%omegaermin)<TOL16) then
2028 : 25 : Ep%omegaermin=zero
2029 : 25 : write(msg,'(a,f8.4)')' Min real frequency set to default [Ha] = ',Ep%omegaermin
2030 : 25 : ABI_WARNING(msg)
2031 : : end if
2032 [ - + ]: 26 : if (Ep%omegaermin>Ep%omegaermax) then
2033 : 0 : ABI_ERROR('freqremin > freqremax !')
2034 : : end if
2035 [ + + ]: 26 : if (Ep%omegaermax<TOL16) then
2036 : 4 : Ep%omegaermax=OMEGAERMAX
2037 : 4 : write(msg,'(a,f8.4)')' Max real frequency set to default [Ha] = ',OMEGAERMAX
2038 : 4 : ABI_WARNING(msg)
2039 : : end if
2040 : : ! Check if a subset of the frequencies is to be used
2041 [ - + ]: 26 : if (Dtset%cd_subset_freq(1)/=0) then
2042 : 0 : istart = Dtset%cd_subset_freq(1)
2043 : 0 : iend = Dtset%cd_subset_freq(2)
2044 [ # # # # ]: 0 : if (istart>iend.or.istart<0.or.iend<0) then
2045 : 0 : ABI_ERROR(' check indices of cd_subset_freq!')
2046 : : end if
2047 : 0 : write(msg,'(2(a,i0))')' Using cd_subset_freq to only do freq. from ',istart,' to ',iend
2048 : 0 : ABI_WARNING(msg)
2049 : : ! Reset the numbers
2050 [ # # ]: 0 : if (Dtset%gw_frqre_tangrid/=1) then ! Normal equidistant grid
2051 : 0 : Ep%nomegaer = iend-istart+1
2052 : 0 : domegareal=(Ep%omegaermax-Ep%omegaermin)/(Ep%nomegaer-1)
2053 : 0 : Ep%omegaermin = Ep%omegaermin+(istart-1)*domegareal
2054 : 0 : Ep%omegaermax = Ep%omegaermin+(iend-1)*domegareal
2055 : : else
2056 : 0 : Ep%nomegaer = iend-istart+1
2057 : : end if
2058 : : end if
2059 : : end if
2060 : :
2061 : : ! Check full grid calculations
2062 [ + + ]: 150 : if (Dtset%cd_full_grid/=0) then
2063 : 1 : ABI_WARNING("FULL GRID IN COMPLEX PLANE CALCULATED. YOU MIGHT NOT BE ABLE TO USE SCREENING FILES!")
2064 [ - + ]: 1 : if (Dtset%cd_subset_freq(1)/=0) then
2065 : 0 : ABI_ERROR('cd_subset_freq cannot be used with cd_full_grid!')
2066 : : end if
2067 : 1 : Ep%nomegaec = Ep%nomegaei*(Ep%nomegaer-1)
2068 : : end if
2069 : :
2070 : 150 : Ep%nomega=Ep%nomegaer+Ep%nomegaei+Ep%nomegaec ! Total number of frequencies.
2071 : :
2072 : : ! ==== Setup of the spectral method ====
2073 : 150 : Ep%spmeth =Dtset%spmeth; Ep%nomegasf=Dtset%nomegasf; Ep%spsmear =Dtset%spbroad
2074 : :
2075 [ + + ]: 150 : if (Ep%spmeth/=0) then
2076 : 6 : write(msg,'(2a,i3,2a,i8)')ch10,&
2077 : 6 : ' setup_screening: using spectral method: ',Ep%spmeth,ch10,&
2078 : 12 : ' Number of frequencies for imaginary part: ',Ep%nomegasf
2079 : 6 : call wrtout(std_out, msg)
2080 [ + + ]: 6 : if (Ep%spmeth==2) then
2081 : 1 : write(msg,'(a,f8.5,a)')' Gaussian broadening = ',Ep%spsmear*Ha_eV,' [eV]'
2082 : 1 : call wrtout(std_out, msg)
2083 : : end if
2084 : : end if
2085 : :
2086 : 150 : Ep%nI=1; Ep%nJ=1
2087 [ + + ]: 150 : if (Dtset%nspinor==2) then
2088 : : !if (Dtset%usepaw==1.and.Dtset%pawspnorb>0) then
2089 : : ! Ep%nI=1; Ep%nJ=4
2090 : : !end if
2091 : : ! For spin-spin interaction
2092 : : ! Ep%nI=4; Ep%nJ=4
2093 [ - + ]: 2 : ABI_CHECK(Ep%npwepG0 == Ep%npwe, "npwepG0 must be == npwe if nspinor==2")
2094 : : !ABI_CHECK(Ep%symchi == 0, "symchi/=0 and nspinor=2 not available")
2095 : : end if
2096 : :
2097 : : ! === Enable the calculations of chi0 on user-specified q-points ===
2098 : 150 : Ep%nqibz=Qmesh%nibz
2099 [ - + + - : 450 : ABI_MALLOC(Ep%qibz,(3,Ep%nqibz))
+ - + - -
+ - + -
+ ]
2100 [ + + + + ]: 3922 : Ep%qibz(:,:)=Qmesh%ibz(:,:)
2101 : :
2102 : 150 : Ep%nqcalc=Ep%nqibz
2103 [ + + ]: 150 : if (Dtset%nqptdm>0) Ep%nqcalc=Dtset%nqptdm
2104 : :
2105 [ - + + - : 450 : ABI_MALLOC(Ep%qcalc,(3,Ep%nqcalc))
+ - + - -
+ - + -
+ ]
2106 [ + + ]: 150 : if (Ep%nqcalc/=Ep%nqibz) then
2107 : 3 : write(msg,'(6a)')ch10,&
2108 : 3 : ' Dielectric matrix will be calculated only for some ',ch10,&
2109 : 3 : ' selected q points provided by the user through the input variables ',ch10,&
2110 : 6 : ' nqptdm and qptdm'
2111 : 3 : call wrtout(units, msg)
2112 : 3 : ltest= Ep%nqcalc <= Qmesh%nibz
2113 [ - + ]: 3 : ABI_CHECK(ltest, 'nqptdm should not exceed the number of q points in the IBZ')
2114 [ + + + + ]: 27 : Ep%qcalc(:,:)=Dtset%qptdm(:,1:Ep%nqcalc)
2115 : : ! Check whether the q-points provided are correct.
2116 [ + + ]: 9 : do iq=1,Ep%nqcalc
2117 : 6 : found=.FALSE.
2118 [ + - ]: 21 : do iqp=1,Qmesh%nibz
2119 [ + + ]: 84 : qtmp(:)=Ep%qcalc(:,iq)-Qmesh%ibz(:,iqp)
2120 : 21 : found=(normv(qtmp,gmet,'G')<GW_TOLQ)
2121 [ + + ]: 21 : if (found) EXIT
2122 : : end do
2123 [ - + ]: 9 : ABI_CHECK(found, 'One or more points specified by Dtset%qptdm do not satisfy q=k1-k2')
2124 : : end do
2125 : : else
2126 [ + + + + ]: 3847 : Ep%qcalc(:,:)=Ep%qibz(:,:)
2127 : : end if
2128 : :
2129 : : ! To write the SCR header correctly, with heads and wings, we have
2130 : : ! to make sure that q==0, if present, is the first q-point in the list.
2131 : 150 : has_q0=.FALSE.
2132 [ + + ]: 155 : do iq=1,Ep%nqcalc
2133 [ + + ]: 155 : if (normv(Ep%qcalc(:,iq),gmet,'G')<GW_TOLQ0) then
2134 : 150 : has_q0=.TRUE.; EXIT
2135 : : end if
2136 : : end do
2137 : :
2138 [ + + - + ]: 150 : if (has_q0.and.normv(Ep%qcalc(:,1),gmet,'G')>=GW_TOLQ0) then
2139 : : write(msg,'(5a)')&
2140 : 0 : 'The list of q-points to be calculated contains the Gamma point, ',ch10,&
2141 : 0 : 'however Gamma is not the first point in the list. ' ,ch10,&
2142 : 0 : 'Please, change your input file accordingly. '
2143 : 0 : ABI_ERROR(msg)
2144 : : end if
2145 : :
2146 : : ! === Initialize the band structure datatype ===
2147 : : ! * Copy KSS energies and occupations up to Ep%nbnds==Dtset%nband(:)
2148 : : ! TODO Recheck symmorphy and inversion
2149 [ + + ]: 1016 : bantot = SUM(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol))
2150 : :
2151 [ - + + - : 22800 : ABI_CALLOC(doccde, (bantot))
+ - + - -
+ - + +
+ ]
2152 [ - + + - : 22650 : ABI_CALLOC(eigen, (bantot))
- + - + +
+ ]
2153 [ - + + - : 22650 : ABI_CALLOC(occfact, (bantot))
- + - + +
+ ]
2154 : :
2155 : : jj=0; ibtot=0
2156 [ + + ]: 302 : do isppol=1,Dtset%nsppol
2157 [ + + ]: 1168 : do ikibz=1,Dtset%nkpt
2158 [ + + ]: 37846 : do ib=1,Hdr_wfk%nband(ikibz+Dtset%nkpt*(isppol-1))
2159 : 36828 : ibtot=ibtot+1
2160 [ + + ]: 37694 : if (ib<=Ep%nbnds) then
2161 : 22350 : jj=jj+1
2162 : 22350 : occfact(jj)=Hdr_wfk%occ(ibtot)
2163 : 22350 : eigen (jj)=energies_p(ib,ikibz,isppol)
2164 : : end if
2165 : : end do
2166 : : end do
2167 : : end do
2168 [ - + ]: 150 : ABI_FREE(energies_p)
2169 : :
2170 : : ! Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of
2171 : : ! the symmetry operations in the old GW code (symmorphy and inversion)
2172 [ + + + - ]: 1008 : ltest = (ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz)) < tol6))
2173 [ - + ]: 150 : if (.not. ltest) then
2174 [ # # ]: 0 : do jj=1,Kmesh%nibz
2175 : 0 : write(std_out, *)"wtk dtset vs kmesh:", dtset%wtk(jj), kmesh%wt(jj)
2176 : : end do
2177 : : end if
2178 : 0 : ABI_CHECK(ltest, 'Mismatch between Dtset%wtk and Kmesh%wt')
2179 : :
2180 [ - + + - : 450 : ABI_MALLOC(npwarr,(Hdr_wfk%nkpt))
+ - + - -
+ - + ]
2181 [ + + ]: 1008 : npwarr(:)=Ep%npwwfn
2182 : :
2183 : : call ks_ebands%init(bantot, Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,&
2184 : : doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,&
2185 : : Kmesh%nibz,npwarr,Dtset%nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,&
2186 : : dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
2187 : 150 : dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
2188 : :
2189 : : ! TODO modify outkss in order to calculate the eigenvalues also if NSCF calculation.
2190 : : ! this fails simply because in case of NSCF occ are zero
2191 : : !ltest=(ALL(ABS(occfact-ks_ebands%occ)<1.d-2))
2192 : : !call assert(ltest,'difference in occfact')
2193 : : !write(std_out,*)MAXVAL(ABS(occfact(:)-ks_ebands%occ(:)))
2194 : :
2195 : : !TODO call ebands_update_occ here
2196 : : !call ks_ebands%update_occ(spinmagntarget,Dtset%prtvol)
2197 : :
2198 [ - + ]: 150 : ABI_FREE(doccde)
2199 [ - + ]: 150 : ABI_FREE(eigen)
2200 [ - + ]: 150 : ABI_FREE(npwarr)
2201 : :
2202 : : ! Initialize abinit header for the screening part
2203 : 150 : call Hdr_out%init(ks_ebands,codvsn,Dtset,Pawtab,pertcase0,Psps,wvl)
2204 : :
2205 : : ! Get Pawrhoij from the header.
2206 [ + + + - : 483 : ABI_MALLOC(Pawrhoij,(Cryst%natom*Dtset%usepaw))
+ - + + -
+ - + - +
+ + ]
2207 [ + + ]: 150 : if (Dtset%usepaw==1) then
2208 : 6 : call pawrhoij_alloc(Pawrhoij,1,Dtset%nspden,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab)
2209 : 6 : call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij)
2210 : : end if
2211 : 150 : call Hdr_out%update(bantot,1.0d20,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1))
2212 : :
2213 [ - + ]: 150 : ABI_FREE(occfact)
2214 : 150 : call pawrhoij_free(Pawrhoij)
2215 [ + - + + : 183 : ABI_FREE(Pawrhoij)
- + - + -
+ - + - +
- + - + ]
2216 : :
2217 : : ! ==== Setup of extrapolar technique ====
2218 : 150 : Ep%gwcomp = Dtset%gwcomp; Ep%gwencomp = Dtset%gwencomp
2219 [ + + ]: 150 : if (Ep%gwcomp == 1) then
2220 : 19 : write(msg,'(a,f8.2,a)')' Using the completeness correction with gwencomp ',Ep%gwencomp*Ha_eV,' [eV] '
2221 : 19 : call wrtout(std_out, msg)
2222 : : end if
2223 : :
2224 : : ! Final compatibility tests
2225 [ + + + - : 1016 : ltest = (ks_ebands%mband == Ep%nbnds .and. ALL(ks_ebands%nband == Ep%nbnds))
+ - - + ]
2226 : 0 : ABI_CHECK(ltest, 'BUG in definition of ks_ebands%nband')
2227 : :
2228 [ + + - + ]: 150 : if (Ep%gwcomp==1 .and. Ep%spmeth>0) then
2229 : 0 : ABI_ERROR("Hilbert transform and extrapolar method are not compatible")
2230 : : end if
2231 : :
2232 [ - + - + : 450 : end subroutine setup_screening
- + - - -
- - - - -
- - - - -
- ]
2233 : : !!***
2234 : :
2235 : : !----------------------------------------------------------------------
2236 : :
2237 : : !!****f* m_screening_driver/chi0_bksmask
2238 : : !! NAME
2239 : : !! chi0_bksmask
2240 : : !!
2241 : : !! FUNCTION
2242 : : !! Compute tables for the distribution and the storage of the wavefunctions in the SCREENING code.
2243 : : !!
2244 : : !! INPUTS
2245 : : !! Dtset<type(dataset_type)>=all input variables for this dataset
2246 : : !! Ep<em1params_t>=Parameters for the screening calculation.
2247 : : !! Kmesh <kmesh_t>=Structure describing the k-point sampling.
2248 : : !! nbvw = Max. number of fully/partially occupied states over spin
2249 : : !! nbcw = Max. number of unoccupied states considering the spin
2250 : : !! nprocs=Total number of MPI processors
2251 : : !! my_rank=Rank of this this processor.
2252 : : !!
2253 : : !! OUTPUT
2254 : : !! bks_mask(Ep%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will treat this state.
2255 : : !! keep_ur(Ep%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will store this state in real space.
2256 : : !! ierr=Exit status.
2257 : : !!
2258 : : !! SOURCE
2259 : :
2260 : 150 : subroutine chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,nprocs,bks_mask,keep_ur,ierr)
2261 : :
2262 : : !Arguments ------------------------------------
2263 : : !scalars
2264 : : integer,intent(in) :: my_rank,nprocs,nbvw,nbcw
2265 : : integer,intent(out) :: ierr
2266 : : type(Dataset_type),intent(in) :: Dtset
2267 : : type(em1params_t),intent(in) :: Ep
2268 : : type(kmesh_t),intent(in) :: Kmesh
2269 : : !arrays
2270 : : logical,intent(out) :: bks_mask(Ep%nbnds,Kmesh%nibz,Dtset%nsppol)
2271 : : logical,intent(out) :: keep_ur(Ep%nbnds,Kmesh%nibz,Dtset%nsppol)
2272 : :
2273 : : !Local variables-------------------------------
2274 : : !scalars
2275 : : integer :: my_nspins,my_maxb,my_minb,isp,spin,nsppol,band,rank_spin,ib
2276 : : character(len=500) :: msg
2277 : : logical :: store_ur
2278 : : !arrays
2279 : 300 : integer :: my_spins(Dtset%nsppol),nprocs_spin(Dtset%nsppol)
2280 : 150 : integer,allocatable :: istart(:),istop(:)
2281 : :
2282 : : ! *************************************************************************
2283 : :
2284 : 150 : ierr=0; nsppol=Dtset%nsppol
2285 : :
2286 [ + - + + : 454 : my_nspins=Dtset%nsppol; my_spins= [(isp,isp=1,nsppol)]
+ + ]
2287 : :
2288 : : ! List of spins for each node, number of processors per each spin
2289 : : ! and the MPI rank in the "spin" communicator.
2290 [ + + ]: 302 : nprocs_spin = nprocs; rank_spin = my_rank
2291 : :
2292 [ + + - + ]: 150 : if (nsppol==2.and.nprocs>1) then
2293 : : ! Distribute spins (optimal distribution if nprocs is even)
2294 : 0 : nprocs_spin(1) = nprocs/2
2295 : 0 : nprocs_spin(2) = nprocs - nprocs/2
2296 : 0 : my_nspins=1; my_spins(1)=1
2297 [ # # ]: 0 : if (my_rank+1>nprocs/2) then
2298 : 0 : my_spins(1)=2
2299 : 0 : rank_spin = my_rank - nprocs/2
2300 : : end if
2301 : : end if
2302 : :
2303 : 150 : store_ur = (MODULO(Dtset%gwmem,10)==1)
2304 [ + + + + : 46886 : bks_mask=.FALSE.; keep_ur=.FALSE.
+ + + + +
+ + + ]
2305 : :
2306 : 182 : select case (Dtset%gwpara)
2307 : : case (1)
2308 : : ! Parallelization over transitions **without** memory distributions (Except for the spin).
2309 : 32 : my_minb=1; my_maxb=Ep%nbnds
2310 [ + + ]: 64 : do isp=1,my_nspins
2311 : 32 : spin = my_spins(isp)
2312 [ + + + + ]: 3888 : bks_mask(my_minb:my_maxb,:,spin) = .TRUE.
2313 [ + - + + : 3920 : if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE.
+ + ]
2314 : : end do
2315 : :
2316 : : case (2)
2317 : : ! Distribute bands and spin.
2318 [ + + ]: 238 : do isp=1,my_nspins
2319 : 120 : spin = my_spins(isp)
2320 : :
2321 [ + - ]: 120 : if (nprocs_spin(spin) <= nbcw) then
2322 : : ! Distribute nbcw empty bands among nprocs_spin (block of bands without replicas).
2323 : : ! Bands are distributed in contiguous blocks because
2324 : : ! this distribution is well suited for the Hilber transform
2325 : : ! since each node will allocate only a smaller frequency interval
2326 : : ! for the spectral function whose size scales with the number of MPI nodes.
2327 : : ! Note it is now meaningless to distinguish gwcomp=0 or 1 since the workload is well balanced later on
2328 [ - + + - : 360 : ABI_MALLOC(istart,(nprocs_spin(spin)))
+ - + - -
+ - + ]
2329 [ - + + - : 240 : ABI_MALLOC(istop,(nprocs_spin(spin)))
- + - + ]
2330 : :
2331 : 120 : call xmpi_split_work2_i4b(nbcw,nprocs_spin(spin),istart,istop)
2332 : :
2333 : 120 : my_minb = nbvw + istart(rank_spin+1)
2334 : 120 : my_maxb = nbvw + istop (rank_spin+1)
2335 : :
2336 : 120 : ABI_FREE(istart)
2337 : 120 : ABI_FREE(istop)
2338 : :
2339 [ - + ]: 120 : if (my_maxb - my_minb + 1 <= 0) then
2340 : : write(msg,'(3a,2(i0,a),2a)')&
2341 : 0 : 'One or more processors has zero number of bands ',ch10,&
2342 : 0 : 'my_minb= ',my_minb,' my_maxb= ',my_maxb,ch10,&
2343 : 0 : 'This is a waste, decrease the number of processors.'
2344 : 0 : ABI_ERROR(msg)
2345 : : end if
2346 : :
2347 [ + + + + ]: 14023 : bks_mask(my_minb:my_maxb,:,spin)=.TRUE.
2348 [ + + + + : 13671 : if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE.
+ + ]
2349 : :
2350 : : else
2351 : : ! New version (alternate bands with replicas if nprocs > nbcw)
2352 : : ! FIXME: Fix segmentation fault with Hilbert transform.
2353 [ # # ]: 0 : do ib=1,nbcw
2354 [ # # ]: 0 : if (xmpi_distrib_with_replicas(ib,nbcw,rank_spin,nprocs_spin(spin))) then
2355 : 0 : band = ib + nbvw
2356 [ # # ]: 0 : bks_mask(band,:,spin)=.TRUE.
2357 [ # # # # ]: 0 : if (store_ur) keep_ur(band,:,spin)=.TRUE.
2358 : : end if
2359 : : end do
2360 : : end if
2361 : :
2362 : : ! This is needed to have all the occupied states on each node.
2363 [ + + + + ]: 3553 : bks_mask(1:nbvw,:,spin) = .TRUE.
2364 [ + + + + : 3563 : if (store_ur) keep_ur(1:nbvw,:,spin)=.TRUE.
+ + ]
2365 : : end do ! isp
2366 : :
2367 : : case default
2368 : 0 : ierr = 1
2369 [ + + - ]: 150 : ABI_WARNING("Wrong value for gwpara")
2370 : : end select
2371 : :
2372 : 150 : end subroutine chi0_bksmask
2373 : : !!***
2374 : :
2375 : : !!****f* m_screening_driver/random_stopping_power
2376 : : !! NAME
2377 : : !! random_stopping_power
2378 : : !!
2379 : : !! FUNCTION
2380 : : !! Calculate the electronic random stopping power
2381 : : !!
2382 : : !! INPUTS
2383 : : !!
2384 : : !! OUTPUT
2385 : : !!
2386 : : !! SIDE EFFECTS
2387 : : !!
2388 : : !! SOURCE
2389 : :
2390 : 12 : subroutine random_stopping_power(iqibz,npvel,pvelmax,Ep,Gsph_epsG0,Qmesh,Vcp,Cryst,Dtfil,epsm1,rspower)
2391 : :
2392 : : use m_splines
2393 : :
2394 : : !Arguments ------------------------------------
2395 : : !scalars
2396 : : integer,intent(in) :: iqibz,npvel
2397 : : real(dp),intent(in) :: pvelmax(3)
2398 : : type(em1params_t),intent(in) :: Ep
2399 : : type(gsphere_t),intent(in) :: Gsph_epsG0
2400 : : type(kmesh_t),intent(in) :: Qmesh
2401 : : type(vcoul_t),intent(in) :: Vcp
2402 : : type(crystal_t),intent(in) :: Cryst
2403 : : type(Datafiles_type),intent(in) :: Dtfil
2404 : : complex(gwp),intent(in) :: epsm1(Ep%npwe,Ep%npwe,Ep%nomega)
2405 : : real(dp),intent(inout) :: rspower(npvel)
2406 : :
2407 : : !Local variables ------------------------------
2408 : : integer :: ipvel,ig, iq_bz,iq_ibz,isym_q,itim_q
2409 : : integer :: iomega,iomegap,nomega_re,unt_rsp, units(2)
2410 : 12 : integer,allocatable :: iomega_re(:)
2411 : :
2412 : : real(dp),parameter :: zp=1.0_dp ! Hard-coded charge of the impinging particle
2413 : : real(dp) :: omega_p
2414 : : real(dp) :: im_epsm1_int(1), qbz(3),qpgcart(3),qpg_red(3)
2415 : 24 : real(dp) :: pvel(3,npvel),pvel_norm(npvel),ypp_i(Ep%nomega),vcoul(Ep%npwe)
2416 : 12 : real(dp),allocatable :: im_epsm1_diag_qbz(:,:),tmp_data(:),omega_re(:)
2417 : : character(len=500) :: msg
2418 : : character(len=fnlen+4) :: fname
2419 : : !************************************************************************
2420 : :
2421 [ + + ]: 36 : units = [std_out, ab_out]
2422 : :
2423 : : ! First set up the velocities array from the input variables npvel and pvelmax(3)
2424 : : ! Remember pvelmax is in Cartesian coordinates and so is pvel
2425 [ + + ]: 492 : do ipvel=1,npvel
2426 [ + + ]: 1920 : pvel(:,ipvel) = REAL(ipvel,dp) / REAL(npvel,dp) * pvelmax(:)
2427 [ + + ]: 1932 : pvel_norm(ipvel) = SQRT( SUM( pvel(:,ipvel)**2 ) )
2428 : : enddo
2429 : : !
2430 : : ! Select the purely real frequency in Ep%omega
2431 : 12 : nomega_re=0
2432 [ + + ]: 492 : do iomega=1,Ep%nomega
2433 [ + - ]: 492 : if( AIMAG(Ep%omega(iomega)) < 1.0e-4_dp ) then
2434 : 480 : nomega_re=nomega_re+1
2435 : : endif
2436 : : enddo
2437 [ - + + - : 36 : ABI_MALLOC(omega_re,(nomega_re))
+ - + - -
+ - + ]
2438 [ - + + - : 36 : ABI_MALLOC(iomega_re,(nomega_re))
+ - + - -
+ - + ]
2439 [ - + + - : 48 : ABI_MALLOC(im_epsm1_diag_qbz,(Ep%npwe,Ep%nomega))
- + + - +
- + - - +
- + ]
2440 [ - + + - : 36 : ABI_MALLOC(tmp_data,(Ep%nomega))
+ - + - -
+ - + ]
2441 : :
2442 : 492 : iomegap=0
2443 [ + + ]: 492 : do iomega=1,Ep%nomega
2444 [ + - ]: 492 : if( AIMAG(Ep%omega(iomega)) < 1.0e-4_dp ) then
2445 : 480 : iomegap=iomegap+1
2446 : 480 : iomega_re(iomegap)=iomega
2447 : 480 : omega_re(iomegap)=REAL(Ep%omega(iomega),dp)
2448 : : endif
2449 : : enddo
2450 : :
2451 : : ! Loop over all the q-points in the full Brillouin zone and select only the
2452 : : ! ones that corresponds to the correct q-point in the irreducible wedge we are
2453 : : ! currently treating (index iqibz)
2454 [ + + ]: 396 : do iq_bz=1,Qmesh%nbz
2455 : :
2456 : : ! Perform the check and obtain the symmetry information
2457 : 384 : call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q)
2458 [ + + ]: 384 : if( iqibz /= iq_ibz ) cycle
2459 : :
2460 : : ! Apply the symmetry operation to the diagonal of epsm1
2461 [ + + ]: 2624 : do iomega=1,nomega_re
2462 [ + + ]: 41024 : do ig=1,Ep%npwe
2463 : 40960 : im_epsm1_diag_qbz(Gsph_epsG0%rottb(ig,itim_q,isym_q),iomega)= AIMAG( epsm1(ig,ig,iomega_re(iomega)) )
2464 : : enddo
2465 : : enddo
2466 : : ! Apply the symmetry operation to the Coulomb interaction
2467 [ + + ]: 1024 : do ig=1,Ep%npwe
2468 : 1024 : vcoul(Gsph_epsG0%rottb(ig,itim_q,isym_q))=Vcp%vc_sqrt(ig,iqibz)**2
2469 : : enddo
2470 : :
2471 : : ! Sum over G vectors
2472 [ + + ]: 1036 : do ig=1,Ep%npwe
2473 : : !
2474 : : ! Loop over velocities
2475 [ + + ]: 39744 : do ipvel=1,npvel
2476 : :
2477 [ + + ]: 153600 : qpg_red(:) = qbz(:) + Gsph_epsG0%gvec(:,ig)
2478 : : ! Transform q + G from reduced to cartesian with the symmetry operation
2479 : : qpgcart(:) = two_pi * Cryst%gprimd(:,1) * qpg_red(1) &
2480 : : + two_pi * Cryst%gprimd(:,2) * qpg_red(2) &
2481 [ + + ]: 153600 : + two_pi * Cryst%gprimd(:,3) * qpg_red(3)
2482 : :
2483 : : ! omega_p = ( q + G ) . v
2484 [ + + ]: 153600 : omega_p = DOT_PRODUCT( qpgcart(:) , pvel(:,ipvel) )
2485 : :
2486 : : ! Check that the calculated frequency omega_p is within the omega
2487 : : ! range of epsm1 and thus that the interpolation will go fine
2488 [ - + ]: 38400 : if ( ABS(omega_p) > omega_re(nomega_re) ) then
2489 : 0 : write(msg,'(a,e16.4,2a,e16.4)') ' freqremax is currently ',omega_re(nomega_re),ch10,&
2490 : 0 : ' increase it to at least ',omega_p
2491 : 0 : ABI_WARNING(msg)
2492 : : endif
2493 : :
2494 : : ! Perform the spline interpolation to obtain epsm1 at the desired omega = omega_p
2495 [ - + + + ]: 1612800 : tmp_data = im_epsm1_diag_qbz(ig,:)
2496 : :
2497 : 38400 : call spline( omega_re, tmp_data, nomega_re, 1.0e+32_dp, 1.0e+32_dp, ypp_i)
2498 [ + + ]: 76800 : call splint( nomega_re, omega_re, tmp_data, ypp_i, 1, (/ ABS(omega_p) /), im_epsm1_int )
2499 : :
2500 : : ! Apply the odd parity of Im epsm1 in omega to recover the causal response function
2501 [ + + ]: 38400 : if (omega_p<zero) im_epsm1_int(1)=-im_epsm1_int(1)
2502 : :
2503 : : ! Calculate 4 * pi / |q+G|**2 * omega_p * Im{ epsm1_GG(q,omega_p) }
2504 : 38400 : im_epsm1_int(1) = omega_p * vcoul(ig) * im_epsm1_int(1)
2505 : :
2506 : : ! Accumulate the final result without the prefactor
2507 : : ! (It will be included at the very end)
2508 : 39360 : rspower(ipvel) = rspower(ipvel) + im_epsm1_int(1)
2509 : :
2510 : : end do ! end of velocity loop
2511 : : end do ! end G-loop
2512 : :
2513 : : enddo ! end of q loop in the full BZ
2514 : :
2515 : : ! If it is the last q, write down the result in the main output file and in a
2516 : : ! separate file _RSP (for Random Stopping Power)
2517 [ + + ]: 12 : if (iqibz == Qmesh%nibz ) then
2518 : :
2519 : : ! Multiply by the prefactors
2520 : : ! Note that this expression differs from Eq. (3.11) in Campillo PRB 58, 10307 (1998) [[cite:Campillo1998]].
2521 : : ! A factor one half is missing in the paper.
2522 [ + + ]: 82 : rspower(:) = - zp**2 / ( Cryst%ucvol * Qmesh%nbz * pvel_norm(:) ) * rspower(:)
2523 : :
2524 : 2 : write(msg,'(2a)') ch10,' ==== Random stopping power along Cartesian direction === '
2525 : 2 : call wrtout(units, msg)
2526 : 2 : write(msg,'(a,3(f12.4,2x),a)') ' ==== ',pvelmax(:),'===='
2527 : 2 : call wrtout(units, msg)
2528 : 2 : write(msg,'(a)') '# |v| (a.u.) , RSP (a.u.) '
2529 : 2 : call wrtout(units, msg)
2530 [ + + ]: 82 : do ipvel=1,npvel
2531 : 80 : write(msg,'(f16.8,4x,f16.8)') pvel_norm(ipvel),rspower(ipvel)
2532 : 82 : call wrtout(units, msg)
2533 : : enddo
2534 : 2 : write(msg,'(2a)') ' ========================================================= ',ch10
2535 : 2 : call wrtout(units, msg)
2536 : :
2537 [ + - + - ]: 2 : fname=TRIM(Dtfil%filnam_ds(4))//'_RSP'
2538 [ - + ]: 2 : if (open_file(fname,msg,newunit=unt_rsp,status='unknown',form='formatted') /= 0) then
2539 : 0 : ABI_ERROR(msg)
2540 : : end if
2541 : :
2542 : 2 : write(msg,'(a)') '# ==== Random stopping power along Cartesian direction === '
2543 : 2 : call wrtout(unt_rsp, msg)
2544 : 2 : write(msg,'(a,3(f12.4,2x))') '# ==== ',pvelmax(:)
2545 : 2 : call wrtout(unt_rsp, msg)
2546 : 2 : write(msg,'(a)') '# |v| (a.u.) , RSP (a.u.) '
2547 : 2 : call wrtout(unt_rsp, msg)
2548 [ + + ]: 82 : do ipvel=1,npvel
2549 : 80 : write(msg,'(f16.8,4x,f16.8)') pvel_norm(ipvel),rspower(ipvel)
2550 : 82 : call wrtout(unt_rsp,msg)
2551 : : enddo
2552 : 2 : close(unt_rsp)
2553 : : end if
2554 : :
2555 [ - + ]: 12 : ABI_FREE(omega_re)
2556 [ - + ]: 12 : ABI_FREE(iomega_re)
2557 [ - + ]: 12 : ABI_FREE(im_epsm1_diag_qbz)
2558 [ - + ]: 12 : ABI_FREE(tmp_data)
2559 : :
2560 : 12 : end subroutine random_stopping_power
2561 : : !!***
2562 : :
2563 : : !!****f* m_screening_driver/calc_rpa_functional
2564 : : !! NAME
2565 : : !! calc_rpa_functional
2566 : : !!
2567 : : !! FUNCTION
2568 : : !! Routine used to calculate the Galitskii-Migdal and RPA approximations to the correlation energy
2569 : : !! from the irreducible polarizability.
2570 : : !!
2571 : : !! INPUTS
2572 : : !! iq=index of the q-point in the array Qmesh%ibz where epsilon^-1 has to be calculated
2573 : : !! Ep<em1params_t>=Structure with parameters and dimensions related to the inverse dielectric matrix.
2574 : : !! Pvc<vcoul_t>=Structure gathering data on the Coulombian interaction
2575 : : !! Qmesh<kmesh_t>=Data type with information on the q-sampling
2576 : : !! Dtfil<Datafiles_type)>=variables related to files
2577 : : !! comm=MPI communicator.
2578 : : !!
2579 : : !! OUTPUT
2580 : : !!
2581 : : !! SOURCE
2582 : :
2583 : 18 : subroutine calc_rpa_functional(gwrpacorr,gwgmcorr,iqcalc,iq,Ep,Pvc,Qmesh,Dtfil,gmet,chi0,comm,ec_rpa,ec_gm)
2584 : :
2585 : : use m_hide_lapack, only : xginv, xheev
2586 : :
2587 : : !Arguments ------------------------------------
2588 : : !scalars
2589 : : integer,intent(in) :: iqcalc,iq,gwrpacorr,gwgmcorr,comm
2590 : : real(dp),intent(inout) :: ec_gm
2591 : : type(kmesh_t),intent(in) :: Qmesh
2592 : : type(vcoul_t),intent(in) :: Pvc
2593 : : type(Datafiles_type),intent(in) :: Dtfil
2594 : : type(em1params_t),intent(in) :: Ep
2595 : : !arrays
2596 : : real(dp),intent(in) :: gmet(3,3)
2597 : : real(dp),intent(inout) :: ec_rpa(gwrpacorr)
2598 : : complex(gwp),intent(inout) :: chi0(Ep%npwe,Ep%npwe,Ep%nomega)
2599 : :
2600 : : !Local variables-------------------------------
2601 : : !scalars
2602 : : integer,parameter :: master = 0
2603 : : integer :: ig1,ig2,ilambda,io,rank,nprocs,unt,ierr, units(2)
2604 : : real(dp) :: ecorr,ecorr_gm,lambda
2605 : : logical :: q_is_gamma
2606 : : character(len=500) :: msg
2607 : : !arrays
2608 : 18 : real(dp),allocatable :: z(:),zl(:),zlw(:),zw(:)
2609 : 18 : complex(gwp),allocatable :: chi0_diag(:),chitmp(:,:),chi0_diag_gm(:),chitmp_gm(:,:)
2610 : 18 : real(gwp),allocatable :: eig(:)
2611 : : ! *************************************************************************
2612 : :
2613 [ + + ]: 54 : units = [std_out, ab_out]
2614 : :
2615 : : ! initialize MPI data
2616 : 18 : rank = xmpi_comm_rank(comm)
2617 : 18 : nprocs = xmpi_comm_size(comm)
2618 : :
2619 [ - + ]: 18 : ABI_CHECK(ep%iw_mesh_type == "gauss_legendre", "only gauss legendre mesh is supported")
2620 : :
2621 : : !if (rank==master) then ! presently only master has chi0 in screening
2622 : :
2623 : : ! vc_sqrt contains vc^{1/2}(q,G), complex-valued to allow for a possible cutoff
2624 : 18 : q_is_gamma = normv(Qmesh%ibz(:,iq),gmet,'G')<GW_TOLQ0
2625 : :
2626 : : ! Calculate Gauss-Legendre quadrature knots and weights for the omega integration
2627 [ - + + - : 54 : ABI_MALLOC(zw, (Ep%nomegaei))
+ - + - -
+ - + ]
2628 [ - + + - : 36 : ABI_MALLOC(z, (Ep%nomegaei))
- + - + ]
2629 : 18 : call coeffs_gausslegint(zero, one, z, zw, Ep%nomegaei)
2630 : :
2631 : : ! Calculate Gauss-Legendre quadrature knots and weights for the lambda integration
2632 [ - + + - : 54 : ABI_MALLOC(zlw, (gwrpacorr))
+ - + - -
+ - + ]
2633 [ - + + - : 36 : ABI_MALLOC(zl, (gwrpacorr))
- + - + ]
2634 : 18 : call coeffs_gausslegint(zero,one,zl,zlw,gwrpacorr)
2635 : :
2636 [ - + + - : 54 : ABI_MALLOC(chi0_diag,(Ep%npwe))
+ - + - -
+ - + ]
2637 [ - + + - : 72 : ABI_MALLOC_OR_DIE(chitmp,(Ep%npwe,Ep%npwe), ierr)
- + + - +
- + - + -
- + - + ]
2638 [ + + ]: 18 : if(gwgmcorr==1) then
2639 [ - + + - : 12 : ABI_MALLOC(chi0_diag_gm,(Ep%npwe))
- + - + ]
2640 [ - + + - : 18 : ABI_MALLOC_OR_DIE(chitmp_gm,(Ep%npwe,Ep%npwe), ierr)
- + + - +
- - + -
+ ]
2641 : : end if
2642 : :
2643 [ + + ]: 810 : do io=2,Ep%nomega
2644 : : !if (q_is_gamma) then
2645 : : ! call wrtout(units, "RPA: Ignoring q==0"); cycle
2646 : : !end if
2647 : :
2648 [ + + ]: 810 : if(gwrpacorr==1) then ! exact integration over the coupling constant
2649 : :
2650 [ + - ]: 774 : if(modulo(io-2,nprocs)/=rank) cycle ! distributing the workload
2651 : :
2652 [ + + ]: 18972 : do ig2=1,Ep%npwe
2653 [ + + ]: 1194498 : do ig1=1,Ep%npwe
2654 : 1193724 : chitmp(ig1,ig2) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig2,iq) * chi0(ig1,ig2,io)
2655 : : end do !ig1
2656 : : end do !ig2
2657 : :
2658 [ - + + - : 2322 : ABI_MALLOC(eig,(Ep%npwe))
+ - + - -
+ - + ]
2659 : 774 : call xheev('N','U',Ep%npwe,chitmp,eig)
2660 : :
2661 [ + + ]: 18972 : do ig1=1,Ep%npwe
2662 : : ec_rpa(:) = ec_rpa(:) &
2663 : : & - zw(io-1) / ( z(io-1) * z(io-1) ) &
2664 [ + + ]: 36396 : & * Qmesh%wt(iq) * (-log( 1.0_dp-eig(ig1) ) - eig(ig1) ) / (2.0_dp * pi )
2665 : : ec_gm = ec_gm &
2666 : : & - zw(io-1) / ( z(io-1) * z(io-1) ) &
2667 : 18972 : & * Qmesh%wt(iq) * ( eig(ig1) / ( 1.0_dp-eig(ig1) ) - eig(ig1) ) / (2.0_dp * pi )
2668 : : end do
2669 [ - + ]: 774 : ABI_FREE(eig)
2670 : :
2671 : : else ! numerical integration over the coupling constant
2672 : :
2673 : : !if(modulo( (ilambda-1)+gwrpacorr*(io-2),nprocs)/=rank) cycle ! distributing the workload
2674 : :
2675 [ + + ]: 162 : do ilambda=1,gwrpacorr
2676 [ + - ]: 144 : if(modulo( (ilambda-1)+gwrpacorr*(io-2),nprocs)/=rank) cycle ! distributing the workload
2677 : 144 : lambda=zl(ilambda)
2678 [ + + ]: 19872 : do ig1=1,Ep%npwe
2679 : 19872 : chi0_diag(ig1) = Pvc%vc_sqrt(ig1,iq)**2 * chi0(ig1,ig1,io)
2680 : : end do
2681 : :
2682 [ + + + - ]: 144 : if(ilambda==1 .and. gwgmcorr==1) then ! Copy v^1/2*Chi0*v^1/2 for Galitskii-Migdal
2683 [ # # ]: 0 : chi0_diag_gm(:) = chi0_diag(:)
2684 : : end if
2685 : :
2686 [ + + ]: 19872 : do ig2=1,Ep%npwe
2687 [ + + ]: 2722464 : do ig1=1,Ep%npwe
2688 : 2702736 : chitmp(ig1,ig2) = - lambda * Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chi0(ig1,ig2,io)
2689 : :
2690 [ + + - + ]: 2722464 : if(ilambda==1 .and. gwgmcorr==1) then ! Use lambda=1 for Galitskii-Migdal
2691 : 0 : chitmp_gm(ig1,ig2) = - Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chi0(ig1,ig2,io)
2692 : : end if
2693 : :
2694 : : end do !ig1
2695 : 19728 : chitmp(ig2,ig2) = chitmp(ig2,ig2) + 1.0_dp
2696 : :
2697 [ + + - + ]: 19872 : if(ilambda==1 .and. gwgmcorr==1) then ! Prepare (1-v^1/2*Chi0*v^1/2) for Galitskii-Migdal
2698 : 0 : chitmp_gm(ig2,ig2) = chitmp_gm(ig2,ig2) + 1.0_dp
2699 : : end if
2700 : :
2701 : : end do !ig2
2702 : 144 : call xginv(chitmp(:,:),Ep%npwe)
2703 [ + + + + : 2722752 : chitmp(:,:) = matmul( chi0(:,:,io) , chitmp(:,:) )
- + - - -
- - - - -
- - + - ]
2704 : :
2705 [ + + - + ]: 144 : if(ilambda==1 .and. gwgmcorr==1) then ! Prepare Chi = [(1-v^1/2*Chi0*v^1/2)]^-1 * Chi0 for Galitskii-Migdal
2706 : 0 : call xginv(chitmp_gm(:,:),Ep%npwe)
2707 [ # # # # : 0 : chitmp_gm(:,:) = matmul( chi0(:,:,io) , chitmp_gm(:,:) )
# # # # #
# # # # #
# # # # ]
2708 : : end if
2709 : :
2710 [ + + ]: 19872 : do ig1=1,Ep%npwe
2711 : 19728 : chi0_diag(ig1) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chitmp(ig1,ig1) - chi0_diag(ig1)
2712 : :
2713 [ + + - + ]: 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
2714 : 0 : chi0_diag_gm(ig1) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chitmp_gm(ig1,ig1) - chi0_diag_gm(ig1)
2715 : : end if
2716 : :
2717 : : end do
2718 : :
2719 [ + + ]: 19890 : do ig1=1,Ep%npwe
2720 : : ec_rpa(ilambda) = ec_rpa(ilambda) &
2721 : 19728 : & - zw(io-1) / ( z(io-1) * z(io-1) ) * Qmesh%wt(iq) * real( chi0_diag(ig1) ) / (2.0_dp * pi )
2722 : :
2723 [ + + - + ]: 19872 : if(ilambda==1 .and. gwgmcorr==1) then ! Integrate [v*Chi-v*Chi0](iw) dw for Galitskii-Migdal
2724 : : ec_gm = ec_gm &
2725 : 0 : & - zw(io-1) / ( z(io-1) * z(io-1) ) * Qmesh%wt(iq) * real( chi0_diag_gm(ig1) ) / (2.0_dp * pi )
2726 : : end if
2727 : :
2728 : : end do
2729 : :
2730 : : end do ! ilambda
2731 : :
2732 : : end if ! exact or numerical integration over the coupling constant
2733 : :
2734 : : end do ! io
2735 : :
2736 : :
2737 : : ! Output the correlation energy when the last q-point to be calculated is reached
2738 : : ! This would allow for a manual parallelization over q-points
2739 [ + + ]: 18 : if(iqcalc==Ep%nqcalc) then
2740 : :
2741 : 5 : call xmpi_sum_master(ec_rpa,master,comm,ierr)
2742 : 5 : call xmpi_sum_master(ec_gm,master,comm,ierr)
2743 : :
2744 [ + - ]: 5 : if(rank==master) then
2745 [ + + ]: 17 : ecorr = sum( zlw(:)*ec_rpa(:) )
2746 : 5 : ecorr_gm = ec_gm
2747 [ - + ]: 5 : if (open_file(dtfil%fnameabo_rpa, msg, newunit=unt) /=0) then
2748 : 0 : ABI_ERROR(msg)
2749 : : end if
2750 : 5 : write(unt,'(a,(2x,f14.8))') '#RPA',ecorr
2751 : 5 : write(msg,'(2a,(2x,f14.8))') ch10,' RPA energy [Ha] :',ecorr
2752 : 5 : call wrtout(units, msg)
2753 [ + + ]: 5 : if(gwrpacorr>1) then
2754 [ + + ]: 9 : do ilambda=1,gwrpacorr
2755 : 8 : write(unt,'(i6,2x,f10.6,2x,e13.6)') ilambda,zl(ilambda),ec_rpa(ilambda)
2756 : 8 : write(msg,'(i6,2x,f10.6,2x,e13.6)') ilambda,zl(ilambda),ec_rpa(ilambda)
2757 : 9 : call wrtout(units, msg)
2758 : : end do
2759 : : end if
2760 [ + + ]: 5 : if(gwgmcorr==1) then ! Only exact integration over the coupling constant
2761 : 1 : write(unt,'(a,(2x,f14.8))') '#GM',ecorr_gm
2762 : 1 : write(msg,'(2a,(2x,f14.8))') ch10,' Galitskii-Migdal energy [Ha] :',ecorr_gm
2763 : 1 : call wrtout(units, msg)
2764 : 1 : write(unt,'(a1)') ' '
2765 : 1 : write(msg,'(a1)') ' '
2766 : 1 : call wrtout(units, msg)
2767 : : end if
2768 : 5 : close(unt)
2769 : : end if
2770 : :
2771 : : end if
2772 : :
2773 [ + + ]: 18 : if(gwgmcorr==1) then
2774 [ - + ]: 6 : ABI_FREE(chitmp_gm)
2775 [ - + ]: 6 : ABI_FREE(chi0_diag_gm)
2776 : : end if
2777 [ - + ]: 18 : ABI_FREE(chi0_diag)
2778 [ - + ]: 18 : ABI_FREE(chitmp)
2779 [ - + ]: 18 : ABI_FREE(zl)
2780 [ - + ]: 18 : ABI_FREE(zlw)
2781 [ - + ]: 18 : ABI_FREE(z)
2782 [ - + ]: 18 : ABI_FREE(zw)
2783 : :
2784 [ - + - + ]: 18 : end subroutine calc_rpa_functional
2785 : : !!***
2786 : :
2787 : 6976 : end module m_screening_driver
2788 : : !!***
|
