Branch data Line data Source code
1 : : !!****m* ABINIT/m_electronpositron
2 : : !! NAME
3 : : !! m_electronpositron
4 : : !!
5 : : !! FUNCTION
6 : : !! This module provides the definition of the electronpositron_type used
7 : : !! used to store data for the electron-positron two-component DFT
8 : : !! as methods to operate on it.
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2008-2026 ABINIT group (MT, GJ)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public License, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !!
16 : : !! SOURCE
17 : :
18 : : #if defined HAVE_CONFIG_H
19 : : #include "config.h"
20 : : #endif
21 : : #include "abi_common.h"
22 : :
23 : : MODULE m_electronpositron
24 : :
25 : : use defs_basis
26 : : use m_abicore
27 : : use m_errors
28 : : use m_energies
29 : : use m_xmpi
30 : : use m_cgtools
31 : : use m_dtset
32 : :
33 : : use defs_abitypes, only : MPI_type
34 : : use m_pawtab, only : pawtab_type
35 : : use m_paw_an, only : paw_an_type
36 : : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy
37 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy
38 : : use m_mpinfo, only : proc_distrb_cycle
39 : : use m_xcpositron, only : xcpositron
40 : : use m_drivexc, only : mkdenpos
41 : : use m_xctk, only : xcden
42 : : use m_fft, only : fourdp
43 : :
44 : : implicit none
45 : :
46 : : private
47 : :
48 : : ! public constants
49 : : integer,public,parameter :: EP_NOTHING =-1
50 : : integer,public,parameter :: EP_ELECTRON = 0
51 : : integer,public,parameter :: EP_POSITRON = 1
52 : : !!***
53 : :
54 : : !!****t* m_electronpositron/electronpositron_type
55 : : !! NAME
56 : : !!
57 : : !! FUNCTION
58 : : !!
59 : : !! NOTES
60 : : !!
61 : : !! SOURCE
62 : :
63 : : type, public :: electronpositron_type
64 : :
65 : : ! Integer scalars
66 : : integer :: calctype ! type of electron-positron calculation:
67 : : ! 0: no calculation
68 : : ! 1: positron in the electrons potential
69 : : ! 2: electrons in the positron potential
70 : : integer :: particle ! current particle stored in electronpositron%xxx_ep arrays
71 : : ! -1: no particle, 0: electron, 1: positron
72 : : integer :: dimcg ! Dimension of cg array dimcg=mcg (=dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol,
73 : : ! or 0 for a process without wavefunctions to treat)
74 : : integer :: dimcprj ! Dimension of cprj array dimcprj=dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol*usecprj
75 : : integer :: dimeigen ! Dimension of eigen array dimeigen=dtset%mband*dtset%nkpt*dtset%nsppol
76 : : integer :: dimocc ! Dimension of occ array dimocc=dtset%mband*dtset%nkpt*dtset%nsppol
77 : : integer :: has_pawrhoij_ep ! flag for pawrhoij_ep (0: not allocated, 1: allocated, 2: computed)
78 : : integer :: has_pos_ham ! flag: 1 if current Hamiltonian in memory (vtrial, vpsp, vhartr, vxc, paw_ij%dij)
79 : : ! is the positronic hamiltonian, 0 is it is the electronic one
80 : : integer :: ixcpositron ! XC type for electron-positron correlation
81 : : integer :: istep ! Current index of TC-DFT SCF step
82 : : integer :: istep_scf ! Current index of DFT SCF step in current electron/positron minimization
83 : : integer :: lmmax ! Max. number of (l,m) moments over all types of atom
84 : : integer :: natom ! Number of atoms
85 : : integer :: nfft ! Number of points in FFT grid
86 : : integer :: nspden ! Number of spin density components
87 : : integer :: nstep ! Max. number of steps for the TC-DFT SCF cycle
88 : :
89 : : ! Logical scalars
90 : : logical :: posdensity0_limit ! True if we are in the zero positron density limit
91 : : logical :: scf_converged ! True if the SCF cycle is converged for a positronic/electronic GS calculation
92 : :
93 : : ! Real(dp) scalars
94 : : real(dp) :: e_hartree ! Hartree electron-positron interaction energy
95 : : real(dp) :: e_xc ! XC electron-positron interaction energy
96 : : real(dp) :: e_xcdc ! Double-counting XC electron-positron interaction energy
97 : : real(dp) :: e_paw ! PAW electron-positron interaction energy
98 : : real(dp) :: e_pawdc ! Double-counting PAW electron-positron interaction energy
99 : : real(dp) :: e0 ! Energy only due to particle(s) currently evolving
100 : : ! calctype=1, energy due to positron only
101 : : ! calctype=2, energy due to electrons only
102 : : real(dp) :: etotal_prev ! Total energy of the previous GS calculation
103 : : real(dp) :: lambda ! Electron-positron annihilation rate
104 : : real(dp) :: lifetime ! Positron lifetime
105 : : real(dp) :: maxfor_prev ! Max. force of the previous GS calculation
106 : : real(dp) :: posocc ! Occupation number for the positron
107 : : real(dp) :: postoldfe ! Tolerance on total energy for the TC-DFT SCF cycle
108 : : real(dp) :: postoldff ! Tolerance on max. force for the TC-DFT SCF cycle
109 : :
110 : : ! Other scalars
111 : : type(energies_type) :: energies_ep ! Energies of the previous electronic/positronic SCF step
112 : :
113 : : ! Logical pointers
114 : : logical, allocatable :: lmselect_ep(:,:)
115 : : ! lmselect_ep(lmmax,my_natom)
116 : : ! flags selecting the non-zero LM-moments of on-site densities
117 : :
118 : : ! Real(dp) pointers
119 : : real(dp), allocatable :: cg_ep(:,:)
120 : : ! cg_ep(2,dimcg)
121 : : ! if typecalc=1: electronic wavefunctions
122 : : ! if typecalc=2: positronic wavefunctions
123 : :
124 : : real(dp), allocatable :: eigen_ep(:)
125 : : ! eigen(dimeigen)
126 : : ! if typecalc=1: electronic eigen energies
127 : : ! if typecalc=2: positronic eigen energies
128 : :
129 : : real(dp), allocatable :: gred_ep(:,:)
130 : : ! gred_ep(3,natom)
131 : : ! if typecalc=1: forces only due to electrons
132 : : ! if typecalc=2: forces only due to positron
133 : :
134 : : real(dp), allocatable :: nhat_ep(:,:)
135 : : ! nhat_ep(nfft,nspden)
136 : : ! if typecalc=1: electronic compensation charge density in real space
137 : : ! if typecalc=2: positronic compensation charge density in real space
138 : :
139 : : real(dp), allocatable :: occ_ep(:)
140 : : ! occ(dimocc)
141 : : ! if typecalc=1: electronic occupations
142 : : ! if typecalc=2: positronic occupations
143 : :
144 : : real(dp), allocatable :: rhor_ep(:,:)
145 : : ! rhor_ep(nfft,nspden)
146 : : ! if typecalc=1: electronic density in real space
147 : : ! if typecalc=2: positronic density in real space
148 : :
149 : : real(dp), allocatable :: stress_ep(:)
150 : : ! stress_ep(6)
151 : : ! if typecalc=1: stresses only due to electrons
152 : : ! if typecalc=2: stresses only due to positron
153 : :
154 : : real(dp), allocatable :: vha_ep(:)
155 : : ! vha_ep(nfft)
156 : : ! if typecalc=1: electronic Hartree potential
157 : : ! if typecalc=2: positronic Hartree potential
158 : :
159 : : ! Other pointers
160 : : type(pawrhoij_type), allocatable :: pawrhoij_ep(:)
161 : : ! pawrhoij_ep(natom)
162 : : ! Relevant only if PAW
163 : : ! if typecalc=1: electronic PAW occupation matrix associated with rhor_ep
164 : : ! if typecalc=2: positronic PAW occupation matrix associated with rhor_ep
165 : :
166 : : type(pawcprj_type), allocatable :: cprj_ep(:,:)
167 : : ! cprj_ep(natom,dimcprj)
168 : : ! Relevant only if PAW
169 : : ! if typecalc=1: electronic WF projected on nl projectors <p_i|Cnk>
170 : : ! if typecalc=2: positronic WF projected on nl projectors <p_i|Cnk>
171 : :
172 : : end type electronpositron_type
173 : :
174 : : ! public procedures
175 : : public :: init_electronpositron
176 : : public :: destroy_electronpositron
177 : : public :: exchange_electronpositron
178 : : public :: electronpositron_calctype
179 : : public :: rhohxcpositron
180 : :
181 : : CONTAINS
182 : :
183 : : !===========================================================
184 : : !!***
185 : :
186 : : !!****f* m_electronpositron/init_electronpositron
187 : : !! NAME
188 : : !! init_electronpositron
189 : : !!
190 : : !! FUNCTION
191 : : !! Init all scalars and pointers in the structure.
192 : : !!
193 : : !! INPUTS
194 : : !! ireadwf=if 1, read the wavefunction
195 : : !! dtset <type(dataset_type)>=all input variables for this dataset
196 : : !! mcg=size of the wave-functions array (cg) on this process
197 : : !! mpi_enreg=information about MPI parallelization
198 : : !! nfft=(effective) number of FFT grid points (for this processor)
199 : : !! pawrhoij(natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies
200 : : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
201 : : !!
202 : : !! SIDE EFFECTS
203 : : !! electronpositron=<type(electronpositron_type)>=electronpositron datastructure
204 : : !!
205 : : !! SOURCE
206 : :
207 : 14 : subroutine init_electronpositron(ireadwf,dtset,electronpositron,mcg,mpi_enreg,nfft,pawrhoij,pawtab)
208 : :
209 : : !Arguments ------------------------------------
210 : : !scalars
211 : : integer,intent(in) :: ireadwf,mcg,nfft
212 : : type(dataset_type),intent(in) :: dtset
213 : : type(electronpositron_type),pointer :: electronpositron
214 : : type(MPI_type),intent(in) :: mpi_enreg
215 : : !arrays
216 : : type(pawrhoij_type), intent(in) :: pawrhoij(mpi_enreg%my_natom*dtset%usepaw)
217 : : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
218 : :
219 : : !Local variables-------------------------------
220 : : !scalars
221 : : integer :: ii,my_nspinor,ncpgr,optfor,optstr
222 : : logical,parameter :: include_nhat_in_gamma=.false.
223 : : !arrays
224 : 14 : integer,allocatable :: nlmn(:)
225 : :
226 : : !************************************************************************
227 : :
228 : : !@electronpositron_type
229 : :
230 [ + - ]: 14 : if (dtset%positron/=0) then
231 : :
232 [ - + ]: 14 : ABI_MALLOC(electronpositron,)
233 : :
234 : 14 : electronpositron%calctype=0
235 : 14 : electronpositron%particle=-1
236 : :
237 : 14 : electronpositron%ixcpositron=dtset%ixcpositron
238 : 14 : electronpositron%natom=dtset%natom
239 : 14 : electronpositron%nfft=nfft
240 : 14 : electronpositron%nspden=dtset%nspden
241 : 14 : electronpositron%istep=0
242 : 14 : electronpositron%istep_scf=0
243 : :
244 : 14 : electronpositron%posocc=dtset%posocc
245 : 14 : electronpositron%nstep=dtset%posnstep
246 : 14 : electronpositron%postoldfe=dtset%postoldfe
247 : 14 : electronpositron%postoldff=dtset%postoldff
248 : 14 : electronpositron%posdensity0_limit=(dtset%ixcpositron/=2)
249 : 14 : electronpositron%scf_converged=.false.
250 : 14 : electronpositron%has_pos_ham=0
251 : :
252 : 14 : call electronpositron%energies_ep%init()
253 : :
254 : 14 : electronpositron%e_hartree =zero
255 : 14 : electronpositron%e_xc =zero
256 : 14 : electronpositron%e_xcdc =zero
257 : 14 : electronpositron%e_paw =zero
258 : 14 : electronpositron%e_pawdc =zero
259 : 14 : electronpositron%e0 =zero
260 : 14 : electronpositron%etotal_prev=zero
261 : 14 : electronpositron%maxfor_prev=zero
262 : :
263 : 14 : electronpositron%lambda=zero
264 : 14 : electronpositron%lifetime=zero
265 : :
266 [ - + + - : 56 : ABI_MALLOC(electronpositron%rhor_ep,(nfft,dtset%nspden))
- + + - +
- + - - +
- + - + ]
267 [ - + + - : 42 : ABI_MALLOC(electronpositron%vha_ep,(nfft))
+ - + - -
+ - + -
+ ]
268 [ + + + - : 119 : ABI_MALLOC(electronpositron%pawrhoij_ep,(mpi_enreg%my_natom*dtset%usepaw))
+ - + + -
+ - + - +
+ + ]
269 : :
270 [ + + ]: 14 : if (dtset%usepaw==1) then
271 : 13 : electronpositron%has_pawrhoij_ep=1
272 [ + - ]: 13 : if (mpi_enreg%my_natom>0) then
273 : : call pawrhoij_alloc(electronpositron%pawrhoij_ep,pawrhoij(1)%cplex_rhoij,pawrhoij(1)%nspden,&
274 : : & pawrhoij(1)%nspinor,pawrhoij(1)%nsppol,dtset%typat,&
275 : : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,&
276 : : & pawtab=pawtab,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,&
277 : : & qphase=pawrhoij(1)%qphase,use_rhoij_=pawrhoij(1)%use_rhoij_,&
278 : 13 : & use_rhoijres=pawrhoij(1)%use_rhoijres)
279 : : end if
280 : 13 : electronpositron%lmmax=0
281 [ + + ]: 26 : do ii=1,dtset%ntypat
282 : 26 : electronpositron%lmmax=max(electronpositron%lmmax,pawtab(ii)%lcut_size**2)
283 : : end do
284 [ - + + - : 52 : ABI_MALLOC(electronpositron%lmselect_ep,(electronpositron%lmmax,mpi_enreg%my_natom))
- + + - +
- + - - +
- + - + ]
285 : : if (maxval(pawtab(1:dtset%ntypat)%usexcnhat)==0.or.(.not.include_nhat_in_gamma)) then
286 [ - + + - : 39 : ABI_MALLOC(electronpositron%nhat_ep,(nfft,dtset%nspden))
- + + - -
+ - + -
+ ]
287 : : end if
288 : : else
289 : 1 : electronpositron%has_pawrhoij_ep=0
290 : 1 : electronpositron%lmmax=0
291 : : end if
292 : :
293 [ + + - + ]: 14 : if (dtset%positron<=-10.or.dtset%posdoppler>0) then
294 : 7 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
295 : 7 : electronpositron%dimcg=mcg
296 : 7 : electronpositron%dimocc=dtset%mband*dtset%nkpt*dtset%nsppol
297 : 7 : electronpositron%dimeigen=dtset%mband*dtset%nkpt*dtset%nsppol
298 [ - + + - : 21 : ABI_MALLOC(electronpositron%cg_ep,(2,electronpositron%dimcg))
+ - + - -
+ - + -
+ ]
299 [ - + + - : 21 : ABI_MALLOC(electronpositron%eigen_ep,(electronpositron%dimeigen))
+ - + - -
+ - + -
+ ]
300 [ - + + - : 21 : ABI_MALLOC(electronpositron%occ_ep,(electronpositron%dimocc))
+ - + - -
+ - + -
+ ]
301 : 7 : electronpositron%dimcprj=0
302 : : ! if (.false.) then !TEMPORARY: will be activated later
303 [ + - + - : 7 : if (dtset%usepaw==1.and.dtset%pawusecp>0.and.dtset%posdoppler>0) then
+ + ]
304 : 3 : electronpositron%dimcprj=my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
305 [ - + ]: 3 : if (mpi_enreg%paral_kgb/=0) electronpositron%dimcprj=electronpositron%dimcprj/mpi_enreg%nproc_band
306 [ - + + - : 840 : ABI_MALLOC(electronpositron%cprj_ep,(dtset%natom,electronpositron%dimcprj))
- + + - +
- + - - +
- + - + +
+ + + ]
307 [ - + + - : 9 : ABI_MALLOC(nlmn,(dtset%natom))
+ - + - -
+ - + ]
308 : 3 : ncpgr=0
309 [ + + ]: 11 : do ii=1,dtset%natom;nlmn(ii)=pawtab(dtset%typat(ii))%lmn_size;end do
310 : 3 : call pawcprj_alloc(electronpositron%cprj_ep,ncpgr,nlmn)
311 [ - + ]: 3 : ABI_FREE(nlmn)
312 : : else
313 [ - + + - : 8 : ABI_MALLOC(electronpositron%cprj_ep,(dtset%natom,electronpositron%dimcprj))
- + - + -
+ - + -
- ]
314 : : end if
315 : : else
316 : 7 : electronpositron%dimcg =0
317 : 7 : electronpositron%dimcprj =0
318 : 7 : electronpositron%dimeigen=0
319 : 7 : electronpositron%dimocc =0
320 : : end if
321 : :
322 : 14 : optfor=0;optstr=0
323 [ + + + - : 14 : if ((dtset%optforces>0.or.dtset%ionmov/=0.or.abs(dtset%toldff)>tiny(0._dp))) optfor=1
- + ]
324 [ + + - + : 14 : if (dtset%optstress>0.and.dtset%iscf>0.and.(dtset%nstep>0.or.ireadwf==1)) optstr=1
- + - - ]
325 : :
326 [ + + ]: 14 : if (optfor>0) then
327 [ - + + - : 9 : ABI_MALLOC(electronpositron%gred_ep,(3,dtset%natom))
+ - + - -
+ - + -
+ ]
328 [ + + + + ]: 95 : electronpositron%gred_ep(:,:)=zero
329 : : end if
330 : :
331 [ + + ]: 14 : if (optstr>0) then
332 [ - + - + ]: 1 : ABI_MALLOC(electronpositron%stress_ep,(6))
333 [ + + ]: 7 : electronpositron%stress_ep(:)=zero
334 : : end if
335 : :
336 : : else !dtset%positron==0
337 : 0 : nullify(electronpositron)
338 : : end if
339 : :
340 [ - + ]: 14 : end subroutine init_electronpositron
341 : : !!***
342 : :
343 : : !----------------------------------------------------------------------
344 : :
345 : : !!****f* m_electronpositron/destroy_electronpositron
346 : : !! NAME
347 : : !! destroy_electronpositron
348 : : !!
349 : : !! FUNCTION
350 : : !! Clean and destroy electronpositron datastructure
351 : : !!
352 : : !! SIDE EFFECTS
353 : : !! electronpositron=<type(electronpositron_type)>=electronpositron datastructure
354 : : !!
355 : : !! SOURCE
356 : :
357 : 14 : subroutine destroy_electronpositron(electronpositron)
358 : :
359 : : !Arguments ------------------------------------
360 : : !scalars
361 : : type(electronpositron_type),pointer :: electronpositron
362 : :
363 : : !************************************************************************
364 : :
365 : : !@electronpositron_type
366 : :
367 [ + - ]: 14 : if (associated(electronpositron)) then
368 : :
369 [ + + ]: 14 : ABI_SFREE(electronpositron%cg_ep)
370 [ + + ]: 14 : ABI_SFREE(electronpositron%eigen_ep)
371 [ + + ]: 14 : ABI_SFREE(electronpositron%occ_ep)
372 [ + - ]: 14 : ABI_SFREE(electronpositron%rhor_ep)
373 [ + + ]: 14 : ABI_SFREE(electronpositron%nhat_ep)
374 [ + - ]: 14 : ABI_SFREE(electronpositron%vha_ep)
375 [ + + ]: 14 : ABI_SFREE(electronpositron%lmselect_ep)
376 [ + + ]: 14 : ABI_SFREE(electronpositron%gred_ep)
377 [ + + ]: 14 : ABI_SFREE(electronpositron%stress_ep)
378 : :
379 [ + + ]: 14 : if (electronpositron%has_pawrhoij_ep/=0) then
380 : 13 : call pawrhoij_free(electronpositron%pawrhoij_ep)
381 : : end if
382 [ + - + + : 91 : ABI_SFREE(electronpositron%pawrhoij_ep)
- + - + -
+ - + - +
- + - + ]
383 : :
384 [ + + ]: 14 : if (electronpositron%dimcprj/=0) then
385 : 3 : call pawcprj_free(electronpositron%cprj_ep)
386 : : end if
387 [ + + + + : 606 : ABI_SFREE(electronpositron%cprj_ep)
- + - + -
+ ]
388 : :
389 : 14 : electronpositron%calctype =0
390 : 14 : electronpositron%particle =-1
391 : 14 : electronpositron%dimcg =0
392 : 14 : electronpositron%dimcprj =0
393 : 14 : electronpositron%dimeigen =0
394 : 14 : electronpositron%dimocc =0
395 : 14 : electronpositron%has_pawrhoij_ep=0
396 : 14 : electronpositron%has_pos_ham =0
397 : 14 : electronpositron%istep =0
398 : 14 : electronpositron%istep_scf =0
399 : :
400 : 14 : electronpositron%posdensity0_limit=.false.
401 : 14 : electronpositron%scf_converged=.false.
402 : :
403 [ - + - + : 14 : ABI_FREE(electronpositron)
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - - -
- - - - -
- + - - -
- - - ]
404 : :
405 : : end if
406 : :
407 : 14 : end subroutine destroy_electronpositron
408 : : !!***
409 : :
410 : : !----------------------------------------------------------------------
411 : :
412 : : !!****f* m_electronpositron/exchange_electronpositron
413 : : !! NAME
414 : : !! exchange_electronpositron
415 : : !!
416 : : !! FUNCTION
417 : : !! Invert electron and positron quantities between an electronpositron datastructure
418 : : !! and current evoving variables
419 : : !! Example: exchange electronpositron%rhor_ep and rhor
420 : : !!
421 : : !! INPUTS
422 : : !! dtset <type(dataset_type)>=all input variables for this dataset
423 : : !! mpi_enreg=information about MPI parallelization
424 : : !! my_natom=number of atoms treated by current proc
425 : : !! nfft=(effective) number of FFT grid points (for this processor)
426 : : !! ngfft(18)=contain all needed information about 3D FFT
427 : : !! npwarr(nkpt)=number of planewaves in basis at this k point
428 : : !! usecprj= 1 if cprj array is stored in memory
429 : : !!
430 : : !! SIDE EFFECTS
431 : : !! cg(2,mcg)=wavefunctions
432 : : !! cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors:
433 : : !! cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector.
434 : : !! electronpositron=<type(electronpositron_type)>=electronpositron datastructure
435 : : !! energies <type(energies_type)>=all part of total energy.
436 : : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
437 : : !! gred(3,natom)=gradients wrt nuclear positions in reduced coordinates
438 : : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
439 : : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
440 : : !! occ(mband*nkpt*nsppol)=occupation number for each band at each k point
441 : : !! paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh
442 : : !! pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
443 : : !! rhog(2,nfft)=Fourier transform of total electron/positron density
444 : : !! rhor(nfft,nspden)=total electron/positron density (el/bohr**3)
445 : : !! stress(6)=components of the stress tensor (hartree/bohr^3) for the
446 : : !! vhartr(nfftf)=array for holding Hartree potential
447 : : !!
448 : : !! SOURCE
449 : :
450 : 6 : subroutine exchange_electronpositron(cg,cprj,dtset,eigen,electronpositron,energies,gred,mcg,mcprj,&
451 : 6 : & mpi_enreg,my_natom,nfft,ngfft,nhat,npwarr,occ,paw_an,pawrhoij,&
452 : 6 : & rhog,rhor,stress,usecprj,vhartr)
453 : :
454 : : !Arguments ------------------------------------
455 : : !scalars
456 : : integer,intent(in) :: mcg,mcprj,my_natom,nfft,usecprj
457 : : type(dataset_type),intent(in) :: dtset
458 : : type(electronpositron_type),pointer :: electronpositron
459 : : type(energies_type),intent(inout) :: energies
460 : : type(MPI_type),intent(in) :: mpi_enreg
461 : : !arrays
462 : : integer,intent(in) :: ngfft(18),npwarr(dtset%nkpt)
463 : : real(dp),intent(inout) :: cg(2,mcg)
464 : : real(dp),intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
465 : : real(dp),intent(inout) :: gred(3,dtset%natom),nhat(nfft,dtset%nspden)
466 : : real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
467 : : real(dp), intent(inout) :: rhog(2,nfft),rhor(nfft,dtset%nspden)
468 : : real(dp),intent(inout) :: stress(6),vhartr(nfft)
469 : : type(pawcprj_type) :: cprj(dtset%natom,mcprj*usecprj)
470 : : type(paw_an_type),intent(inout) :: paw_an(my_natom*dtset%usepaw)
471 : : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*dtset%usepaw)
472 : :
473 : : !Local variables-------------------------------
474 : : !scalars
475 : : integer :: comm,iatom,ib,ibsp,icg,icgb,ifft,ii,ilm,ikpt
476 : : integer :: ispden,isppol,ispinor,me,my_nspinor,nband_k,npw_k,sz1,sz2,sz3
477 : : logical :: ltmp
478 : : real(dp) :: rtmp
479 : : type(energies_type) :: energies_tmp
480 : : !arrays
481 : 6 : integer,allocatable :: nlmn(:),typ(:)
482 : : real(dp) :: ctmp(2)
483 : 6 : type(pawcprj_type),allocatable :: cprj_tmp(:,:)
484 : 6 : type(pawrhoij_type),allocatable :: pawrhoij_tmp(:)
485 : :
486 : : !*********************************************************************
487 : :
488 : 6 : if (associated(electronpositron)) then
489 [ + - ]: 6 : if (electronpositron%particle/=EP_NOTHING) then
490 : :
491 : : ! Type of particle stored
492 [ + - ]: 6 : if (electronpositron%particle==EP_ELECTRON) then
493 : 6 : electronpositron%particle=EP_POSITRON
494 [ # # ]: 0 : else if (electronpositron%particle==EP_POSITRON) then
495 : 0 : electronpositron%particle=EP_ELECTRON
496 : : end if
497 : :
498 : : ! Energies
499 : 6 : ctmp(1)=energies%e_electronpositron
500 : : ! ctmp(2)=energies%edc_electronpositron
501 : 6 : call electronpositron%energies_ep%copy(energies_tmp)
502 : 6 : call energies%copy(electronpositron%energies_ep)
503 : 6 : call energies_tmp%copy(energies)
504 : 6 : energies%e_electronpositron=ctmp(1)
505 : : ! energies%edc_electronpositron=ctmp(2)
506 : 6 : energies%e0_electronpositron=electronpositron%e0
507 : 6 : electronpositron%e0=electronpositron%energies_ep%e0_electronpositron
508 : :
509 : : ! Density and PAW occupation matrix
510 [ + + ]: 12 : do ispden=1,dtset%nspden
511 [ + + ]: 147372 : do ifft=1,nfft
512 : 147366 : rtmp=rhor(ifft,ispden)
513 : 147366 : rhor(ifft,ispden)=electronpositron%rhor_ep(ifft,ispden)
514 : 147372 : electronpositron%rhor_ep(ifft,ispden)=rtmp
515 : : end do
516 [ + - ]: 12 : if (allocated(electronpositron%nhat_ep).and.size(nhat,2)>0) then
517 [ + + ]: 147372 : do ifft=1,nfft
518 : 147366 : rtmp=nhat(ifft,ispden)
519 : 147366 : nhat(ifft,ispden)=electronpositron%nhat_ep(ifft,ispden)
520 : 147372 : electronpositron%nhat_ep(ifft,ispden)=rtmp
521 : : end do
522 : : end if
523 : : end do
524 : 6 : call fourdp(1,rhog,rhor,-1,mpi_enreg,nfft,1,ngfft,0)
525 [ + - + - ]: 6 : if (dtset%usepaw==1.and.my_natom>0) then
526 [ + - ]: 6 : if (electronpositron%has_pawrhoij_ep==1) then
527 [ - + + - : 38 : ABI_MALLOC(pawrhoij_tmp,(my_natom))
+ - - + -
+ - + +
+ ]
528 [ - + + - : 18 : ABI_MALLOC(typ,(my_natom))
+ - - + -
+ ]
529 [ - + + - : 12 : ABI_MALLOC(nlmn,(my_natom))
- + - + -
+ ]
530 [ + + ]: 26 : do iatom=1,my_natom
531 : 20 : typ(iatom)=iatom
532 : 26 : nlmn(iatom)=pawrhoij(iatom)%lmn_size
533 : : end do
534 : : ! Be careful: parallelism over atoms is ignored...
535 : : call pawrhoij_alloc(pawrhoij_tmp,pawrhoij(1)%cplex_rhoij,pawrhoij(1)%nspden,&
536 : : & pawrhoij(1)%nspinor,pawrhoij(1)%nsppol,typ, &
537 : : & lmnsize=nlmn,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,&
538 : : & qphase=pawrhoij(1)%qphase,use_rhoij_=pawrhoij(1)%use_rhoij_,&
539 : 6 : & use_rhoijres=pawrhoij(1)%use_rhoijres)
540 [ - + ]: 6 : ABI_FREE(typ)
541 [ - + ]: 6 : ABI_FREE(nlmn)
542 : 6 : call pawrhoij_copy(pawrhoij,pawrhoij_tmp)
543 : 6 : call pawrhoij_copy(electronpositron%pawrhoij_ep,pawrhoij)
544 : 6 : call pawrhoij_copy(pawrhoij_tmp,electronpositron%pawrhoij_ep)
545 [ - + - - ]: 6 : if (pawrhoij_tmp(1)%ngrhoij>0.and.pawrhoij(1)%ngrhoij==0) then
546 [ # # ]: 0 : do iatom=1,my_natom
547 : 0 : sz1=pawrhoij_tmp(iatom)%ngrhoij
548 : 0 : sz2=pawrhoij_tmp(iatom)%cplex_rhoij*pawrhoij_tmp(iatom)%qphase*pawrhoij_tmp(iatom)%lmn2_size
549 : 0 : sz3=pawrhoij_tmp(iatom)%nspden
550 [ # # # # : 0 : ABI_MALLOC(pawrhoij(iatom)%grhoij,(sz1,sz2,sz3))
# # # # #
# # # # #
# # # # #
# # # ]
551 [ # # # # : 0 : pawrhoij(iatom)%grhoij(:,:,:)=pawrhoij_tmp(iatom)%grhoij(:,:,:)
# # ]
552 : : end do
553 : : end if
554 [ + - - + ]: 6 : if (pawrhoij_tmp(1)%use_rhoijres>0.and.pawrhoij(1)%use_rhoijres==0) then
555 [ # # ]: 0 : do iatom=1,my_natom
556 : 0 : sz1=pawrhoij_tmp(iatom)%cplex_rhoij*pawrhoij_tmp(iatom)%qphase*pawrhoij_tmp(iatom)%lmn2_size
557 : 0 : sz2=pawrhoij_tmp(iatom)%nspden
558 [ # # # # : 0 : ABI_MALLOC(pawrhoij(iatom)%rhoijres,(sz1,sz2))
# # # # #
# # # # #
# # # # ]
559 [ # # # # ]: 0 : pawrhoij(iatom)%rhoijres(:,:)=pawrhoij_tmp(iatom)%rhoijres(:,:)
560 : : end do
561 : : end if
562 [ - + - - ]: 6 : if (pawrhoij_tmp(1)%use_rhoij_>0.and.pawrhoij(1)%use_rhoij_==0) then
563 [ # # ]: 0 : do iatom=1,my_natom
564 : 0 : sz1=pawrhoij_tmp(iatom)%cplex_rhoij*pawrhoij_tmp(iatom)%qphase*pawrhoij_tmp(iatom)%lmn2_size
565 : 0 : sz2=pawrhoij_tmp(iatom)%nspden
566 [ # # # # : 0 : ABI_MALLOC(pawrhoij(iatom)%rhoij_,(sz1,sz2))
# # # # #
# # # # #
# # # # ]
567 [ # # # # ]: 0 : pawrhoij(iatom)%rhoij_(:,:)=pawrhoij_tmp(iatom)%rhoij_(:,:)
568 : : end do
569 : : end if
570 [ + - - + ]: 6 : if (pawrhoij_tmp(1)%lmnmix_sz>0.and.pawrhoij(1)%lmnmix_sz==0) then
571 [ # # ]: 0 : do iatom=1,my_natom
572 [ # # # # : 0 : ABI_MALLOC(pawrhoij(iatom)%kpawmix,(pawrhoij_tmp(iatom)%lmnmix_sz))
# # # # #
# # # #
# ]
573 [ # # ]: 0 : pawrhoij(iatom)%kpawmix(:)=pawrhoij_tmp(iatom)%kpawmix(:)
574 : : end do
575 : : end if
576 : 6 : call pawrhoij_free(pawrhoij_tmp)
577 [ + - + + : 26 : ABI_FREE(pawrhoij_tmp)
- + - + -
+ - + - +
- + - + ]
578 : : else
579 [ # # ]: 0 : do iatom=1,my_natom
580 [ # # # # ]: 0 : pawrhoij(iatom)%rhoijp=zero
581 : : end do
582 : : end if
583 : : end if
584 : :
585 : : ! Hartree potential
586 [ + + ]: 147372 : do ifft=1,nfft
587 : 147366 : rtmp=vhartr(ifft)
588 : 147366 : vhartr(ifft)=electronpositron%vha_ep(ifft)
589 : 147372 : electronpositron%vha_ep(ifft)=rtmp
590 : : end do
591 : :
592 : : ! PAW LM-moment selection flags
593 [ + - + - ]: 6 : if (dtset%usepaw==1.and.my_natom>0) then
594 [ + - ]: 6 : if (electronpositron%lmmax>0) then
595 [ + + ]: 26 : do iatom=1,my_natom
596 [ + + ]: 206 : do ilm=1,paw_an(iatom)%lm_size
597 : 180 : ltmp=electronpositron%lmselect_ep(ilm,iatom)
598 : 180 : electronpositron%lmselect_ep(ilm,iatom)=paw_an(iatom)%lmselect(ilm)
599 : 200 : paw_an(iatom)%lmselect(ilm)=ltmp
600 : : end do
601 : : end do
602 : : else
603 [ # # ]: 0 : do iatom=1,my_natom
604 [ # # ]: 0 : paw_an(iatom)%lmselect(:)=.true.
605 : : end do
606 : : end if
607 : : end if
608 : :
609 : : ! Wave-functions
610 [ + + ]: 6 : if (electronpositron%dimcg>0) then
611 [ + + ]: 186193 : do ii=1,electronpositron%dimcg
612 [ + + ]: 558564 : ctmp(1:2)=electronpositron%cg_ep(1:2,ii)
613 [ + + ]: 558564 : electronpositron%cg_ep(1:2,ii)=cg(1:2,ii)
614 [ + + ]: 558569 : cg(1:2,ii)=ctmp(1:2)
615 : : end do
616 : : else
617 : 1 : icg=0
618 : 1 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
619 : 1 : comm=mpi_enreg%comm_cell
620 : 1 : me=xmpi_comm_rank(comm)
621 [ + + ]: 2 : do isppol=1,dtset%nsppol
622 [ + + ]: 6 : do ikpt=1,dtset%nkpt
623 : 4 : npw_k=npwarr(ikpt);nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
624 [ - + ]: 4 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
625 : : icgb=icg;ibsp=0
626 [ + + ]: 64 : do ib=1,nband_k
627 [ + + + + ]: 86910 : cg(:,icgb+1:icgb+my_nspinor*npw_k)=zero
628 [ + + ]: 120 : do ispinor=1,my_nspinor
629 [ + - ]: 120 : ibsp=ibsp+1;if (ibsp<my_nspinor*npw_k) cg(1,icgb+ibsp)=one
630 : : end do
631 : 64 : icgb=icgb+my_nspinor*npw_k
632 : : end do
633 [ - + ]: 5 : if (dtset%mkmem/=0) icg=icg+my_nspinor*npw_k*nband_k
634 : : end do
635 : : end do
636 : : end if
637 [ + - ]: 6 : if (dtset%usepaw==1) then
638 [ + + ]: 6 : if(electronpositron%dimcprj>0) then
639 [ - + + - : 9 : ABI_MALLOC(nlmn,(dtset%natom))
+ - + - -
+ - + -
+ ]
640 [ - + + - : 840 : ABI_MALLOC(cprj_tmp,(dtset%natom,electronpositron%dimcprj))
- + + - +
- + - - +
- + + + +
+ ]
641 [ + + ]: 11 : do iatom=1,dtset%natom;nlmn(iatom)=cprj(iatom,1)%nlmn;end do
642 : 3 : call pawcprj_alloc(cprj_tmp,cprj(1,1)%ncpgr,nlmn)
643 [ - + ]: 3 : ABI_FREE(nlmn)
644 : 3 : call pawcprj_copy(electronpositron%cprj_ep,cprj_tmp)
645 : 3 : call pawcprj_copy(cprj,electronpositron%cprj_ep)
646 : 3 : call pawcprj_copy(cprj_tmp,cprj)
647 : 3 : call pawcprj_free(cprj_tmp)
648 [ + - + + : 595 : ABI_FREE(cprj_tmp)
- + - + -
+ ]
649 : : else
650 : : !TO BE ACTIVATED WHEN cprj IS PRESENT
651 : : ! call pawcprj_set_zero(cprj)
652 : : end if
653 : : end if
654 : :
655 : : ! Eigenvalues
656 [ + + ]: 6 : if (electronpositron%dimeigen>0) then
657 [ + + ]: 391 : do ii=1,electronpositron%dimeigen
658 : 386 : rtmp=eigen(ii)
659 : 386 : eigen(ii)=electronpositron%eigen_ep(ii)
660 : 391 : electronpositron%eigen_ep(ii)=rtmp
661 : : end do
662 : : else
663 [ + + ]: 61 : eigen(:)=9.99999_dp
664 : : end if
665 : :
666 : : ! Occupations
667 [ + + ]: 6 : if (electronpositron%dimocc>0) then
668 [ + + ]: 391 : do ii=1,electronpositron%dimocc
669 : 386 : rtmp=occ(ii)
670 : 386 : occ(ii)=electronpositron%occ_ep(ii)
671 : 391 : electronpositron%occ_ep(ii)=rtmp
672 : : end do
673 : : else
674 [ + + ]: 61 : occ(:)=9.99999_dp
675 : : end if
676 : :
677 : : ! Forces
678 [ + + ]: 6 : if (allocated(electronpositron%gred_ep)) then
679 [ + + ]: 5 : do iatom=1,dtset%natom
680 [ + + ]: 17 : electronpositron%gred_ep(1:3,iatom)=gred(1:3,iatom)-electronpositron%gred_ep(1:3,iatom)
681 : : end do
682 : : end if
683 : :
684 : : ! Stresses
685 [ - + ]: 6 : if (allocated(electronpositron%stress_ep)) then
686 [ # # ]: 0 : electronpositron%stress_ep(1:6)=stress(1:6)-electronpositron%stress_ep(1:6)
687 : : end if
688 : :
689 : : end if
690 : : end if
691 : :
692 [ + - - + : 12 : end subroutine exchange_electronpositron
- + - - -
- - - - -
- - - - -
- - + - +
- - - - -
- ]
693 : : !!***
694 : :
695 : : !----------------------------------------------------------------------
696 : :
697 : : !!****f* m_electronpositron/electronpositron_calctype
698 : : !! NAME
699 : : !! electronpositron_calctype
700 : : !!
701 : : !! FUNCTION
702 : : !! Returns the value of the calculation type from an electronpositron
703 : : !! structure (can be eventually unassociated)
704 : : !!
705 : : !! INPUTS
706 : : !! electronpositron=<type(electronpositron_type)>=electronpositron datastructure
707 : : !!
708 : : !! SOURCE
709 : :
710 : 182466 : integer function electronpositron_calctype(electronpositron)
711 : :
712 : : !Arguments ------------------------------------
713 : : !scalars
714 : : type(electronpositron_type),pointer :: electronpositron
715 : :
716 : : !************************************************************************
717 : :
718 [ + + ]: 182199 : if (associated(electronpositron)) then
719 : 3034 : electronpositron_calctype=electronpositron%calctype
720 : : else
721 : : electronpositron_calctype=0
722 : : end if
723 : :
724 : :
725 : 182199 : end function electronpositron_calctype
726 : : !!***
727 : :
728 : : !----------------------------------------------------------------------
729 : :
730 : : !!****f* ABINIT/rhohxcpositron
731 : : !! NAME
732 : : !! rhohxcpositron
733 : : !!
734 : : !! FUNCTION
735 : : !! Calculate the electrons/positron correlation term for the positron
736 : : !!
737 : : !! INPUTS
738 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
739 : : !! mpi_enreg=information about MPI parallelization
740 : : !! nfft=(effective) number of FFT grid points (for this processor)
741 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
742 : : !! nhat(nfft,nspden*usepaw)= -PAW only- compensation density
743 : : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
744 : : !! nspden=number of spin density components
745 : : !! n3xccc=dimension of the xccc3d array (0 or nfft).
746 : : !! paral_kgb=flag for (k,band,FFT) parallelism
747 : : !! rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
748 : : !! ucvol = unit cell volume (Bohr**3)
749 : : !! usexcnhat= -PAW only- flag controlling use of compensation density in Vxc
750 : : !! usepaw=flag for PAW
751 : : !! xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
752 : : !! xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
753 : : !!
754 : : !! OUTPUT
755 : : !! electronpositron%e_xc=electron-positron XC energy
756 : : !! electronpositron%e_xcdc=Double-counting electron-positron XC energy
757 : : !! strsxc(6)= contribution of xc to stress tensor (hartree/bohr^3),
758 : : !! vhartr(nfft)=Hartree potential (returned if option/=0 and option/=10)
759 : : !! vxcapn=XC electron-positron XC potential for the positron
760 : : !! vxcavg=unit cell average of Vxc = (1/ucvol) Int [Vxc(r) d^3 r].
761 : : !! kxcapn(nfft,nkxc)=electron-positron XC kernel (returned only if nkxc/=0)
762 : : !!
763 : : !! SIDE EFFECTS
764 : : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
765 : : !!
766 : : !! SOURCE
767 : :
768 : 267 : subroutine rhohxcpositron(electronpositron,gprimd,kxcapn,mpi_enreg,nfft,ngfft,nhat,nkxc,nspden,n3xccc,&
769 : 267 : & paral_kgb,rhor,strsxc,ucvol,usexcnhat,usepaw,vhartr,vxcapn,vxcavg,xccc3d,xc_denpos)
770 : :
771 : : !Arguments ------------------------------------
772 : : !scalars
773 : : integer,intent(in) :: nfft,nkxc,nspden,n3xccc,paral_kgb,usexcnhat,usepaw
774 : : real(dp),intent(in) :: ucvol,xc_denpos
775 : : real(dp),intent(out) :: vxcavg
776 : : type(electronpositron_type),pointer :: electronpositron
777 : : !arrays
778 : : integer,intent(in) :: ngfft(18)
779 : : real(dp),intent(in) :: gprimd(3,3)
780 : : real(dp),intent(in) :: nhat(nfft,nspden*usepaw),rhor(nfft,nspden),xccc3d(n3xccc)
781 : : real(dp),intent(out) :: kxcapn(nfft,nkxc),strsxc(6),vhartr(nfft),vxcapn(nfft,nspden)
782 : : type(MPI_type),intent(in) :: mpi_enreg
783 : :
784 : : !Local variables-------------------------------
785 : : !scalars
786 : : integer :: cplex,ierr,ifft,ishift,iwarn,iwarnp,nfftot,ngr,ngrad,nspden_ep
787 : : real(dp) :: exc,excdc,strdiag
788 : : character(len=500) :: message
789 : : !arrays
790 : : real(dp),parameter :: qphon(3)=(/0._dp,0._dp,0._dp/)
791 : : real(dp) :: vxcavg_tmp(1)
792 [ + - ]: 267 : real(dp),allocatable :: fxcapn(:),grho2apn(:),rhoe(:,:,:),rhop(:,:),rhotote(:),vxc_ep(:),vxcgr_ep(:)
793 : :
794 : : ! *************************************************************************
795 : :
796 [ - + ]: 267 : if (electronpositron_calctype(electronpositron)/=1) then
797 : 0 : message = 'Only electronpositron%calctype=1 allowed !'
798 : 0 : ABI_BUG(message)
799 : : end if
800 : :
801 [ - + ]: 267 : if (nkxc>3) then
802 : 0 : message = 'nkxc>3 (Kxc for GGA) not yet implemented !'
803 : 0 : ABI_ERROR(message)
804 : : end if
805 : :
806 : : !Hartree potential of the positron is zero
807 [ + + ]: 28747977 : vhartr=zero
808 : :
809 : : !Some allocations/inits
810 [ - + ]: 267 : ngrad=1;if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2
811 [ - + ]: 267 : ngr=0;if (ngrad==2) ngr=nfft
812 [ - + + - : 801 : ABI_MALLOC(fxcapn,(nfft))
+ - + - -
+ - + ]
813 [ + - - - : 801 : ABI_MALLOC(grho2apn,(ngr))
+ - - + -
+ - + ]
814 : 267 : nspden_ep=1;cplex=1;ishift=0
815 : 267 : iwarn=0;iwarnp=1
816 : :
817 : : !Compute total electronic density
818 [ - + + - : 534 : ABI_MALLOC(rhotote,(nfft))
- + - + ]
819 [ + + ]: 28747977 : rhotote(:)=electronpositron%rhor_ep(:,1)
820 [ + + + + ]: 28433300 : if (n3xccc>0) rhotote(:)=rhotote(:)+xccc3d(:)
821 [ + + + + : 26867562 : if (usepaw==1.and.usexcnhat==0) rhotote(:)=rhotote(:)-electronpositron%nhat_ep(:,1)
+ + ]
822 : :
823 : : !Extra total electron/positron densities; compute gradients for GGA
824 [ - + + - : 1068 : ABI_MALLOC(rhoe,(nfft,nspden_ep,ngrad**2))
- + + - +
- + - - +
- + ]
825 [ - + + - : 534 : ABI_MALLOC(rhop,(nfft,nspden_ep))
+ - - + -
+ ]
826 : 267 : call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_ep,qphon,rhotote,rhoe)
827 [ - + - - ]: 267 : if (ngrad==2) grho2apn(:)=rhoe(:,1,2)**2+rhoe(:,1,3)**2+rhoe(:,1,4)**2
828 [ + + + + : 55615468 : rhop(:,1)=rhor(:,1);if (usepaw==1.and.usexcnhat==0) rhop(:,1)=rhop(:,1)-nhat(:,1)
+ + + + ]
829 [ - + ]: 267 : ABI_FREE(rhotote)
830 : :
831 : : !Make the densities positive
832 : 267 : call mkdenpos(iwarn ,nfft,nspden_ep,1,rhoe(:,1,1),xc_denpos)
833 [ + + ]: 267 : if (.not.electronpositron%posdensity0_limit) then
834 : 7 : call mkdenpos(iwarnp,nfft,nspden_ep,1,rhop,xc_denpos)
835 : : end if
836 : :
837 : : !Compute electron-positron Vxc_pos, Vxc_el, Fxc, Kxc, ...
838 [ - + + - : 534 : ABI_MALLOC(vxc_ep,(nfft))
- + - + ]
839 [ + - - - : 534 : ABI_MALLOC(vxcgr_ep,(ngr))
- + - + ]
840 [ + + ]: 267 : if (nkxc==0) then
841 : : call xcpositron(fxcapn,grho2apn,electronpositron%ixcpositron,ngr,nfft,electronpositron%posdensity0_limit,&
842 : 143 : & rhoe(:,1,1),rhop(:,1),vxc_ep,vxcgr_ep,vxcapn)
843 : : else
844 : : call xcpositron(fxcapn,grho2apn,electronpositron%ixcpositron,ngr,nfft,electronpositron%posdensity0_limit,&
845 : 124 : & rhoe(:,1,1),rhop(:,1),vxc_ep,vxcgr_ep,vxcapn,dvxce=kxcapn)
846 : : end if
847 [ - + ]: 267 : ABI_FREE(rhoe)
848 [ - + ]: 267 : ABI_FREE(vxc_ep)
849 [ - + ]: 267 : ABI_FREE(vxcgr_ep)
850 [ - + ]: 267 : ABI_FREE(grho2apn)
851 : :
852 : : !Store Vxc and Kxc according to spin components
853 [ - + - - ]: 267 : if (nspden>=2) vxcapn(:,2)=vxcapn(:,1)
854 [ - + - - : 267 : if (nspden==4) vxcapn(:,3:4)=zero
- - ]
855 [ - + ]: 267 : if (nkxc==3) then
856 [ # # ]: 0 : kxcapn(:,1)=two*kxcapn(:,1)
857 [ # # ]: 0 : kxcapn(:,2)=kxcapn(:,1)
858 [ # # ]: 0 : kxcapn(:,3)=kxcapn(:,1)
859 : : end if
860 : :
861 : : !Compute XC energies and contribution to stress tensor
862 : 267 : electronpositron%e_xc =zero
863 : 267 : electronpositron%e_xcdc=zero
864 : 267 : strdiag=zero
865 [ + + ]: 1068 : nfftot=PRODUCT(ngfft(1:3))
866 [ + + ]: 28747977 : do ifft=1,nfft
867 : 28747710 : electronpositron%e_xc =electronpositron%e_xc +fxcapn(ifft)
868 : 28747710 : electronpositron%e_xcdc=electronpositron%e_xcdc+vxcapn(ifft,1)*rhor(ifft,1)
869 : : ! strdiag=strdiag+fxcapn(ifft) ! Already stored in rhotoxc !
870 : 28747977 : strdiag=strdiag-vxcapn(ifft,1)*rhop(ifft,1)
871 : : end do
872 [ + + + + ]: 267 : if (usepaw==1.and.usexcnhat==0) then
873 [ + + ]: 26867491 : do ifft=1,nfft
874 : 26867491 : electronpositron%e_xcdc=electronpositron%e_xcdc-vxcapn(ifft,1)*nhat(ifft,1)
875 : : end do
876 : : end if
877 : 267 : electronpositron%e_xc =electronpositron%e_xc *ucvol/dble(nfftot)
878 : 267 : electronpositron%e_xcdc=electronpositron%e_xcdc*ucvol/dble(nfftot)
879 : 267 : strdiag=strdiag/dble(nfftot)
880 [ - + ]: 267 : ABI_FREE(fxcapn)
881 [ - + ]: 267 : ABI_FREE(rhop)
882 : :
883 : : !Store stress tensor
884 [ + + ]: 1068 : strsxc(1:3)=strdiag
885 [ + + ]: 1068 : strsxc(4:6)=zero
886 : :
887 : : !Reduction in case of parallelism
888 [ - + ]: 267 : if(mpi_enreg%paral_kgb==1)then
889 [ # # ]: 0 : if(paral_kgb/=0)then
890 : 0 : exc=electronpositron%e_xc;excdc=electronpositron%e_xcdc
891 : 0 : call xmpi_sum(exc ,mpi_enreg%comm_fft,ierr)
892 : 0 : call xmpi_sum(excdc,mpi_enreg%comm_fft,ierr)
893 : 0 : electronpositron%e_xc=exc;electronpositron%e_xcdc=excdc
894 : 0 : call xmpi_sum(strsxc,mpi_enreg%comm_fft,ierr)
895 : : end if
896 : : end if
897 : :
898 : : !Compute vxcavg
899 : 267 : call mean_fftr(vxcapn(:,1),vxcavg_tmp,nfft,nfftot,1)
900 : 267 : vxcavg=vxcavg_tmp(1)
901 : :
902 : 267 : end subroutine rhohxcpositron
903 : : !!***
904 : :
905 [ # # # # : 0 : END MODULE m_electronpositron
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
906 : : !!***
|
