Branch data 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 : : !!***
|
