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