Line data Source code
1 : !!****m* ABINIT/m_positron
2 : !! NAME
3 : !! m_positron
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (GJ, MT, JW)
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_positron
23 :
24 : use defs_basis
25 : use m_efield
26 : use m_errors
27 : use m_abicore
28 : use m_energies
29 : use m_wffile
30 : use m_electronpositron
31 : use m_hdr
32 : use m_xmpi
33 : use m_bandfft_kpt
34 : use m_dtset
35 : use m_dtfil
36 : use m_extfpmd
37 :
38 : use defs_datatypes, only : pseudopotential_type
39 : use defs_abitypes, only : MPI_type
40 : use m_special_funcs, only : sbf8
41 : use m_ioarr, only : ioarr, read_rhor
42 : use m_paw_atomorb, only : atomorb_type,destroy_atomorb
43 : use m_pawang, only : pawang_type
44 : use m_paw_sphharm, only : realgaunt
45 : use m_pawrad, only : pawrad_type, simp_gen, nderiv_gen
46 : use m_pawtab, only : pawtab_type
47 : use m_paw_ij, only : paw_ij_type
48 : use m_pawfgrtab,only : pawfgrtab_type
49 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_alloc, pawrhoij_free,&
50 : pawrhoij_nullify, pawrhoij_gather, pawrhoij_inquire_dim, pawrhoij_symrhoij
51 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_mpi_send, &
52 : pawcprj_mpi_recv, pawcprj_free, pawcprj_copy, pawcprj_bcast
53 : use m_pawfgr, only : pawfgr_type
54 : use m_paw_nhat, only : pawmknhat
55 : use m_fock, only : fock_type
56 : use m_kg, only : getcut
57 : use defs_wvltypes, only : wvl_data
58 : use m_spacepar, only : hartre
59 : use m_mkrho, only : initro
60 : use m_paw_occupancies, only : initrhoij, pawaccrhoij
61 : use m_gammapositron, only : gammapositron, gammapositron_fft
62 : use m_forstr, only : forstr
63 : use m_pawxc, only : pawxcsum
64 : use m_paw_denpot, only : pawdensities
65 : use m_drivexc, only : mkdenpos
66 :
67 : use m_paw_sphharm, only : initylmr
68 : use m_pawpsp, only : pawpsp_init_core
69 : use m_crystal, only : crystal_t
70 : use m_mpinfo, only : ptabs_fourdp,set_mpi_enreg_fft,unset_mpi_enreg_fft,destroy_mpi_enreg, initmpi_seq, proc_distrb_cycle
71 : use m_io_tools, only : open_file,close_unit,get_unit
72 : use m_fftcore, only : sphereboundary
73 : use m_prep_kgb, only : prep_fourwf
74 : use m_fft, only : fourwf, fourdp
75 : use m_cgprj, only : ctocprj
76 : use m_xg_nonlop, only : xg_nonlop_t
77 :
78 : implicit none
79 :
80 : private
81 : !!***
82 :
83 : public :: setup_positron
84 : public :: poslifetime
85 : public :: posdoppler
86 : !!***
87 :
88 : contains
89 : !!***
90 :
91 : !!****f* ABINIT/setup_positron
92 : !! NAME
93 : !! setup_positron
94 : !!
95 : !! FUNCTION
96 : !! Do various initializations for the positron lifetime calculation
97 : !!
98 : !! INPUTS
99 : !! atindx(natom)=index table for atoms (see gstate.f)
100 : !! atindx1(natom)=index table for atoms, inverse of atindx
101 : !! dtefield <type(efield_type)> = variables related to Berry phase
102 : !! dtfil <type(datafiles_type)>=variables related to files
103 : !! dtset <type(dataset_type)>=all input variables for this dataset
104 : !! ecore=core psp energy (part of total energy) (hartree)
105 : !! etotal=current value of total energy
106 : !! extfpmd <type(extfpmd_type)>=extended first-principles molecular dynamics type
107 : !! fock <type(fock_type)>= quantities to calculate Fock exact exchange
108 : !! forces_needed=if >0 forces are needed
109 : !! gred(3,natom)=gradients wrt nuclear positions in reduced coordinates
110 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space
111 : !! gmet(3,3)=reciprocal space metric
112 : !! grchempottn(3,natom)=d(E_chemical_potential)/d(xred) (hartree)
113 : !! grcondft(3,natom)=d(E_constrainedDFT)/d(xred) (hartree)
114 : !! grewtn(3,natom)=d(Ewald)/d(xred) (hartree)
115 : !! grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D dispersion (hartree)
116 : !! gsqcut=cutoff value on G**2 for sphere inside fft box
117 : !! hdr <type(hdr_type)>=the header of wf, den and pot files
118 : !! ifirst_gs= 0 if we are in a single ground-state calculation
119 : !! or in the first ground-state calculation of a structural minimization/dynamics
120 : !! indsym(4,nsym,natom)=index showing transformation of atom labels
121 : !! under symmetry operations (computed in symatm)
122 : !! istep=index of the number of steps in the routine scfcv
123 : !! istep_mix=index of the number of steps for the SCF mixing (can be <istep)
124 : !! kg(3,mpw*mkmem)=reduced (integer) coordinates of G vecs in basis sphere
125 : !! kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0
126 : !! maxfor=maximum absolute value of fcart (forces)
127 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
128 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
129 : !! mgfft=maximum size of 1D FFTs
130 : !! mpi_enreg=information about MPI parallelization
131 : !! my_natom=number of atoms treated by current processor
132 : !! n3xccc=dimension of the xccc3d array (0 or nfftf).
133 : !! nattyp(ntypat)= # atoms of each type.
134 : !! nfft=(effective) number of FFT grid points (for this processor)
135 : !! ngfft(18)=contain all needed information about 3D FFT
136 : !! ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc
137 : !! nhat(nfftf,nspden*usepaw)= -PAW only- compensation density
138 : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
139 : !! npwarr(nkpt)=number of planewaves in basis and on boundary for each k
140 : !! nvresid(nfftf,nspden)=array for the residual of the density/potential
141 : !! optres=0 if the potential residual has to be used for forces corrections
142 : !! =1 if the density residual has to be used for forces corrections
143 : !! paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
144 : !! pawang <type(pawang_type)>=paw angular mesh and related data
145 : !! pawfgr(my_natom*usepaw) <type(pawfgr_type)>=fine grid parameters and related data
146 : !! pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
147 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
148 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
149 : !! ph1d(2,3*(2*mgfftf+1)*natom)=one-dimensional structure factor information (fine FFT grid)
150 : !! ph1dc(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases (coarse FFT grid)
151 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
152 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
153 : !! stress_needed=if >0 stresses are needed
154 : !! strscondft(6)=cDFT correction to stress
155 : !! strsxc(6)=xc correction to stress
156 : !! symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates
157 : !! ucvol=unit cell volume in bohr**3.
158 : !! usecprj= 1 if cprj array is stored in memory
159 : !! usevxctau=1 if kinetic energy density contribution has to be included (mGGA)
160 : !! vhartr(nfftf)=array for holding Hartree potential
161 : !! vpsp(nfftf)=array for holding local psp
162 : !! vxc(nfftf,nspden)=exchange-correlation potential (hartree) in real space
163 : !! vxctau(nfft,nspden,4*usevxctau)=(only for meta-GGA) derivative of XC energy density
164 : !! wrt kinetic energy density (depsxcdtau)
165 : !! xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
166 : !! xcctau3d(n3xccc*usekden)=(only for meta-GGA): 3D core electron kinetic energy density for XC core correction
167 : !! xred(3,natom)=reduced dimensionless atomic coordinates
168 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
169 : !! ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
170 : !!
171 : !! SIDE EFFECTS
172 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
173 : !! energies <type(energies_type)>=all part of total energy.
174 : !! cg(2,mcg)=wavefunctions
175 : !! cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors:
176 : !! cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector.
177 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
178 : !! occ(mband*nkpt*nsppol)=occupation number for each band at each k point
179 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
180 : !! rhog(2,nfft)=Fourier transform of total electron/positron density
181 : !! rhor(nfft,nspden)=total electron/positron density (el/bohr**3)
182 : !!
183 : !! SOURCE
184 :
185 525 : subroutine setup_positron(atindx,atindx1,cg,cprj,dtefield,dtfil,dtset,ecore,eigen,etotal,electronpositron,&
186 525 : & energies,fock,forces_needed,gred,gmet,gprimd,grchempottn,&
187 525 : & grcondft,grewtn,grvdw,gsqcut,hdr,extfpmd,ifirst_gs,indsym,istep,istep_mix,kg,&
188 525 : & kxc,maxfor,mcg,mcprj,mgfft,mpi_enreg,my_natom,n3xccc,nattyp,nfft,ngfft,ngrvdw,nhat,nkxc,npwarr,nvresid,occ,optres,&
189 525 : & paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1dc,psps,rhog,rhor,&
190 525 : & rmet,rprimd,stress_needed,strscondft,strsxc,symrec,ucvol,usecprj,usevxctau,vhartr,vpsp,vxc,vxctau,&
191 525 : & xccc3d,xcctau3d,xred,ylm,ylmgr,xg_nonlop)
192 :
193 : !Arguments ------------------------------------
194 : !scalars
195 : integer,intent(in) :: forces_needed,ifirst_gs,istep,mcg,mcprj,mgfft,my_natom,n3xccc,nfft
196 : integer,intent(in) :: ngrvdw,nkxc,optres,stress_needed,usecprj,usevxctau
197 : integer,intent(inout) :: istep_mix
198 : real(dp),intent(in) :: ecore,etotal,gsqcut,maxfor,ucvol
199 : type(efield_type),intent(in) :: dtefield
200 : type(datafiles_type),intent(in) :: dtfil
201 : type(dataset_type),intent(in) :: dtset
202 : type(electronpositron_type),pointer :: electronpositron
203 : type(energies_type),intent(inout) :: energies
204 : type(hdr_type),intent(inout) :: hdr
205 : type(extfpmd_type),pointer,intent(inout) :: extfpmd
206 : type(MPI_type),intent(inout) :: mpi_enreg
207 : type(pawang_type),intent(in) :: pawang
208 : type(pawfgr_type),intent(in) :: pawfgr
209 : type(pseudopotential_type), intent(in) :: psps
210 : type(fock_type),pointer, intent(inout) :: fock
211 : type(xg_nonlop_t), intent(inout) :: xg_nonlop
212 : !arrays
213 : integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom)
214 : integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),nattyp(dtset%natom),ngfft(18)
215 : integer,intent(in) :: npwarr(dtset%nkpt),symrec(3,3,dtset%nsym)
216 : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),grchempottn(3,dtset%natom),grcondft(3,dtset%natom)
217 : real(dp),intent(in) :: grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc)
218 : real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),ph1dc(2,(3*(2*dtset%mgfft+1)*dtset%natom)*dtset%usepaw)
219 : real(dp),intent(in) :: rmet(3,3),strscondft(6),strsxc(6),vhartr(nfft),vpsp(nfft),vxc(nfft,dtset%nspden)
220 : real(dp),intent(in) :: vxctau(nfft,dtset%nspden,4*usevxctau)
221 : real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
222 : real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
223 : real(dp),intent(inout) :: cg(2,mcg)
224 : real(dp),intent(inout) :: nhat(nfft,dtset%nspden*dtset%usepaw)
225 : real(dp),intent(inout) :: nvresid(nfft,dtset%nspden)
226 : real(dp),intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol),gred(3,dtset%natom)
227 : real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
228 : real(dp),intent(inout) :: rhog(2,nfft),rhor(nfft,dtset%nspden),rprimd(3,3)
229 : real(dp),intent(inout) :: xccc3d(n3xccc),xcctau3d(n3xccc*dtset%usekden),xred(3,dtset%natom)
230 : type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj*usecprj)
231 : type(paw_ij_type),intent(in) :: paw_ij(my_natom*dtset%usepaw)
232 : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*dtset%usepaw)
233 : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw)
234 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
235 : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*dtset%usepaw)
236 :
237 : !Local variables-------------------------------
238 : !scalars
239 : integer,parameter :: cplex1=1
240 : integer :: history_level,iatom,iband,icalctype,icalctype0,icg,ifft,ikpt
241 : integer :: iocc,ireadwf,ispden,isppol,n3xccc0,nocc,optfor,optstr,rdwrpaw,comm_cell
242 : logical,parameter :: always_restart=.false. ! Set to true to restart by a pure electronic step at each new atomic structure
243 : logical :: need_scocc,new_calctype
244 : real(dp) :: boxcut_dum,diffor_dum,ecut_eff,eigtmp,etotal_read,gsqcut_eff,maxfor_dum
245 : real(dp) :: maxocc,nelect,occlast,occtmp,rhotmp
246 : character(len=69) :: TypeCalcStrg
247 : character(len=500) :: message
248 : character(len=fnlen) :: fname
249 : type(energies_type) :: energies_tmp
250 525 : type(wvl_data) :: wvl
251 525 : type(hdr_type) :: hdr_den
252 : !arrays
253 525 : integer,allocatable :: nlmn(:)
254 : real(dp) :: cgtmp(2)
255 : real(dp),parameter :: qphon(3)=(/zero,zero,zero/)
256 525 : real(dp),allocatable :: favg_dum(:),fcart_dum(:,:),forold_dum(:,:),gred_tmp(:,:)
257 525 : real(dp),allocatable :: gresid_dum(:,:),grhf_dum(:,:),grxc_dum(:,:)
258 525 : real(dp),allocatable :: rhog_ep(:,:),scocc(:),str_tmp(:),synlgr_dum(:,:)
259 : real(dp) :: nhatgr(0,0,0)
260 525 : type(pawcprj_type),allocatable :: cprj_tmp(:,:)
261 525 : type(pawrhoij_type),allocatable :: pawrhoij_tmp(:)
262 :
263 : ! *************************************************************************
264 :
265 : !Compatibility tests
266 525 : if (dtset%positron==0) then
267 0 : ABI_BUG('Not valid for dtset%positron=0!')
268 : end if
269 :
270 525 : if (istep>1.and.nfft/=electronpositron%nfft) then
271 0 : ABI_BUG('Invalid value for nfft!')
272 : end if
273 :
274 525 : if (dtset%usewvl==1) then
275 0 : ABI_BUG('Not valid for wavelets!')
276 : end if
277 :
278 525 : if (dtset%positron==1) then
279 10 : do isppol=1,dtset%nsppol
280 15 : do ikpt=1,dtset%nkpt
281 10 : if (dtset%nband(ikpt+dtset%nkpt*(isppol-1))/=dtset%nband(1)) then
282 0 : message = "dtset%positron needs nband to be the same at each k-point !"
283 0 : ABI_ERROR(message)
284 : end if
285 : end do
286 : end do
287 : end if
288 :
289 525 : comm_cell = mpi_enreg%comm_cell
290 :
291 : !-----------------------------------------------------------------------
292 : !Compute new value for calctype (type of electron-positron calculation)
293 : !-----------------------------------------------------------------------
294 525 : icalctype0=electronpositron%calctype
295 :
296 525 : new_calctype=.false.
297 525 : if (dtset%positron==1.or.dtset%positron==2) then
298 6 : if (istep==1) new_calctype=.true.
299 519 : else if (dtset%positron<0) then
300 519 : if (ifirst_gs/=0.and.istep==1.and.(.not.always_restart)) new_calctype=.true.
301 525 : if (electronpositron%scf_converged) new_calctype=.true.
302 : end if
303 :
304 : !Comment:
305 : !history_level=-1: not used
306 : !history_level= 0: rhor from scratch, rhor_ep from scratch or read
307 : !history_level= 1: rhor in memory, rhor_ep from scratch or read
308 : !history_level= 2: rhor_ep <-rhor, rhor from scratch
309 : !history_level= 3: rhor_ep <-> rhor
310 : !history_level= 4: rhor in memory, rhor_ep in memory
311 525 : history_level=-1
312 525 : if (dtset%positron==1.or.dtset%positron==2) then
313 6 : if (ifirst_gs==0.and.istep==1) history_level=0
314 6 : if (ifirst_gs/=0.and.istep==1) history_level=4
315 519 : else if (dtset%positron<0) then
316 519 : if (.not.electronpositron%scf_converged) then
317 449 : if (ifirst_gs/=0.and.istep==1.and.(.not.always_restart)) history_level=4
318 : else if (electronpositron%scf_converged) then
319 70 : if (icalctype0==0) history_level=2
320 62 : if (icalctype0> 0) history_level=3
321 : end if
322 : end if
323 :
324 : electronpositron%calctype=icalctype0
325 525 : if (dtset%positron==1.or.dtset%positron==2) then
326 6 : electronpositron%calctype=dtset%positron
327 519 : else if (dtset%positron<0) then
328 519 : if (electronpositron%scf_converged) then
329 70 : if (icalctype0==0) electronpositron%calctype=1
330 70 : if (icalctype0>0 ) electronpositron%calctype=3-electronpositron%calctype
331 449 : else if (ifirst_gs/=0.and.istep==1) then
332 3 : if (always_restart) then
333 : electronpositron%calctype=0
334 : else
335 : electronpositron%calctype=2
336 : ! if (electronpositron%particle==EP_POSITRON) electronpositron%calctype=1
337 : end if
338 : end if
339 : end if
340 :
341 525 : electronpositron%scf_converged=.false.
342 525 : if (new_calctype) electronpositron%istep=electronpositron%istep+1
343 525 : if (istep==1) electronpositron%istep=1
344 525 : ireadwf=dtfil%ireadwf;if (electronpositron%istep>1) ireadwf=0
345 :
346 : !============================================
347 : !The following lines occur only when the type
348 : !of electron-positron calculation changes
349 : !============================================
350 525 : if (new_calctype) then
351 :
352 : ! Reset some indexes
353 79 : if (electronpositron%calctype==0) then
354 0 : electronpositron%particle=EP_NOTHING
355 79 : else if (electronpositron%calctype==1) then
356 43 : electronpositron%particle=EP_ELECTRON
357 36 : else if (electronpositron%calctype==2) then
358 36 : electronpositron%particle=EP_POSITRON
359 : end if
360 79 : electronpositron%has_pos_ham=mod(electronpositron%calctype,2)
361 79 : electronpositron%istep_scf=1
362 79 : istep_mix=1
363 :
364 : ! -----------------------------------------------------------------------------------------
365 : ! Update forces and stresses
366 : ! If electronpositron%calctype==1: gred_ep/stress_ep are the electronic gred/stress
367 : ! If electronpositron%calctype==2: gred_ep/stress_ep are the positronic gred/stress
368 : ! -----------------------------------------------------------------------------------------
369 79 : if (history_level==2.or.history_level==3) then
370 70 : optstr=0;optfor=0
371 70 : if (allocated(electronpositron%stress_ep)) optstr=stress_needed
372 70 : if (allocated(electronpositron%gred_ep).and.forces_needed==2) optfor=1
373 70 : if (optfor>0.or.optstr>0) then
374 0 : ABI_MALLOC(favg_dum,(3))
375 0 : ABI_MALLOC(fcart_dum,(3,dtset%natom))
376 0 : ABI_MALLOC(forold_dum,(3,dtset%natom))
377 0 : ABI_MALLOC(gresid_dum,(3,dtset%natom))
378 0 : ABI_MALLOC(grhf_dum,(3,dtset%natom))
379 0 : ABI_MALLOC(grxc_dum,(3,dtset%natom))
380 0 : ABI_MALLOC(synlgr_dum,(3,dtset%natom))
381 0 : ABI_MALLOC(gred_tmp,(3,dtset%natom))
382 0 : ABI_MALLOC(str_tmp,(6))
383 0 : forold_dum=zero;n3xccc0=n3xccc
384 0 : icalctype=electronpositron%calctype;electronpositron%calctype=-icalctype0
385 0 : if (electronpositron%calctype==0) electronpositron%calctype=-100
386 0 : if (electronpositron%calctype==-1) n3xccc0=0 ! Note: if calctype=-1, previous calculation was positron
387 : call forstr(atindx1,cg,cprj,diffor_dum,dtefield,dtset,eigen,electronpositron,energies,&
388 : & favg_dum,fcart_dum,fock,forold_dum,gred_tmp,grchempottn,grcondft,gresid_dum,grewtn,grhf_dum,grvdw,grxc_dum,gsqcut,&
389 : & extfpmd,indsym,kg,kxc,maxfor_dum,mcg,mcprj,mgfft,mpi_enreg,my_natom,n3xccc0,nattyp,nfft,ngfft,&
390 : & ngrvdw,nhat,nkxc,npwarr,dtset%ntypat,nvresid,occ,optfor,optres,paw_ij,pawang,pawfgr,&
391 : & pawfgrtab,pawrad,pawrhoij,pawtab,ph1dc,ph1d,psps,rhog,rhor,rprimd,optstr,strscondft,strsxc,str_tmp,symrec,&
392 0 : & synlgr_dum,ucvol,usecprj,usevxctau,vhartr,vpsp,vxc,vxctau,wvl,xccc3d,xcctau3d,xred,ylm,ylmgr,0.0_dp,xg_nonlop)
393 0 : electronpositron%calctype=icalctype
394 0 : if (optfor>0) electronpositron%gred_ep(:,:)=gred_tmp(:,:)
395 0 : if (optstr>0) electronpositron%stress_ep(:)=str_tmp(:)
396 0 : ABI_FREE(favg_dum)
397 0 : ABI_FREE(fcart_dum)
398 0 : ABI_FREE(forold_dum)
399 0 : ABI_FREE(gresid_dum)
400 0 : ABI_FREE(grhf_dum)
401 0 : ABI_FREE(grxc_dum)
402 0 : ABI_FREE(synlgr_dum)
403 0 : ABI_FREE(gred_tmp)
404 0 : ABI_FREE(str_tmp)
405 : end if
406 70 : if (optfor==0.and.forces_needed>0.and.allocated(electronpositron%gred_ep)) then
407 2369 : electronpositron%gred_ep(:,:)=gred(:,:)-electronpositron%gred_ep(:,:)
408 : end if
409 : end if
410 :
411 : ! ----------------------------------------------------------------------------------------------------
412 : ! Initialize/Update densities
413 : ! If electronpositron%calctype==1: rhor is the positronic density, rhor_ep is the electronic density
414 : ! If electronpositron%calctype==2: rhor is the electronic density, rhor_ep is the positronic density
415 : ! ---------------------------------------------------------------------------------------------------
416 237 : ABI_MALLOC(rhog_ep,(2,nfft))
417 :
418 : ! ===== PREVIOUS DENSITY RHOR_EP:
419 79 : if (history_level==0.or.history_level==1) then
420 : ! ----- Read from disk
421 6 : if (dtset%positron>0) then
422 6 : rdwrpaw=dtset%usepaw
423 6 : fname=trim(dtfil%fildensin);if (dtset%positron==2) fname=trim(dtfil%fildensin)//'_POSITRON'
424 : call read_rhor(trim(fname), cplex1, dtset%nspden, nfft, ngfft, rdwrpaw, mpi_enreg, electronpositron%rhor_ep, &
425 6 : hdr_den, electronpositron%pawrhoij_ep, comm_cell, check_hdr=hdr)
426 6 : etotal_read = hdr_den%etot; call hdr_den%free()
427 6 : call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,1,ngfft,0)
428 6 : if (dtset%usepaw==1.and.allocated(electronpositron%nhat_ep)) then
429 : call pawmknhat(occtmp,1,0,0,0,0,gprimd,my_natom,dtset%natom,nfft,ngfft,0,&
430 : & dtset%nspden,dtset%ntypat,pawang,pawfgrtab,nhatgr,electronpositron%nhat_ep,&
431 : & electronpositron%pawrhoij_ep,electronpositron%pawrhoij_ep,pawtab,&
432 : & qphon,rprimd,ucvol,dtset%usewvl,xred,distribfft=mpi_enreg%distribfft,&
433 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
434 5 : & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0)
435 : end if
436 : end if
437 : ! ----- Electronic from scratch
438 6 : if (dtset%positron<0.and.electronpositron%calctype==1) then
439 0 : ecut_eff=dtset%pawecutdg*(dtset%dilatmx)**2
440 0 : call getcut(boxcut_dum,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,std_out,qphon,ngfft)
441 : call initro(atindx,dtset%densty,gmet,gsqcut_eff,dtset%usepaw,mgfft,mpi_enreg,&
442 : & psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,dtset%ntypat,&
443 : & psps,pawtab,ph1d,psps%qgrid_vl,rhog_ep,electronpositron%rhor_ep,&
444 0 : & dtset%spinat,ucvol,dtset%usepaw,dtset%ziontypat,dtset%znucl)
445 0 : if (dtset%usepaw==1) then
446 0 : if (size(electronpositron%pawrhoij_ep)>0) then
447 0 : ABI_MALLOC(pawrhoij_tmp,(my_natom))
448 : call initrhoij(electronpositron%pawrhoij_ep(1)%cplex_rhoij,dtset%lexexch,&
449 : & dtset%lpawu,my_natom,dtset%natom,dtset%nspden,&
450 : & electronpositron%pawrhoij_ep(1)%nspinor,dtset%nsppol,&
451 : & dtset%ntypat,pawrhoij_tmp,dtset%pawspnorb,pawtab,cplex1,dtset%spinat,dtset%typat,&
452 : & ngrhoij=electronpositron%pawrhoij_ep(1)%ngrhoij,&
453 : & nlmnmix=electronpositron%pawrhoij_ep(1)%lmnmix_sz,&
454 : & use_rhoij_=electronpositron%pawrhoij_ep(1)%use_rhoij_,&
455 : & use_rhoijres=electronpositron%pawrhoij_ep(1)%use_rhoijres,&
456 0 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
457 0 : if (electronpositron%pawrhoij_ep(1)%lmnmix_sz>0) then
458 0 : do iatom=1,my_natom
459 0 : pawrhoij_tmp(iatom)%kpawmix(:)=electronpositron%pawrhoij_ep(iatom)%kpawmix(:)
460 : end do
461 : end if
462 0 : call pawrhoij_copy(pawrhoij_tmp,electronpositron%pawrhoij_ep)
463 0 : call pawrhoij_free(pawrhoij_tmp)
464 0 : ABI_FREE(pawrhoij_tmp)
465 : end if
466 0 : if (allocated(electronpositron%nhat_ep)) then
467 : call pawmknhat(occtmp,1,0,0,0,0,gprimd,my_natom,dtset%natom,nfft,ngfft,0,&
468 : & dtset%nspden,dtset%ntypat,pawang,pawfgrtab,nhatgr,electronpositron%nhat_ep,&
469 : & electronpositron%pawrhoij_ep,electronpositron%pawrhoij_ep,pawtab,qphon,rprimd,&
470 : & ucvol,dtset%usewvl,xred,distribfft=mpi_enreg%distribfft,&
471 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
472 0 : & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0)
473 : end if
474 : end if
475 : end if
476 : ! ----- Positronic from scratch
477 6 : if (dtset%positron<0.and.electronpositron%calctype==2) then
478 0 : electronpositron%rhor_ep(:,1)=one/ucvol
479 0 : if (dtset%nspden>=2) electronpositron%rhor_ep(:,2)=half/ucvol
480 0 : if (dtset%nspden==4) electronpositron%rhor_ep(:,3:4)=zero
481 0 : rhog_ep=zero;rhog_ep(1,1)=one/ucvol
482 0 : if (dtset%usepaw==1) then
483 0 : do iatom=1,dtset%natom
484 0 : electronpositron%pawrhoij_ep(iatom)%rhoijp(:,:)=zero
485 0 : electronpositron%pawrhoij_ep(iatom)%nrhoijsel=0
486 : end do
487 0 : if (allocated(electronpositron%nhat_ep)) electronpositron%nhat_ep=zero
488 : end if
489 : end if
490 : end if
491 : ! ----- Deduced from rhor in memory
492 79 : if (history_level==2) then
493 462310 : electronpositron%rhor_ep(:,:)=rhor(:,:)
494 1386890 : rhog_ep(:,:)=rhog(:,:)
495 8 : if (dtset%usepaw==1) then
496 8 : call pawrhoij_copy(pawrhoij,electronpositron%pawrhoij_ep)
497 462310 : if (allocated(electronpositron%nhat_ep)) electronpositron%nhat_ep(:,:)=nhat(:,:)
498 : end if
499 : end if
500 :
501 : ! ===== CURRENT DENSITY RHOR:
502 79 : if (history_level==0.or.history_level==2) then
503 14 : if (ireadwf==0) then
504 : ! ----- Positronic from scratch
505 14 : if (electronpositron%calctype==1) then
506 705820 : rhor(:,1)=one/ucvol
507 13 : if (dtset%nspden>=2) rhor(:,2)=half/ucvol
508 13 : if (dtset%nspden==4) rhor(:,3:4)=zero
509 2117434 : rhog=zero;rhog(1,1)=one/ucvol
510 13 : if (dtset%usepaw==1) then
511 85 : do iatom=1,my_natom
512 2942 : pawrhoij(iatom)%rhoijp(:,:)=zero
513 85 : pawrhoij(iatom)%nrhoijsel=0
514 : end do
515 686148 : nhat(:,:)=zero
516 : end if
517 : end if
518 : ! ----- Electronic from scratch
519 14 : if (electronpositron%calctype==2) then
520 1 : ecut_eff=dtset%pawecutdg*(dtset%dilatmx)**2
521 1 : call getcut(boxcut_dum,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,std_out,qphon,ngfft)
522 : call initro(atindx,dtset%densty,gmet,gsqcut_eff,dtset%usepaw,mgfft,mpi_enreg,&
523 : & psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,dtset%ntypat,&
524 : & psps,pawtab,ph1d,psps%qgrid_vl,rhog,rhor,dtset%spinat,ucvol,&
525 1 : & dtset%usepaw,dtset%ziontypat,dtset%znucl)
526 :
527 1 : if (dtset%usepaw==1) then
528 1 : if (size(pawrhoij)>0) then
529 7 : ABI_MALLOC(pawrhoij_tmp,(my_natom))
530 : call initrhoij(pawrhoij(1)%cplex_rhoij,dtset%lexexch,dtset%lpawu,&
531 : & my_natom,dtset%natom,dtset%nspden,pawrhoij(1)%nspinor,dtset%nsppol,&
532 : & dtset%ntypat,pawrhoij_tmp,dtset%pawspnorb,pawtab,pawrhoij(1)%qphase,dtset%spinat,&
533 : & dtset%typat,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,&
534 : & use_rhoij_=pawrhoij(1)%use_rhoij_,use_rhoijres=pawrhoij(1)%use_rhoijres,&
535 1 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
536 5 : do iatom=1,my_natom
537 149 : pawrhoij_tmp(iatom)%kpawmix(:)=pawrhoij(iatom)%kpawmix(:)
538 : end do
539 1 : call pawrhoij_copy(pawrhoij_tmp,pawrhoij)
540 1 : call pawrhoij_free(pawrhoij_tmp)
541 5 : ABI_FREE(pawrhoij_tmp)
542 : end if
543 : call pawmknhat(occtmp,1,0,0,0,0,gprimd,my_natom,dtset%natom,nfft,ngfft,0,&
544 : & dtset%nspden,dtset%ntypat,pawang,pawfgrtab,nhatgr,nhat,&
545 : & pawrhoij,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred, &
546 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
547 : & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,&
548 1 : & me_g0=mpi_enreg%me_g0,distribfft=mpi_enreg%distribfft)
549 : end if
550 :
551 : end if
552 : end if
553 : end if
554 :
555 : ! ===== EXCHANGE POSITRONIC AND ELECTRONIC DENSITY (CURRENT AND PREVIOUS)
556 79 : if (history_level==3) then
557 124 : do ispden=1,dtset%nspden
558 8719180 : do ifft=1,nfft
559 8719056 : rhotmp=rhor(ifft,ispden)
560 8719056 : rhor(ifft,ispden)=electronpositron%rhor_ep(ifft,ispden)
561 8719118 : electronpositron%rhor_ep(ifft,ispden)=rhotmp
562 : end do
563 : end do
564 26157230 : rhog_ep(:,:)=rhog
565 62 : call fourdp(1,rhog,rhor,-1,mpi_enreg,nfft,1,ngfft,0)
566 : ! If PAW, exchange "positronic" and "electronic" rhoij
567 62 : if (dtset%usepaw==1) then
568 62 : if (size(pawrhoij)>0.and.size(electronpositron%pawrhoij_ep)>0) then
569 1028 : ABI_MALLOC(pawrhoij_tmp,(my_natom))
570 : call pawrhoij_alloc(pawrhoij_tmp,pawrhoij(1)%cplex_rhoij,pawrhoij(1)%nspden,&
571 : & pawrhoij(1)%nspinor,pawrhoij(1)%nsppol,dtset%typat,&
572 : & pawtab=pawtab,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,&
573 : & use_rhoij_=pawrhoij(1)%use_rhoij_,use_rhoijres=pawrhoij(1)%use_rhoijres, &
574 62 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
575 62 : call pawrhoij_copy(pawrhoij,pawrhoij_tmp)
576 62 : call pawrhoij_copy(electronpositron%pawrhoij_ep,pawrhoij)
577 62 : call pawrhoij_copy(pawrhoij_tmp,electronpositron%pawrhoij_ep)
578 62 : call pawrhoij_free(pawrhoij_tmp)
579 904 : ABI_FREE(pawrhoij_tmp)
580 : end if
581 62 : if (allocated(electronpositron%nhat_ep)) then
582 124 : do ispden=1,dtset%nspden
583 8719180 : do ifft=1,nfft
584 8719056 : rhotmp=nhat(ifft,ispden)
585 8719056 : nhat(ifft,ispden)=electronpositron%nhat_ep(ifft,ispden)
586 8719118 : electronpositron%nhat_ep(ifft,ispden)=rhotmp
587 : end do
588 : end do
589 : end if
590 : end if
591 : end if
592 :
593 : ! ===== COMPUTE HARTREE POTENTIAL ASSOCIATED TO RHOR_EP
594 79 : if (history_level==4) then
595 3 : call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,1,ngfft,0)
596 : end if
597 79 : if (history_level/=-1) then
598 : call hartre(1,gsqcut,dtset%icutcoul,dtset%usepaw,mpi_enreg,nfft,ngfft,&
599 79 : & dtset%nkpt,dtset%rcut,rhog_ep,rprimd,dtset%vcutgeo,electronpositron%vha_ep)
600 9924334 : electronpositron%vha_ep=-electronpositron%vha_ep
601 : else
602 0 : electronpositron%vha_ep=zero
603 : end if
604 79 : ABI_FREE(rhog_ep)
605 :
606 : ! ----------------------------------------------------------------------
607 : ! Initialize/Update energies
608 : ! ----------------------------------------------------------------------
609 79 : electronpositron%etotal_prev=etotal
610 79 : electronpositron%maxfor_prev=maxfor
611 :
612 : ! Inits/exchange news energies
613 : ! Retrieve energy of non-evolving particle(s)
614 79 : if (history_level== 0) then
615 6 : call energies%init()
616 6 : call electronpositron%energies_ep%init()
617 6 : if (dtset%positron>0) energies%e0_electronpositron=etotal_read
618 6 : if (dtset%positron<0) energies%e0_electronpositron=zero
619 : else if (history_level== 1) then
620 : call electronpositron%energies_ep%init()
621 : if (dtset%positron>0) energies%e0_electronpositron=etotal_read
622 73 : else if (history_level== 2) then
623 8 : call energies%copy(electronpositron%energies_ep)
624 8 : call energies%init()
625 8 : energies%e0_electronpositron=electronpositron%e0
626 65 : else if (history_level== 3) then
627 62 : call electronpositron%energies_ep%copy(energies_tmp)
628 62 : call energies%copy(electronpositron%energies_ep)
629 62 : call energies_tmp%copy(energies)
630 62 : energies%e0_electronpositron=electronpositron%e0
631 : ! else if (history_level== 4) then
632 : end if
633 :
634 : ! Adjust core psps energy
635 79 : if (electronpositron%calctype/=1) energies%e_corepsp=ecore/ucvol
636 :
637 : ! -----------------------------------------------------------------------------------------
638 : ! Update wavefunctions
639 : ! If electronpositron%calctype==1: cg are the positronic WFs, cg_ep are the electronic WFs
640 : ! If electronpositron%calctype==2: cg are the electronic WFs, cg_ep are the positronic WFs
641 : ! -----------------------------------------------------------------------------------------
642 79 : if (electronpositron%dimcg>0.or.electronpositron%dimcprj>0) then
643 :
644 70 : if (history_level==0.or.history_level==1) then
645 0 : electronpositron%cg_ep=zero
646 : end if
647 :
648 70 : if (history_level==2) then
649 7 : if (electronpositron%dimcg>0) then
650 439563 : do icg=1,electronpositron%dimcg
651 1318675 : electronpositron%cg_ep(1:2,icg)=cg(1:2,icg)
652 : end do
653 : end if
654 7 : if (dtset%usepaw==1.and.electronpositron%dimcprj>0) then
655 3 : call pawcprj_copy(cprj,electronpositron%cprj_ep)
656 : end if
657 : end if
658 :
659 70 : if (history_level==3) then
660 60 : if (electronpositron%dimcg>0) then
661 7045116 : do icg=1,electronpositron%dimcg
662 21135168 : cgtmp(1:2)=electronpositron%cg_ep(1:2,icg)
663 21135168 : electronpositron%cg_ep(1:2,icg)=cg(1:2,icg)
664 21135228 : cg(1:2,icg)=cgtmp(1:2)
665 : end do
666 : end if
667 60 : if (dtset%usepaw==1.and.electronpositron%dimcprj>0) then
668 6 : ABI_MALLOC(nlmn,(dtset%natom))
669 608 : ABI_MALLOC(cprj_tmp,(dtset%natom,electronpositron%dimcprj))
670 10 : do iatom=1,dtset%natom
671 10 : nlmn(iatom)=cprj(iatom,1)%nlmn
672 : end do
673 2 : call pawcprj_alloc(cprj_tmp,cprj(1,1)%ncpgr,nlmn)
674 2 : ABI_FREE(nlmn)
675 2 : call pawcprj_copy(electronpositron%cprj_ep,cprj_tmp)
676 2 : call pawcprj_copy(cprj,electronpositron%cprj_ep)
677 2 : call pawcprj_copy(cprj_tmp,cprj)
678 2 : call pawcprj_free(cprj_tmp)
679 482 : ABI_FREE(cprj_tmp)
680 : end if
681 : end if
682 :
683 : end if ! dimcg>0 or dimcprj>0
684 :
685 : ! -----------------------------------------------------------------------------------------------------------
686 : ! Initialize/Update occupations
687 : ! If electronpositron%calctype==1: occ are the positronic occupations, occ_ep are the electronic occupations
688 : ! If electronpositron%calctype==2: occ are the electronic occupations, occ_ep are the positronic occupations
689 : ! -----------------------------------------------------------------------------------------------------------
690 : ! When needed, precompute electronic occupations with semiconductor occupancies
691 79 : need_scocc=.false.
692 79 : if (electronpositron%dimocc>0.and.electronpositron%calctype==1.and. &
693 79 : & (history_level==0.or.history_level==1)) need_scocc=.true.
694 79 : if (electronpositron%calctype==2.and.ireadwf==0.and. &
695 : & (history_level==0.or.history_level==2.or. &
696 : & (history_level==3.and.electronpositron%dimocc==0))) need_scocc=.true.
697 77 : if (need_scocc) then
698 2 : nelect=-dtset%cellcharge(1)
699 10 : do iatom=1,dtset%natom
700 10 : nelect=nelect+dtset%ziontypat(dtset%typat(iatom))
701 : end do
702 2 : maxocc=two/real(dtset%nsppol*dtset%nspinor,dp)
703 2 : nocc=int((nelect-tol8)/maxocc) + 1
704 2 : nocc=min(nocc,dtset%nband(1)*dtset%nsppol)
705 2 : occlast=nelect-maxocc*(nocc-1)
706 6 : ABI_MALLOC(scocc,(dtset%nband(1)*dtset%nsppol))
707 32 : scocc=zero
708 12 : if (1<nocc) scocc(1:nocc-1)=maxocc
709 2 : if (1<=nocc) scocc(nocc)=occlast
710 : end if
711 :
712 : ! ===== PREVIOUS OCCUPATIONS OCC_EP:
713 79 : if (electronpositron%dimocc>0) then
714 70 : if (history_level==0.or.history_level==1) then
715 : ! ----- Electronic from scratch
716 0 : if (electronpositron%calctype==1) then
717 : ! Initialize electronic occupations with semiconductor occupancies
718 0 : do ikpt=1,dtset%nkpt
719 0 : do iband=1,dtset%nband(1)
720 0 : do isppol=1,dtset%nsppol
721 : electronpositron%occ_ep(iband+dtset%nband(1)*(ikpt-1+dtset%nkpt*(isppol-1)))=&
722 0 : & scocc(isppol+dtset%nsppol*(iband-1))
723 : end do
724 : end do
725 : end do
726 : end if
727 : ! ----- Positronic from scratch
728 0 : if (electronpositron%calctype==1) then
729 : ! Initialize positronic occupations with only one positron (or less)
730 0 : electronpositron%occ_ep(:)=zero
731 0 : isppol=1;iocc=1
732 0 : do ikpt=1,dtset%nkpt
733 0 : electronpositron%occ_ep(iocc)=electronpositron%posocc
734 0 : iocc=iocc+dtset%nband(ikpt+dtset%nkpt*(isppol-1))
735 : end do
736 : end if
737 : end if
738 : ! ----- Deduced from occ in memory
739 70 : if (history_level==2) then
740 609 : electronpositron%occ_ep(:)=occ(:)
741 : end if
742 : end if ! dimocc>0
743 :
744 : ! ===== CURRENT OCCUPATIONS OCC:
745 79 : if (history_level==0.or.history_level==2.or.(history_level==3.and.electronpositron%dimocc==0)) then
746 16 : if (ireadwf==0) then
747 : ! ----- Positronic from scratch
748 16 : if (electronpositron%calctype==1) then
749 : ! Initialize positronic occupations with only one positron (or less)
750 808 : occ(:)=zero
751 14 : isppol=1;iocc=1
752 63 : do ikpt=1,dtset%nkpt
753 49 : occ(iocc)=electronpositron%posocc
754 63 : iocc=iocc+dtset%nband(ikpt+dtset%nkpt*(isppol-1))
755 : end do
756 : end if
757 : ! ----- Electronic from scratch
758 16 : if (electronpositron%calctype==2) then
759 : ! Initialize electronic occupations with semiconductor occupancies
760 10 : do ikpt=1,dtset%nkpt
761 130 : do iband=1,dtset%nband(1)
762 248 : do isppol=1,dtset%nsppol
763 : occ(iband+dtset%nband(1)*(ikpt-1+dtset%nkpt*(isppol-1)))=&
764 240 : & scocc(isppol+dtset%nsppol*(iband-1))
765 : end do
766 : end do
767 : end do
768 : end if
769 : end if
770 : end if
771 :
772 : ! ===== EXCHANGE POSITRONIC AND ELECTRONIC OCCUPATIONS (CURRENT AND PREVIOUS)
773 76 : if (history_level==3.and.electronpositron%dimocc>0) then
774 6312 : do iocc=1,electronpositron%dimocc
775 6252 : occtmp=occ(iocc)
776 6252 : occ(iocc)=electronpositron%occ_ep(iocc)
777 6312 : electronpositron%occ_ep(iocc)=occtmp
778 : end do
779 : end if
780 :
781 79 : if (need_scocc) then
782 2 : ABI_FREE(scocc)
783 : end if
784 :
785 : ! -----------------------------------------------------------------------------------------------------------
786 : ! Initialize/Update eigen energies
787 : ! If electronpositron%calctype==1: eigen are the positronic eigen E, eigen_ep are the electronic eigen E
788 : ! If electronpositron%calctype==2: eigen are the electronic eigen E, eigen_ep are the positronic eigen E
789 : ! -----------------------------------------------------------------------------------------------------------
790 :
791 : ! ===== PREVIOUS EIGEN ENERGIES EIGEN_EP:
792 79 : if (electronpositron%dimeigen>0) then
793 70 : if (history_level==0.or.history_level==1) then
794 : ! ----- Electronic or positronic from scratch
795 0 : electronpositron%eigen_ep(:)=zero
796 : end if
797 : ! ----- Deduced from eigen in memory
798 70 : if (history_level==2) then
799 609 : electronpositron%eigen_ep(:)=eigen(:)
800 : end if
801 : end if ! dimeigen>0
802 :
803 : ! ===== CURRENT EIGEN ENERGIES EIGEN:
804 79 : if (history_level==0.or.history_level==2.or.(history_level==3.and.electronpositron%dimeigen==0)) then
805 16 : if (ireadwf==0) then
806 : ! ----- Electronic or positronic from scratch
807 930 : eigen(:)=zero
808 : end if
809 : end if
810 :
811 : ! ===== EXCHANGE POSITRONIC AND ELECTRONIC EIGEN ENERGIES (CURRENT AND PREVIOUS)
812 76 : if (history_level==3.and.electronpositron%dimeigen>0) then
813 6312 : do iocc=1,electronpositron%dimeigen
814 6252 : eigtmp=eigen(iocc)
815 6252 : eigen(iocc)=electronpositron%eigen_ep(iocc)
816 6312 : electronpositron%eigen_ep(iocc)=eigtmp
817 : end do
818 : end if
819 :
820 : ! In some cases cprj are kept in memory, so we have to update them before the call of vtorho
821 79 : if (dtset%cprj_in_memory==2) then
822 0 : iatom=0
823 0 : call wrtout(std_out,' Computing cprj from wavefunctions (positron)')
824 : call ctocprj(atindx,cg,1,cprj,gmet,gprimd,iatom,0,&
825 : & 0,dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,&
826 : & dtset%mpw,dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,&
827 : & dtset%nsppol,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,dtfil%unpaw,&
828 0 : & xred,ylm,ylmgr)
829 0 : call wrtout(std_out,' cprj is computed')
830 : end if
831 :
832 : ! =============================================
833 : end if ! the type of e-p calculation changes
834 : !=============================================
835 :
836 : !------------------------------------------------------------------
837 : !Write messages
838 : !------------------------------------------------------------------
839 525 : if (istep_mix==1.and.dtset%positron/=0) then
840 : ! Log message
841 87 : if (electronpositron%calctype==0) then
842 8 : message = 'Were are now performing an electronic ground-state calculation...'
843 79 : else if (electronpositron%calctype==1) then
844 43 : message = 'Were are now performing a positronic ground-state calculation...'
845 36 : else if (electronpositron%calctype==2) then
846 36 : message = 'Were are now performing an electronic ground-state calculation in presence of a positron...'
847 : end if
848 87 : ABI_COMMENT(message)
849 : ! Output message
850 87 : if (dtset%positron<0) then
851 81 : if (electronpositron%calctype==0) then
852 8 : TypeCalcStrg='ELECTRONIC GROUND-STATE CALCULATION'
853 73 : else if (electronpositron%calctype==1) then
854 38 : TypeCalcStrg='POSITRONIC GROUND-STATE CALCULATION IN PRESENCE OF ELECTRONS AND IONS'
855 35 : else if (electronpositron%calctype==2) then
856 35 : TypeCalcStrg='ELECTRONIC GROUND-STATE CALCULATION IN PRESENCE OF A POSITRON'
857 : end if
858 81 : if (istep>1) then
859 70 : write(message,'(2a,i3,2a)') ch10,'TC-DFT STEP ',electronpositron%istep,' - ',trim(TypeCalcStrg)
860 : else
861 11 : write(message,'(a,i3,2a)') 'TC-DFT STEP ',electronpositron%istep,' - ',trim(TypeCalcStrg)
862 : end if
863 81 : call wrtout(ab_out,message,'COLL')
864 : end if
865 : end if
866 :
867 1050 : end subroutine setup_positron
868 : !!***
869 :
870 : !!****f* ABINIT/poslifetime
871 : !! NAME
872 : !! poslifetime
873 : !!
874 : !! FUNCTION
875 : !! Calculate the positron lifetime
876 : !!
877 : !! INPUTS
878 : !! dtset <type(dataset_type)>=all input variables for this dataset
879 : !! | nspden=number of spin-density components
880 : !! | ntypat=number of atom types
881 : !! | paral_kgb=flag controlling (k,g,bands) parallelization
882 : !! | pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
883 : !! | usepaw=flag for PAW
884 : !! gprimd(3,3)= dimensional reciprocal space primitive translations
885 : !! mpi_enreg= information about MPI parallelization
886 : !! my_natom=number of atoms treated by current processor
887 : !! n3xccc= dimension of the xccc3d array (0 or nfft).
888 : !! nfft= number of FFT grid points
889 : !! ngfft(18)= contain all needed information about 3D FFT
890 : !! nhat(nfft,nspden)=charge compensation density (content depends on electronpositron%particle)
891 : !! option= if 1, calculate positron lifetime for whole density
892 : !! if 2, calculate positron lifetime for given state
893 : !! if 3, calculate positron lifetime for given state with IPM
894 : !! pawang <type(pawang)>=paw angular mesh and related data
895 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
896 : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies
897 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
898 : !! rhor(nfft,nspden)=total electron/positron density (content depends on electronpositron%particle)
899 : !! ucvol=unit cell volume in bohr**3.
900 : !! xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
901 : !! ===== Optional arguments, used only if option>1 =====
902 : !! pawrhoij_dop_el(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies of one state
903 : !! rhor_dop_el(nfft)=electron density of given state for the state dependent scheme
904 : !! ===== Optional argument =====
905 : !! pawrhoij_ep(my_natom*usepaw) <type(pawrhoij_type)>= atomic occupancies to be used in place of
906 : !! electronpositron%pawrhoij_ep
907 : !!
908 : !! OUTPUT
909 : !! rate= annihilation rate of a given state needed for state dependent scheme for doppler broadening
910 : !!
911 : !! SIDE EFFECTS
912 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
913 : !!
914 : !! SOURCE
915 :
916 658 : subroutine poslifetime(dtset,electronpositron,gprimd,my_natom,mpi_enreg,n3xccc,nfft,ngfft,nhat,&
917 329 : & option,pawang,pawrad,pawrhoij,pawtab,rate,rate_paw,rhor,ucvol,xccc3d,&
918 312 : & rhor_dop_el,pawrhoij_dop_el,pawrhoij_ep) ! optional arguments
919 :
920 : !Arguments ------------------------------------
921 : !scalars
922 : integer,intent(in) :: my_natom,n3xccc,nfft,option
923 : real(dp),intent(in) :: ucvol
924 : real(dp),intent(out) :: rate,rate_paw
925 : type(dataset_type), intent(in) :: dtset
926 : type(electronpositron_type),pointer :: electronpositron
927 : type(MPI_type),intent(in) :: mpi_enreg
928 : type(pawang_type), intent(in) :: pawang
929 : !arrays
930 : integer,intent(in) :: ngfft(18)
931 : real(dp),intent(in) :: gprimd(3,3),nhat(nfft,dtset%nspden*dtset%usepaw),xccc3d(n3xccc)
932 : real(dp),intent(in),target :: rhor(nfft,dtset%nspden)
933 : real(dp),optional,intent(in) :: rhor_dop_el(nfft)
934 : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw)
935 : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*dtset%usepaw)
936 : type(pawrhoij_type),optional,intent(in) :: pawrhoij_dop_el(my_natom*dtset%usepaw)
937 : type(pawrhoij_type),optional,target,intent(in) :: pawrhoij_ep(my_natom*dtset%usepaw)
938 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
939 :
940 : !Local variables-------------------------------
941 : !scalars
942 : integer :: cplex,iatom,ierr,ifft,igam,ii,ilm,ilm1,ilm2,iloop,ipt,ir,isel
943 : integer :: itypat,iwarn,iwarnj,iwarnp,lm_size,lmn2_size,mesh_size
944 : integer :: nfftot,ngamma,ngr,ngrad,nspden_ep,opt_dens,usecore
945 : logical,parameter :: include_nhat_in_gamma=.false.
946 : real(dp),parameter :: delta=1.d-4
947 : real(dp) :: fact,fact2,intg
948 : real(dp) :: lambda_core ,lambda_core_ipm ,lambda ,lambda_ipm
949 : real(dp) :: lambda_core_paw,lambda_core_paw_ipm,lambda_paw,lambda_paw_ipm
950 : real(dp) :: lifetime,lifetime_ipm,nbec,nbev,nbp,rdum,sqfpi,units
951 : character(len=500) :: msg
952 : !arrays
953 329 : integer,allocatable :: igamma(:)
954 329 : logical,allocatable :: lmselect(:),lmselect_ep(:),lmselect_dum(:)
955 : real(dp) :: mpibuf(4)
956 : real(dp),parameter :: qphon(3)=(/zero,zero,zero/),lsign(2)=(/one,-one/)
957 329 : real(dp),allocatable :: d1gam(:,:,:),d2gam(:,:,:),ff(:),gam_(:,:,:),gamma(:,:),gammam(:,:,:),gg(:,:)
958 329 : real(dp),allocatable :: grhocore2(:),grhocor2_(:),grhoe2(:),grho2_(:)
959 329 : real(dp),allocatable :: nhat1(:,:,:),nhat1_ep(:,:,:),nhat1_j(:,:,:)
960 329 : real(dp),allocatable :: rho_(:),rho_ep_(:),rho1(:,:,:),rho1_ep(:,:,:),rho1_j(:,:,:)
961 329 : real(dp),allocatable :: rhoarr1(:),rhoarr1_ep(:),rhoarr1_j(:),rhoarr2(:)
962 329 : real(dp),allocatable :: rhocore(:),rhocor_(:),rhoe(:,:),rhop(:,:),rhor_dop_el_(:)
963 329 : real(dp),allocatable :: rhosph(:),rhosph_ep(:),rhosph_j(:),rhotot(:,:),rhotot_ep(:,:)
964 329 : real(dp),allocatable :: rhotot_j(:,:),trho1(:,:,:),trho1_ep(:,:,:),trho1_j(:,:,:)
965 329 : real(dp),allocatable :: v1sum(:,:),v2sum(:,:,:)
966 329 : real(dp),pointer :: rhor_(:,:),rhor_ep_(:,:)
967 329 : type(pawrhoij_type),pointer :: pawrhoij_ep_(:)
968 :
969 : ! *************************************************************************
970 :
971 : DBG_ENTER("COLL")
972 :
973 : !Tests for developers
974 329 : if (.not.associated(electronpositron)) then
975 0 : msg='electronpositron variable must be associated!'
976 0 : ABI_BUG(msg)
977 : end if
978 329 : if (option/=1) then
979 312 : if ((.not.present(rhor_dop_el)).or.(.not.present(pawrhoij_dop_el))) then
980 0 : msg='when option/=1, rhor_dop_el and pawrhoij_dop_el must be present!'
981 0 : ABI_BUG(msg)
982 : end if
983 : end if
984 :
985 : ! This to avoid using uninitialized variables.
986 329 : lambda_core = zero; lambda_paw = zero; lambda_core_paw = zero
987 :
988 : !Constants
989 329 : fact=0.0
990 329 : cplex=1;nspden_ep=1
991 329 : usecore=n3xccc/nfft
992 329 : nfftot=ngfft(1)*ngfft(2)*ngfft(3)
993 329 : ngrad=1;if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2
994 329 : iwarn=0;iwarnj=0;iwarnp=1
995 329 : sqfpi=sqrt(four_pi)
996 :
997 : !Compatibility tests
998 329 : if (electronpositron%particle==EP_NOTHING) then
999 0 : msg='Not valid for electronpositron%particle=NOTHING!'
1000 0 : ABI_BUG(msg)
1001 : end if
1002 329 : if (electronpositron%nfft/=nfft) then
1003 0 : msg='nfft/=electronpositron%nfft!'
1004 0 : ABI_BUG(msg)
1005 : end if
1006 329 : if (dtset%usepaw==1) then
1007 328 : if(dtset%pawxcdev==0.and.ngrad==2) then
1008 0 : msg='GGA is not implemented for pawxcdev=0 (use dtset%pawxcdev/=0)!'
1009 0 : ABI_BUG(msg)
1010 : end if
1011 : end if
1012 :
1013 : !Select type(s) of enhancement factor
1014 329 : if ((electronpositron%ixcpositron==1.or.electronpositron%ixcpositron==3).and.option==1) then
1015 11 : ngamma=2
1016 11 : ABI_MALLOC(igamma,(ngamma))
1017 11 : igamma(1)=1;igamma(2)=2
1018 : else
1019 318 : ngamma=1
1020 318 : ABI_MALLOC(igamma,(ngamma))
1021 318 : if (electronpositron%ixcpositron==-1) igamma(1)=0
1022 318 : if (electronpositron%ixcpositron== 2) igamma(1)=4
1023 318 : if (electronpositron%ixcpositron==11.or.electronpositron%ixcpositron==31) igamma(1)=3
1024 318 : if (electronpositron%ixcpositron==1.or.electronpositron%ixcpositron==3) igamma(1)=2
1025 : end if
1026 :
1027 : !Select density according to nhat choice
1028 329 : if (dtset%usepaw==0.or.include_nhat_in_gamma) then
1029 1 : rhor_ => rhor
1030 1 : rhor_ep_ => electronpositron%rhor_ep
1031 : else
1032 1312 : ABI_MALLOC(rhor_,(nfft,dtset%nspden))
1033 984 : ABI_MALLOC(rhor_ep_,(nfft,dtset%nspden))
1034 15474040 : rhor_=rhor-nhat
1035 15474040 : rhor_ep_=electronpositron%rhor_ep-electronpositron%nhat_ep
1036 : end if
1037 :
1038 : !Eventually overwrite electronpositron%pawrhoij_ep
1039 329 : if (present(pawrhoij_ep)) then
1040 312 : pawrhoij_ep_ => pawrhoij_ep
1041 : else
1042 17 : pawrhoij_ep_ => electronpositron%pawrhoij_ep
1043 : end if
1044 :
1045 : !Loop on different enhancement factors
1046 669 : do igam=1,ngamma
1047 :
1048 : ! Compute electron-positron annihilation rate using pseudo densities (plane waves)
1049 : ! ----------------------------------------------------------------------------------------
1050 :
1051 : ! Select the densities and make them positive
1052 1360 : ABI_MALLOC(rhoe,(nfft,nspden_ep))
1053 1020 : ABI_MALLOC(rhop,(nfft,nspden_ep))
1054 340 : if (electronpositron%particle==EP_ELECTRON) then
1055 920061 : rhoe(:,1)=rhor_ep_(:,1);rhop(:,1)=rhor_(:,1)
1056 331 : else if (electronpositron%particle==EP_POSITRON) then
1057 16678771 : rhoe(:,1)=rhor_(:,1);rhop(:,1)=rhor_ep_(:,1)
1058 : end if
1059 340 : call mkdenpos(iwarn ,nfft,nspden_ep,1,rhoe,dtset%xc_denpos)
1060 340 : call mkdenpos(iwarnp,nfft,nspden_ep,1,rhop,dtset%xc_denpos)
1061 340 : if (option/=1) then
1062 936 : ABI_MALLOC(rhor_dop_el_,(nfft))
1063 6551160 : rhor_dop_el_(:)=rhor_dop_el(:)
1064 312 : call mkdenpos(iwarnp,nfft,1,1,rhor_dop_el_,dtset%xc_denpos)
1065 : end if
1066 :
1067 : ! Compute enhancement factor at each FFT grid point
1068 : ! gamma(:,1): using total electronic density
1069 : ! gamma(:,2): using valence electronic density
1070 1020 : ABI_MALLOC(gamma,(nfft,2))
1071 340 : if (option==1.or.option==2) then
1072 : call gammapositron_fft(electronpositron,gamma,gprimd,igamma(igam),mpi_enreg,&
1073 184 : & n3xccc,nfft,ngfft,rhoe,rhop,xccc3d)
1074 : else
1075 6551316 : gamma=one
1076 : end if
1077 :
1078 : ! Compute positron annihilation rates
1079 340 : lambda =zero;lambda_ipm =zero
1080 340 : lambda_core=zero;lambda_core_ipm=zero
1081 340 : if (option==1) then
1082 2248426 : do ifft=1,nfft
1083 2248398 : lambda =lambda +rhop(ifft,1)*rhoe(ifft,1)*gamma(ifft,1)
1084 2248426 : lambda_ipm=lambda_ipm+rhop(ifft,1)*rhoe(ifft,1)*gamma(ifft,2)
1085 : end do
1086 : else
1087 6551160 : do ifft=1,nfft
1088 6550848 : lambda =lambda +rhop(ifft,1)*rhor_dop_el_(ifft)*gamma(ifft,1)
1089 6551160 : lambda_ipm=lambda_ipm+rhop(ifft,1)*rhor_dop_el_(ifft)*gamma(ifft,2)
1090 : end do
1091 : end if
1092 340 : if (usecore==1) then
1093 4941522 : do ifft=1,nfft
1094 4941378 : lambda_core =lambda_core +rhop(ifft,1)*xccc3d(ifft)*gamma(ifft,1)
1095 4941522 : lambda_core_ipm=lambda_core_ipm+rhop(ifft,1)*xccc3d(ifft)
1096 : end do
1097 : end if
1098 340 : lambda =lambda *ucvol/dble(nfftot)
1099 340 : lambda_ipm =lambda_ipm *ucvol/dble(nfftot)
1100 340 : lambda_core =lambda_core *ucvol/dble(nfftot)
1101 340 : lambda_core_ipm=lambda_core_ipm*ucvol/dble(nfftot)
1102 340 : ABI_FREE(gamma)
1103 340 : ABI_FREE(rhoe)
1104 340 : ABI_FREE(rhop)
1105 340 : if (option/=1) then
1106 312 : ABI_FREE(rhor_dop_el_)
1107 : end if
1108 : ! NC pseudopotential: check electrons/positron number
1109 340 : if (dtset%usepaw==0.and.igam==ngamma) then
1110 1 : nbec=zero;nbev=zero;nbp=zero
1111 1 : if (electronpositron%particle==EP_ELECTRON) then
1112 19684 : do ifft=1,nfft
1113 19683 : nbec=nbec+xccc3d(ifft)
1114 19683 : nbev=nbev+electronpositron%rhor_ep(ifft,1)
1115 19684 : nbp =nbp +rhor(ifft,1)
1116 : end do
1117 : else
1118 0 : do ifft=1,nfft
1119 0 : nbec=nbec+xccc3d(ifft)
1120 0 : nbev=nbev+rhor(ifft,1)
1121 0 : nbp =nbp +electronpositron%rhor_ep(ifft,1)
1122 : end do
1123 : end if
1124 1 : nbec=nbec*ucvol/dble(nfftot)
1125 1 : nbev=nbev*ucvol/dble(nfftot)
1126 1 : nbp =nbp *ucvol/dble(nfftot)
1127 : end if
1128 :
1129 : ! MPI parallelization
1130 340 : if(mpi_enreg%nproc_fft>1)then
1131 0 : call xmpi_sum(lambda ,mpi_enreg%comm_fft,ierr)
1132 0 : call xmpi_sum(lambda_ipm,mpi_enreg%comm_fft,ierr)
1133 0 : call xmpi_sum(lambda_core ,mpi_enreg%comm_fft,ierr)
1134 0 : call xmpi_sum(lambda_core_ipm,mpi_enreg%comm_fft,ierr)
1135 0 : if (dtset%usepaw==0.and.igam==ngamma) then
1136 0 : call xmpi_sum(nbec,mpi_enreg%comm_fft,ierr)
1137 0 : call xmpi_sum(nbev,mpi_enreg%comm_fft,ierr)
1138 0 : call xmpi_sum(nbp ,mpi_enreg%comm_fft,ierr)
1139 : end if
1140 : end if
1141 :
1142 :
1143 : ! PAW: add on-site contributions to electron-positron annihilation rate
1144 : ! ----------------------------------------------------------------------------------------
1145 340 : if (dtset%usepaw==1) then
1146 :
1147 338 : lambda_paw =zero;lambda_paw_ipm =zero
1148 338 : lambda_core_paw=zero;lambda_core_paw_ipm=zero
1149 :
1150 : ! Loop on atoms
1151 1294 : do iatom=1,my_natom
1152 :
1153 956 : itypat=pawrhoij(iatom)%itypat
1154 956 : lmn2_size=pawtab(itypat)%lmn2_size
1155 956 : mesh_size=pawtab(itypat)%mesh_size
1156 956 : lm_size=pawtab(itypat)%lcut_size**2
1157 956 : cplex=1
1158 956 : ngr=0;if (ngrad==2) ngr=mesh_size
1159 :
1160 : ! Allocations of "on-site" densities
1161 4780 : ABI_MALLOC(rho1 ,(cplex*mesh_size,lm_size,nspden_ep))
1162 3824 : ABI_MALLOC(trho1,(cplex*mesh_size,lm_size,nspden_ep))
1163 3824 : ABI_MALLOC(rho1_ep ,(cplex*mesh_size,lm_size,nspden_ep))
1164 3824 : ABI_MALLOC(trho1_ep,(cplex*mesh_size,lm_size,nspden_ep))
1165 956 : if (option/=1) then
1166 2944 : ABI_MALLOC(rho1_j ,(cplex*mesh_size,lm_size,nspden_ep))
1167 2944 : ABI_MALLOC(trho1_j,(cplex*mesh_size,lm_size,nspden_ep))
1168 : else
1169 220 : ABI_MALLOC(rho1_j ,(0,0,0))
1170 220 : ABI_MALLOC(trho1_j,(0,0,0))
1171 : end if
1172 : if (include_nhat_in_gamma) then
1173 : ABI_MALLOC(nhat1,(cplex*mesh_size,lm_size,nspden_ep))
1174 : ABI_MALLOC(nhat1_ep,(cplex*mesh_size,lm_size,nspden_ep))
1175 : else
1176 956 : ABI_MALLOC(nhat1,(0,0,0))
1177 956 : ABI_MALLOC(nhat1_ep,(0,0,0))
1178 : end if
1179 : if (include_nhat_in_gamma.and.option/=1) then
1180 : ABI_MALLOC(nhat1_j,(cplex*mesh_size,lm_size,nspden_ep))
1181 : else
1182 956 : ABI_MALLOC(nhat1_j,(0,0,0))
1183 : end if
1184 2868 : ABI_MALLOC(lmselect,(lm_size))
1185 1912 : ABI_MALLOC(lmselect_ep,(lm_size))
1186 1912 : ABI_MALLOC(lmselect_dum,(lm_size))
1187 :
1188 : ! Compute "on-site" densities (n1, ntild1, nhat1) for electron and positron =====
1189 9560 : lmselect(:)=.true.
1190 956 : opt_dens=1;if (include_nhat_in_gamma) opt_dens=0
1191 : call pawdensities(rdum,cplex,iatom,lmselect,lmselect_dum,lm_size,nhat1,nspden_ep,1,&
1192 : & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij(iatom),&
1193 956 : & pawtab(itypat),rho1,trho1)
1194 9560 : lmselect_ep(:)=.true.
1195 : call pawdensities(rdum,cplex,iatom,lmselect_ep,lmselect_dum,lm_size,nhat1_ep,nspden_ep,1,&
1196 : & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij_ep_(iatom),&
1197 956 : & pawtab(itypat),rho1_ep,trho1_ep)
1198 :
1199 : ! For state dependent scheme in Doppler =====
1200 : ! Compute "on-site" densities (n1, ntild1, nhat1) for a given electron state j=====
1201 956 : if (option/=1) then
1202 : opt_dens=1;if (include_nhat_in_gamma) opt_dens=0
1203 : call pawdensities(rdum,cplex,iatom,lmselect,lmselect_dum,lm_size,nhat1_j,nspden_ep,1,&
1204 : & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij_dop_el(iatom),&
1205 736 : & pawtab(itypat),rho1_j,trho1_j)
1206 : end if
1207 :
1208 : ! Compute contribution to annihilation rate:
1209 : ! Loop: first step: compute all-electron contribution (from n^1, n_c)
1210 : ! 2nd step: compute pseudo contribution (from tild_n^1, hat_n^1, tild_n_c)
1211 2868 : do iloop=1,2
1212 1912 : if (iloop==1) usecore=1
1213 956 : if (iloop==2) usecore=pawtab(itypat)%usetcore
1214 5736 : ABI_MALLOC(rhocore,(mesh_size))
1215 :
1216 : ! First formalism: use densities on r,theta,phi
1217 1912 : if (dtset%pawxcdev==0) then
1218 :
1219 0 : ABI_MALLOC(gamma,(mesh_size,2))
1220 0 : ABI_MALLOC(rhoarr1,(mesh_size))
1221 0 : ABI_MALLOC(rhoarr1_ep,(mesh_size))
1222 0 : if (option/=1) then
1223 0 : ABI_MALLOC(rhoarr1_j,(mesh_size))
1224 : end if
1225 : ! Loop on the angular part
1226 0 : do ipt=1,pawang%angl_size
1227 : ! Build densities
1228 0 : rhoarr1=zero;rhoarr1_ep=zero;rhocore=zero
1229 0 : if (option/=1) rhoarr1_j=zero
1230 0 : if (iloop==1) then
1231 0 : do ilm=1,lm_size
1232 0 : if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+rho1(:,ilm,1)*pawang%ylmr(ilm,ipt)
1233 : end do
1234 0 : if (option/=1) then
1235 0 : do ilm=1,lm_size
1236 0 : if (lmselect(ilm)) rhoarr1_j(:)=rhoarr1_j(:)+rho1_j(:,ilm,1)*pawang%ylmr(ilm,ipt)
1237 : end do
1238 : end if
1239 0 : do ilm=1,lm_size
1240 0 : if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+rho1_ep(:,ilm,1)*pawang%ylmr(ilm,ipt)
1241 : end do
1242 0 : if (usecore==1) rhocore(:)=pawtab(itypat)%coredens(:)
1243 : else
1244 : if (include_nhat_in_gamma) then
1245 : do ilm=1,lm_size
1246 : if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+(trho1(:,ilm,1)+nhat1(:,ilm,1))*pawang%ylmr(ilm,ipt)
1247 : end do
1248 : if (option/=1) then
1249 : do ilm=1,lm_size
1250 : if (lmselect(ilm)) rhoarr1_j(:)=rhoarr1_j(:)+(trho1_j(:,ilm,1)+nhat1_j(:,ilm,1))*pawang%ylmr(ilm,ipt)
1251 : end do
1252 : end if
1253 : do ilm=1,lm_size
1254 : if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+(trho1_ep(:,ilm,1)+nhat1_ep(:,ilm,1))*pawang%ylmr(ilm,ipt)
1255 : end do
1256 : else
1257 0 : do ilm=1,lm_size
1258 0 : if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+trho1(:,ilm,1)*pawang%ylmr(ilm,ipt)
1259 : end do
1260 0 : if (option/=1) then
1261 0 : do ilm=1,lm_size
1262 0 : if (lmselect(ilm)) rhoarr1_j(:)=rhoarr1(:)+trho1_j(:,ilm,1)*pawang%ylmr(ilm,ipt)
1263 : end do
1264 : end if
1265 0 : do ilm=1,lm_size
1266 0 : if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+trho1_ep(:,ilm,1)*pawang%ylmr(ilm,ipt)
1267 : end do
1268 : end if
1269 0 : if (usecore==1) rhocore(:)=pawtab(itypat)%tcoredens(:,1)
1270 : end if
1271 : ! Make the densities positive
1272 0 : if (electronpositron%particle==EP_ELECTRON) then
1273 0 : call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1 ,dtset%xc_denpos)
1274 0 : call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos)
1275 0 : else if (electronpositron%particle==EP_POSITRON) then
1276 0 : call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1 ,dtset%xc_denpos)
1277 0 : call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos)
1278 0 : if (option/=1) then
1279 0 : call mkdenpos(iwarnj,mesh_size,1,1,rhoarr1_j,dtset%xc_denpos)
1280 : end if
1281 : end if
1282 : ! Compute Gamma
1283 0 : ABI_MALLOC(grhoe2,(ngr))
1284 0 : ABI_MALLOC(grhocore2,(ngr))
1285 0 : if (option==1.or.option==2) then
1286 0 : if (electronpositron%particle==EP_ELECTRON) then
1287 : call gammapositron(gamma,grhocore2,grhoe2,igamma(igam),ngr,mesh_size,&
1288 0 : & rhocore,rhoarr1_ep,rhoarr1,usecore)
1289 0 : else if (electronpositron%particle==EP_POSITRON) then
1290 : call gammapositron(gamma,grhocore2,grhoe2,igamma(igam),ngr,mesh_size,&
1291 0 : & rhocore,rhoarr1,rhoarr1_ep,usecore)
1292 : end if
1293 : else
1294 0 : gamma(:,:)=one
1295 : end if
1296 0 : ABI_FREE(grhoe2)
1297 0 : ABI_FREE(grhocore2)
1298 : ! Compute contribution to annihilation rates
1299 0 : ABI_MALLOC(ff,(mesh_size))
1300 0 : if (option/=1) rhoarr1(:)=rhoarr1_j(:)
1301 0 : do ii=1,4
1302 0 : if (ii==1) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhoarr1_ep(1:mesh_size) &
1303 0 : & *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2
1304 0 : if (ii==2) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhoarr1_ep(1:mesh_size) &
1305 0 : & *gamma(1:mesh_size,2)*pawrad(itypat)%rad(1:mesh_size)**2
1306 0 : if (electronpositron%particle==EP_ELECTRON) then
1307 0 : if (ii==3) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhocore(1:mesh_size) &
1308 0 : & *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2
1309 0 : if (ii==4) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhocore(1:mesh_size) &
1310 0 : & *pawrad(itypat)%rad(1:mesh_size)**2
1311 : else
1312 0 : if (ii==3) ff(1:mesh_size)=rhoarr1_ep(1:mesh_size)*rhocore(1:mesh_size) &
1313 0 : & *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2
1314 0 : if (ii==4) ff(1:mesh_size)=rhoarr1_ep(1:mesh_size)*rhocore(1:mesh_size) &
1315 0 : & *pawrad(itypat)%rad(1:mesh_size)**2
1316 : end if
1317 0 : call simp_gen(intg,ff,pawrad(itypat))
1318 0 : intg=intg*pawang%angwgth(ipt)*four_pi
1319 0 : if (ii==1) lambda_paw =lambda_paw +lsign(iloop)*intg
1320 0 : if (ii==2) lambda_paw_ipm =lambda_paw_ipm +lsign(iloop)*intg
1321 0 : if (ii==3) lambda_core_paw =lambda_core_paw +lsign(iloop)*intg
1322 0 : if (ii==4) lambda_core_paw_ipm=lambda_core_paw_ipm+lsign(iloop)*intg
1323 : end do
1324 0 : ABI_FREE(ff)
1325 : end do ! ipt
1326 0 : ABI_FREE(gamma)
1327 0 : ABI_FREE(rhoarr1)
1328 0 : ABI_FREE(rhoarr1_ep)
1329 0 : if (option/=1) then
1330 0 : ABI_FREE(rhoarr1_j)
1331 : end if
1332 :
1333 : ! Second formalism: use (l,m) moments for densities
1334 : else if (dtset%pawxcdev/=0) then
1335 :
1336 : ! Build densities
1337 7648 : ABI_MALLOC(gammam,(mesh_size,2,lm_size))
1338 7648 : ABI_MALLOC(rhotot,(mesh_size,lm_size))
1339 5736 : ABI_MALLOC(rhotot_ep,(mesh_size,lm_size))
1340 3824 : ABI_MALLOC(rhosph,(mesh_size))
1341 3824 : ABI_MALLOC(rhosph_ep,(mesh_size))
1342 1912 : if (option/=1) then
1343 4416 : ABI_MALLOC(rhotot_j,(mesh_size,lm_size))
1344 2944 : ABI_MALLOC(rhosph_j,(mesh_size))
1345 : end if
1346 535816 : if (usecore==0) rhocore(:)=zero
1347 1912 : if (iloop==1) then
1348 10112456 : rhotot (:,:)=rho1 (:,:,1)
1349 10112456 : rhotot_ep(:,:)=rho1_ep(:,:,1)
1350 7530140 : if (option/=1) rhotot_j (:,:)=rho1_j (:,:,1)
1351 1123500 : if (usecore==1) rhocore(:)=pawtab(itypat)%coredens(:)
1352 : else
1353 : if (include_nhat_in_gamma) then
1354 : rhotot (:,:)=trho1 (:,:,1)+nhat1 (:,:,1)
1355 : rhotot_ep(:,:)=trho1_ep(:,:,1)+nhat1_ep(:,:,1)
1356 : if (option/=1) rhotot_j (:,:)=trho1_j (:,:,1)+nhat1_j (:,:,1)
1357 : else
1358 10112456 : rhotot (:,:)=trho1 (:,:,1)
1359 10112456 : rhotot_ep(:,:)=trho1_ep(:,:,1)
1360 7530140 : if (option/=1) rhotot_j (:,:)=trho1_j (:,:,1)
1361 : end if
1362 589988 : if (usecore==1) rhocore(:)=pawtab(itypat)%tcoredens(:,1)
1363 : end if
1364 2247000 : rhosph (:)=rhotot (:,1)/sqfpi
1365 2247000 : rhosph_ep(:)=rhotot_ep(:,1)/sqfpi
1366 1673592 : if (option/=1) rhosph_j (:)=rhotot_j(:,1)/sqfpi
1367 : ! Make spherical densities positive
1368 1912 : if (electronpositron%particle==EP_ELECTRON) then
1369 84 : call mkdenpos(iwarnp,mesh_size,1,1,rhosph ,dtset%xc_denpos)
1370 84 : call mkdenpos(iwarn ,mesh_size,1,1,rhosph_ep,dtset%xc_denpos)
1371 1828 : else if (electronpositron%particle==EP_POSITRON) then
1372 1828 : call mkdenpos(iwarn ,mesh_size,1,1,rhosph ,dtset%xc_denpos)
1373 1828 : call mkdenpos(iwarnp,mesh_size,1,1,rhosph_ep,dtset%xc_denpos)
1374 1828 : if (option/=1) then
1375 1472 : call mkdenpos(iwarnp,mesh_size,1,1,rhosph_j,dtset%xc_denpos)
1376 : end if
1377 : end if
1378 : ! Need gradients of electronic densities for GGA
1379 5736 : ABI_MALLOC(grhoe2,(ngr))
1380 3824 : ABI_MALLOC(grhocore2,(ngr))
1381 1912 : if (ngr>0) then
1382 0 : if (electronpositron%particle==EP_ELECTRON) then
1383 0 : call nderiv_gen(grhoe2,rhosph_ep,pawrad(itypat))
1384 0 : else if (electronpositron%particle==EP_POSITRON) then
1385 0 : call nderiv_gen(grhoe2,rhosph,pawrad(itypat))
1386 : end if
1387 0 : grhoe2(:)=grhoe2(:)**2
1388 0 : if (usecore==1) then
1389 0 : call nderiv_gen(grhocore2,rhocore,pawrad(itypat))
1390 0 : grhocore2(:)=grhocore2(:)**2
1391 : end if
1392 : end if
1393 : ! Compute Gamma for (rho-,rho+),
1394 : ! (rho- +drho-,rho+), (rho- -drho-,rho+),
1395 : ! (rho-,rho+ +drho+), (rho-,rho+ -drho+),
1396 : ! (rho- +drho-,rho+ +drho+), (rho- -drho-,rho+ -drho+)
1397 : ! Do a seven steps loop
1398 7648 : ABI_MALLOC(gam_,(mesh_size,2,7))
1399 3824 : ABI_MALLOC(rho_,(mesh_size))
1400 3824 : ABI_MALLOC(rho_ep_,(mesh_size))
1401 3824 : ABI_MALLOC(rhocor_,(mesh_size))
1402 3824 : ABI_MALLOC(grho2_,(ngr))
1403 3824 : ABI_MALLOC(grhocor2_,(ngr))
1404 15296 : do ii=1,7
1405 : ! Apply delta to get perturbed densities
1406 43449672 : rho_(:)=rhosph(:);rho_ep_(:)=rhosph_ep(:);if (usecore==1) rhocor_(:)=rhocore(:)
1407 13384 : if (ngr>0) grho2_(:)=grhoe2(:)
1408 13384 : if (ngr>0) grhocor2_(:)=grhocore2(:)
1409 13384 : if (ii==2.or.ii==4.or.ii==6) fact=(one+delta)
1410 13384 : if (ii==3.or.ii==5.or.ii==7) fact=(one-delta)
1411 13384 : fact2=fact**2
1412 13384 : if (ii==2.or.ii==3.or.ii==6.or.ii==7) then
1413 8988000 : rho_(:)=fact*rho_(:)
1414 7648 : if (electronpositron%particle==EP_POSITRON) then
1415 7312 : if (ngr>0) grho2_(:)=fact2*grho2_(:)
1416 6440192 : if (usecore==1)rhocor_(:)=fact*rhocor_(:)
1417 7312 : if (ngr>0.and.usecore==1) grhocor2_(:)=fact2*grhocor2_(:)
1418 : end if
1419 : end if
1420 13384 : if (ii==4.or.ii==5.or.ii==6.or.ii==7) then
1421 8988000 : rho_ep_(:)=fact*rho_ep_(:)
1422 7648 : if (electronpositron%particle==EP_ELECTRON) then
1423 336 : if (ngr>0) grho2_(:)=fact2*grho2_(:)
1424 419840 : if (usecore==1)rhocor_(:)=fact*rhocor_(:)
1425 336 : if (ngr>0.and.usecore==1) grhocor2_(:)=fact2*grhocor2_(:)
1426 : end if
1427 : end if
1428 : ! Compute gamma for these perturbed densities
1429 15296 : if (option==1.or.option==2) then
1430 8232 : if (electronpositron%particle==EP_ELECTRON) then
1431 588 : call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma(igam),ngr,mesh_size,rhocor_,rho_ep_,rho_,usecore)
1432 7644 : else if (electronpositron%particle==EP_POSITRON) then
1433 7644 : call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma(igam),ngr,mesh_size,rhocor_,rho_,rho_ep_,usecore)
1434 : end if
1435 : else
1436 82025664 : gam_(:,:,:)=one
1437 : end if
1438 : end do ! end loop ii=1,7
1439 :
1440 1912 : ABI_FREE(rhocor_)
1441 1912 : ABI_FREE(grho2_)
1442 1912 : ABI_FREE(grhocor2_)
1443 1912 : ABI_FREE(grhoe2)
1444 1912 : ABI_FREE(grhocore2)
1445 3855220 : rho_ (:)=rhosph (:);if (electronpositron%particle==EP_POSITRON.and.usecore==1) rho_ (:)=rho_ (:)+rhocore(:)
1446 2351876 : rho_ep_(:)=rhosph_ep(:);if (electronpositron%particle==EP_ELECTRON.and.usecore==1) rho_ep_(:)=rho_ep_(:)+rhocore(:)
1447 : ! Compute numerical first and second derivatives of Gamma
1448 : ! d1gam(1) = dgam/drho+ (particle=ELECTRON), dgam/drho- (particle=POSITRON)
1449 : ! d1gam(2) = dgam/drho- (particle=ELECTRON), dgam/drho+ (particle=POSITRON)
1450 7648 : ABI_MALLOC(d1gam,(mesh_size,2,2))
1451 8993736 : d1gam(:,:,:)=zero
1452 2247000 : do ir=1,mesh_size
1453 2245088 : if (rho_ (ir)>tol14) d1gam(ir,1,1)=(gam_(ir,1,2)-gam_(ir,1,3))*half/(delta*rho_ (ir))
1454 2245088 : if (rhosph (ir)>tol14) d1gam(ir,2,1)=(gam_(ir,2,2)-gam_(ir,2,3))*half/(delta*rhosph (ir))
1455 2245088 : if (rho_ep_ (ir)>tol14) d1gam(ir,1,2)=(gam_(ir,1,4)-gam_(ir,1,5))*half/(delta*rho_ep_ (ir))
1456 2247000 : if (rhosph_ep(ir)>tol14) d1gam(ir,2,2)=(gam_(ir,2,4)-gam_(ir,2,5))*half/(delta*rhosph_ep(ir))
1457 : end do
1458 :
1459 : ! d2gam(1) = d2gam/drho+_drho+ (particle=ELECTRON), dgam/drho-_drho- (particle=POSITRON)
1460 : ! d2gam(2) = d2gam/drho-_drho+ (particle=ELECTRON), dgam/drho+_drho- (particle=POSITRON)
1461 : ! d2gam(3) = d2gam/drho-_drho- (particle=ELECTRON), dgam/drho+_drho+ (particle=POSITRON)
1462 7648 : ABI_MALLOC(d2gam,(mesh_size,2,3))
1463 13489648 : d2gam(:,:,:)=zero
1464 2247000 : do ir=1,mesh_size
1465 2245088 : if (rho_ (ir)>tol14) d2gam(ir,1,1)=(gam_(ir,1,2)+gam_(ir,1,3)-two*gam_(ir,1,1))/(delta*rho_ (ir))**2
1466 2245088 : if (rhosph(ir)>tol14) d2gam(ir,2,1)=(gam_(ir,2,2)+gam_(ir,2,3)-two*gam_(ir,2,1))/(delta*rhosph(ir))**2
1467 2245088 : if (rho_ep_(ir)>tol14) then
1468 2245088 : d2gam(ir,1,3)=(gam_(ir,1,4)+gam_(ir,1,5)-two*gam_(ir,1,1))/(delta*rho_ep_(ir))**2
1469 2245088 : if (rho_(ir)>tol14) then
1470 : d2gam(ir,1,2)=(gam_(ir,1,6)+gam_(ir,1,7)+two*gam_(ir,1,1) &
1471 : & -gam_(ir,1,2)-gam_(ir,1,3)-gam_(ir,1,4)-gam_(ir,1,5)) &
1472 2245088 : & *half/(delta*rho_(ir))/(delta*rho_ep_(ir))
1473 : end if
1474 : end if
1475 2247000 : if (rhosph_ep(ir)>tol14) then
1476 2245088 : d2gam(ir,2,3)=(gam_(ir,2,4)+gam_(ir,2,5)-two*gam_(ir,2,1))/(delta*rhosph_ep(ir))**2
1477 2245088 : if (rhosph(ir)>tol14) then
1478 : d2gam(ir,2,2)=(gam_(ir,2,6)+gam_(ir,2,7)+two*gam_(ir,2,1) &
1479 : & -gam_(ir,2,2)-gam_(ir,2,3)-gam_(ir,2,4)-gam_(ir,2,5)) &
1480 2245088 : & *half/(delta*rhosph(ir))/(delta*rhosph_ep(ir))
1481 : end if
1482 : end if
1483 : end do
1484 1912 : ABI_FREE(rho_)
1485 1912 : ABI_FREE(rho_ep_)
1486 : ! Compute useful sums of densities
1487 5736 : ABI_MALLOC(v1sum,(mesh_size,3))
1488 1912 : if ( dtset%pawxcdev>=2) then
1489 0 : ABI_MALLOC(v2sum,(mesh_size,lm_size,3))
1490 : else
1491 1912 : ABI_MALLOC(v2sum,(0,0,0))
1492 : end if
1493 4494000 : rhotot(:,1)=sqfpi*rhosph(:);rhotot_ep(:,1)=sqfpi*rhosph_ep(:)
1494 : call pawxcsum(1,1,1,lmselect,lmselect_ep,lm_size,mesh_size,3,dtset%pawxcdev,&
1495 1912 : & pawang,rhotot,rhotot_ep,v1sum,v2sum)
1496 : ! Compute final development of gamma moments
1497 40465120 : gammam(:,:,:)=zero
1498 4495912 : gammam(:,:,1)=gam_(:,:,1)*sqfpi
1499 : gammam(:,1,1)=gammam(:,1,1)+(d2gam(:,1,2)*v1sum(:,2) &
1500 2247000 : & +half*(d2gam(:,1,1)*v1sum(:,1)+d2gam(:,1,3)*v1sum(:,3)))/sqfpi
1501 : gammam(:,2,1)=gammam(:,2,1)+(d2gam(:,2,2)*v1sum(:,2) &
1502 2247000 : & +half*(d2gam(:,2,1)*v1sum(:,1)+d2gam(:,2,3)*v1sum(:,3)))/sqfpi
1503 17208 : do ilm=2,lm_size
1504 15296 : if (lmselect(ilm)) then
1505 17976000 : gammam(:,1,ilm)=gammam(:,1,ilm)+d1gam(:,1,1)*rhotot(:,ilm)
1506 17976000 : gammam(:,2,ilm)=gammam(:,2,ilm)+d1gam(:,2,1)*rhotot(:,ilm)
1507 : end if
1508 17208 : if (lmselect_ep(ilm)) then
1509 17976000 : gammam(:,1,ilm)=gammam(:,1,ilm)+d1gam(:,1,2)*rhotot_ep(:,ilm)
1510 17976000 : gammam(:,2,ilm)=gammam(:,2,ilm)+d1gam(:,2,2)*rhotot_ep(:,ilm)
1511 : end if
1512 : end do
1513 1912 : if (dtset%pawxcdev>1) then
1514 0 : do ilm=2,lm_size
1515 : gammam(:,1,ilm)=gammam(:,1,ilm)+d2gam(:,1,2)*v2sum(:,ilm,2) &
1516 0 : & +half*(d2gam(:,1,1)*v2sum(:,ilm,1)+d2gam(:,1,3)*v2sum(:,ilm,3))
1517 : gammam(:,2,ilm)=gammam(:,2,ilm)+d2gam(:,2,2)*v2sum(:,ilm,2) &
1518 0 : & +half*(d2gam(:,2,1)*v2sum(:,ilm,1)+d2gam(:,2,3)*v2sum(:,ilm,3))
1519 : end do
1520 : end if
1521 1912 : ABI_FREE(gam_)
1522 1912 : ABI_FREE(d1gam)
1523 1912 : ABI_FREE(d2gam)
1524 1912 : ABI_FREE(v1sum)
1525 1912 : ABI_FREE(v2sum)
1526 :
1527 : ! Compute contribution to annihilation rate
1528 :
1529 : ! In state dependent scheme for Doppler, replace electronic density with one state
1530 1912 : if (option/=1) then
1531 15059840 : rhotot(:,:) = rhotot_j(:,:)
1532 1673152 : rhosph (:) = rhosph_j (:)
1533 : end if
1534 :
1535 3824 : ABI_MALLOC(gg,(mesh_size,4))
1536 8989912 : gg=zero
1537 3824 : ABI_MALLOC(rhoarr1,(mesh_size))
1538 3824 : ABI_MALLOC(rhoarr2,(mesh_size))
1539 19120 : do ilm=1,lm_size
1540 173992 : do ilm1=1,lm_size
1541 172080 : if (lmselect(ilm1)) then
1542 20377872 : if (ilm1==1) rhoarr1(:)=sqfpi*rhosph(:)
1543 161784000 : if (ilm1/=1) rhoarr1(:)=rhotot(:,ilm1)
1544 1548720 : do ilm2=1,lm_size
1545 1548720 : if (lmselect_ep(ilm2)) then
1546 183400848 : if (ilm2==1) rhoarr2(:)=sqfpi*rhosph_ep(:)
1547 1456056000 : if (ilm2/=1) rhoarr2(:)=rhotot_ep(:,ilm2)
1548 1393848 : if (ilm1>=ilm2) then
1549 774360 : isel=pawang%gntselect(ilm,ilm2+ilm1*(ilm1-1)/2)
1550 : else
1551 619488 : isel=pawang%gntselect(ilm,ilm1+ilm2*(ilm2-1)/2)
1552 : end if
1553 1393848 : if (isel>0) then
1554 158696 : fact=pawang%realgnt(isel)
1555 186501000 : gg(:,1)=gg(:,1)+fact*rhoarr1(:)*rhoarr2(:)*gammam(:,1,ilm)
1556 186501000 : gg(:,2)=gg(:,2)+fact*rhoarr1(:)*rhoarr2(:)*gammam(:,2,ilm)
1557 : end if
1558 : end if
1559 : end do
1560 : end if
1561 : end do
1562 : end do
1563 1912 : ABI_FREE(rhoarr1)
1564 1912 : ABI_FREE(rhoarr2)
1565 1912 : if (electronpositron%particle==EP_ELECTRON) then
1566 840 : do ilm=1,lm_size
1567 993000 : if (lmselect(ilm)) gg(:,3)=gg(:,3)+rhotot(:,ilm)*rhocore(:)*gammam(:,1,ilm)
1568 : end do
1569 110324 : gg(:,4)=sqfpi*rhotot(:,1)*rhocore(:)
1570 1828 : else if (electronpositron%particle==EP_POSITRON) then
1571 18280 : do ilm=1,lm_size
1572 19231912 : if (lmselect_ep(ilm)) gg(:,3)=gg(:,3)+rhotot_ep(:,ilm)*rhocore(:)*gammam(:,1,ilm)
1573 : end do
1574 2136676 : gg(:,4)=sqfpi*rhotot_ep(:,1)*rhocore(:)
1575 : end if
1576 9560 : do ii=1,4
1577 8988000 : gg(1:mesh_size,ii)=gg(1:mesh_size,ii)*pawrad(itypat)%rad(1:mesh_size)**2
1578 7648 : call simp_gen(intg,gg(:,ii),pawrad(itypat))
1579 7648 : if (ii==1) lambda_paw =lambda_paw +lsign(iloop)*intg
1580 7648 : if (ii==2) lambda_paw_ipm =lambda_paw_ipm +lsign(iloop)*intg
1581 7648 : if (ii==3) lambda_core_paw =lambda_core_paw +lsign(iloop)*intg
1582 9560 : if (ii==4) lambda_core_paw_ipm=lambda_core_paw_ipm+lsign(iloop)*intg
1583 : end do
1584 1912 : ABI_FREE(gg)
1585 1912 : ABI_FREE(gammam)
1586 1912 : ABI_FREE(rhotot)
1587 1912 : ABI_FREE(rhotot_ep)
1588 1912 : ABI_FREE(rhosph)
1589 1912 : ABI_FREE(rhosph_ep)
1590 1912 : if (option/=1) then
1591 1472 : ABI_FREE(rhotot_j)
1592 1472 : ABI_FREE(rhosph_j)
1593 : end if
1594 :
1595 : end if ! dtset%pawxcdev
1596 :
1597 2868 : ABI_FREE(rhocore)
1598 :
1599 : end do ! iloop
1600 :
1601 956 : ABI_FREE(rho1)
1602 956 : ABI_FREE(trho1)
1603 956 : ABI_FREE(rho1_ep)
1604 956 : ABI_FREE(trho1_ep)
1605 956 : ABI_FREE(rho1_j)
1606 956 : ABI_FREE(trho1_j)
1607 956 : ABI_FREE(nhat1)
1608 956 : ABI_FREE(nhat1_ep)
1609 956 : ABI_FREE(nhat1_j)
1610 956 : ABI_FREE(lmselect)
1611 956 : ABI_FREE(lmselect_ep)
1612 1294 : ABI_FREE(lmselect_dum)
1613 :
1614 : end do ! iatom
1615 :
1616 : ! Reduction in case of distribution over atomic sites
1617 338 : if (mpi_enreg%nproc_atom>1) then
1618 0 : mpibuf(1)=lambda_paw ;mpibuf(2)=lambda_paw_ipm
1619 0 : mpibuf(3)=lambda_core_paw;mpibuf(4)=lambda_core_paw_ipm
1620 0 : call xmpi_sum(mpibuf,mpi_enreg%comm_atom,ierr)
1621 0 : lambda_paw=mpibuf(1) ;lambda_paw_ipm=mpibuf(2)
1622 0 : lambda_core_paw=mpibuf(3);lambda_core_paw_ipm=mpibuf(4)
1623 : end if
1624 :
1625 : ! Add plane-wave and PAW contributions to annihilation rates
1626 338 : lambda =lambda +lambda_paw
1627 338 : lambda_ipm =lambda_ipm +lambda_paw_ipm
1628 338 : lambda_core =lambda_core +lambda_core_paw
1629 338 : lambda_core_ipm=lambda_core_ipm+lambda_core_paw_ipm
1630 : end if ! dtset%usepaw
1631 :
1632 :
1633 : ! Convert into proper units and print
1634 : ! ---------------------------------------------------------------------------------------
1635 :
1636 : ! Sum valence and core contributions to annihilation rates
1637 340 : lambda =lambda +lambda_core
1638 340 : lambda_ipm =lambda_ipm +lambda_core_ipm
1639 340 : if (dtset%usepaw==1) then
1640 338 : lambda_paw =lambda_paw +lambda_core_paw
1641 338 : lambda_paw_ipm=lambda_paw_ipm+lambda_core_paw_ipm
1642 : end if
1643 :
1644 : ! Set annihilation rate in proper unit (picosec.)
1645 340 : units=pi*(one/InvFineStruct)**3/Time_Sec/1.e12_dp/electronpositron%posocc
1646 :
1647 340 : lambda =lambda *units
1648 340 : lambda_ipm =lambda_ipm *units
1649 340 : lambda_core =lambda_core *units
1650 340 : lambda_core_ipm=lambda_core_ipm*units
1651 340 : lifetime =one/lambda
1652 340 : lifetime_ipm=one/lambda_ipm
1653 340 : electronpositron%lambda=lambda
1654 340 : electronpositron%lifetime=lifetime
1655 340 : if (dtset%usepaw==1) then
1656 338 : lambda_paw =lambda_paw *units
1657 338 : lambda_paw_ipm =lambda_paw_ipm *units
1658 338 : lambda_core_paw =lambda_core_paw *units
1659 338 : lambda_core_paw_ipm=lambda_core_paw_ipm*units
1660 : end if
1661 340 : rate=lambda-lambda_core-lambda_paw+lambda_core_paw
1662 340 : rate_paw=lambda_paw-lambda_core_paw
1663 : ! Print life time and additional information
1664 669 : if (option==1) then
1665 28 : if (igam==1) then
1666 17 : write(msg,'(a,80("-"),2a)') ch10,ch10,' Results for electron-positron annihilation:'
1667 17 : call wrtout(ab_out,msg,'COLL')
1668 17 : call wrtout(std_out,msg,'COLL')
1669 : end if
1670 28 : if (ngamma>1.and.igam==1) then
1671 11 : write(msg,'(a,i1,a)') ch10,ngamma,&
1672 22 : & ' computations of positron lifetime have been performed (with different enhancement factors).'
1673 11 : call wrtout(ab_out,msg,'COLL')
1674 11 : call wrtout(std_out,msg,'COLL')
1675 : end if
1676 28 : if (ngamma>1) then
1677 22 : write(msg,'(2a,i1)') ch10,"########## Lifetime computation ",igam
1678 22 : call wrtout(ab_out,msg,'COLL')
1679 22 : call wrtout(std_out,msg,'COLL')
1680 : end if
1681 28 : if (abs(electronpositron%ixcpositron)==1) then
1682 22 : write(msg,'(4a)') ch10,' # Zero-positron density limit of Arponen and Pajanne provided by Boronski & Nieminen',&
1683 44 : & ch10,' Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)' ! [[cite:Boronski1986]]
1684 6 : else if (electronpositron%ixcpositron==11) then
1685 4 : write(msg,'(4a)') ch10,' # Zero-positron density limit of Arponen and Pajanne fitted by Sterne & Kaiser',&
1686 8 : & ch10,' Ref.: P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991)' ! [[cite:Sterne1991]]
1687 2 : else if (electronpositron%ixcpositron==2) then
1688 2 : write(msg,'(4a)') ch10,' # Electron-positron correlation provided by Puska, Seitsonen, and Nieminen',&
1689 4 : & ch10,' Ref: M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994)' ! [[cite:Puska1994]]
1690 0 : else if (electronpositron%ixcpositron==3) then
1691 0 : write(msg,'(8a)') ch10,' # Zero-positron density limit of Arponen and Pajanne provided by Boronski & Nieminen',&
1692 0 : & ch10,' + GGA corrections',&
1693 0 : & ch10,' Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)',& ! [[cite:Boronski1986]]
1694 0 : & ch10,' B. Barbiellini, M.J. Puska, T. Torsti and R.M.Nieminen, Phys. Rev. B 51, 7341 (1995)' ! [[cite:Barbiellini1995]]
1695 0 : else if (electronpositron%ixcpositron==31) then
1696 0 : write(msg,'(8a)') ch10,' # Zero-positron density limit of Arponen and Pajanne fitted by Sterne & Kaiser',&
1697 0 : & ch10,' + GGA corrections',&
1698 0 : & ch10,' Ref.: P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991)',& ! [[cite:Sterne1991]]
1699 0 : & ch10,' B. Barbiellini, M.J. Puska, T. Torsti and R.M. Nieminen, Phys. Rev. B 51, 7341 (1995)' ! [[cite:Barbiellini1995]]
1700 : end if
1701 28 : call wrtout(ab_out,msg,'COLL')
1702 28 : call wrtout(std_out, msg,'COLL')
1703 28 : if (igamma(igam)==0) then
1704 0 : write(msg,'(a)') ' # Enhancement factor set to one (test)'
1705 : else if (igamma(igam)==1) then
1706 11 : write(msg,'(3a)') ' # Enhancement factor of Boronski & Nieminen',&
1707 22 : & ch10,' Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)' ! [[cite:Boronski1986]]
1708 : else if (igamma(igam)==2) then
1709 11 : write(msg,'(3a)') ' # Enhancement factor of Boronski & Nieminen IN THE RPA LIMIT',&
1710 22 : & ch10,' Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)' ! [[cite:Boronski1986]]
1711 : else if (igamma(igam)==3) then
1712 4 : write(msg,'(3a)') ' # Enhancement factor of Sterne & Kaiser',&
1713 8 : & ch10,' Ref.: P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991)' ! [[cite:Sterne1991]]
1714 : else if (igamma(igam)==4) then
1715 2 : write(msg,'(3a)') ' # Enhancement factor of Puska, Seitsonen, and Nieminen',&
1716 4 : & ch10,' Ref.: M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994)' ! [[cite:Puska1994]]
1717 : end if
1718 28 : call wrtout(ab_out,msg,'COLL')
1719 28 : call wrtout(std_out,msg,'COLL')
1720 28 : write(msg, '(4(2a,es16.8))' ) ch10,&
1721 28 : & ' Positron lifetime (ps) =',lifetime ,ch10,&
1722 28 : & ' Positron lifetime with IPM for core elec. (ps) =',lifetime_ipm,ch10,&
1723 28 : & ' Annihilation rate (ns-1) =',lambda *1000._dp,ch10,&
1724 56 : & ' Annihilation rate with IPM for core elec. (ns-1) =',lambda_ipm*1000._dp
1725 28 : call wrtout(ab_out,msg,'COLL')
1726 28 : call wrtout(std_out,msg,'COLL')
1727 28 : write(msg,'(2a,5(2a,es16.8))' ) ch10,&
1728 28 : & ' Annihilation rate core/valence decomposition:',ch10,&
1729 28 : & ' Core contribution to ann.rate (ns-1) =', lambda_core *1000._dp,ch10,&
1730 28 : & ' Valence contribution to ann.rate (ns-1) =',(lambda-lambda_core) *1000._dp,ch10,&
1731 28 : & ' Core contribution to ann.rate with IPM (ns-1) =', lambda_core_ipm *1000._dp,ch10,&
1732 56 : & ' Valence contribution to ann.rate with IPM (ns-1) =',(lambda_ipm-lambda_core_ipm) *1000._dp
1733 28 : call wrtout(ab_out,msg,'COLL')
1734 28 : call wrtout(std_out,msg,'COLL')
1735 28 : if (dtset%usepaw==1) then
1736 26 : write(msg, '(2a,6(2a,es16.8))' ) ch10,&
1737 26 : & ' Annihilation rate PAW decomposition:',ch10,&
1738 26 : & ' Plane-wave contribution to ann.rate (ns-1) =',(lambda-lambda_paw)*1000._dp,ch10,&
1739 26 : & ' Plane-wave valence contribution to ann.rate (ns-1) =',(lambda-lambda_paw-lambda_core+lambda_core_paw)*1000._dp,ch10,&
1740 26 : & ' On-site core contribution to ann.rate (ns-1) =', lambda_core_paw*1000._dp,ch10,&
1741 26 : & ' On-site valence contribution to ann.rate (ns-1) =',(lambda_paw-lambda_core_paw)*1000._dp,ch10,&
1742 26 : & ' Plane-wave contribution to ann.rate with IPM (ns-1) =',(lambda_ipm-lambda_paw_ipm)*1000._dp,ch10,&
1743 52 : & ' Plane-wave core contrb. to ann.rate with IPM (ns-1) =',(lambda_core_ipm-lambda_core_paw_ipm)*1000._dp
1744 26 : call wrtout(ab_out,msg,'COLL')
1745 26 : call wrtout(std_out,msg,'COLL')
1746 : end if
1747 28 : if (dtset%usepaw==0.and.igam==ngamma) then ! These tests are not relevant with PAW
1748 1 : write(msg, '(2a,3(2a,es16.8))' ) ch10,&
1749 1 : & ' ########## Some checks, for testing purpose:',ch10,&
1750 1 : & ' Number of core electrons =',nbec,ch10,&
1751 1 : & ' Number of valence electrons =',nbev,ch10,&
1752 2 : & ' Number of positrons =',nbp
1753 1 : call wrtout(ab_out,msg,'COLL')
1754 1 : call wrtout(std_out,msg,'COLL')
1755 : end if
1756 : end if !end if option
1757 : end do ! Big loop on igam
1758 :
1759 329 : if (option==1) then
1760 17 : write(msg, '(3a)' ) ch10,' (*) IPM=Independent particle Model',ch10
1761 17 : call wrtout(ab_out,msg,'COLL')
1762 17 : call wrtout(std_out,msg,'COLL')
1763 : end if !end if option
1764 :
1765 : !Deallocate memory
1766 329 : ABI_FREE(igamma)
1767 329 : if (dtset%usepaw==1.and.(.not.include_nhat_in_gamma)) then
1768 328 : ABI_FREE(rhor_)
1769 328 : ABI_FREE(rhor_ep_)
1770 : end if
1771 :
1772 : DBG_EXIT("COLL")
1773 :
1774 970 : end subroutine poslifetime
1775 : !!***
1776 :
1777 : !!****f* ABINIT/posdoppler
1778 : !! NAME
1779 : !! posdoppler
1780 : !!
1781 : !! FUNCTION
1782 : !! Calculate the momentum distribution annihilating electrons-positron (Doppler broadening)
1783 : !!
1784 : !! INPUTS
1785 : !! cg(2,mcg)=planewave coefficients of wavefunctions.
1786 : !! cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk>
1787 : !! and each |p_lmn> non-local projector
1788 : !! Cryst<Crystal_structure> = Info on unit cell and its symmetries
1789 : !! dimcprj(natom)=array of dimensions of array cprj (not ordered)
1790 : !! dtfil <type(datafiles_type)>=variables related to files
1791 : !! | unpaw=unit number for temporary PAW files
1792 : !! dtset <type(dataset_type)>=all input variables for this dataset
1793 : !! | istwfk=input option=1 parameter that describes the storage of wfs
1794 : !! | mband=maximum number of bands
1795 : !! | mgfft=maximum size of 1D FFTs for the "coarse" grid
1796 : !! | mkmem=number of k points treated by this node.
1797 : !! | mpw=maximum dimensioned size of npw
1798 : !! | natom=number of atoms
1799 : !! | nband=number of bands at each k point
1800 : !! | ngfft=contain all needed information about 3D FFT (coarse grid)
1801 : !! | nkpt=number of k points
1802 : !! | nspden=number of spin-density components
1803 : !! | nspinor=number of spinorial components of the wavefunctions
1804 : !! | nsppol=1 for unpolarized, 2 for spin-polarized
1805 : !! | usepaw=flag for PAW
1806 : !! | gpu_option=GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
1807 : !! | wtk(=weights associated with various k points
1808 : !! filpsp(ntypat)=name(s) of the pseudopotential file(s)
1809 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
1810 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
1811 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
1812 : !! mpi_enreg= information about MPI parallelization
1813 : !! my_natom=number of atoms treated by current processor
1814 : !! n3xccc= dimension of the xccc3d array (0 or nfft).
1815 : !! nfft= number of FFT grid points
1816 : !! ngfft(18)= contain all needed information about 3D FFT
1817 : !! nhat(nfft,nspden)=charge compensation density (content depends on electronpositron%particle)
1818 : !! npwarr(nkpt)=number of planewaves in basis at this k point
1819 : !! occ(mband*nkpt*nsppol)=occupancies for each band and k point
1820 : !! pawang <type(pawang)>=paw angular mesh and related data
1821 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
1822 : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies
1823 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
1824 : !! rhor(nfft,nspden)=total electron/positron density (content depends on electronpositron%particle)
1825 : !! xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
1826 : !!
1827 : !! OUTPUT
1828 : !!
1829 : !! SIDE EFFECTS
1830 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
1831 : !!
1832 : !! TODO
1833 : !! print a warning if the core wave function is not localized in the PAW sphere
1834 : !! implement PAW on-site contribution for state-independent scheme
1835 : !!
1836 : !! SOURCE
1837 :
1838 : !Macro to go from row-column indexing to combined indexing
1839 : #define RCC(glmn,hlmn) max(glmn,hlmn)*(max(glmn,hlmn)-1)/2+min(glmn,hlmn)
1840 : !Macro to go from l,m angular momentum indexing to combined indexing
1841 : #define LMC(lval,mval) lval*lval+lval+mval+1
1842 :
1843 3 : subroutine posdoppler(cg,cprj,Crystal,dimcprj,dtfil,dtset,electronpositron,&
1844 3 : & filpsp,kg,mcg,mcprj,mpi_enreg,my_natom,&
1845 3 : & n3xccc,nfft,ngfft,nhat,npwarr,occ,pawang,pawrad,&
1846 3 : & pawrhoij,pawtab,rhor,xccc3d)
1847 :
1848 : !Arguments ------------------------------------
1849 : !scalars
1850 : integer,intent(in) :: mcg,mcprj,my_natom,n3xccc,nfft
1851 : type(crystal_t) :: Crystal
1852 : type(datafiles_type),intent(in) :: dtfil
1853 : type(dataset_type),intent(in) :: dtset
1854 : type(electronpositron_type),pointer :: electronpositron
1855 : type(MPI_type),intent(inout) :: mpi_enreg
1856 : type(pawang_type),intent(in) :: pawang
1857 : !arrays
1858 : integer,intent(in) :: dimcprj(dtset%natom)
1859 : integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),ngfft(18),npwarr(dtset%nkpt)
1860 : real(dp),intent(in) :: nhat(nfft,dtset%nspden*dtset%usepaw),xccc3d(n3xccc)
1861 : real(dp),intent(in),target :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
1862 : real(dp),intent(in),target :: rhor(nfft,dtset%nspden)
1863 : real(dp),intent(inout),target :: cg(2,mcg)
1864 : character(len=fnlen),intent(in) :: filpsp(dtset%ntypat)
1865 : type(pawcprj_type),target :: cprj(dtset%natom,mcprj)
1866 : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw)
1867 : type(pawrhoij_type),intent(in),target :: pawrhoij(my_natom*dtset%usepaw)
1868 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
1869 :
1870 : !Local variables-------------------------------
1871 : !scalars
1872 : integer :: accessfil,basis_size,bandpp,bdtot_index,bdtot_index_pos,blocksize,cplex,cplex_rhoij
1873 : integer :: glmij,i0lmn,i1,i2,i3,iat,iatm,iatom
1874 : integer :: ib,ib_cprj,ib_cprj_pos,ib_pos,ibg,ibg_pos
1875 : integer :: iblock,iblock_pos,ibpp,ibpp_pos
1876 : integer :: icg,icg_pos,id1,id2,id3,ierr,ig1,ig2,ig3,igamma,ii,ikg,ikg_pos,ikpt
1877 : integer :: ikpt_pos,il,ilm,ilmn,iln,indx,indx0,iorder_cprj,iproc,ir,isppol,isppol_pos,istwf_k
1878 : integer :: istwf_k_pos,itypat,iwarn,iwavef,iwavef_pos,j2,j3,jj,jkpt,jl,jlm,jlmn,jln
1879 : integer :: klm,kln,klmn,l_size,l_size_max,ll,llmax,llmin,lm,lmn_size,lmn_size_c,lmn2_size
1880 : integer :: mband_cprj,mband_cprj_pos,mcg_pos
1881 : integer :: mcprj_k,mcprj_k_pos,me_band,me_fft,me_kpt,me_kptband
1882 : integer :: mesh_size,meshsz,mm,my_ngrid,my_nspinor,my_nsppol,my_n2,n1,n2,n3,n4,n5,n6
1883 : integer :: nband_cprj_eff_pos,nband_cprj_k,nband_cprj_k_pos
1884 : integer :: nband_eff_pos,nband_k,nband_k_pos
1885 : integer :: nblock_band,nblock_band_eff_pos,nkpt
1886 : integer :: nproc_band,nproc_fft,nproc_spkpt,nproc_kptband,npw_k,npw_k_pos
1887 : integer :: nspden_rhoij,option,tag,unit_doppler
1888 : integer :: tim_fourdp=0,tim_fourwf=-36
1889 : integer :: ylmr_normchoice,ylmr_npts,ylmr_option
1890 : logical,parameter :: include_nhat_in_gamma=.false.,state_dependent=.true.
1891 : logical,parameter :: kgamma_only_positron=.true.,wf_conjugate=.false.
1892 : logical :: cprj_paral_band,mykpt,mykpt_pos,use_timerev,use_zeromag
1893 : real(dp) :: arg,bessarg,cpi,cpr,cp11,cp12,cp21,cp22,gammastate,intg
1894 : real(dp) :: lambda_v1,lambda_v2,lambda_core,lambda_pw,occ_el,occ_pos
1895 : real(dp) :: pnorm,pr,rate,rate_ipm,ratec,ratec_ipm,rate_paw,rate_paw_ipm
1896 : real(dp) :: scale_,units_,weight,weight_pos,wf_fact,wtk_k,wtk_k_pos,vec
1897 : character(len=fnlen) :: filename_dop
1898 : character(len=1000) :: msg
1899 : type(bandfft_kpt_type),pointer :: bandfft_kpt_el,bandfft_kpt_pos
1900 3 : type(MPI_type) :: mpi_enreg_seq
1901 : type(wffile_type) :: wff
1902 : type(atomorb_type), pointer :: atm=> null()
1903 : !arrays
1904 3 : integer,allocatable :: gbound(:,:),gbound_pos(:,:),kg_k(:,:),kg_k_pos(:,:)
1905 3 : integer,allocatable :: lmncmax(:),my_ffttab(:),my_gridtab(:),nphicor(:)
1906 3 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
1907 3 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1908 3 : logical,allocatable :: have_intc(:,:,:),have_rad(:,:)
1909 : real(dp) :: buf(4),contrib(2),cp(2),cp_pos(2),expipr(2),pbn(3),pcart(3)
1910 : real(dp) :: radsumnfftc(2),ylmgr(1,1,0),ylmr_nrm(1)
1911 3 : real(dp),allocatable :: cwaveg(:,:),cwaveg_pos(:,:),cwaver(:),cwaver_pos(:),cwaver_pos_block(:)
1912 3 : real(dp),allocatable :: cg_k_pos(:,:),cwaveaug(:,:,:,:),cwaveaug_pos(:,:,:,:)
1913 3 : real(dp),allocatable :: denpot_dum(:,:,:),ff(:),fofgout_dum(:,:)
1914 3 : real(dp),allocatable :: gamma(:,:),intc(:,:,:),j_bessel(:,:),jbes(:),mpibuf(:,:)
1915 3 : real(dp),allocatable :: occ_k(:),occ_k_pos(:),pcart_k(:,:)
1916 3 : real(dp),allocatable :: radint1(:,:),radint2(:,:),radint3(:,:)
1917 3 : real(dp),allocatable :: radsumnfft1(:,:),radsumnfft2(:,:),radsumnfft3(:,:)
1918 3 : real(dp),allocatable :: rho_contrib(:),rho_contrib_g(:,:)
1919 3 : real(dp),allocatable :: rho_contrib_paw1(:,:),rho_contrib_paw2(:,:),rho_contrib_paw3(:,:)
1920 3 : real(dp),allocatable :: rho_moment_v1(:,:),rho_moment_v2(:,:)
1921 3 : real(dp),allocatable :: rho_moment_core(:,:),rho_moment_k(:),rho_moment_k2(:)
1922 3 : real(dp),allocatable :: rho_pw(:,:),rhor_dop_el(:)
1923 3 : real(dp),allocatable :: rhocorej(:),rhoe(:,:),rhop(:,:),ylmp(:)
1924 3 : real(dp),pointer :: cg_pos_ptr(:,:),cg_ptr(:,:),occ_ptr(:),occ_pos_ptr(:)
1925 3 : real(dp),pointer :: rhor_(:,:),rhor_ep_(:,:)
1926 : complex(dp) :: ifac ! (-i)^L mod 4
1927 : complex(dp),dimension(0:3) :: ilfac(0:3)=(/(1.0,0.0),(0.0,-1.0),(-1.0,0.0),(0.0,1.0)/)
1928 3 : type(coeff1_type),allocatable :: gammastate_c(:)
1929 3 : type(coeffi2_type),allocatable :: indlmncor(:)
1930 3 : type(coeff2_type),allocatable :: phicor(:)
1931 3 : type(coeff6_type),allocatable :: radsum1(:),radsum2(:),radsum3(:)
1932 3 : type(coeff7_type),allocatable :: radsumc(:)
1933 3 : type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_k_pos(:,:),cprj_pos(:,:)
1934 3 : type(pawcprj_type),pointer :: cprj_pos_ptr(:,:),cprj_ptr(:,:)
1935 3 : type(pawrhoij_type),allocatable :: pawrhoij_dop_el(:)
1936 3 : type(pawrhoij_type),pointer :: pawrhoij_ptr(:),pawrhoij_all(:),pawrhoij_ep_all(:)
1937 :
1938 : ! *************************************************************************
1939 :
1940 : DBG_ENTER("COLL")
1941 :
1942 : !Compatibility tests
1943 3 : if (.not.associated(electronpositron)) then
1944 0 : ABI_BUG('electronpositron variable must be associated!')
1945 : end if
1946 3 : if (allocated(mpi_enreg%proc_distrb)) then
1947 6 : do isppol=1,dtset%nsppol
1948 26 : do ikpt=1,dtset%nkpt
1949 20 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1950 259 : if (any(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol)/=mpi_enreg%proc_distrb(ikpt,1,isppol))) then
1951 0 : msg='proc_distrib cannot be distributed over bands!'
1952 0 : ABI_BUG(msg)
1953 : end if
1954 : end do
1955 : end do
1956 : end if
1957 3 : if (dtset%nspinor==2) then
1958 0 : msg='Doppler broadening not available for spinorial wave functions (nspinor=2)!'
1959 0 : ABI_BUG(msg)
1960 : end if
1961 3 : if (mcprj==0) then
1962 0 : msg='<p|Psi> (cprj) datastructure must be kept in memory (see pawusecp input keyword)!'
1963 0 : ABI_BUG(msg)
1964 : end if
1965 3 : if (dtset%usepaw==0) then
1966 0 : write(msg,'(5a)') 'Momentum distribution of annihilating electron-positron pairs',ch10,&
1967 0 : & 'in the Norm-conserving Pseudopotential formalism is incomplete!',ch10,&
1968 0 : & 'No core contribution is included.'
1969 0 : ABI_WARNING(msg)
1970 : end if
1971 23 : if (any(dtset%nband(:)/=dtset%nband(1))) then
1972 0 : write(msg,'(a)') 'Number of bands has to be the same for all k-points!'
1973 0 : ABI_BUG(msg)
1974 : end if
1975 3 : if (dtset%usepaw==1) then
1976 3 : if (size(pawrhoij)/=mpi_enreg%my_natom) then
1977 0 : write(msg,'(a)') 'wrong size for pawrhoij! '
1978 0 : ABI_BUG(msg)
1979 : end if
1980 : end if
1981 :
1982 : !Various initializations
1983 3 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1984 3 : n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6)
1985 3 : id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2
1986 3 : iorder_cprj=0 ; cplex=2 ; iwarn=1
1987 3 : wf_fact=one;if (wf_conjugate) wf_fact=-one
1988 3 : nkpt=dtset%nkpt
1989 :
1990 : !Manage kpt/spin parallelism
1991 9 : ABI_MALLOC(my_gridtab,(nkpt))
1992 23 : my_gridtab=0
1993 23 : do ii=1,nkpt
1994 23 : if (any(mpi_enreg%my_isppoltab(:)==1)) my_gridtab(ii)=mpi_enreg%my_kpttab(ii)
1995 : end do
1996 23 : my_ngrid=count(my_gridtab(:)/=0)
1997 9 : my_nsppol=sum(mpi_enreg%my_isppoltab(:))
1998 3 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
1999 :
2000 : !Parallel settings
2001 3 : if (mpi_enreg%paral_kgb/=0) then
2002 0 : nproc_spkpt=mpi_enreg%nproc_spkpt
2003 0 : nproc_band=mpi_enreg%nproc_band
2004 0 : nproc_fft=mpi_enreg%nproc_fft
2005 0 : nproc_kptband=xmpi_comm_size(mpi_enreg%comm_kptband)
2006 0 : me_kpt=mpi_enreg%me_kpt
2007 0 : me_band=mpi_enreg%me_band
2008 0 : me_fft=mpi_enreg%me_fft
2009 0 : me_kptband=xmpi_comm_rank(mpi_enreg%comm_kptband)
2010 0 : bandpp=mpi_enreg%bandpp
2011 0 : my_n2=n2/nproc_fft
2012 0 : accessfil=IO_MODE_FORTRAN;if(nproc_fft>1)accessfil=IO_MODE_MPI
2013 : else
2014 3 : nproc_spkpt=mpi_enreg%nproc_spkpt
2015 3 : nproc_band=1;nproc_fft=1
2016 3 : nproc_kptband=nproc_spkpt
2017 3 : me_band=0;me_fft=0
2018 3 : me_kpt=mpi_enreg%me_kpt
2019 3 : me_kptband=me_kpt
2020 3 : bandpp=1 ; my_n2=n2
2021 3 : accessfil=IO_MODE_FORTRAN
2022 : end if
2023 3 : blocksize=nproc_band*bandpp
2024 3 : nblock_band=dtset%nband(1)/blocksize
2025 :
2026 : !Select density according to nhat choice0
2027 3 : if (dtset%usepaw==0.or.include_nhat_in_gamma) then
2028 0 : rhor_ => rhor
2029 0 : rhor_ep_ => electronpositron%rhor_ep
2030 : else
2031 12 : ABI_MALLOC(rhor_,(nfft,dtset%nspden))
2032 9 : ABI_MALLOC(rhor_ep_,(nfft,dtset%nspden))
2033 132741 : rhor_=rhor-nhat
2034 132741 : rhor_ep_=electronpositron%rhor_ep-electronpositron%nhat_ep
2035 : end if
2036 :
2037 : !Select type(s) of enhancement factor
2038 3 : igamma=0
2039 3 : if (electronpositron%ixcpositron==-1) igamma=0
2040 : if (electronpositron%ixcpositron== 1) igamma=2
2041 : if (electronpositron%ixcpositron== 2) igamma=4
2042 : if (electronpositron%ixcpositron== 3) igamma=2
2043 : if (electronpositron%ixcpositron==11) igamma=3
2044 3 : if (electronpositron%ixcpositron==31) igamma=3
2045 :
2046 : !Select electronic and positronic states
2047 3 : if (electronpositron%particle==EP_ELECTRON) then !we should not be in this case
2048 0 : cg_ptr => electronpositron%cg_ep
2049 0 : cprj_ptr => electronpositron%cprj_ep
2050 0 : occ_ptr => electronpositron%occ_ep
2051 0 : pawrhoij_ptr => electronpositron%pawrhoij_ep
2052 0 : cg_pos_ptr => cg
2053 0 : cprj_pos_ptr => cprj
2054 0 : occ_pos_ptr => occ
2055 : end if
2056 3 : if (electronpositron%particle==EP_POSITRON) then
2057 3 : cg_ptr => cg
2058 3 : cprj_ptr => cprj
2059 3 : occ_ptr => occ
2060 3 : pawrhoij_ptr => pawrhoij
2061 3 : cg_pos_ptr => electronpositron%cg_ep
2062 3 : cprj_pos_ptr => electronpositron%cprj_ep
2063 3 : occ_pos_ptr => electronpositron%occ_ep
2064 : end if
2065 :
2066 : !Determine if cprj datastructures are distributed over bands
2067 3 : mband_cprj=size(cprj_ptr,2)/(my_nspinor*dtset%mkmem*dtset%nsppol)
2068 3 : mband_cprj_pos=size(cprj_pos_ptr,2)/(my_nspinor*dtset%mkmem*dtset%nsppol)
2069 3 : cprj_paral_band=(mband_cprj<dtset%mband)
2070 :
2071 : !Get the distrib associated with the fft_grid
2072 3 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2073 :
2074 : !===============================================================================
2075 : !================ Calculate the PAW on-site constants ==========================
2076 :
2077 3 : if (dtset%usepaw==1) then
2078 :
2079 3 : ylmr_normchoice = 0 ! input to initylmr are normalized
2080 3 : ylmr_npts = 1 ! only 1 point to compute in initylmr
2081 3 : ylmr_nrm(1) = one ! weight of normed point for initylmr
2082 3 : ylmr_option = 1 ! compute only ylm's in initylmr
2083 :
2084 : !Prepare radial integral for PAW correction for each atom type
2085 12 : ABI_MALLOC(radsum1,(dtset%ntypat))
2086 9 : ABI_MALLOC(radsum2,(dtset%ntypat))
2087 9 : ABI_MALLOC(radsum3,(dtset%ntypat))
2088 12 : ABI_MALLOC(radsumc,(dtset%ntypat))
2089 :
2090 12 : ABI_MALLOC(indlmncor,(dtset%ntypat))
2091 12 : ABI_MALLOC(phicor,(dtset%ntypat))
2092 17 : ABI_MALLOC(gammastate_c,(dtset%natom))
2093 9 : ABI_MALLOC(nphicor,(dtset%ntypat))
2094 6 : ABI_MALLOC(lmncmax,(dtset%ntypat))
2095 :
2096 : ! Reading of core wave functions
2097 3 : if (mpi_enreg%me_cell==0) then
2098 3 : ABI_MALLOC(atm,)
2099 6 : do itypat=1,dtset%ntypat
2100 3 : call pawpsp_init_core(atm,psp_filename=trim(filpsp(itypat)),radmesh=pawrad(itypat))
2101 12 : ABI_MALLOC(indlmncor(itypat)%value,(size(atm%indlmn(:,1)),atm%lmn_size))
2102 12 : ABI_MALLOC(phicor(itypat)%value,(atm%mesh_size,atm%ln_size))
2103 83 : indlmncor(itypat)%value=atm%indlmn
2104 3 : lmncmax(itypat)=atm%lmn_size
2105 3 : nphicor(itypat)=atm%ln_size
2106 9829 : phicor(itypat)%value=atm%phi(:,:,1)
2107 6 : call destroy_atomorb(atm)
2108 : end do
2109 3 : ABI_FREE(atm)
2110 : end if
2111 3 : if (mpi_enreg%nproc_cell>1) then
2112 0 : call xmpi_bcast(indlmncor,0,mpi_enreg%comm_cell,ierr)
2113 0 : call xmpi_bcast(phicor,0,mpi_enreg%comm_cell,ierr)
2114 0 : call xmpi_bcast(nphicor,0,mpi_enreg%comm_cell,ierr)
2115 0 : call xmpi_bcast(lmncmax,0,mpi_enreg%comm_cell,ierr)
2116 : end if
2117 :
2118 6 : do itypat=1,dtset%ntypat
2119 :
2120 3 : mesh_size = pawtab(itypat)%mesh_size
2121 3 : l_size = pawtab(itypat)%l_size
2122 3 : lmn_size = pawtab(itypat)%lmn_size
2123 3 : lmn2_size = pawtab(itypat)%lmn2_size
2124 3 : basis_size = pawtab(itypat)%basis_size
2125 3 : lmn_size_c=lmncmax(itypat)
2126 14 : llmax=maxval(indlmncor(itypat)%value(1,1:lmn_size_c))
2127 3 : l_size_max=max(l_size,2*llmax+1)
2128 :
2129 12 : ABI_MALLOC(j_bessel,(mesh_size,l_size_max))
2130 9 : ABI_MALLOC(ylmp,(l_size_max*l_size_max))
2131 15 : ABI_MALLOC(have_intc,(l_size_max,basis_size,nphicor(itypat)))
2132 15 : ABI_MALLOC(intc,(l_size_max,basis_size,nphicor(itypat)))
2133 12 : ABI_MALLOC(have_rad,(l_size,pawtab(itypat)%ij_size))
2134 12 : ABI_MALLOC(radint1,(l_size,pawtab(itypat)%ij_size))
2135 9 : ABI_MALLOC(radint2,(l_size,pawtab(itypat)%ij_size))
2136 9 : ABI_MALLOC(radint3,(l_size,pawtab(itypat)%ij_size))
2137 :
2138 24 : ABI_MALLOC(radsumc(itypat)%value,(2,lmn_size,lmn_size_c,n1,my_n2,n3,my_ngrid))
2139 21 : ABI_MALLOC(radsum1(itypat)%value,(2,lmn2_size,n1,my_n2,n3,my_ngrid))
2140 18 : ABI_MALLOC(radsum2(itypat)%value,(2,lmn2_size,n1,my_n2,n3,my_ngrid))
2141 18 : ABI_MALLOC(radsum3(itypat)%value,(2,lmn2_size,n1,my_n2,n3,my_ngrid))
2142 39447935 : radsumc(itypat)%value=zero
2143 65955455 : radsum1(itypat)%value=zero
2144 65955455 : radsum2(itypat)%value=zero
2145 65955455 : radsum3(itypat)%value=zero
2146 :
2147 9 : ABI_MALLOC(jbes,(l_size_max))
2148 9 : ABI_MALLOC(ff,(mesh_size))
2149 3 : meshsz=pawrad(itypat)%int_meshsz
2150 3 : if (meshsz>mesh_size) ff(meshsz+1:mesh_size)=zero
2151 :
2152 3 : indx=0;jkpt=0
2153 23 : do ikpt=1,nkpt
2154 20 : if (my_gridtab(ikpt)==0) cycle
2155 20 : jkpt=jkpt+1
2156 575 : do i3=1,n3
2157 552 : ig3=i3-(i3/id3)*n3-1
2158 15836 : do i2=1,n2
2159 15264 : if (me_fft/=fftn2_distrib(i2)) cycle
2160 15264 : j2=ffti2_local(i2)
2161 15264 : indx=n1*(my_n2*(i3-1)+(j2-1))
2162 15264 : ig2=i2-(i2/id2)*n2-1
2163 438744 : do i1=1,n1
2164 422928 : ig1=i1-(i1/id1)*n1-1
2165 422928 : indx=indx+1;if (mod(indx-1,nproc_band)/=me_band) cycle
2166 :
2167 : pcart(:)=Crystal%gprimd(:,1)*real(ig1+dtset%kpt(1,ikpt))+&
2168 : & Crystal%gprimd(:,2)*real(ig2+dtset%kpt(2,ikpt))+&
2169 1691712 : & Crystal%gprimd(:,3)*real(ig3+dtset%kpt(3,ikpt))
2170 1691712 : pnorm=dsqrt(dot_product(pcart,pcart))
2171 :
2172 422928 : if (pnorm < tol12) then
2173 2 : pbn(:) = zero
2174 20 : ylmp(:) = zero
2175 2 : ylmp(1) = 1.d0/sqrt(four_pi)
2176 : else
2177 1691704 : pbn(:) = pcart(:)/pnorm ! unit vector
2178 422926 : call initylmr(l_size,ylmr_normchoice,ylmr_npts,ylmr_nrm,ylmr_option,pbn(:),ylmp(:),ylmgr)
2179 : end if
2180 :
2181 422928 : pnorm=two_pi*pnorm ! re-normed for call to bessel
2182 498838392 : do ir = 1, mesh_size
2183 498415464 : bessarg = pnorm*pawrad(itypat)%rad(ir)
2184 498415464 : call sbf8(l_size,bessarg,jbes)
2185 1994084784 : j_bessel(ir,:)=jbes(:)
2186 : end do
2187 :
2188 : ! ===== Core part =====
2189 : ! Need intc=\int phi phi_core jl (pr) dr
2190 :
2191 35373456 : have_intc(:,:,:)=.FALSE. ; intc(:,:,:)=zero
2192 :
2193 1907712 : do jlmn = 1,lmn_size_c
2194 1484784 : jln = indlmncor(itypat)%value(5,jlmn)
2195 1484784 : jlm = indlmncor(itypat)%value(4,jlmn)
2196 1484784 : jl = indlmncor(itypat)%value(1,jlmn)
2197 14415840 : do ilmn = 1,lmn_size
2198 12508128 : iln = pawtab(itypat)%indlmn(5,ilmn)
2199 12508128 : ilm = pawtab(itypat)%indlmn(4,ilmn)
2200 12508128 : il = pawtab(itypat)%indlmn(1,ilmn)
2201 :
2202 12508128 : llmin = abs(il-jl)
2203 12508128 : llmax = il+jl
2204 12508128 : klm = RCC(ilm,jlm)
2205 31279392 : do ll=llmin,llmax,2
2206 17286480 : ifac=ilfac(mod(ll,4))
2207 :
2208 21947760 : if (.not.have_intc(ll+1,iln,jln)) then
2209 : ff(1:mesh_size)=(pawtab(itypat)%phi(1:mesh_size,iln)*phicor(itypat)%value(1:mesh_size,jln))&
2210 5268009744 : & *j_bessel(1:mesh_size,ll+1)
2211 4661280 : call simp_gen(intg,ff,pawrad(itypat))
2212 4661280 : intc(ll+1,iln,jln)=intg
2213 4661280 : have_intc(ll+1,iln,jln)=.true.
2214 : end if
2215 :
2216 78585552 : do mm=-ll,ll
2217 48790944 : lm = LMC(ll,mm)
2218 48790944 : glmij=pawang%gntselect(lm,klm)
2219 66077424 : if (glmij>0) then
2220 15162768 : arg=ylmp(lm)*pawang%realgnt(glmij)*intc(ll+1,iln,jln)
2221 : radsumc(itypat)%value(1,ilmn,jlmn,i1,j2,i3,jkpt) = &
2222 32449248 : & radsumc(itypat)%value(1,ilmn,jlmn,i1,j2,i3,jkpt)+arg*real(ifac)
2223 : radsumc(itypat)%value(2,ilmn,jlmn,i1,j2,i3,jkpt) = &
2224 32449248 : & radsumc(itypat)%value(2,ilmn,jlmn,i1,j2,i3,jkpt)+arg*aimag(ifac)
2225 : end if
2226 : end do !mm
2227 : end do !ll
2228 : end do !ilmn
2229 : end do !jlmn
2230 :
2231 : ! ===== Valence part =====
2232 : ! Need int1=\int phi_i phi_j jl (pr) dr
2233 : ! and int2=\int tphi_i tphi_j jl (pr) dr
2234 :
2235 95805072 : have_rad(:,:)= .FALSE.;radint1=zero;radint2=zero;radint3=zero
2236 :
2237 22277088 : do klmn=1,pawtab(itypat)%lmn2_size
2238 21838896 : klm=pawtab(itypat)%indklmn(1,klmn);kln=pawtab(itypat)%indklmn(2,klmn)
2239 21838896 : llmin=pawtab(itypat)%indklmn(3,klmn);llmax=pawtab(itypat)%indklmn(4,klmn)
2240 :
2241 56761344 : do ll=llmin,llmax,2
2242 34499520 : ifac=ilfac(mod(ll,4))
2243 :
2244 57607200 : if (.not.have_rad(ll+1,kln)) then
2245 9487422648 : ff(1:mesh_size)=pawtab(itypat)%phiphj(1:mesh_size,kln)*j_bessel(1:mesh_size,ll+1)
2246 7702560 : call simp_gen(intg,ff,pawrad(itypat))
2247 7702560 : radint1(ll+1,kln)=intg
2248 9487422648 : ff(1:mesh_size)=pawtab(itypat)%tphitphj(1:mesh_size,kln)*j_bessel(1:mesh_size,ll+1)
2249 7702560 : call simp_gen(intg,ff,pawrad(itypat))
2250 7702560 : radint2(ll+1,kln)=intg
2251 : ff(1:mesh_size)=(pawtab(itypat)%phiphj (1:mesh_size,kln) &
2252 : & -pawtab(itypat)%tphitphj(1:mesh_size,kln))&
2253 9487422648 : & *j_bessel(1:mesh_size,ll+1)
2254 7702560 : call simp_gen(intg,ff,pawrad(itypat))
2255 7702560 : radint3(ll+1,kln)=intg
2256 7702560 : have_rad(ll+1,kln)=.true.
2257 : end if
2258 :
2259 156354624 : do mm=-ll,ll
2260 100016208 : lm = LMC(ll,mm)
2261 100016208 : glmij=pawang%gntselect(lm,klm)
2262 134515728 : if (glmij>0) then
2263 30544776 : arg=ylmp(lm)*pawang%realgnt(glmij)
2264 : radsum1(itypat)%value(1,klmn,i1,j2,i3,jkpt) = &
2265 30544776 : & radsum1(itypat)%value(1,klmn,i1,j2,i3,jkpt)+real(ifac) *arg*radint1(ll+1,kln)
2266 : radsum1(itypat)%value(2,klmn,i1,j2,i3,jkpt) = &
2267 30544776 : & radsum1(itypat)%value(2,klmn,i1,j2,i3,jkpt)+aimag(ifac)*arg*radint1(ll+1,kln)
2268 : radsum2(itypat)%value(1,klmn,i1,j2,i3,jkpt) = &
2269 30544776 : & radsum2(itypat)%value(1,klmn,i1,j2,i3,jkpt)+real(ifac) *arg*radint2(ll+1,kln)
2270 : radsum2(itypat)%value(2,klmn,i1,j2,i3,jkpt) = &
2271 30544776 : & radsum2(itypat)%value(2,klmn,i1,j2,i3,jkpt)+aimag(ifac)*arg*radint2(ll+1,kln)
2272 : radsum3(itypat)%value(1,klmn,i1,j2,i3,jkpt) = &
2273 65044296 : & radsum3(itypat)%value(1,klmn,i1,j2,i3,jkpt)+real(ifac) *arg*radint3(ll+1,kln)
2274 : radsum3(itypat)%value(2,klmn,i1,j2,i3,jkpt) = &
2275 65044296 : & radsum3(itypat)%value(2,klmn,i1,j2,i3,jkpt)+aimag(ifac)*arg*radint3(ll+1,kln)
2276 : end if
2277 : end do !mm
2278 : end do !ll
2279 : end do !klmn
2280 :
2281 : end do ! end loop over i1
2282 : end do ! end loop over i2
2283 : end do ! end loop over i3
2284 : end do ! end loop over ikpt
2285 :
2286 3 : ABI_FREE(ff)
2287 3 : ABI_FREE(jbes)
2288 :
2289 3 : ABI_FREE(j_bessel)
2290 3 : ABI_FREE(ylmp)
2291 :
2292 3 : ABI_FREE(intc)
2293 3 : ABI_FREE(have_intc)
2294 :
2295 3 : ABI_FREE(radint1)
2296 3 : ABI_FREE(radint2)
2297 3 : ABI_FREE(radint3)
2298 3 : ABI_FREE(have_rad)
2299 :
2300 3 : call xmpi_sum(radsumc(itypat)%value,mpi_enreg%comm_band,ierr)
2301 3 : call xmpi_sum(radsum1(itypat)%value,mpi_enreg%comm_band,ierr)
2302 3 : call xmpi_sum(radsum2(itypat)%value,mpi_enreg%comm_band,ierr)
2303 6 : call xmpi_sum(radsum3(itypat)%value,mpi_enreg%comm_band,ierr)
2304 :
2305 : end do ! end loop over atom types
2306 : end if ! PAW
2307 :
2308 : !Allocate main memory
2309 9 : ABI_MALLOC(rho_contrib,(cplex*nfft))
2310 12 : ABI_MALLOC(rho_contrib_g,(cplex,nfft))
2311 9 : ABI_MALLOC(rho_contrib_paw1,(cplex,nfft))
2312 9 : ABI_MALLOC(rho_contrib_paw2,(cplex,nfft))
2313 9 : ABI_MALLOC(rho_contrib_paw3,(cplex,nfft))
2314 :
2315 12 : ABI_MALLOC(rho_moment_v1,(nfft,my_ngrid))
2316 9 : ABI_MALLOC(rho_moment_v2,(nfft,my_ngrid))
2317 9 : ABI_MALLOC(rho_moment_core,(nfft,my_ngrid))
2318 9 : ABI_MALLOC(rho_pw,(nfft,my_ngrid))
2319 845899 : rho_moment_v1=zero;rho_moment_v2=zero
2320 845899 : rho_pw=zero;rho_moment_core=zero
2321 :
2322 : !Prepare gamma(r) for the state independent scheme
2323 9 : ABI_MALLOC(gamma,(nfft,2))
2324 : if (.not.state_dependent) then
2325 : ABI_MALLOC(rhoe,(nfft,1))
2326 : ABI_MALLOC(rhop,(nfft,1))
2327 : if (electronpositron%particle==EP_ELECTRON) then
2328 : rhoe(:,1)=rhor_ep_(:,1);rhop(:,1)=rhor_(:,1)
2329 : else if (electronpositron%particle==EP_POSITRON) then
2330 : rhoe(:,1)=rhor_(:,1);rhop(:,1)=rhor_ep_(:,1)
2331 : end if
2332 : call mkdenpos(iwarn,nfft,1,1,rhoe(:,1),dtset%xc_denpos)
2333 : call mkdenpos(iwarn,nfft,1,1,rhop(:,1),dtset%xc_denpos)
2334 : call gammapositron_fft(electronpositron,gamma,Crystal%gprimd,igamma,mpi_enreg,&
2335 : & n3xccc,nfft,ngfft,rhoe(:,1),rhop(:,1),xccc3d)
2336 : ABI_FREE(rhoe)
2337 : ABI_FREE(rhop)
2338 : else
2339 132741 : gamma=one
2340 : end if
2341 :
2342 : !Some allocations for state-dependent scheme
2343 : if (state_dependent) then
2344 : ! Fake MPI data to be used in poslifetime; allow only FFT parallelism
2345 3 : call initmpi_seq(mpi_enreg_seq)
2346 3 : mpi_enreg_seq%my_natom=dtset%natom
2347 : call set_mpi_enreg_fft(mpi_enreg_seq,mpi_enreg%comm_fft,mpi_enreg%distribfft,&
2348 3 : & mpi_enreg%me_g0,mpi_enreg%paral_kgb)
2349 : ! Allocate memory for state-dependent scheme
2350 9 : ABI_MALLOC(rhor_dop_el,(nfft))
2351 3 : if (dtset%usepaw==1) then
2352 17 : ABI_MALLOC(pawrhoij_dop_el,(dtset%natom))
2353 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
2354 3 : & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cpxocc=dtset%pawcpxocc)
2355 : call pawrhoij_alloc(pawrhoij_dop_el,cplex_rhoij,nspden_rhoij,&
2356 : dtset%nspinor,dtset%nsppol,dtset%typat,&
2357 3 : pawtab=pawtab,use_rhoij_=1,use_rhoijp=1)
2358 : ! Cancel distribution of PAW data over atomic sites
2359 : ! We use here pawrhoij because polifetime routine
2360 : ! detects by itself the particle described by pawrhoij
2361 6 : if (mpi_enreg%my_natom<dtset%natom) then
2362 0 : ABI_MALLOC(pawrhoij_all,(dtset%natom))
2363 0 : call pawrhoij_nullify(pawrhoij_all)
2364 : call pawrhoij_gather(pawrhoij,pawrhoij_all,-1,mpi_enreg%comm_atom, &
2365 0 : & with_rhoijres=.false.,with_rhoij_=.false.,with_lmnmix=.false.)
2366 0 : ABI_MALLOC(pawrhoij_ep_all,(dtset%natom))
2367 0 : call pawrhoij_nullify(pawrhoij_ep_all)
2368 : call pawrhoij_gather(electronpositron%pawrhoij_ep,pawrhoij_ep_all,-1,mpi_enreg%comm_atom, &
2369 0 : & with_rhoijres=.false.,with_rhoij_=.false.,with_lmnmix=.false.)
2370 : else
2371 3 : pawrhoij_all => pawrhoij
2372 3 : pawrhoij_ep_all => electronpositron%pawrhoij_ep
2373 : end if
2374 : end if
2375 : end if
2376 :
2377 : !==============================================================================
2378 : !================ Loop over positronic states =================================
2379 :
2380 : !LOOP OVER k POINTS
2381 3 : ibg_pos=0;icg_pos=0;ikg_pos=0;bdtot_index_pos=0;isppol_pos=1
2382 6 : do ikpt_pos=1,merge(1,nkpt,kgamma_only_positron)
2383 :
2384 : ! Extract data for this kpt_pos
2385 3 : npw_k_pos=npwarr(ikpt_pos)
2386 3 : wtk_k_pos=dtset%wtk(ikpt_pos); if (kgamma_only_positron) wtk_k_pos=one
2387 3 : istwf_k_pos=dtset%istwfk(ikpt_pos)
2388 3 : nband_k_pos=dtset%nband(ikpt_pos+(isppol_pos-1)*nkpt)
2389 3 : nband_cprj_k_pos=nband_k_pos/nproc_band
2390 : mykpt_pos=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_pos,1,nband_k_pos,&
2391 3 : & isppol_pos,mpi_enreg%me_kpt))
2392 :
2393 : ! Retrieve additional data for this kpt_pos
2394 9 : ABI_MALLOC(occ_k_pos,(nband_k_pos))
2395 40 : occ_k_pos(:)=occ_pos_ptr(1+bdtot_index_pos:nband_k_pos+bdtot_index_pos)
2396 : nband_eff_pos=1
2397 40 : do ib_pos=1,nband_k_pos
2398 40 : if (occ_k_pos(ib_pos)>tol8) nband_eff_pos=ib_pos
2399 : end do
2400 3 : if (mod(nband_eff_pos,blocksize)/=0) nband_eff_pos=((nband_eff_pos/blocksize)+1)*blocksize
2401 :
2402 3 : nblock_band_eff_pos=nband_eff_pos/blocksize
2403 :
2404 3 : mcg_pos=npw_k_pos*my_nspinor*nband_eff_pos
2405 9 : ABI_MALLOC(cg_k_pos,(2,mcg_pos))
2406 :
2407 3 : mcprj_k_pos=0
2408 3 : if (dtset%usepaw==1) then
2409 3 : nband_cprj_eff_pos=nband_eff_pos/nproc_band
2410 3 : mcprj_k_pos=my_nspinor*nband_cprj_eff_pos
2411 23 : ABI_MALLOC(cprj_k_pos,(dtset%natom,mcprj_k_pos))
2412 3 : call pawcprj_alloc(cprj_k_pos,0,dimcprj)
2413 : end if
2414 :
2415 3 : if (mpi_enreg%paral_kgb==0) then
2416 12 : ABI_MALLOC(gbound_pos,(2*dtset%mgfft+8,2))
2417 9 : ABI_MALLOC(kg_k_pos,(3,npw_k_pos))
2418 0 : else if (mykpt_pos) then
2419 : nullify(bandfft_kpt_pos)
2420 : else
2421 0 : ABI_MALLOC(bandfft_kpt_pos,)
2422 0 : call bandfft_kpt_reset(bandfft_kpt_pos)
2423 : end if
2424 :
2425 : ! Exchange data (WF components) between procs
2426 3 : if (mykpt_pos) then
2427 3882 : cg_k_pos(:,1:mcg_pos)=cg_pos_ptr(:,icg_pos+1:icg_pos+mcg_pos)
2428 5175 : if (mpi_enreg%paral_kgb==0) kg_k_pos(:,1:npw_k_pos)=kg(:,1+ikg_pos:npw_k_pos+ikg_pos)
2429 3 : if (dtset%usepaw==1) then
2430 : call pawcprj_get(Crystal%atindx1,cprj_k_pos,cprj_pos_ptr,dtset%natom,1,ibg_pos,ikpt_pos,iorder_cprj,&
2431 : & isppol_pos,mband_cprj_pos,dtset%mkmem,dtset%natom,nband_cprj_eff_pos,nband_k_pos,my_nspinor,&
2432 3 : & dtset%nsppol,dtfil%unpaw,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
2433 : end if
2434 3 : if (mpi_enreg%paral_kgb/=0) then
2435 0 : jj=mpi_enreg%my_kpttab(ikpt_pos)
2436 0 : bandfft_kpt_pos => bandfft_kpt(jj)
2437 : end if
2438 6 : do ii=0,mpi_enreg%nproc_spkpt-1
2439 6 : if (ii/=mpi_enreg%me_kpt) then
2440 0 : tag=ikpt_pos+(isppol_pos-1)*nkpt+2*nkpt*ii
2441 0 : call xmpi_send(cg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr)
2442 0 : tag=tag+nkpt*(1+2*mpi_enreg%nproc_spkpt)
2443 0 : if (mpi_enreg%paral_kgb==0) then
2444 0 : call xmpi_send(kg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr)
2445 : else
2446 0 : call bandfft_kpt_mpi_send(bandfft_kpt_pos,ii,tag,mpi_enreg%comm_kpt,ierr,profile='fourwf')
2447 : end if
2448 0 : if (dtset%usepaw==1) then
2449 0 : call pawcprj_mpi_send(dtset%natom,mcprj_k_pos,dimcprj,0,cprj_k_pos,ii,mpi_enreg%comm_kpt,ierr)
2450 : end if
2451 : end if
2452 : end do
2453 : else
2454 0 : ii=0;if (allocated(mpi_enreg%proc_distrb)) ii=mpi_enreg%proc_distrb(ikpt_pos,1,isppol_pos)
2455 0 : tag=ikpt_pos+(isppol_pos-1)*nkpt+2*nkpt*mpi_enreg%me_kpt
2456 0 : call xmpi_recv(cg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr)
2457 0 : tag=tag+nkpt*(1+2*mpi_enreg%nproc_spkpt)
2458 0 : if (mpi_enreg%paral_kgb==0) then
2459 0 : call xmpi_recv(kg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr)
2460 : else
2461 0 : call bandfft_kpt_mpi_recv(bandfft_kpt_pos,ii,tag,mpi_enreg%comm_kpt,ierr)
2462 : end if
2463 0 : if (dtset%usepaw==1) then
2464 0 : call pawcprj_mpi_recv(dtset%natom,mcprj_k_pos,dimcprj,0,cprj_k_pos,ii,mpi_enreg%comm_kpt,ierr)
2465 : end if
2466 : end if
2467 :
2468 3 : if (mpi_enreg%paral_kgb==0) then
2469 3 : call sphereboundary(gbound_pos,istwf_k_pos,kg_k_pos,dtset%mgfft,npw_k_pos)
2470 : end if
2471 :
2472 9 : ABI_MALLOC(cwaver_pos,(cplex*nfft))
2473 9 : ABI_MALLOC(cwaver_pos_block,(cplex*nfft*bandpp))
2474 3 : if (dtset%usepaw==1) then
2475 23 : ABI_MALLOC(cprj_pos,(dtset%natom,my_nspinor))
2476 3 : call pawcprj_alloc(cprj_pos,0,dimcprj)
2477 : end if
2478 :
2479 : ! ============================================================================
2480 : ! Loops on positronic bands
2481 :
2482 6 : do iblock_pos=1,nblock_band_eff_pos
2483 3 : ib_pos=1+(iblock_pos-1)*blocksize
2484 6 : if (any(abs(occ_k_pos(ib_pos:ib_pos+blocksize-1))>tol8)) then
2485 :
2486 9 : ABI_MALLOC(cwaveg_pos,(2,npw_k_pos*blocksize))
2487 15 : ABI_MALLOC(cwaveaug_pos,(2,n4,n5,n6*bandpp))
2488 15 : ABI_MALLOC(denpot_dum,(n4,n5,n6))
2489 6 : ABI_MALLOC(fofgout_dum,(2,npw_k_pos*blocksize))
2490 3 : iwavef_pos=(iblock_pos-1)*npw_k_pos*blocksize
2491 : cwaveg_pos(:,1:npw_k_pos*blocksize)= &
2492 3882 : & cg_k_pos(:,iwavef_pos+1:iwavef_pos+npw_k_pos*blocksize)
2493 :
2494 : ! Get positronic wave function in real space
2495 3 : option=0
2496 3 : if (mpi_enreg%paral_kgb==0) then
2497 3 : weight_pos=occ_k_pos(ib_pos)*wtk_k_pos
2498 : call fourwf(1,denpot_dum,cwaveg_pos,fofgout_dum,cwaveaug_pos,&
2499 : & gbound_pos,gbound_pos,istwf_k_pos,kg_k_pos,kg_k_pos,&
2500 : & dtset%mgfft,mpi_enreg,1,ngfft,npw_k_pos,npw_k_pos,&
2501 : & n4,n5,n6,option,tim_fourwf,weight_pos,weight_pos,&
2502 3 : & gpu_option=dtset%gpu_option)
2503 : else
2504 : call prep_fourwf(denpot_dum,blocksize,cwaveg_pos,cwaveaug_pos,&
2505 : & iblock_pos,istwf_k_pos,dtset%mgfft,mpi_enreg,nband_k_pos,&
2506 : & bandpp,ngfft,npw_k_pos,n4,n5,n6,occ_k_pos,option,Crystal%ucvol,wtk_k_pos,1,&
2507 0 : & bandfft_kpt_tab=bandfft_kpt_pos,gpu_option=dtset%gpu_option)
2508 : end if
2509 :
2510 132735 : cwaver_pos_block=zero
2511 6 : do ii=1,bandpp
2512 3 : j3=(ii-1)*n3
2513 3 : indx0=1+(ii-1)*cplex*nfft
2514 90 : do i3=1,n3
2515 87 : if (me_fft==fftn3_distrib(i3)) then
2516 84 : indx=indx0+cplex*n1*n2*(ffti3_local(i3)-1)
2517 2442 : do i2=1,n2
2518 68808 : do i1=1,n1
2519 66366 : cwaver_pos_block(indx )=cwaveaug_pos(1,i1,i2,i3+j3)
2520 66366 : cwaver_pos_block(indx+1)=cwaveaug_pos(2,i1,i2,i3+j3)
2521 68724 : indx=indx+2
2522 : end do
2523 : end do
2524 : end if
2525 : end do
2526 : end do
2527 3 : ABI_FREE(fofgout_dum)
2528 3 : ABI_FREE(denpot_dum)
2529 3 : ABI_FREE(cwaveaug_pos)
2530 3 : ABI_FREE(cwaveg_pos)
2531 :
2532 : ! At this stage, each band proc has bandpp bands in real space
2533 : ! (distributed on FFT procs)
2534 :
2535 : ! ========================================================================
2536 : ! Compute core contribution for this positronic band (PAW only)
2537 :
2538 3 : if (dtset%usepaw==1) then
2539 6 : do ibpp_pos=1,bandpp
2540 3 : ib_cprj_pos=(iblock_pos-1)*bandpp+ibpp_pos
2541 3 : weight_pos=occ_k_pos(ib_pos+ibpp_pos-1+me_band*bandpp)*wtk_k_pos
2542 : ! Calculate the annihilation rate for each core state for state dependent scheme
2543 3 : iatm=0
2544 6 : do itypat=1,dtset%ntypat
2545 3 : mesh_size = pawtab(itypat)%mesh_size
2546 14 : do iat=1,Crystal%nattyp(itypat)
2547 8 : iatm=iatm+1;iatom=Crystal%atindx1(iatm)
2548 24 : ABI_MALLOC(gammastate_c(iatom)%value,(lmncmax(itypat)))
2549 43 : do jlmn=1,lmncmax(itypat)
2550 32 : jln = indlmncor(itypat)%value(5,jlmn)
2551 32 : contrib(:)=zero
2552 96 : ABI_MALLOC(rhocorej,(mesh_size))
2553 28734 : rhocorej(1:mesh_size)=2*phicor(itypat)%value(1:mesh_size,jln)**2
2554 : call posratecore(dtset,electronpositron,iatom,dtset%natom,mesh_size,mpi_enreg_seq,&
2555 32 : & 1,pawang,pawrad,pawrhoij_all,pawrhoij_ep_all,pawtab,ratec,rhocorej)
2556 :
2557 : call posratecore(dtset,electronpositron,iatom,dtset%natom,mesh_size,mpi_enreg_seq,&
2558 32 : & 2,pawang,pawrad,pawrhoij_all,pawrhoij_ep_all,pawtab,ratec_ipm,rhocorej)
2559 :
2560 32 : gammastate_c(iatom)%value(jlmn)=ratec/ratec_ipm
2561 40 : ABI_FREE(rhocorej)
2562 : end do
2563 : end do
2564 : end do
2565 3 : jkpt=0
2566 26 : do ikpt=1,nkpt
2567 20 : if (my_gridtab(ikpt)==0) cycle
2568 20 : jkpt=jkpt+1
2569 575 : do i3=1,n3
2570 552 : ig3=i3-(i3/id3)*n3-1
2571 15836 : do i2=1,n2
2572 15816 : if (me_fft==fftn2_distrib(i2)) then
2573 15264 : j2=ffti2_local(i2)
2574 15264 : ig2=i2-(i2/id2)*n2-1
2575 15264 : indx=n1*(my_n2*(i3-1)+(j2-1))
2576 438192 : do i1=1,n1
2577 422928 : ig1=i1-(i1/id1)*n1-1
2578 422928 : indx=indx+1
2579 :
2580 : ! Loop on atoms (type sorted)
2581 422928 : iatm=0
2582 861120 : do itypat=1,dtset%ntypat
2583 422928 : lmn_size = pawtab(itypat)%lmn_size
2584 :
2585 1907712 : do iat=1,Crystal%nattyp(itypat)
2586 1061856 : iatm=iatm+1;iatom=Crystal%atindx1(iatm)
2587 :
2588 : pcart(:)=Crystal%gprimd(:,1)*real(ig1+dtset%kpt(1,ikpt))+&
2589 : & Crystal%gprimd(:,2)*real(ig2+dtset%kpt(2,ikpt))+&
2590 4247424 : & Crystal%gprimd(:,3)*real(ig3+dtset%kpt(3,ikpt))
2591 1061856 : pnorm=dsqrt(dot_product(pcart,pcart))
2592 4247424 : pr=dot_product(pcart,Crystal%xcart(:,iatom))
2593 1061856 : expipr(1)= cos(two_pi*pr)
2594 1061856 : expipr(2)=-sin(two_pi*pr)
2595 :
2596 : ! Loop on ij states
2597 5534352 : do jlmn = 1,lmncmax(itypat)
2598 : contrib(:)=zero
2599 37705824 : do ilmn = 1,lmn_size
2600 : radsumnfftc(1)=expipr(1)*radsumc(itypat)%value(1,ilmn,jlmn,i1,j2,i3,jkpt)&
2601 33656256 : & -expipr(2)*radsumc(itypat)%value(2,ilmn,jlmn,i1,j2,i3,jkpt)
2602 : radsumnfftc(2)=expipr(1)*radsumc(itypat)%value(2,ilmn,jlmn,i1,j2,i3,jkpt)&
2603 33656256 : & +expipr(2)*radsumc(itypat)%value(1,ilmn,jlmn,i1,j2,i3,jkpt)
2604 100968768 : cp_pos(:)=cprj_k_pos(iatom,ib_cprj_pos)%cp(:,ilmn)
2605 : contrib(1)=contrib(1)+four_pi*(cp_pos(1)*radsumnfftc(1) &
2606 33656256 : & -cp_pos(2)*radsumnfftc(2))
2607 : contrib(2)=contrib(2)+four_pi*(cp_pos(1)*radsumnfftc(2) &
2608 37705824 : & +cp_pos(2)*radsumnfftc(1))
2609 : end do ! end loop over ilmn
2610 : ! 2 - electron state weight for 2 spins
2611 : rho_moment_core(indx,jkpt) = rho_moment_core(indx,jkpt) &
2612 5111424 : & +gammastate_c(iatom)%value(jlmn)*2*weight_pos*(contrib(1)**2+contrib(2)**2)
2613 : end do ! end loop over jlmn
2614 :
2615 : end do !end loop over atoms
2616 : end do !end loop over atom types
2617 :
2618 : end do ! end loop over i1
2619 : end if ! end loop over i2
2620 : end do
2621 : end do ! end loop over i3
2622 : end do ! jkpt
2623 : end do ! ibpp_pos
2624 : end if
2625 :
2626 : ! We now loop over positronic bands inside a block
2627 : ! and select occupied ones
2628 6 : do ibpp_pos=1,blocksize
2629 3 : ib_pos=(iblock_pos-1)*blocksize+ibpp_pos
2630 3 : occ_pos=occ_k_pos(ib_pos)
2631 6 : if (abs(occ_pos)>tol8) then
2632 :
2633 : ! Parallelism: dirty trick (broadcast bands) but there should be few positronic bands (~1)
2634 3 : if (nproc_band>1) then
2635 0 : iproc=(ibpp_pos-1)/bandpp
2636 0 : if (me_band==iproc) then
2637 0 : indx=mod((ibpp_pos-1),bandpp)*cplex*nfft
2638 0 : cwaver_pos(1:cplex*nfft)=cwaver_pos_block(indx+1:indx+cplex*nfft)
2639 : end if
2640 0 : call xmpi_bcast(cwaver_pos,iproc,mpi_enreg%comm_band,ierr)
2641 0 : if (dtset%usepaw==1) then
2642 0 : if (me_band==iproc) then
2643 0 : indx=mod((ibpp_pos-1),bandpp)*my_nspinor
2644 0 : call pawcprj_copy(cprj_k_pos(:,indx+1:indx+my_nspinor),cprj_pos)
2645 : end if
2646 : call pawcprj_bcast(cprj_pos,dtset%natom,my_nspinor,dimcprj,0,iproc,&
2647 0 : & mpi_enreg%comm_band,ierr)
2648 : end if
2649 : else
2650 132735 : cwaver_pos(1:cplex*nfft)=cwaver_pos_block(1:cplex*nfft)
2651 3 : if (dtset%usepaw==1) then
2652 3 : call pawcprj_copy(cprj_k_pos(:,(ib_pos-1)*my_nspinor+1:ib_pos*my_nspinor),cprj_pos)
2653 : end if
2654 : end if
2655 :
2656 : ! ========================================================================
2657 : ! ================ Loop over electronic states ===========================
2658 :
2659 : ! Loop over spins
2660 3 : ibg=0;icg=0;ikg=0;bdtot_index=0
2661 6 : do isppol=1,dtset%nsppol
2662 : ! Loop over k points
2663 3 : ikg=0;jkpt=0
2664 26 : do ikpt=1,nkpt
2665 :
2666 : ! Extract data for this kpt_pos
2667 20 : npw_k=npwarr(ikpt)
2668 20 : wtk_k=dtset%wtk(ikpt)
2669 20 : istwf_k=dtset%istwfk(ikpt)
2670 20 : nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
2671 20 : nband_cprj_k=nband_k/nproc_band
2672 : mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,&
2673 20 : & isppol,mpi_enreg%me_kpt))
2674 :
2675 : ! Select k-points for current proc
2676 20 : if (mykpt) then
2677 :
2678 : ! Retrieve additional data for this kpt_pos
2679 20 : jkpt=jkpt+1
2680 60 : ABI_MALLOC(occ_k,(nband_k))
2681 256 : occ_k(:)=occ_ptr(1+bdtot_index:nband_k+bdtot_index)
2682 :
2683 20 : mcprj_k=0
2684 20 : if (dtset%usepaw==1) then
2685 20 : mcprj_k=my_nspinor*nband_cprj_k
2686 908 : ABI_MALLOC(cprj_k,(dtset%natom,mcprj_k))
2687 20 : call pawcprj_alloc(cprj_k,0,dimcprj)
2688 : call pawcprj_get(Crystal%atindx1,cprj_k,cprj_ptr,dtset%natom,1,ibg,ikpt,iorder_cprj,&
2689 : & isppol,mband_cprj,dtset%mkmem,dtset%natom,nband_cprj_k,nband_cprj_k,my_nspinor,&
2690 20 : & dtset%nsppol,dtfil%unpaw,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
2691 : end if
2692 :
2693 20 : if (mpi_enreg%paral_kgb==0) then
2694 80 : ABI_MALLOC(gbound,(2*dtset%mgfft+8,2))
2695 60 : ABI_MALLOC(kg_k,(3,npw_k))
2696 33988 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
2697 20 : call sphereboundary(gbound,istwf_k,kg_k,dtset%mgfft,npw_k)
2698 : else
2699 0 : jj=mpi_enreg%my_kpttab(ikpt)
2700 0 : bandfft_kpt_el => bandfft_kpt(jj)
2701 : end if
2702 :
2703 60 : ABI_MALLOC(cwaver,(cplex*nfft*bandpp))
2704 :
2705 : ! ==================================================================
2706 : ! Loops on electronic bands
2707 :
2708 256 : do iblock=1,nblock_band
2709 236 : ib=1+(iblock-1)*blocksize
2710 :
2711 336 : if (any(abs(occ_k(ib:ib+blocksize-1))>tol8)) then
2712 :
2713 : ! Retrieve electronic wave function
2714 468 : ABI_MALLOC(cwaveg,(2,npw_k*blocksize))
2715 780 : ABI_MALLOC(cwaveaug,(2,n4,n5,n6*bandpp))
2716 780 : ABI_MALLOC(denpot_dum,(n4,n5,n6))
2717 312 : ABI_MALLOC(fofgout_dum,(2,npw_k*blocksize))
2718 156 : iwavef=(iblock-1)*npw_k*blocksize
2719 : cwaveg(:,1:npw_k*blocksize)= &
2720 222273 : & cg_ptr(:,icg+iwavef+1:icg+iwavef+npw_k*blocksize)
2721 :
2722 : ! Get electronic wave function in real space
2723 156 : option=0
2724 156 : if (mpi_enreg%paral_kgb==0) then
2725 156 : weight=occ_k(ib)*wtk_k
2726 : call fourwf(1,denpot_dum,cwaveg,fofgout_dum,cwaveaug,&
2727 : & gbound,gbound,istwf_k,kg_k,kg_k,&
2728 : & dtset%mgfft,mpi_enreg,1,ngfft,npw_k,npw_k,&
2729 : & n4,n5,n6,option,tim_fourwf,weight,weight,&
2730 156 : & gpu_option=dtset%gpu_option)
2731 : else
2732 : call prep_fourwf(denpot_dum,blocksize,cwaveg,cwaveaug,&
2733 : & iblock,istwf_k,dtset%mgfft,mpi_enreg,nband_k,&
2734 : & bandpp,ngfft,npw_k,n4,n5,n6,occ_k,option,Crystal%ucvol,wtk_k,1,&
2735 0 : & bandfft_kpt_tab=bandfft_kpt_el,gpu_option=dtset%gpu_option)
2736 : end if
2737 :
2738 6551004 : cwaver=zero
2739 312 : do ii=1,bandpp
2740 156 : j3=(ii-1)*n3
2741 156 : indx0=1+(ii-1)*cplex*nfft
2742 4608 : do i3=1,n3
2743 4452 : if (me_fft==fftn3_distrib(i3)) then
2744 4296 : indx=indx0+cplex*n1*n2*(ffti3_local(i3)-1)
2745 122808 : do i2=1,n2
2746 3398232 : do i1=1,n1
2747 3275424 : cwaver(indx )=cwaveaug(1,i1,i2,i3+j3)
2748 3275424 : cwaver(indx+1)=cwaveaug(2,i1,i2,i3+j3)
2749 3393936 : indx=indx+2
2750 : end do
2751 : end do
2752 : end if
2753 : end do
2754 : end do
2755 156 : ABI_FREE(fofgout_dum)
2756 156 : ABI_FREE(denpot_dum)
2757 156 : ABI_FREE(cwaveaug)
2758 156 : ABI_FREE(cwaveg)
2759 : ! At this stage, each band proc has bandpp bands in real space
2760 : ! (distributed on FFT procs)
2761 :
2762 : ! We now loop on the bandpp bands
2763 : ! and select occupied ones
2764 312 : do ibpp=1,bandpp
2765 156 : occ_el=occ_k(ib+ibpp-1+me_band*bandpp)
2766 312 : if (abs(occ_el)>tol8) then
2767 :
2768 : ! ==============================================================
2769 : ! Compute state-dependent annihilation rate
2770 : ! Avoid parallelism over kpt/bands/atoms
2771 156 : gammastate=one;rate_paw=one
2772 : if (state_dependent) then
2773 156 : weight=occ_el*wtk_k
2774 156 : ib_cprj=(iblock-1)*bandpp+ibpp
2775 156 : indx=1+(ibpp-1)*cplex*nfft
2776 3275580 : do ii=1,nfft
2777 3275424 : rhor_dop_el(ii)=weight*(cwaver(indx)*cwaver(indx)+cwaver(indx+1)*cwaver(indx+1))
2778 3275580 : indx=indx+2
2779 : end do
2780 156 : if (dtset%usepaw==1) then
2781 524 : do iatom=1,dtset%natom
2782 22204 : pawrhoij_dop_el(iatom)%rhoij_=zero
2783 : end do
2784 156 : cplex_rhoij=2;if (istwf_k>1) cplex_rhoij=1
2785 156 : use_timerev=(dtset%kptopt>0.and.dtset%kptopt<3)
2786 156 : use_zeromag=(pawrhoij_dop_el(1)%nspden==4.and.dtset%nspden==1)
2787 : call pawaccrhoij(Crystal%atindx,cplex_rhoij,cprj_k(:,ib_cprj),&
2788 : & cprj_k(:,ib_cprj),0,isppol,dtset%natom,dtset%natom,dtset%nspinor,&
2789 156 : & occ_el,1,pawrhoij_dop_el,use_timerev,use_zeromag,wtk_k)
2790 : ! Is it correct to apply symetries here (on a single band)?
2791 : ! If not, call pawrhoij_symrhoij with nsym=1
2792 : call pawrhoij_symrhoij(pawrhoij_dop_el,pawrhoij_dop_el,1,Crystal%gprimd,&
2793 : & Crystal%indsym,0,dtset%natom,Crystal%nsym,dtset%ntypat,1,pawang,-10001,&
2794 156 : & pawtab,Crystal%rprimd,Crystal%symafm,Crystal%symrec,dtset%typat)
2795 : end if
2796 : ! Has to call poslifetime in sequential because we are in a parallel section
2797 : ! Only FFT parallelism is allowed
2798 : call poslifetime(dtset,electronpositron,Crystal%gprimd,dtset%natom,mpi_enreg_seq,n3xccc,&
2799 : & nfft,ngfft,nhat,2,pawang,pawrad,pawrhoij_all,pawtab,rate,rate_paw,rhor,Crystal%ucvol,xccc3d,&
2800 156 : & rhor_dop_el=rhor_dop_el,pawrhoij_dop_el=pawrhoij_dop_el,pawrhoij_ep=pawrhoij_ep_all)
2801 : call poslifetime(dtset,electronpositron,Crystal%gprimd,dtset%natom,mpi_enreg_seq,n3xccc,&
2802 : & nfft,ngfft,nhat,3,pawang,pawrad,pawrhoij_all,pawtab,rate_ipm,rate_paw_ipm,rhor,Crystal%ucvol,xccc3d,&
2803 156 : & rhor_dop_el=rhor_dop_el,pawrhoij_dop_el=pawrhoij_dop_el,pawrhoij_ep=pawrhoij_ep_all)
2804 156 : gammastate=rate/rate_ipm
2805 156 : rate_paw=rate_paw/rate_paw_ipm
2806 : end if
2807 :
2808 : ! ==============================================================
2809 : ! Compute plane-wave contribution to momentum distribution
2810 :
2811 : ! Compute Psi^+(r) * Psi^-(r) * gamma(r) in real space
2812 6551004 : rho_contrib(:)=zero
2813 156 : indx=(ibpp-1)*cplex*nfft
2814 156 : if (cplex==2) then
2815 3275580 : do jj=1,nfft
2816 3275424 : ii=2*jj-1
2817 : rho_contrib(ii) =sqrt(gamma(jj,2))*(cwaver_pos(ii)*cwaver(indx+ii)&
2818 3275424 : & -wf_fact*cwaver_pos(ii+1)*cwaver(indx+ii+1))
2819 : rho_contrib(ii+1)=sqrt(gamma(jj,2))*(cwaver_pos(ii)*cwaver(indx+ii+1) &
2820 3275580 : & +wf_fact*cwaver_pos(ii+1)*cwaver(indx+ii))
2821 : end do
2822 : else
2823 0 : do ii=1,nfft
2824 0 : rho_contrib(ii)=sqrt(gamma(ii,2))*cwaver_pos(ii)*cwaver(indx+ii)
2825 : end do
2826 : end if
2827 :
2828 : ! FFT of (Psi+.Psi-.gamma) to get Intg[(Psi+.Psi-.gamma).exp(-igr)]
2829 : call fourdp(cplex,rho_contrib_g,rho_contrib,-1,mpi_enreg,nfft,1,ngfft,&
2830 156 : & tim_fourdp)
2831 :
2832 : rho_pw(1:nfft,jkpt)=rho_pw(1:nfft,jkpt) +gammastate*occ_el*occ_pos &
2833 3275580 : & *(rho_contrib_g(1,1:nfft)**2+rho_contrib_g(2,1:nfft)**2)
2834 :
2835 : ! ==============================================================
2836 : ! Compute PAW on-site contribution to momentum distribution
2837 :
2838 156 : if (dtset%usepaw==1) then
2839 :
2840 9826428 : rho_contrib_paw1(:,:)= zero
2841 9826428 : rho_contrib_paw2(:,:)= zero
2842 9826428 : rho_contrib_paw3(:,:)= zero
2843 :
2844 156 : ib_cprj=(iblock-1)*bandpp+ibpp
2845 :
2846 : ! Loop on moments
2847 156 : indx=0
2848 4452 : do i3=1,n3
2849 4296 : ig3=i3-(i3/id3)*n3-1
2850 122964 : do i2=1,n2
2851 122808 : if (me_fft==fftn2_distrib(i2)) then
2852 118512 : j2=ffti2_local(i2)
2853 118512 : ig2=i2-(i2/id2)*n2-1
2854 118512 : indx=n1*(my_n2*(i3-1)+(j2-1))
2855 3393936 : do i1=1,n1
2856 3275424 : ig1=i1-(i1/id1)*n1-1
2857 3275424 : indx=indx+1
2858 :
2859 : pcart(:)=Crystal%gprimd(:,1)*real(ig1+dtset%kpt(1,ikpt))+&
2860 : & Crystal%gprimd(:,2)*real(ig2+dtset%kpt(2,ikpt))+&
2861 13101696 : & Crystal%gprimd(:,3)*real(ig3+dtset%kpt(3,ikpt))
2862 3275424 : pnorm=dsqrt(dot_product(pcart,pcart))
2863 :
2864 : ! Loop on atoms (type-sorted)
2865 3275424 : iatm=0
2866 6550848 : do itypat=1,dtset%ntypat
2867 3275424 : lmn_size=pawtab(itypat)%lmn_size
2868 3275424 : lmn2_size=pawtab(itypat)%lmn2_size
2869 9826272 : ABI_MALLOC(radsumnfft1,(2,lmn2_size))
2870 6550848 : ABI_MALLOC(radsumnfft2,(2,lmn2_size))
2871 6550848 : ABI_MALLOC(radsumnfft3,(2,lmn2_size))
2872 :
2873 11338272 : do iat=1,Crystal%nattyp(itypat)
2874 8062848 : iatm=iatm+1;iatom=Crystal%atindx1(iatm)
2875 :
2876 32251392 : pr=dot_product(pcart,Crystal%xcart(:,iatom))
2877 8062848 : expipr(1)= cos(two_pi*pr)
2878 8062848 : expipr(2)=-sin(two_pi*pr)
2879 :
2880 457049088 : do klmn=1,lmn2_size
2881 : radsumnfft1(1,klmn)=expipr(1)*radsum1(itypat)%value(1,klmn,i1,j2,i3,jkpt)&
2882 448986240 : & -expipr(2)*radsum1(itypat)%value(2,klmn,i1,j2,i3,jkpt)
2883 : radsumnfft1(2,klmn)=expipr(1)*radsum1(itypat)%value(2,klmn,i1,j2,i3,jkpt)&
2884 448986240 : & +expipr(2)*radsum1(itypat)%value(1,klmn,i1,j2,i3,jkpt)
2885 : radsumnfft2(1,klmn)=expipr(1)*radsum2(itypat)%value(1,klmn,i1,j2,i3,jkpt)&
2886 448986240 : & -expipr(2)*radsum2(itypat)%value(2,klmn,i1,j2,i3,jkpt)
2887 : radsumnfft2(2,klmn)=expipr(1)*radsum2(itypat)%value(2,klmn,i1,j2,i3,jkpt)&
2888 448986240 : & +expipr(2)*radsum2(itypat)%value(1,klmn,i1,j2,i3,jkpt)
2889 : radsumnfft3(1,klmn)=expipr(1)*radsum3(itypat)%value(1,klmn,i1,j2,i3,jkpt)&
2890 448986240 : & -expipr(2)*radsum3(itypat)%value(2,klmn,i1,j2,i3,jkpt)
2891 : radsumnfft3(2,klmn)=expipr(1)*radsum3(itypat)%value(2,klmn,i1,j2,i3,jkpt)&
2892 457049088 : & +expipr(2)*radsum3(itypat)%value(1,klmn,i1,j2,i3,jkpt)
2893 : end do
2894 :
2895 : ! Loop on ij states
2896 90957600 : do ilmn = 1, lmn_size
2897 79619328 : i0lmn = ilmn*(ilmn-1)/2
2898 906035328 : do jlmn = 1, lmn_size
2899 818353152 : klmn = i0lmn+jlmn
2900 818353152 : if (jlmn>ilmn) then
2901 369366912 : i0lmn=jlmn*(jlmn-1)/2; klmn=i0lmn+ilmn
2902 : end if
2903 : ! Transform 3-dimentional radsum to 1-dimentional radsumnfft
2904 2455059456 : cp(:)=cprj_k(iatom,ib_cprj)%cp(:,ilmn)
2905 2455059456 : cp_pos(:)=cprj_pos(iatom,1)%cp(:,jlmn)
2906 818353152 : cp11= cp(1)*cp_pos(1)
2907 818353152 : cp22= cp(2)*cp_pos(2)*wf_fact
2908 818353152 : cp12= cp(1)*cp_pos(2)*wf_fact
2909 818353152 : cp21= cp(2)*cp_pos(1)
2910 818353152 : cpr=cp11-cp22 ; cpi=cp12+cp21
2911 : rho_contrib_paw1(1,indx) = rho_contrib_paw1(1,indx) &
2912 818353152 : & + four_pi*(cpr*radsumnfft1(1,klmn)-cpi*radsumnfft1(2,klmn))
2913 : rho_contrib_paw1(2,indx) = rho_contrib_paw1(2,indx) &
2914 818353152 : & + four_pi*(cpr*radsumnfft1(2,klmn)+cpi*radsumnfft1(1,klmn))
2915 : rho_contrib_paw2(1,indx) = rho_contrib_paw2(1,indx) &
2916 818353152 : & + four_pi*(cpr*radsumnfft2(1,klmn)-cpi*radsumnfft2(2,klmn))
2917 : rho_contrib_paw2(2,indx) = rho_contrib_paw2(2,indx) &
2918 818353152 : & + four_pi*(cpr*radsumnfft2(2,klmn)+cpi*radsumnfft2(1,klmn))
2919 : rho_contrib_paw3(1,indx) = rho_contrib_paw3(1,indx) &
2920 818353152 : & + four_pi*(cpr*radsumnfft3(1,klmn)-cpi*radsumnfft3(2,klmn))
2921 : rho_contrib_paw3(2,indx) = rho_contrib_paw3(2,indx) &
2922 897972480 : & + four_pi*(cpr*radsumnfft3(2,klmn)+cpi*radsumnfft3(1,klmn))
2923 : end do ! end loop over jlmn
2924 : end do ! end loop over ilmn
2925 :
2926 : end do !end loop over atoms
2927 :
2928 3275424 : ABI_FREE(radsumnfft1)
2929 3275424 : ABI_FREE(radsumnfft2)
2930 6550848 : ABI_FREE(radsumnfft3)
2931 : end do !end loop over atom types
2932 :
2933 : rho_moment_v1(indx,jkpt) = rho_moment_v1(indx,jkpt) &
2934 : & +occ_el*occ_pos &
2935 : & *(gammastate*(rho_contrib_g(1,indx)**2+rho_contrib_g(2,indx)**2) &
2936 : & +rate_paw*(rho_contrib_paw1(1,indx)**2+rho_contrib_paw1(2,indx)**2 &
2937 3275424 : & -rho_contrib_paw2(1,indx)**2-rho_contrib_paw2(2,indx)**2))
2938 : rho_moment_v2(indx,jkpt) = rho_moment_v2(indx,jkpt) &
2939 : & +occ_el*occ_pos*gammastate &
2940 : & *((rho_contrib_g(1,indx)+rho_contrib_paw3(1,indx))**2+&
2941 3393936 : & (rho_contrib_g(2,indx)+rho_contrib_paw3(2,indx))**2)
2942 :
2943 : end do ! end loop over i1
2944 :
2945 : end if ! end loop over i2
2946 : end do
2947 : end do ! end loop over i3
2948 :
2949 : end if ! PAW
2950 :
2951 : ! ================================================================
2952 : ! End loops on electronic bands
2953 :
2954 : end if ! occ>1.e-8
2955 : end do ! ibpp
2956 : end if ! occ_block>1.e-8
2957 : end do ! iblock
2958 :
2959 : ! End loops over k points and spins (electrons)
2960 20 : icg = icg + npw_k*my_nspinor*nband_k
2961 20 : ibg = ibg + my_nspinor*nband_cprj_k
2962 20 : ikg = ikg + npw_k
2963 :
2964 20 : ABI_FREE(cwaver)
2965 20 : ABI_FREE(occ_k)
2966 20 : if (mpi_enreg%paral_kgb==0) then
2967 20 : ABI_FREE(kg_k)
2968 20 : ABI_FREE(gbound)
2969 : else
2970 : nullify(bandfft_kpt_el)
2971 : end if
2972 20 : if (dtset%usepaw==1) then
2973 20 : call pawcprj_free(cprj_k)
2974 612 : ABI_FREE(cprj_k)
2975 : end if
2976 :
2977 : end if ! mykpt
2978 23 : bdtot_index=bdtot_index+nband_k
2979 : end do ! ikpt
2980 : end do ! isppol
2981 :
2982 : ! ================================================================
2983 : ! End loops on positronic bands
2984 :
2985 : end if ! occ>1.e-8
2986 : end do ! ibpp_pos
2987 : end if ! occ(block)>1.e-8
2988 : end do ! iblock_pos
2989 :
2990 : ! End loop over k points (positron)
2991 3 : if (mykpt_pos) then
2992 3 : icg_pos = icg_pos + npw_k_pos*my_nspinor*nband_k_pos
2993 3 : ibg_pos = ibg_pos + my_nspinor*nband_cprj_k_pos
2994 3 : ikg_pos = ikg_pos + npw_k_pos
2995 : end if
2996 3 : bdtot_index_pos=bdtot_index_pos+nband_k_pos
2997 :
2998 3 : ABI_FREE(cwaver_pos)
2999 3 : ABI_FREE(cwaver_pos_block)
3000 3 : ABI_FREE(cg_k_pos)
3001 3 : ABI_FREE(occ_k_pos)
3002 3 : if (mpi_enreg%paral_kgb==0) then
3003 3 : ABI_FREE(kg_k_pos)
3004 3 : ABI_FREE(gbound_pos)
3005 0 : else if (mykpt_pos) then
3006 : nullify(bandfft_kpt_pos)
3007 : else
3008 0 : call bandfft_kpt_destroy(bandfft_kpt_pos)
3009 0 : ABI_FREE(bandfft_kpt_pos)
3010 : end if
3011 6 : if (dtset%usepaw==1) then
3012 3 : call pawcprj_free(cprj_pos)
3013 11 : ABI_FREE(cprj_pos)
3014 3 : call pawcprj_free(cprj_k_pos)
3015 11 : ABI_FREE(cprj_k_pos)
3016 : end if
3017 :
3018 : end do ! ikpt_pos
3019 :
3020 : !================================================================
3021 : !Final computations and printing
3022 :
3023 : !In case of parallelism, sum over the communicator(s)
3024 3 : if (nproc_band>1) then
3025 0 : ABI_MALLOC(mpibuf,(3*nfft,my_ngrid))
3026 0 : do jkpt=1,my_ngrid
3027 0 : mpibuf( 1: nfft,jkpt)=rho_moment_v1(1:nfft,jkpt)
3028 0 : mpibuf( nfft+1:2*nfft,jkpt)=rho_moment_v2(1:nfft,jkpt)
3029 0 : mpibuf(2*nfft+1:3*nfft,jkpt)=rho_pw (1:nfft,jkpt)
3030 : end do
3031 0 : call xmpi_sum(mpibuf,mpi_enreg%comm_band,ierr)
3032 0 : do jkpt=1,my_ngrid
3033 0 : rho_moment_v1(1:nfft,jkpt)=mpibuf( 1: nfft,jkpt)
3034 0 : rho_moment_v2(1:nfft,jkpt)=mpibuf( nfft+1:2*nfft,jkpt)
3035 0 : rho_pw(1:nfft,jkpt) =mpibuf(2*nfft+1:3*nfft,jkpt)
3036 : end do
3037 0 : ABI_FREE(mpibuf)
3038 : end if
3039 3 : if (dtset%usepaw==1) then
3040 3 : call xmpi_sum(rho_moment_core,mpi_enreg%comm_band,ierr)
3041 : end if
3042 :
3043 : !Add valence and core contributions
3044 3 : if (dtset%usepaw==1) then
3045 3 : if (dtset%nsppol==2.and.my_nsppol==1) rho_moment_core(:,:)=half*rho_moment_core(:,:)
3046 422951 : rho_moment_v1(:,:)=rho_moment_v1(:,:)+rho_moment_core(:,:)
3047 422951 : rho_moment_v2(:,:)=rho_moment_v2(:,:)+rho_moment_core(:,:)
3048 : end if
3049 :
3050 3 : units_=pi*(one/InvFineStruct)**3/Time_Sec/1.e12_dp/electronpositron%posocc
3051 3 : scale_=(two_pi**2)/(Crystal%ucvol**two_thirds)
3052 :
3053 : !Integrate rho_moment over p
3054 422951 : buf(1)=sum(rho_moment_v1(1:nfft,1:my_ngrid))
3055 422951 : buf(2)=sum(rho_moment_v2(1:nfft,1:my_ngrid))
3056 422951 : buf(3)=sum(rho_moment_core(1:nfft,1:my_ngrid))
3057 422951 : buf(4)=sum(rho_pw(1:nfft,1:my_ngrid))
3058 3 : call xmpi_sum(buf,mpi_enreg%comm_kpt,ierr)
3059 3 : call xmpi_sum(buf,mpi_enreg%comm_fft,ierr)
3060 3 : lambda_v1=buf(1)*units_/Crystal%ucvol/nkpt
3061 3 : lambda_v2=buf(2)*units_/Crystal%ucvol/nkpt
3062 3 : lambda_core=buf(3)*units_/Crystal%ucvol/nkpt
3063 3 : lambda_pw=buf(4)*units_/Crystal%ucvol/nkpt
3064 :
3065 : !Write result in _DOPPLER file
3066 : !Requires MPI-IO if nproc_fft>1
3067 3 : if (me_band==0) then
3068 3 : if (me_kpt==0) then
3069 3 : filename_dop=trim(dtfil%filnam_ds(4))//'_DOPPLER'
3070 : vec=sqrt(dot_product(Crystal%gprimd(:,3),Crystal%gprimd(:,3)))
3071 9 : ABI_MALLOC(pcart_k,(3,nfft))
3072 9 : ABI_MALLOC(rho_moment_k,(nfft))
3073 3 : if (dtset%nsppol==2) then
3074 0 : ABI_MALLOC(rho_moment_k2,(nfft))
3075 : end if
3076 3 : if (accessfil==IO_MODE_FORTRAN) then ! >>>>> Fortran access
3077 : ! Open file and write first line
3078 3 : ierr=open_file(filename_dop,msg,newunit=unit_doppler,form='unformatted')
3079 3 : write(unit_doppler) nfft,nkpt,Crystal%ucvol,Crystal%rprimd(:,:)
3080 : else ! >>>>> MPI-IO access
3081 0 : unit_doppler=get_unit()
3082 : ! Open file and write first line
3083 0 : call WffOpen(IO_MODE_MPI,mpi_enreg%comm_fft,filename_dop,ierr,wff,0,me_fft,unit_doppler)
3084 0 : if (me_fft==0) then
3085 0 : call xderiveWRecInit(wff,ierr)
3086 0 : call xderiveWrite(wff,n1*n2*n3,ierr)
3087 0 : call xderiveWrite(wff,nkpt,ierr)
3088 0 : call xderiveWrite(wff,Crystal%ucvol,ierr)
3089 0 : call xderiveWrite(wff,Crystal%rprimd(:,:),ierr)
3090 0 : call xderiveWRecEnd(wff,ierr)
3091 : else
3092 0 : call xmoveOff(wff,n_int=2,n_dp=10,n_mark=2)
3093 : end if
3094 : ! Store table of FFT points treated by current proc
3095 0 : ABI_MALLOC(my_ffttab,(nfft))
3096 0 : my_ffttab=0
3097 0 : do i3=1,n3
3098 0 : do i2=1,n2
3099 0 : if (me_fft==fftn2_distrib(i2)) then
3100 0 : indx0=n1*(n2*(i3-1)+(i2-1))
3101 0 : indx=n1*(my_n2*(i3-1)+(ffti2_local(i2)-1))
3102 0 : my_ffttab(indx+1:indx+n1)=(/(indx0+ii,ii=1,n1)/)
3103 : end if
3104 : end do
3105 : end do
3106 0 : ABI_MALLOC(mpibuf,(1,nfft))
3107 : end if
3108 : end if
3109 :
3110 3 : jkpt=0
3111 23 : do ikpt=1,nkpt
3112 20 : if (nproc_spkpt==1) then
3113 422948 : rho_moment_k(1:nfft)=rho_moment_v2(1:nfft,ikpt)
3114 : else
3115 0 : if (my_gridtab(ikpt)/=0) jkpt=jkpt+1
3116 0 : if (me_kpt==0) then
3117 0 : if (my_gridtab(ikpt)==0) then
3118 0 : tag=ikpt;iproc=mpi_enreg%proc_distrb(ikpt,1,1)
3119 0 : call xmpi_recv(rho_moment_k,iproc,tag,mpi_enreg%comm_kpt,ierr)
3120 0 : if (dtset%nsppol==2) then
3121 0 : tag=2*ikpt;iproc=mpi_enreg%proc_distrb(ikpt,1,2)
3122 0 : call xmpi_recv(rho_moment_k2,iproc,tag,mpi_enreg%comm_kpt,ierr)
3123 0 : rho_moment_k(1:nfft)=rho_moment_k(1:nfft)+rho_moment_k2(1:nfft)
3124 : end if
3125 0 : else if (any(mpi_enreg%my_isppoltab(:)==1)) then
3126 0 : rho_moment_k(1:nfft)=rho_moment_v2(1:nfft,jkpt)
3127 0 : if (dtset%nsppol==2) then
3128 0 : ii=2;if (mpi_enreg%my_isppoltab(2)==1) ii=1
3129 0 : tag=ii*ikpt;iproc=mpi_enreg%proc_distrb(ikpt,1,ii)
3130 0 : call xmpi_recv(rho_moment_k2,iproc,tag,mpi_enreg%comm_kpt,ierr)
3131 0 : rho_moment_k(1:nfft)=rho_moment_k(1:nfft)+rho_moment_k2(1:nfft)
3132 : end if
3133 : end if
3134 0 : else if (my_gridtab(ikpt)/=0) then
3135 0 : if (mpi_enreg%my_isppoltab(1)==1) then
3136 0 : tag=ikpt
3137 0 : call xmpi_send(rho_moment_v2(1:nfft,jkpt),0,tag,mpi_enreg%comm_kpt,ierr)
3138 : end if
3139 0 : if (dtset%nsppol==2.and.mpi_enreg%my_isppoltab(2)==1) then
3140 0 : tag=2*ikpt
3141 0 : call xmpi_send(rho_moment_v2(1:nfft,jkpt),0,tag,mpi_enreg%comm_kpt,ierr)
3142 : end if
3143 : end if
3144 : end if ! nproc_spkpt>1
3145 23 : if (me_kpt==0) then
3146 20 : indx=0
3147 572 : do i3=1,n3
3148 552 : ig3=i3-(i3/id3)*n3-1
3149 15836 : do i2=1,n2
3150 15264 : if (me_fft/=fftn2_distrib(i2)) cycle
3151 15264 : ig2=i2-(i2/id2)*n2-1
3152 438744 : do i1=1,n1
3153 422928 : ig1=i1-(i1/id1)*n1-1
3154 422928 : indx=indx+1
3155 : pcart_k(:,indx)=Crystal%gprimd(:,1)*real(ig1+dtset%kpt(1,ikpt))+&
3156 : & Crystal%gprimd(:,2)*real(ig2+dtset%kpt(2,ikpt))+&
3157 1706976 : & Crystal%gprimd(:,3)*real(ig3+dtset%kpt(3,ikpt))
3158 : end do
3159 : end do
3160 : end do
3161 20 : if (accessfil==IO_MODE_FORTRAN) then
3162 20 : write(unit_doppler) pcart_k(1:3,1:nfft),rho_moment_k(1:nfft)
3163 : else
3164 0 : mpibuf(1,1:nfft)=rho_moment_k(1:nfft)
3165 0 : call xderiveWRecInit(wff,ierr)
3166 0 : call xderiveWrite(wff,pcart_k,3,nfft,mpi_enreg%comm_fft,my_ffttab,ierr)
3167 0 : call xderiveWrite(wff,mpibuf ,1,nfft,mpi_enreg%comm_fft,my_ffttab,ierr)
3168 0 : call xderiveWRecEnd(wff,ierr)
3169 : end if
3170 : end if
3171 : end do
3172 3 : if (me_kpt==0) then
3173 3 : ABI_FREE(pcart_k)
3174 3 : ABI_FREE(rho_moment_k)
3175 3 : if (dtset%nsppol==2) then
3176 0 : ABI_FREE(rho_moment_k2)
3177 : end if
3178 3 : if (accessfil==IO_MODE_FORTRAN) then
3179 3 : ierr=close_unit(unit_doppler,msg)
3180 : else
3181 0 : call WffClose(wff,ierr)
3182 0 : ABI_FREE(my_ffttab)
3183 0 : ABI_FREE(mpibuf)
3184 : end if
3185 : end if
3186 : end if ! me_band==0
3187 :
3188 : !Write results
3189 : write(msg,'(7a)') &
3190 3 : & ' Computation of electron-positron pairs momentum distribution completed.',ch10,&
3191 3 : & '-File ',trim(filename_dop),' has been created.',ch10,&
3192 6 : & '-Use ~abinit/scripts/post_processing/posdopspectra.F90 to process it.'
3193 3 : call wrtout(ab_out,msg,'COLL')
3194 3 : call wrtout(std_out,msg,'COLL')
3195 3 : msg=' Some e-p annihilation rates (ns-1) obtained by integration of e-p pairs momentum distribution:'
3196 3 : call wrtout(std_out,msg,'COLL')
3197 : write(msg,'(a,es22.12,3(2a,es22.12))') &
3198 3 : & ' Lambda (from module of sum of PAW contrib.) = ',lambda_v2*1000._dp,ch10,&
3199 3 : & ' = lambda_core: ',lambda_core*1000._dp,ch10,&
3200 3 : & ' +lambda_pw : ',lambda_pw*1000._dp,ch10,&
3201 6 : & ' +lambda_paw : ',(lambda_v2-lambda_core-lambda_pw)*1000._dp
3202 3 : call wrtout(std_out,msg,'COLL')
3203 : write(msg,'(4(a,es22.12,a))') &
3204 3 : & ' Lambda (from sum of modules of PAW contrib.) = ',lambda_v1*1000._dp,ch10,&
3205 3 : & ' = lambda_core: ',lambda_core*1000._dp,ch10,&
3206 3 : & ' +lambda_pw : ',lambda_pw*1000._dp,ch10,&
3207 6 : & ' +lambda_paw : ',(lambda_v1-lambda_core-lambda_pw)*1000._dp,ch10
3208 3 : call wrtout(std_out,msg,'COLL')
3209 3 : write(msg,'(4a,es22.12,2a)') ch10,&
3210 3 : & ' Annihilation rate obtained from integration of e-p pairs momentum distribution:',ch10,&
3211 6 : & ' lambda=',lambda_v2*1000._dp,' ns-1',ch10
3212 3 : call wrtout(ab_out,msg,'COLL')
3213 :
3214 : !Deallocate remaining memory
3215 3 : ABI_FREE(my_gridtab)
3216 3 : ABI_FREE(rho_pw)
3217 3 : ABI_FREE(rho_moment_v1)
3218 3 : ABI_FREE(rho_moment_v2)
3219 3 : ABI_FREE(rho_moment_core)
3220 3 : ABI_FREE(rho_contrib)
3221 3 : ABI_FREE(rho_contrib_g)
3222 3 : ABI_FREE(rho_contrib_paw1)
3223 3 : ABI_FREE(rho_contrib_paw2)
3224 3 : ABI_FREE(rho_contrib_paw3)
3225 : if (state_dependent) then
3226 3 : call unset_mpi_enreg_fft(mpi_enreg_seq)
3227 3 : call destroy_mpi_enreg(mpi_enreg_seq)
3228 3 : ABI_FREE(rhor_dop_el)
3229 3 : if (dtset%usepaw==1) then
3230 3 : call pawrhoij_free(pawrhoij_dop_el)
3231 11 : ABI_FREE(pawrhoij_dop_el)
3232 3 : if (mpi_enreg%my_natom<dtset%natom) then
3233 0 : call pawrhoij_free(pawrhoij_all)
3234 0 : call pawrhoij_free(pawrhoij_ep_all)
3235 0 : ABI_FREE(pawrhoij_all)
3236 0 : ABI_FREE(pawrhoij_ep_all)
3237 : end if
3238 : end if
3239 : end if
3240 :
3241 3 : ABI_FREE(gamma)
3242 :
3243 3 : if (dtset%usepaw==1.and.(.not.include_nhat_in_gamma)) then
3244 3 : ABI_FREE(rhor_)
3245 3 : ABI_FREE(rhor_ep_)
3246 : end if
3247 :
3248 3 : if (dtset%usepaw==1) then
3249 3 : ABI_FREE(nphicor)
3250 3 : ABI_FREE(lmncmax)
3251 6 : do itypat=1,dtset%ntypat
3252 3 : if (allocated(phicor(itypat)%value)) then
3253 3 : ABI_FREE(phicor(itypat)%value)
3254 : end if
3255 3 : if (allocated(indlmncor(itypat)%value)) then
3256 3 : ABI_FREE(indlmncor(itypat)%value)
3257 : end if
3258 3 : if (allocated(radsumc(itypat)%value)) then
3259 3 : ABI_FREE(radsumc(itypat)%value)
3260 : end if
3261 3 : if (allocated(radsum1(itypat)%value)) then
3262 3 : ABI_FREE(radsum1(itypat)%value)
3263 : end if
3264 3 : if (allocated(radsum2(itypat)%value)) then
3265 3 : ABI_FREE(radsum2(itypat)%value)
3266 : end if
3267 6 : if (allocated(radsum3(itypat)%value)) then
3268 3 : ABI_FREE(radsum3(itypat)%value)
3269 : end if
3270 : end do
3271 11 : do iatom=1,dtset%natom
3272 11 : if (allocated(gammastate_c(iatom)%value)) then
3273 8 : ABI_FREE(gammastate_c(iatom)%value)
3274 : end if
3275 : end do
3276 6 : ABI_FREE(phicor)
3277 6 : ABI_FREE(indlmncor)
3278 6 : ABI_FREE(radsumc)
3279 6 : ABI_FREE(radsum1)
3280 6 : ABI_FREE(radsum2)
3281 6 : ABI_FREE(radsum3)
3282 11 : ABI_FREE(gammastate_c)
3283 : end if
3284 :
3285 : DBG_EXIT("COLL")
3286 :
3287 6 : end subroutine posdoppler
3288 : !!***
3289 :
3290 : !!****f* ABINIT/posratecore
3291 : !! NAME
3292 : !! posratecore
3293 : !!
3294 : !! FUNCTION
3295 : !! Calculate the annihilataion rate of a given core state
3296 : !!
3297 : !! INPUTS
3298 : !! dtset <type(dataset_type)>=all input variables for this dataset
3299 : !! | nspden=number of spin-density components
3300 : !! | ntypat=number of atom types
3301 : !! | paral_kgb=flag controlling (k,g,bands) parallelization
3302 : !! | pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
3303 : !! | usepaw=flag for PAW
3304 : !! iatom= index of the current atom in posdoppler
3305 : !! mesh_sizej= size of the radial mesh for the current atom in posdoppler
3306 : !! mpi_enreg= information about MPI parallelization
3307 : !! my_natom=number of atoms treated by current processor
3308 : !! option= if 1, use gamma
3309 : !! if 2, use IPM (gamma=1)
3310 : !! pawang <type(pawang)>=paw angular mesh and related data
3311 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
3312 : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies
3313 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
3314 : !!
3315 : !! OUTPUT
3316 : !! rate= annihilation rate of a given core state needed for state dependent scheme for doppler broadening
3317 : !!
3318 : !! SIDE EFFECTS
3319 : !!
3320 : !! SOURCE
3321 :
3322 64 : subroutine posratecore(dtset,electronpositron,iatom,my_natom,mesh_sizej,mpi_enreg,&
3323 64 : & option,pawang,pawrad,pawrhoij,pawrhoij_ep,&
3324 64 : & pawtab,rate,rhocorej)
3325 :
3326 : !Arguments ------------------------------------
3327 : !scalars
3328 : integer,intent(in) :: iatom,my_natom,option,mesh_sizej
3329 : real(dp),intent(out) :: rate
3330 : type(dataset_type), intent(in) :: dtset
3331 : type(electronpositron_type),pointer :: electronpositron
3332 : type(MPI_type),intent(in) :: mpi_enreg
3333 : type(pawang_type), intent(in) :: pawang
3334 : !arrays
3335 : real(dp),intent(in) :: rhocorej(mesh_sizej)
3336 : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw)
3337 : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*dtset%usepaw)
3338 : type(pawrhoij_type),intent(in),target :: pawrhoij_ep(my_natom*dtset%usepaw)
3339 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
3340 :
3341 : !Local variables-------------------------------
3342 : !scalars
3343 : integer :: cplex,ierr,igamma,ii,ilm,ipt,ir
3344 : integer :: itypat,iwarn,iwarnj,iwarnp,lm_size,lmn2_size,mesh_size
3345 : integer :: ngr,ngrad,nspden_ep,opt_dens
3346 : logical,parameter :: include_nhat_in_gamma=.false.
3347 : real(dp),parameter :: delta=1.d-4
3348 : real(dp) :: fact,fact2,intg
3349 : real(dp) :: mpibuf,rdum,sqfpi
3350 : character(len=500) :: msg
3351 : !arrays
3352 64 : logical,allocatable :: lmselect(:),lmselect_ep(:),lmselect_dum(:)
3353 : real(dp),parameter :: qphon(3)=(/zero,zero,zero/),lsign(2)=(/one,-one/)
3354 64 : real(dp),allocatable :: d1gam(:,:),d2gam(:,:),ff(:),gam_(:,:,:),gamma(:,:),gammam(:,:),gg(:,:)
3355 64 : real(dp),allocatable :: grhocore2(:),grhocor2_(:),grhoe2(:),grho2_(:)
3356 64 : real(dp),allocatable :: nhat1(:,:,:),nhat1_ep(:,:,:)
3357 64 : real(dp),allocatable :: rho_(:),rho_ep_(:),rho1(:,:,:),rho1_ep(:,:,:)
3358 64 : real(dp),allocatable :: rhoarr1(:),rhoarr1_ep(:),rhoarr2(:)
3359 64 : real(dp),allocatable :: rhocore(:),rhocor_(:)
3360 64 : real(dp),allocatable :: rhosph(:),rhosph_ep(:),rhotot(:,:),rhotot_ep(:,:)
3361 64 : real(dp),allocatable :: trho1(:,:,:),trho1_ep(:,:,:)
3362 64 : real(dp),allocatable :: v1sum(:,:),v2sum(:,:,:)
3363 64 : type(pawrhoij_type),pointer :: pawrhoij_ep_(:)
3364 :
3365 : ! *************************************************************************
3366 :
3367 : DBG_ENTER("COLL")
3368 :
3369 : !Tests for developers
3370 64 : if (.not.associated(electronpositron)) then
3371 0 : msg='electronpositron variable must be associated!'
3372 0 : ABI_BUG(msg)
3373 : end if
3374 : !Constants
3375 64 : fact=0.0
3376 64 : cplex=1;nspden_ep=1
3377 64 : ngrad=1;if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2
3378 64 : iwarn=0;iwarnj=0;iwarnp=1
3379 64 : sqfpi=sqrt(four_pi)
3380 :
3381 : !Compatibility tests
3382 64 : if (electronpositron%particle==EP_NOTHING) then
3383 0 : msg='Not valid for electronpositron%particle=NOTHING!'
3384 0 : ABI_BUG(msg)
3385 : end if
3386 :
3387 64 : if (dtset%usepaw==1) then
3388 64 : if(dtset%pawxcdev==0.and.ngrad==2) then
3389 0 : msg='GGA is not implemented for pawxcdev=0 (use dtset%pawxcdev/=0)!'
3390 0 : ABI_BUG(msg)
3391 : end if
3392 : end if
3393 :
3394 : !Select type(s) of enhancement factor
3395 64 : if (electronpositron%ixcpositron==-1) igamma=0
3396 64 : if (electronpositron%ixcpositron== 2) igamma=4
3397 64 : if (electronpositron%ixcpositron==11.or.electronpositron%ixcpositron==31) igamma=3
3398 64 : if (electronpositron%ixcpositron==1.or.electronpositron%ixcpositron==3) igamma=2
3399 64 : if (option==2) igamma=0
3400 :
3401 64 : pawrhoij_ep_ => pawrhoij_ep
3402 :
3403 64 : rate=zero
3404 :
3405 64 : itypat=pawrhoij(iatom)%itypat
3406 64 : lmn2_size=pawtab(itypat)%lmn2_size
3407 64 : mesh_size=pawtab(itypat)%mesh_size
3408 64 : lm_size=pawtab(itypat)%lcut_size**2
3409 : cplex=1
3410 64 : ngr=0;if (ngrad==2) ngr=mesh_size
3411 :
3412 : !Allocations of "on-site" densities
3413 256 : ABI_MALLOC(rho1 ,(cplex*mesh_size,lm_size,nspden_ep))
3414 192 : ABI_MALLOC(trho1,(cplex*mesh_size,lm_size,nspden_ep))
3415 192 : ABI_MALLOC(rho1_ep ,(cplex*mesh_size,lm_size,nspden_ep))
3416 192 : ABI_MALLOC(trho1_ep,(cplex*mesh_size,lm_size,nspden_ep))
3417 192 : ABI_MALLOC(lmselect,(lm_size))
3418 128 : ABI_MALLOC(lmselect_ep,(lm_size))
3419 128 : ABI_MALLOC(lmselect_dum,(lm_size))
3420 : if (include_nhat_in_gamma) then
3421 : ABI_MALLOC(nhat1,(cplex*mesh_size,lm_size,nspden_ep))
3422 : ABI_MALLOC(nhat1_ep,(cplex*mesh_size,lm_size,nspden_ep))
3423 : else
3424 64 : ABI_MALLOC(nhat1,(0,0,0))
3425 64 : ABI_MALLOC(nhat1_ep,(0,0,0))
3426 : end if
3427 :
3428 : !Compute "on-site" densities (n1, ntild1, nhat1) for electron and positron =====
3429 640 : lmselect(:)=.true.
3430 64 : opt_dens=1;if (include_nhat_in_gamma) opt_dens=0
3431 : call pawdensities(rdum,cplex,iatom,lmselect,lmselect_dum,lm_size,nhat1,nspden_ep,1,&
3432 : & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij(iatom),&
3433 64 : & pawtab(itypat),rho1,trho1)
3434 640 : lmselect_ep(:)=.true.
3435 : call pawdensities(rdum,cplex,iatom,lmselect_ep,lmselect_dum,lm_size,nhat1_ep,nspden_ep,1,&
3436 : & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij_ep_(iatom),&
3437 64 : & pawtab(itypat),rho1_ep,trho1_ep)
3438 : !Compute contribution to annihilation rate
3439 :
3440 192 : ABI_MALLOC(rhocore,(mesh_size))
3441 :
3442 : !First formalism: use densities on r,theta,phi
3443 64 : if (dtset%pawxcdev==0) then
3444 :
3445 0 : ABI_MALLOC(gamma,(mesh_size,2))
3446 0 : ABI_MALLOC(rhoarr1,(mesh_size))
3447 0 : ABI_MALLOC(rhoarr1_ep,(mesh_size))
3448 :
3449 : ! Loop on the angular part
3450 0 : do ipt=1,pawang%angl_size
3451 : ! Build densities
3452 0 : rhoarr1=zero;rhoarr1_ep=zero;rhocore=zero
3453 0 : do ilm=1,lm_size
3454 0 : if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+rho1(:,ilm,1)*pawang%ylmr(ilm,ipt)
3455 : end do
3456 0 : do ilm=1,lm_size
3457 0 : if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+rho1_ep(:,ilm,1)*pawang%ylmr(ilm,ipt)
3458 : end do
3459 0 : rhocore(:)=pawtab(itypat)%coredens(:)
3460 : ! Make the densities positive
3461 0 : if (electronpositron%particle==EP_ELECTRON) then
3462 0 : call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1 ,dtset%xc_denpos)
3463 0 : call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos)
3464 0 : else if (electronpositron%particle==EP_POSITRON) then
3465 0 : call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1 ,dtset%xc_denpos)
3466 0 : call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos)
3467 : end if
3468 : ! Compute Gamma
3469 0 : ABI_MALLOC(grhoe2,(ngr))
3470 0 : ABI_MALLOC(grhocore2,(ngr))
3471 0 : if (electronpositron%particle==EP_ELECTRON) then
3472 : call gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,mesh_size,&
3473 0 : & rhocore,rhoarr1_ep,rhoarr1,1)
3474 0 : else if (electronpositron%particle==EP_POSITRON) then
3475 : call gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,mesh_size,&
3476 0 : & rhocore,rhoarr1,rhoarr1_ep,1)
3477 : end if
3478 0 : ABI_FREE(grhoe2)
3479 0 : ABI_FREE(grhocore2)
3480 : ! Compute contribution to annihilation rates
3481 :
3482 0 : ABI_MALLOC(ff,(mesh_size))
3483 : ff(1:mesh_size)=rhoarr1_ep(1:mesh_size)*rhocorej(1:mesh_size) &
3484 0 : & *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2
3485 0 : call simp_gen(intg,ff,pawrad(itypat))
3486 0 : intg=intg*pawang%angwgth(ipt)*four_pi
3487 0 : rate =rate +intg
3488 0 : ABI_FREE(ff)
3489 : end do ! ipt
3490 0 : ABI_FREE(gamma)
3491 0 : ABI_FREE(rhoarr1)
3492 0 : ABI_FREE(rhoarr1_ep)
3493 :
3494 : !Second formalism: use (l,m) moments for densities
3495 : else if (dtset%pawxcdev/=0) then
3496 :
3497 : ! Build densities
3498 192 : ABI_MALLOC(gammam,(mesh_size,lm_size))
3499 192 : ABI_MALLOC(rhotot,(mesh_size,lm_size))
3500 192 : ABI_MALLOC(rhotot_ep,(mesh_size,lm_size))
3501 128 : ABI_MALLOC(rhosph,(mesh_size))
3502 128 : ABI_MALLOC(rhosph_ep,(mesh_size))
3503 :
3504 517276 : rhotot (:,:)=rho1 (:,:,1)
3505 517276 : rhotot_ep(:,:)=rho1_ep(:,:,1)
3506 57468 : rhocore(:)=pawtab(itypat)%coredens(:)
3507 57468 : rhosph (:)=rhotot (:,1)/sqfpi
3508 57468 : rhosph_ep(:)=rhotot_ep(:,1)/sqfpi
3509 : ! Make spherical densities positive
3510 64 : if (electronpositron%particle==EP_ELECTRON) then
3511 0 : call mkdenpos(iwarnp,mesh_size,1,1,rhosph ,dtset%xc_denpos)
3512 0 : call mkdenpos(iwarn ,mesh_size,1,1,rhosph_ep,dtset%xc_denpos)
3513 64 : else if (electronpositron%particle==EP_POSITRON) then
3514 64 : call mkdenpos(iwarn ,mesh_size,1,1,rhosph ,dtset%xc_denpos)
3515 64 : call mkdenpos(iwarnp,mesh_size,1,1,rhosph_ep,dtset%xc_denpos)
3516 : end if
3517 :
3518 : ! Need gradients of electronic densities for GGA
3519 192 : ABI_MALLOC(grhoe2,(ngr))
3520 128 : ABI_MALLOC(grhocore2,(ngr))
3521 64 : if (ngr>0) then
3522 0 : if (electronpositron%particle==EP_ELECTRON) then
3523 0 : call nderiv_gen(grhoe2,rhosph_ep,pawrad(itypat))
3524 0 : else if (electronpositron%particle==EP_POSITRON) then
3525 0 : call nderiv_gen(grhoe2,rhosph,pawrad(itypat))
3526 : end if
3527 0 : grhoe2(:)=grhoe2(:)**2
3528 0 : call nderiv_gen(grhocore2,rhocore,pawrad(itypat))
3529 0 : grhocore2(:)=grhocore2(:)**2
3530 : end if
3531 : ! Compute Gamma for (rho-,rho+),
3532 : ! (rho- +drho-,rho+), (rho- -drho-,rho+),
3533 : ! (rho-,rho+ +drho+), (rho-,rho+ -drho+),
3534 : ! (rho- +drho-,rho+ +drho+), (rho- -drho-,rho+ -drho+)
3535 : ! Do a seven steps loop
3536 256 : ABI_MALLOC(gam_,(mesh_size,2,7))
3537 128 : ABI_MALLOC(rho_,(mesh_size))
3538 128 : ABI_MALLOC(rho_ep_,(mesh_size))
3539 128 : ABI_MALLOC(rhocor_,(mesh_size))
3540 128 : ABI_MALLOC(grho2_,(ngr))
3541 128 : ABI_MALLOC(grhocor2_,(ngr))
3542 :
3543 512 : do ii=1,7
3544 : ! Apply delta to get perturbed densities
3545 1206828 : rho_(:)=rhosph(:);rho_ep_(:)=rhosph_ep(:);rhocor_(:)=rhocore(:)
3546 448 : if (ngr>0) grho2_(:)=grhoe2(:)
3547 448 : if (ngr>0) grhocor2_(:)=grhocore2(:)
3548 448 : if (ii==2.or.ii==4.or.ii==6) fact=(one+delta)
3549 448 : if (ii==3.or.ii==5.or.ii==7) fact=(one-delta)
3550 448 : fact2=fact**2
3551 448 : if (ii==2.or.ii==3.or.ii==6.or.ii==7) then
3552 229872 : rho_(:)=fact*rho_(:)
3553 256 : if (electronpositron%particle==EP_POSITRON) then
3554 256 : if (ngr>0) grho2_(:)=fact2*grho2_(:)
3555 229872 : rhocor_(:)=fact*rhocor_(:)
3556 256 : if (ngr>0) grhocor2_(:)=fact2*grhocor2_(:)
3557 : end if
3558 : end if
3559 :
3560 448 : if (ii==4.or.ii==5.or.ii==6.or.ii==7) then
3561 229872 : rho_ep_(:)=fact*rho_ep_(:)
3562 256 : if (electronpositron%particle==EP_ELECTRON) then
3563 0 : if (ngr>0) grho2_(:)=fact2*grho2_(:)
3564 0 : rhocor_(:)=fact*rhocor_(:)
3565 0 : if (ngr>0) grhocor2_(:)=fact2*grhocor2_(:)
3566 : end if
3567 : end if
3568 : ! Compute gamma for these perturbed densities
3569 512 : if (electronpositron%particle==EP_ELECTRON) then
3570 0 : call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma,ngr,mesh_size,rhocor_,rho_ep_,rho_,1)
3571 448 : else if (electronpositron%particle==EP_POSITRON) then
3572 448 : call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma,ngr,mesh_size,rhocor_,rho_,rho_ep_,1)
3573 : end if
3574 :
3575 : end do ! end loop ii=1,7
3576 :
3577 64 : ABI_FREE(rhocor_)
3578 64 : ABI_FREE(grho2_)
3579 64 : ABI_FREE(grhocor2_)
3580 64 : ABI_FREE(grhoe2)
3581 64 : ABI_FREE(grhocore2)
3582 114936 : rho_ (:)=rhosph (:);if (electronpositron%particle==EP_POSITRON) rho_ (:)=rho_ (:)+rhocore(:)
3583 57468 : rho_ep_(:)=rhosph_ep(:);if (electronpositron%particle==EP_ELECTRON) rho_ep_(:)=rho_ep_(:)+rhocore(:)
3584 : ! Compute numerical first and second derivatives of Gamma
3585 : ! d1gam(1) = dgam/drho+ (particle=ELECTRON), dgam/drho- (particle=POSITRON)
3586 : ! d1gam(2) = dgam/drho- (particle=ELECTRON), dgam/drho+ (particle=POSITRON)
3587 192 : ABI_MALLOC(d1gam,(mesh_size,2))
3588 115000 : d1gam(:,:)=zero
3589 57468 : do ir=1,mesh_size
3590 57404 : if (rho_ (ir)>tol14) d1gam(ir,1)=(gam_(ir,1,2)-gam_(ir,1,3))*half/(delta*rho_ (ir))
3591 57468 : if (rho_ep_ (ir)>tol14) d1gam(ir,2)=(gam_(ir,1,4)-gam_(ir,1,5))*half/(delta*rho_ep_ (ir))
3592 : end do
3593 :
3594 : ! d2gam(1) = d2gam/drho+_drho+ (particle=ELECTRON), dgam/drho-_drho- (particle=POSITRON)
3595 : ! d2gam(2) = d2gam/drho-_drho+ (particle=ELECTRON), dgam/drho+_drho- (particle=POSITRON)
3596 : ! d2gam(3) = d2gam/drho-_drho- (particle=ELECTRON), dgam/drho+_drho+ (particle=POSITRON)
3597 192 : ABI_MALLOC(d2gam,(mesh_size,3))
3598 172468 : d2gam(:,:)=zero
3599 57468 : do ir=1,mesh_size
3600 57404 : if (rho_ (ir)>tol14) d2gam(ir,1)=(gam_(ir,1,2)+gam_(ir,1,3)-two*gam_(ir,1,1))/(delta*rho_ (ir))**2
3601 57468 : if (rho_ep_(ir)>tol14) then
3602 57404 : d2gam(ir,3)=(gam_(ir,1,4)+gam_(ir,1,5)-two*gam_(ir,1,1))/(delta*rho_ep_(ir))**2
3603 57404 : if (rho_(ir)>tol14) then
3604 : d2gam(ir,2)=(gam_(ir,1,6)+gam_(ir,1,7)+two*gam_(ir,1,1) &
3605 : & -gam_(ir,1,2)-gam_(ir,1,3)-gam_(ir,1,4)-gam_(ir,1,5)) &
3606 57404 : & *half/(delta*rho_(ir))/(delta*rho_ep_(ir))
3607 : end if
3608 : end if
3609 : end do
3610 :
3611 64 : ABI_FREE(rho_)
3612 64 : ABI_FREE(rho_ep_)
3613 : ! Compute useful sums of densities
3614 128 : ABI_MALLOC(v1sum,(mesh_size,3))
3615 64 : if ( dtset%pawxcdev>=2) then
3616 0 : ABI_MALLOC(v2sum,(mesh_size,lm_size,3))
3617 : else
3618 64 : ABI_MALLOC(v2sum,(0,0,0))
3619 : end if
3620 114936 : rhotot(:,1)=sqfpi*rhosph(:);rhotot_ep(:,1)=sqfpi*rhosph_ep(:)
3621 : call pawxcsum(1,1,1,lmselect,lmselect_ep,lm_size,mesh_size,3,dtset%pawxcdev,&
3622 64 : & pawang,rhotot,rhotot_ep,v1sum,v2sum)
3623 : ! Compute final development of gamma moments
3624 517276 : gammam(:,:)=zero
3625 57468 : gammam(:,1)=gam_(:,1,1)*sqfpi
3626 : gammam(:,1)=gammam(:,1)+(d2gam(:,2)*v1sum(:,2) &
3627 57468 : & +half*(d2gam(:,1)*v1sum(:,1)+d2gam(:,3)*v1sum(:,3)))/sqfpi
3628 576 : do ilm=2,lm_size
3629 512 : if (lmselect(ilm)) then
3630 459744 : gammam(:,ilm)=gammam(:,ilm)+d1gam(:,1)*rhotot(:,ilm)
3631 : end if
3632 576 : if (lmselect_ep(ilm)) then
3633 459744 : gammam(:,ilm)=gammam(:,ilm)+d1gam(:,2)*rhotot_ep(:,ilm)
3634 : end if
3635 : end do
3636 64 : if (dtset%pawxcdev>1) then
3637 0 : do ilm=2,lm_size
3638 : gammam(:,ilm)=gammam(:,ilm)+d2gam(:,2)*v2sum(:,ilm,2) &
3639 0 : & +half*(d2gam(:,1)*v2sum(:,ilm,1)+d2gam(:,3)*v2sum(:,ilm,3))
3640 : end do
3641 : end if
3642 :
3643 64 : ABI_FREE(gam_)
3644 64 : ABI_FREE(d1gam)
3645 64 : ABI_FREE(d2gam)
3646 64 : ABI_FREE(v1sum)
3647 64 : ABI_FREE(v2sum)
3648 : ! Compute contribution to annihilation rate
3649 192 : ABI_MALLOC(gg,(mesh_size,4))
3650 229936 : gg=zero
3651 128 : ABI_MALLOC(rhoarr2,(mesh_size))
3652 640 : do ilm=1,lm_size
3653 517276 : if (lmselect_ep(ilm)) gg(:,1)=gg(:,1)+rhotot_ep(:,ilm)*rhocorej(:)*gammam(:,ilm)
3654 : end do
3655 64 : ABI_FREE(rhoarr2)
3656 57468 : gg(1:mesh_size,1)=gg(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2
3657 64 : call simp_gen(intg,gg(:,1),pawrad(itypat))
3658 64 : rate =rate +intg
3659 64 : ABI_FREE(gg)
3660 64 : ABI_FREE(gammam)
3661 64 : ABI_FREE(rhotot)
3662 64 : ABI_FREE(rhotot_ep)
3663 64 : ABI_FREE(rhosph)
3664 64 : ABI_FREE(rhosph_ep)
3665 :
3666 : end if ! dtset%pawxcdev
3667 64 : ABI_FREE(rhocore)
3668 :
3669 64 : ABI_FREE(rho1)
3670 64 : ABI_FREE(trho1)
3671 64 : ABI_FREE(rho1_ep)
3672 64 : ABI_FREE(trho1_ep)
3673 64 : ABI_FREE(lmselect)
3674 64 : ABI_FREE(lmselect_ep)
3675 64 : ABI_FREE(lmselect_dum)
3676 64 : ABI_FREE(nhat1)
3677 64 : ABI_FREE(nhat1_ep)
3678 :
3679 : !Reduction in case of distribution over atomic sites
3680 64 : if (mpi_enreg%nproc_atom>1) then
3681 0 : mpibuf=rate
3682 0 : call xmpi_sum(mpibuf,mpi_enreg%comm_atom,ierr)
3683 0 : rate=mpibuf
3684 : end if
3685 :
3686 : DBG_EXIT("COLL")
3687 :
3688 128 : end subroutine posratecore
3689 : !!***
3690 :
3691 : end module m_positron
3692 : !!***
|
