Branch data Line data Source code
1 : : !!****m* ABINIT/m_rhotoxc
2 : : !! NAME
3 : : !! m_rhotox
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, MF, GZ, DRH, MT, SPr)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : module m_rhotoxc
22 : :
23 : : use defs_basis
24 : : use m_xmpi
25 : : use m_abicore
26 : : use m_errors
27 : : use m_cgtools
28 : : use m_xcdata
29 : : use m_xc_vdw
30 : : use libxc_functionals
31 : :
32 : : use defs_abitypes, only : MPI_type
33 : : use m_time, only : timab
34 : : use m_geometry, only : metric
35 : : use m_electronpositron, only : electronpositron_type,electronpositron_calctype
36 : : use m_xcpositron, only : xcpositron
37 : : use m_drivexc, only : size_dvxc,drivexc,xcmult,mkdenpos
38 : : use m_xclda, only : xctfw
39 : : use m_xctk, only : xcden, xcpot
40 : :
41 : : implicit none
42 : :
43 : : private
44 : : !!***
45 : :
46 : : public :: rhotoxc
47 : : !!***
48 : :
49 : : contains
50 : : !!***
51 : :
52 : : !!****f* ABINIT/rhotoxc
53 : : !! NAME
54 : : !! rhotoxc
55 : : !!
56 : : !! FUNCTION
57 : : !! Start from the density or spin-density, and
58 : : !! compute xc correlation potential and energies.
59 : : !! Eventually compute xc kernel (if option=-2, 2, 3, 10 or 12).
60 : : !! Cannot be used with wavelets.
61 : : !!
62 : : !! INPUTS
63 : : !! mpi_enreg=information about MPI parallelization
64 : : !! nfft=(effective) number of FFT grid points (for this processor)
65 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
66 : : !! nhat(nfft,xcdata%nspden*nhatdim)= -PAW only- compensation density
67 : : !! nhatdim= -PAW only- 0 if nhat array is not used ; 1 otherwise
68 : : !! nhatgr(nfft,xcdata%nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density
69 : : !! nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise
70 : : !! nkxc=second dimension of the kxc array. If /=0,
71 : : !! the exchange-correlation kernel must be computed.
72 : : !! non_magnetic_xc= if true, handle density/potential as non-magnetic (even if it is)
73 : : !! n3xccc=dimension of the xccc3d array (0 or nfft or cplx*nfft).
74 : : !! option=0 or 1 for xc only (exc, vxc, strsxc),
75 : : !! 2 for xc and kxc (no paramagnetic part if xcdata%nspden=1)
76 : : !! 10 for xc and kxc with only partial derivatives wrt density part (d2Exc/drho^2)
77 : : !! 12 for xc and kxc with only partial derivatives wrt density part (d2Exc/drho^2)
78 : : !! and, in the case of hybrid functionals, substitution of the hybrid functional
79 : : !! by the related auxiliary GGA functional for the computation of the xc kernel (not for other quantities)
80 : : !! 3 for xc, kxc and k3xc
81 : : !! -2 for xc and kxc (with paramagnetic part if xcdata%nspden=1)
82 : : !! rhor(nfft,xcdata%nspden)=electron density in real space in electrons/bohr**3
83 : : !! (total in first half and spin-up in second half if xcdata%nspden=2)
84 : : !! (total in first comp. and magnetization in comp. 2 to 4 if xcdata%nspden=4)
85 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
86 : : !! usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
87 : : !! [vhartr(nfft)=Hartree potential (only needed for Fermi-Amaldi functional)]
88 : : !! xcdata <type(xcdata_type)>=storage for different input variables and derived parameters needed to compute the XC functional
89 : : !! xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
90 : : !!
91 : : !! === optional inputs ===
92 : : !! [add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy
93 : : !! [taur(nfftf,xcdata%nspden*xcdata%usekden)]=array for kinetic energy density
94 : : !! [xc_funcs(2)]= <type(libxc_functional_type)>, optional : libxc XC functionals. Must be coherent with xcdata.
95 : : !! [xccctau3d(n3xccc)]=3D core electron kinetic energy density for XC core correction (bohr^-3)
96 : : !!
97 : : !! OUTPUT
98 : : !! bigexc=returned exchange and correlation energy (hartree).
99 : : !! vxc(nfft,xcdata%nspden)=xc potential
100 : : !! (spin up in first half and spin down in second half if xcdata%nspden=2)
101 : : !! (v^11, v^22, Re[V^12], Im[V^12] if xcdata%nspden=4)
102 : : !! vxcavg=<Vxc>=unit cell average of Vxc = (1/ucvol) Int [Vxc(r) d^3 r].
103 : : !!
104 : : !! === Only if abs(option)=2, -2, 3, 10, 12 (in case 12, for hybrids, substitution of the related GGA) ===
105 : : !! kxc(nfft,nkxc)=exchange and correlation kernel (returned only if nkxc/=0)
106 : : !! Content of Kxc array:
107 : : !! ===== if LDA
108 : : !! if xcdata%nspden==1: kxc(:,1)= d2Exc/drho2
109 : : !! that is 1/2 ( d2Exc/drho_up drho_up + d2Exc/drho_up drho_dn )
110 : : !! kxc(:,2)= d2Exc/drho_up drho_dn
111 : : !! if xcdata%nspden>=2: kxc(:,1)= d2Exc/drho_up drho_up
112 : : !! kxc(:,2)= d2Exc/drho_up drho_dn
113 : : !! kxc(:,3)= d2Exc/drho_dn drho_dn
114 : : !! ===== if GGA or mGGA
115 : : !! if xcdata%nspden==1:
116 : : !! kxc(:,1)= d2Exc/drho2
117 : : !! kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)|
118 : : !! kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho
119 : : !! kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| )
120 : : !! kxc(:,5)= gradx(rho)
121 : : !! kxc(:,6)= grady(rho)
122 : : !! kxc(:,7)= gradz(rho)
123 : : !! if xcdata%nspden>=2:
124 : : !! kxc(:,1)= d2Exc/drho_up drho_up
125 : : !! kxc(:,2)= d2Exc/drho_up drho_dn
126 : : !! kxc(:,3)= d2Exc/drho_dn drho_dn
127 : : !! kxc(:,4)= 1/|grad(rho_up)| dEx/d|grad(rho_up)|
128 : : !! kxc(:,5)= 1/|grad(rho_dn)| dEx/d|grad(rho_dn)|
129 : : !! kxc(:,6)= 1/|grad(rho_up)| d2Ex/d|grad(rho_up)| drho_up
130 : : !! kxc(:,7)= 1/|grad(rho_dn)| d2Ex/d|grad(rho_dn)| drho_dn
131 : : !! kxc(:,8)= 1/|grad(rho_up)| * d/d|grad(rho_up)| ( 1/|grad(rho_up)| dEx/d|grad(rho_up)| )
132 : : !! kxc(:,9)= 1/|grad(rho_dn)| * d/d|grad(rho_dn)| ( 1/|grad(rho_dn)| dEx/d|grad(rho_dn)| )
133 : : !! kxc(:,10)=1/|grad(rho)| dEc/d|grad(rho)|
134 : : !! kxc(:,11)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_up
135 : : !! kxc(:,12)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_dn
136 : : !! kxc(:,13)=1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dEc/d|grad(rho)| )
137 : : !! kxc(:,14)=gradx(rho_up)
138 : : !! kxc(:,15)=gradx(rho_dn)
139 : : !! kxc(:,16)=grady(rho_up)
140 : : !! kxc(:,17)=grady(rho_dn)
141 : : !! kxc(:,18)=gradz(rho_up)
142 : : !! kxc(:,19)=gradz(rho_dn)
143 : : !! Note about mGGA: 2nd derivatives involving Tau or Laplacian are not output
144 : : !!
145 : : !! === Only if abs(option)=3 ===
146 : : !! [k3xc(nfft,nk3xc)]= -- optional -- third derivative of the XC energy functional of the density,
147 : : !! at each point of the real space grid (only in the LDA or LSDA)
148 : : !! Content of K3xc array:
149 : : !! ===== if LDA
150 : : !! if xcdata%nspden==1: k3xc(:,1)= d3Exc/drho3
151 : : !! if xcdata%nspden>=2, k3xc(:,1)= d3Exc/drho_up drho_up drho_up
152 : : !! k3xc(:,2)= d3Exc/drho_up drho_up drho_dn
153 : : !! k3xc(:,3)= d3Exc/drho_up drho_dn drho_dn
154 : : !! k3xc(:,4)= d3Exc/drho_dn drho_dn drho_dn
155 : : !!
156 : : !! === Additional optional output ===
157 : : !! [exc_vdw_out]= vdW-DF contribution to bigexc (hartree)
158 : : !! [vxctau(nfft,xcdata%nspden,4*xcdata%usekden)]=(only for meta-GGA)=
159 : : !! vxctau(:,:,1): derivative of XC energy density with respect to kinetic energy density (depsxcdtau).
160 : : !! vxctau(:,:,2:4): gradient of vxctau (gvxctau)
161 : : !! [strsxc(6)]= contribution of xc to stress tensor (hartree/bohr^3),
162 : : !! given in order (1,1), (2,2), (3,3), (3,2), (3,1), (2,1).
163 : : !! Explicitely : strsxc(mu,nu) = (1/N) Sum(i=1,N)
164 : : !! { delta(mu,nu) * [ exc(i)rhotot(i)
165 : : !! - depsxc_drho(up,i)*rhor(up,i)-depsxc_drho(dn,i)*rhor(dn,i)]
166 : : !! - gradrho(up,mu)*gradrho(up,nu) * depsxc_dgradrho(up,i) / gradrho(up,i)
167 : : !! - gradrho(dn,mu)*gradrho(dn,nu) * depsxc_dgradrho(dn,i) / gradrho(dn,i) }
168 : : !! (note: there are additional terms in case of metaGGA)
169 : : !! === For the TB09 XC functional (modified Becke-Johnson) ===
170 : : !! [grho1_over_rho1]=Integral of |Grad(rho^1)|/rho^1 over the augmentation region
171 : : !! Used to compute the c parameter of the TB09 XC functional
172 : : !!
173 : : !! SIDE EFFECTS
174 : : !! electronpositron <type(electronpositron_type)>= -- optional argument -- quantities for the electron-positron annihilation
175 : : !!
176 : : !! NOTES
177 : : !! Start from the density, and compute Hartree (if option>=1) and xc correlation potential and energies.
178 : : !! Eventually compute xc kernel (if option=-2, 2, 3, 10 or 12 - in the latter case, substitution by the related GGA kernel).
179 : : !! Allows a variety of exchange-correlation functionals
180 : : !! according to ixc. Here is a list of allowed values.
181 : : !! subroutine name
182 : : !! <0 means use of libxc
183 : : !! 0 means no xc applied (usually for testing)
184 : : !! *LDA,LSD
185 : : !! 1 means new Teter (4/93) with spin-pol option xcspol
186 : : !! 2 means Perdew-Zunger-Ceperley-Alder xcpzca
187 : : !! 3 means old Teter (4/91) fit to Ceperley-Alder data xctetr
188 : : !! 4 means Wigner xcwign
189 : : !! 5 means Hedin-Lundqvist xchelu
190 : : !! 6 means "X-alpha" xc xcxalp
191 : : !! 7 mean Perdew-Wang 92 LSD fit to Ceperley-Alder data xcpbe
192 : : !! 8 mean Perdew-Wang 92 LSD , exchange-only xcpbe
193 : : !! 9 mean Perdew-Wang 92 Ex+Ec_RPA energy xcpbe
194 : : !! 10 means RPA LSD energy (only the energy !!) xcpbe
195 : : !! *GGA
196 : : !! 11 means Perdew-Burke-Ernzerhof GGA functional xcpbe
197 : : !! 12 means x-only Perdew-Burke-Ernzerhof GGA functional xcpbe
198 : : !! 13 means LDA (ixc==7), except that the xc potential
199 : : !! is given within the van Leeuwen-Baerends GGA xclb
200 : : !! 14 means revPBE GGA functional xcpbe
201 : : !! 15 means RPBE GGA functional xcpbe
202 : : !! 16 means HCTH GGA functional xchcth
203 : : !! 23 means WC GGA functional xcpbe
204 : : !! 24 means C09x GGA exchange functional xcpbe
205 : : !! *Fermi-Amaldi
206 : : !! 20 means Fermi-Amaldi correction
207 : : !! 21 means Fermi-Amaldi correction with LDA(ixc=1) kernel
208 : : !! 22 means Fermi-Amaldi correction with hybrid BPG kernel
209 : : !! *Hybrid GGA
210 : : !! 41 means PBE0-1/4 xcpbe
211 : : !! 42 means PBE0-1/3 xcpbe
212 : : !! *Temperature-dependant exchange-correlation functionals
213 : : !! 50 means IIT xc (TLDA functional) xciit
214 : : !! 51 means KSDT xc (TLDA functional) xcksdt
215 : : !! 60 means KDT16 xc (PBE TGGA functional) xckdt16
216 : : !!
217 : : !! NOTE: please update echo_xc_name.F90 if you add new functional (apart from libxc)
218 : : !!
219 : : !! Allow for improved xc quadrature (intxc=1) by using the usual FFT grid
220 : : !! as well as another, shifted, grid, and combining both results.
221 : : !! Spin-polarization is allowed only with ixc=0, 1, and GGAs until now.
222 : : !!
223 : : !! To make the variable names easier to understand, a rule notation is tentatively proposed here:
224 : : !! rho ---> means density
225 : : !! tau ---> means kinetic energy density
226 : : !! exc ---> means exchange-correlation energy density per particle
227 : : !! rhoexc ---> means rho*exc == exchange-correlation energy density
228 : : !! vxc ---> means exchange-correlation potential
229 : : !! bigexc ---> means exchange-correlation energy E_xc
230 : : !! bigsxc ---> means exchange-correlation entropy S_xc (zero when standard xc functionals are used)
231 : : !! m_norm ---> means norm of magnetization
232 : : !!
233 : : !! In the case where finite-temperature exchange-correlation functionals are used:
234 : : !! exc_b ---> means exchange-correlation free energy density per particle
235 : : !! rhoexc ---> means rho*exc == exchange-correlation free energy density
236 : : !! bigexc ---> means exchange-correlation internal energy E_xc
237 : : !! bigsxc ---> means exchange-correlation entropy S_xc
238 : : !! tsxc_b ---> means exchange-correlation entropy energy density per particle
239 : : !! rhotsxc --> means rho*tsxc == exchange-correlation entropy energy density
240 : : !!
241 : : !! g... --> means gradient of something (e.g. : grho --> means gradient of electron density)
242 : : !! g...2 -> means square norm of gradient of something (e.g. : grho2 -> means square norm of gradient of electron density)
243 : : !! l... --> means laplacian of something (e.g. : lrho --> means laplacian of electron density)
244 : : !! d...d... --> means derivative of something with regards to something else.
245 : : !! (d2...d...d... ---> means second derivative of ... with regards to ... and to ...) etc...
246 : : !! d... --> without the occurence of the second "d" means that this is an array of
247 : : !! several derivative of the same quantity (e.g. : depsxc)
248 : : !!
249 : : !! ..._b ----> means a block of the quantity "..." (use in mpi loops which treat the data block by block)
250 : : !! ..._updn -> means that spin up and spin down is available in that array
251 : : !! as (..,1) and (..,2). (if xcdata%nspden >=2 of course).
252 : : !! ..._apn --> in case of positrons are concerned.
253 : : !!
254 : : !! for more details about notations please see pdf in /doc/theory/MGGA/
255 : : !!
256 : : !! SOURCE
257 : :
258 : 113660 : subroutine rhotoxc(bigexc,bigsxc,kxc,mpi_enreg,nfft,ngfft, &
259 : 56830 : & nhat,nhatdim,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,option, &
260 : 113660 : & rhor,rprimd,usexcnhat,vxc,vxcavg,xccc3d,xcdata, &
261 [ + + + - : 112151 : & add_tfw,exc_vdw_out,grho1_over_rho1,electronpositron,k3xc,strsxc,taur,vhartr,vxctau,xc_funcs,xcctau3d) ! optional arguments
- + + + +
- - + + +
+ - - + ]
262 : :
263 : : !Arguments ------------------------------------
264 : : !scalars
265 : : integer,intent(in) :: nk3xc,n3xccc,nfft,nhatdim,nhatgrdim,nkxc,option
266 : : integer,intent(in) :: usexcnhat
267 : : logical,intent(in) :: non_magnetic_xc
268 : : logical,intent(in),optional :: add_tfw
269 : : real(dp),intent(out) :: bigexc,vxcavg,bigsxc
270 : : real(dp),intent(out),optional :: exc_vdw_out,grho1_over_rho1
271 : : type(MPI_type),intent(in) :: mpi_enreg
272 : : type(electronpositron_type),pointer,optional :: electronpositron
273 : : type(xcdata_type), intent(in) :: xcdata
274 : : !arrays
275 : : integer,intent(in) :: ngfft(18)
276 : : real(dp),intent(in) :: nhat(nfft,xcdata%nspden*nhatdim)
277 : : real(dp),intent(in) :: nhatgr(nfft,xcdata%nspden,3*nhatgrdim)
278 : : real(dp),intent(in),target :: rhor(nfft,xcdata%nspden)
279 : : real(dp),intent(in) :: rprimd(3,3),xccc3d(n3xccc)
280 : : real(dp),intent(in),optional :: xcctau3d(:)
281 : : real(dp),intent(out) :: kxc(nfft,nkxc),vxc(nfft,xcdata%nspden)
282 : : real(dp),intent(in),optional :: vhartr(nfft)
283 : : real(dp),intent(in),target,optional :: taur(:,:)
284 : : real(dp),intent(out),optional :: strsxc(6),k3xc(1:nfft,1:nk3xc),vxctau(:,:,:)
285 : : type(libxc_functional_type),intent(inout),optional :: xc_funcs(2)
286 : :
287 : : !Local variables-------------------------------
288 : : !scalars
289 : : integer :: auxc_ixc,cplex,ierr,ifft,ii,ixc,ixc_from_lib,indx,ipositron,ipts,ishift,ispden,iwarn,iwarnp
290 : : integer :: jj,mpts,ndvxc,nd2vxc,nfftot,ngr,ngrad,ngrad_apn,nkxc_eff,npts
291 : : integer :: nspden,nspden_apn,nspden_eff,nspden_updn,nspgrad,nvxcgrho,nvxclrho,nvxctau
292 : : integer :: n3xctau,order,nproc_fft,comm_fft,usegradient,usekden,uselaplacian
293 : : logical :: compute_stress,my_add_tfw
294 : : real(dp),parameter :: mot=-one/3.0_dp
295 : : real(dp) :: coeff,divshft,doti,dstrsxc,dvdn,dvdz,rhoexc,rhotsxc,factor,m_norm_min,s1,s2,s3
296 : : real(dp) :: strdiag,strsxc1_tot,strsxc2_tot,strsxc3_tot,strsxc4_tot
297 : : real(dp) :: strsxc5_tot,strsxc6_tot,ucvol
298 : : real(dp) :: deltae_vdw,exc_vdw
299 : : logical :: test_nhat,test_tb09,need_nhat,need_nhatgr,with_vxctau
300 : : character(len=500) :: message
301 : : real(dp) :: hyb_mixing, hyb_mixing_sr, hyb_range
302 : : !arrays
303 : : real(dp) :: d2rho(6),gm_norm(3),grho(3),gmet(3,3),gprimd(3,3),qphon(3),rmet(3,3)
304 : : real(dp) :: tsec(2),vxcmean(4)
305 : 56830 : real(dp),allocatable :: d2rhonow(:,:,:)
306 : 56830 : real(dp),allocatable :: d2vxc_b(:,:),depsxc(:,:),depsxc_apn(:,:),dvxc_apn(:),dvxc_b(:,:)
307 : 56830 : real(dp),allocatable :: exc_b(:),tsxc_b(:),fxc_apn(:),grho2_apn(:),grho2_b_updn(:,:)
308 : 56830 : real(dp),allocatable :: lrhonow(:,:),lrho_b_updn(:,:)
309 : 56830 : real(dp),allocatable :: m_norm(:),nhat_up(:),rho_b_updn(:,:),rho_b(:),rhonow_apn(:,:,:)
310 : 56830 : real(dp),allocatable :: tau_b_updn(:,:),vxc_apn(:,:),vxcgr_apn(:),vxcgrho_b_updn(:,:),vxcrho_b_updn(:,:)
311 : 56830 : real(dp),allocatable :: vxc_b_apn(:),vxc_ep(:),vxctau_b_updn(:,:),vxclrho_b_updn(:,:)
312 : 56830 : real(dp),allocatable,target :: rhonow(:,:,:),taunow(:,:,:)
313 : 56830 : real(dp),pointer :: rhocorval(:,:),rhor_(:,:),taucorval(:,:),taur_(:,:)
314 : 56830 : real(dp),contiguous, pointer :: rhonow_ptr(:,:,:)
315 : 113660 : real(dp) :: decdrho_vdw(nfft,xcdata%nspden),decdgrho_vdw(nfft,3,xcdata%nspden)
316 : : real(dp) :: strsxc_vdw(3,3)
317 [ + + ]: 170490 : type(libxc_functional_type) :: xc_funcs_auxc(2)
318 : : ! *************************************************************************
319 : :
320 : : ! Note: the following cases seem to never be tested (should be fixed)
321 : : ! - ipositron==2 and ngrad_apn==2
322 : : ! - usewvl/=0
323 : : ! - test_nhat and usexcnhat==1 and nspden==4
324 : :
325 : 56830 : call timab(81,1,tsec)
326 : :
327 : : !Optional arguments
328 : 56830 : compute_stress=present(strsxc)
329 [ + + ]: 56830 : my_add_tfw=.false.;if (present(add_tfw)) my_add_tfw=add_tfw
330 : :
331 : : !Useful scalars
332 : 56830 : nspden=xcdata%nspden
333 : 56830 : ixc=xcdata%ixc
334 : 56830 : auxc_ixc=xcdata%auxc_ixc
335 : 56830 : n3xctau=0
336 : :
337 : : !nspden_updn: 1 for non-polarized, 2 for polarized
338 : 56830 : nspden_updn=min(nspden,2)
339 : :
340 : : !The variable order indicates to which derivative of the energy
341 : : !the computation must be done. Computing exc and vxc needs order=1 .
342 : : !Meaningful values are 1, 2, 3. Lower than 1 is the same as 1, and larger
343 : : !than 3 is the same as 3.
344 : : !order=1 or 2 supported for all LSD and GGA ixc
345 : : !order=3 supported only for ixc=3 and ixc=7
346 : 56830 : order=1
347 [ + + ]: 56830 : if(option==2.or.option==10.or.option==12)order=2
348 [ + + ]: 56830 : if(option==-2)order=-2
349 [ + + ]: 56830 : if(option==3)order=3
350 : :
351 : : !Sizes of local arrays
352 [ + + ]: 56830 : if (present(xc_funcs)) then
353 : : call size_dvxc(ixc,order,nspden_updn,&
354 : : & usegradient=usegradient,uselaplacian=uselaplacian,usekden=usekden,&
355 : : & nvxcgrho=nvxcgrho,nvxclrho=nvxclrho,nvxctau=nvxctau,&
356 : 521 : & ndvxc=ndvxc,nd2vxc=nd2vxc,add_tfw=my_add_tfw,xc_funcs=xc_funcs)
357 : : else
358 : : call size_dvxc(ixc,order,nspden_updn,&
359 : : & usegradient=usegradient,uselaplacian=uselaplacian,usekden=usekden,&
360 : : & nvxcgrho=nvxcgrho,nvxclrho=nvxclrho,nvxctau=nvxctau,&
361 : 56309 : & ndvxc=ndvxc,nd2vxc=nd2vxc,add_tfw=my_add_tfw)
362 : : end if
363 : :
364 : : !ngrad=1 is for LDAs or LSDs, ngrad=2 is for GGAs/mGGAs
365 [ + + + + ]: 56830 : ngrad=1;if(xcdata%xclevel==2.or.usegradient==1) ngrad=2
366 : :
367 : : !nspden_eff: effective value of nspden used to compute gradients of density:
368 : : ! 1 for non-polarized system,
369 : : ! 2 for collinear polarized system or LDA (can be reduced to a collinear system)
370 : : ! 4 for non-collinear polarized system and GGA
371 [ + + ]: 56830 : nspden_eff=nspden_updn;if (nspden==4.and.ngrad==2) nspden_eff=4
372 : :
373 : : !Number of kcxc components depends on option (force LDA type if option==10 or 12)
374 [ + + ]: 56830 : nkxc_eff=nkxc;if (option==10.or.option==12) nkxc_eff=min(nkxc,3)
375 : :
376 : : !Check options
377 [ + + + + : 56830 : if(option==3.and.nd2vxc==0.and.ixc/=0)then
- + ]
378 : : write(message, '(3a,i0)' )&
379 : 0 : & 'Third-order xc kernel can only be computed for ixc = 0, 3, 7 or 8,',ch10,&
380 : 0 : & 'while it is found to be ',ixc
381 : 0 : ABI_ERROR(message)
382 : : end if
383 [ + + + + : 56830 : if(nspden==4.and.xcdata%xclevel==2.and.(abs(option)==2))then
- + ]
384 : 0 : ABI_BUG('When nspden==4 and GGA, the absolute value of option cannot be 2 !')
385 : : end if
386 [ + + ]: 56830 : if(ixc<0) then
387 [ + + ]: 13232 : if (present(xc_funcs)) then
388 : 521 : ixc_from_lib=libxc_functionals_ixc(xc_functionals=xc_funcs)
389 : : else
390 : 12711 : ixc_from_lib=libxc_functionals_ixc()
391 : : end if
392 : : ! Check consistency between ixc passed in input and the one used to initialize the library.
393 [ - + ]: 13232 : if (ixc /= ixc_from_lib) then
394 : : write(message, '(a,i0,2a,i0,2a)')&
395 : 0 : & 'The value of ixc specified in input, ixc = ',ixc,ch10,&
396 : 0 : & 'differs from the one used to initialize the functional ',ixc_from_lib,ch10,&
397 : 0 : & 'Action: reinitialize the global structure funcs, see NOTES in m_libxc_functionals'
398 : 0 : ABI_BUG(message)
399 : : end if
400 : : end if
401 : :
402 : : !Handling of mGGA functionals
403 : 56830 : with_vxctau=(present(vxctau))
404 [ + + + + ]: 221110 : if (with_vxctau) with_vxctau=(size(vxctau)>0)
405 [ + + ]: 56830 : if (usekden==1) then
406 [ - + ]: 666 : if (.not.present(taur)) then
407 : 0 : message=' For mGGA functionals, kinetic energy density is needed. Set input variable usekden to 1.'
408 [ # # ]: 0 : message=trim(message)//' Also use NC pseudopotentials without non-linear XC core correction.'
409 : 0 : ABI_BUG(message)
410 [ + + - + ]: 1998 : else if (size(taur)/=nfft*nspden) then
411 : 0 : ABI_BUG('Invalid size for taur!')
412 : : end if
413 [ - + ]: 666 : if (present(xcctau3d)) then
414 : 666 : n3xctau=size(xcctau3d)
415 [ + + + - ]: 666 : if (n3xctau/=0.and.n3xctau/=nfft) then
416 : 0 : ABI_BUG('Invalid size for xccctau3d!')
417 : : end if
418 : : end if
419 [ + - ]: 666 : if (with_vxctau) then
420 [ + + - + ]: 2664 : if (size(vxctau)/=nfft*nspden*4) then
421 : 0 : ABI_BUG('Invalid size for vxctau!')
422 : : end if
423 : : end if
424 : : end if
425 [ + + + + : 56830 : if((usekden==1.or.uselaplacian==1).and.nspden==4)then
- + ]
426 : : !mGGA en NC-magnetism: how do we rotate tau kinetic energy density?
427 : 0 : message=' At present, meta-GGA (usekden=1 or uselaplacian=1) is not compatible with non-collinear magnetism (nspden=4).'
428 : 0 : ABI_ERROR(message)
429 : : end if
430 : :
431 : : !MPI FFT communicator
432 : 56830 : comm_fft = mpi_enreg%comm_fft; nproc_fft = mpi_enreg%nproc_fft
433 : :
434 : : !Compute different geometric tensor, as well as ucvol, from rprimd
435 : 56830 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
436 : :
437 : : !In this routine, hartre, xcden and xcpot are called for real
438 : : !densities and potentials, corresponding to zero wavevector
439 : 56830 : cplex=1
440 : 56830 : qphon(:)=zero
441 : 56830 : iwarn=0
442 : 56830 : nfftot=ngfft(1)*ngfft(2)*ngfft(3)
443 : :
444 : : !Initializations
445 : 56830 : bigexc=zero
446 : 56830 : bigsxc=zero
447 : 56830 : rhoexc=zero
448 : 56830 : rhotsxc=zero
449 [ + + + + ]: 1617965775 : vxc(:,:)=zero
450 : 56830 : vxcavg=zero
451 [ + + ]: 56830 : if (compute_stress) then
452 : 55711 : strsxc(:)=zero
453 : 55711 : strsxc1_tot=zero;strsxc2_tot=zero;strsxc3_tot=zero
454 : 55711 : strsxc4_tot=zero;strsxc5_tot=zero;strsxc6_tot=zero
455 : : strsxc_vdw(:,:)=zero
456 : : end if
457 [ + + + + : 151205330 : if (with_vxctau) vxctau(:,:,:)=zero
+ + + + ]
458 [ + + + + : 760769291 : if (nkxc/=0) kxc(:,:)=zero
+ + ]
459 [ + + + - : 371551 : if(abs(option)==3.and.nk3xc/=0) k3xc(:,:)=zero
+ + + + ]
460 : 56830 : ipositron=0
461 [ + + ]: 56830 : if (present(electronpositron)) then
462 : 288 : ipositron=electronpositron_calctype(electronpositron)
463 [ + - ]: 288 : if (ipositron==2) then
464 : 288 : electronpositron%e_xc =zero
465 : 288 : electronpositron%e_xcdc=zero
466 : : end if
467 : : end if
468 : 56830 : deltae_vdw = zero
469 : 56830 : exc_vdw = zero
470 [ + + + + ]: 1617965775 : decdrho_vdw(:,:) = zero
471 [ + + + + : 4853854184 : decdgrho_vdw(:,:,:) = zero
+ + ]
472 [ + + ]: 56830 : if (present(grho1_over_rho1)) grho1_over_rho1=zero
473 : :
474 [ + + - + : 56830 : if ((xcdata%xclevel==0.or.ixc==0).and.(.not.my_add_tfw)) then
+ - ]
475 : : ! No xc at all is applied (usually for testing)
476 : 735 : ABI_WARNING('Note that no xc is applied (ixc=0).')
477 : :
478 [ + + ]: 56095 : else if (ixc/=20) then
479 : :
480 : : ! Test: has a compensation density to be added/substracted (PAW) ?
481 [ + + + + ]: 56053 : need_nhat=(nhatdim==1.and.usexcnhat==0)
482 [ + + + + : 56053 : need_nhatgr=(nhatdim==1.and.nhatgrdim==1.and.ngrad==2.and.xcdata%intxc==0)
+ - + + ]
483 : 56053 : test_nhat=(need_nhat.or.need_nhatgr)
484 : :
485 : : ! The different components of depsxc will be
486 : : ! for nspden=1, depsxc(:,1)=d(rho.exc)/d(rho) == (depsxcdrho) == (vxcrho)
487 : : ! and if ngrad=2, depsxc(:,2)=1/2*1/|grad rho_up|*d(rho.exc)/d(|grad rho_up|)
488 : : ! +1/|grad rho|*d(rho.exc)/d(|grad rho|)
489 : : ! == (1/2 * 1/|grho_up| * depsxcd|grho_up|) + 1/|grho| * depsxcd|grho|
490 : : ! (vxcgrho=1/|grho| * depsxcd|grho|)
491 : : ! (do not forget : |grad rho| /= |grad rho_up| + |grad rho_down|
492 : : ! and if use_laplacian, depsxc(:,3)=d(rho.exc)/d(lapl rho) == (depsxcdlrho) == (vxclrho)
493 : : !
494 : : ! for nspden>=2, depsxc(:,1)=d(rho.exc)/d(rho_up) == (depsxcdrho_up) == (vxcrho_up)
495 : : ! depsxc(:,2)=d(rho.exc)/d(rho_down) == (depsxcdrho_dn) == (vxcrho_dn)
496 : : ! and if ngrad=2, depsxc(:,3)=1/|grad rho_up|*d(rho.exc)/d(|grad rho_up|) == (1/|grho_up| * depsxcd|grho_up|) == (vxcgrho_up)
497 : : ! depsxc(:,4)=1/|grad rho_down|*d(rho.exc)/d(|grad rho_down|) == (1/|grho_dn| * depsxcd|grho_dn|) == (vxcgrho_dn)
498 : : ! depsxc(:,5)=1/|grad rho|*d(rho.exc)/d(|grad rho|) == (1/|grho| * depsxcd|grho|) == (vxcgrho)
499 : : ! and if use_laplacian, depsxc(:,6)=d(rho.exc)/d(lapl rho_up) == (depsxcdlrho_up) == (vxclrho_up)
500 : : ! depsxc(:,7)=d(rho.exc)/d(lapl rho_dn) == (depsxcdlrho_dn) == (vxclrho_dn)
501 : : ! Note: if nspden=4, rho_up=(rho+|m|)/2, rho_down=(rho-|m|)/2
502 [ + + + + ]: 56053 : nspgrad=nspden_updn*ngrad;if(nspden_updn==2.and.ngrad==2)nspgrad=5
503 [ + + ]: 56053 : if(uselaplacian==1) nspgrad=nspgrad+nspden_updn
504 [ - + + - : 224212 : ABI_MALLOC(depsxc,(nfft,nspgrad))
- + + - +
- + - - +
- + ]
505 [ + + + + ]: 2056458306 : depsxc(:,:)=zero
506 : :
507 : : ! PAW: select the valence density (and magnetization) to use:
508 : : ! link the correct density, according to usexcnhat option
509 [ + + + + ]: 56053 : if ((.not.need_nhat).and.(.not.non_magnetic_xc)) then
510 : 44890 : rhor_ => rhor
511 : : else
512 [ - + + - : 44652 : ABI_MALLOC(rhor_,(nfft,nspden))
- + + - +
- + - - +
- + ]
513 [ + + ]: 11163 : if (need_nhat) then
514 [ + + ]: 23754 : do ispden=1,nspden
515 [ + + ]: 438005024 : do ifft=1,nfft
516 : 437994148 : rhor_(ifft,ispden)=rhor(ifft,ispden)-nhat(ifft,ispden)
517 : : end do
518 : : end do
519 : : else
520 [ + + ]: 1241 : do ispden=1,nspden
521 [ + + ]: 13531913 : do ifft=1,nfft
522 : 13531626 : rhor_(ifft,ispden)=rhor(ifft,ispden)
523 : : end do
524 : : end do
525 : : end if
526 [ + + ]: 11163 : if(non_magnetic_xc) then
527 [ + + + + ]: 3554560 : if(nspden==2) rhor_(:,2)=rhor_(:,1)*half
528 [ + + + + : 12484498 : if(nspden==4) rhor_(:,2:4)=zero
+ + ]
529 : : endif
530 : : end if
531 [ + + ]: 56053 : if (usekden==1) then
532 [ - + ]: 666 : if(non_magnetic_xc) then
533 [ # # # # : 0 : ABI_MALLOC(taur_,(nfft,nspden))
# # # # #
# # # # #
# # ]
534 [ # # # # ]: 0 : if(nspden==2) taur_(:,2)=taur_(:,1)*half
535 [ # # # # : 0 : if(nspden==4) taur_(:,2:4)=zero
# # ]
536 : : else
537 : 666 : taur_ => taur
538 : : end if
539 : : end if
540 : :
541 : : ! Some initializations for the electron-positron correlation
542 [ + + ]: 56053 : if (ipositron==2) then
543 : 288 : nspden_apn=1;ngrad_apn=1;iwarnp=1
544 [ - + ]: 288 : if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad_apn=2
545 [ - + - - ]: 288 : if (ngrad_apn==2.and.xcdata%xclevel<2) then
546 : 0 : message = 'GGA for the positron can only be performed with GGA pseudopotentials for the electron !'
547 : 0 : ABI_ERROR(message)
548 : : end if
549 [ - + - - ]: 288 : if (ngrad_apn>1.and.option/=0.and.option/=1.and.option/=10.and.option/=12) then
550 : 0 : message = 'You cannot compute full GGA XC kernel for electrons-positron systems !'
551 : 0 : ABI_ERROR(message)
552 : : end if
553 [ - + + - : 1152 : ABI_MALLOC(depsxc_apn,(nfft,ngrad_apn))
- + + - +
- + - - +
- + ]
554 : : end if
555 : :
556 : : ! Non-collinear magnetism: store norm of magnetization
557 : : ! m_norm_min= EPSILON(0.0_dp)**2 ! EB: TOO SMALL!!!
558 : 56053 : m_norm_min=tol8 ! EB: tol14 is still too small, tests are underway
559 [ + + ]: 56053 : if (nspden==4) then
560 [ - + + - : 7749 : ABI_MALLOC(m_norm,(nfft))
+ - + - -
+ - + ]
561 [ + + ]: 32013089 : m_norm(:)=sqrt(rhor_(:,2)**2+rhor_(:,3)**2+rhor_(:,4)**2)
562 : : end if
563 : :
564 : : ! rhocorval will contain effective density used to compute gradients:
565 : : ! - with core density (if NLCC)
566 : : ! - without compensation density (if PAW under certain conditions)
567 : : ! - in (up+dn,up) or (n,mx,my,mz) format according to collinearity
568 : : ! of polarization and use of gradients (GGA)
569 [ + + + + ]: 56053 : if (n3xccc>0.or.test_nhat.or.nspden_eff/=nspden) then
570 [ - + + - : 126908 : ABI_MALLOC(rhocorval,(nfft,nspden_eff))
- + + - +
- + - - +
- + ]
571 [ + + ]: 31727 : if (nspden==nspden_eff) then
572 [ + - + + : 1757424894 : rhocorval(:,1:nspden)=rhor_(:,1:nspden)
+ + + + +
+ ]
573 [ + - ]: 2003 : else if (nspden==4) then
574 [ + - + + : 50894017 : rhocorval(:,1)=rhor_(:,1)
+ + ]
575 [ + + + + ]: 50894017 : rhocorval(:,2)=half*(rhor_(:,1)+m_norm(:))
576 : : else
577 [ # # # # ]: 0 : rhocorval=zero
578 : : end if
579 : : else
580 : 24326 : rhocorval => rhor_
581 : : end if
582 [ + + + + ]: 56053 : if (usekden==1.and.(n3xctau>0.or.nspden_eff/=nspden)) then
583 [ - + + - : 2124 : ABI_MALLOC(taucorval,(nfft,nspden_eff))
- + + - +
- + - - +
- + ]
584 [ + - ]: 531 : if (nspden==nspden_eff) then
585 [ + - + + : 57817401 : taucorval(:,1:nspden)=taur_(:,1:nspden)
+ + + + +
+ ]
586 : : else
587 [ # # # # ]: 0 : taucorval=zero
588 : : end if
589 : : else
590 : 55522 : taucorval => taur_
591 : : end if
592 : :
593 : : ! Add core electron density to effective density
594 [ + + ]: 56053 : if (n3xccc>0) then
595 [ + + ]: 632221693 : rhocorval(:,1)=rhocorval(:,1)+xccc3d(:)
596 [ + + ]: 26646 : if(nspden_eff==2) then
597 [ + + ]: 87540113 : rhocorval(:,2)=rhocorval(:,2)+half*xccc3d(:)
598 : : end if
599 : : end if
600 [ + + ]: 56053 : if (n3xctau>0) then
601 [ + + ]: 28908435 : taucorval(:,1)=taucorval(:,1)+xcctau3d(:)
602 [ - + ]: 531 : if(nspden_eff==2) then
603 [ # # ]: 0 : taucorval(:,2)=taucorval(:,2)+half*xcctau3d(:)
604 : : end if
605 : : end if
606 : :
607 : : ! If PAW, substract compensation density from effective density:
608 : : ! - if GGA, because nhat gradients are computed separately
609 [ + + + + ]: 56053 : if (test_nhat.and.usexcnhat==1) then
610 [ + - ]: 730 : if (nspden==nspden_eff) then
611 [ + + + + ]: 28633169 : rhocorval(:,1:nspden)=rhocorval(:,1:nspden)-nhat(:,1:nspden)
612 [ # # ]: 0 : else if (nspden==4) then
613 [ # # # # : 0 : ABI_MALLOC(nhat_up,(nfft))
# # # # #
# # # ]
614 [ # # ]: 0 : do ifft=1,nfft
615 [ # # ]: 0 : if (m_norm(ifft)>m_norm_min) then
616 : : nhat_up(ifft)=half*(nhat(ifft,1) &
617 : : & +(rhor_(ifft,2)*nhat(ifft,2) &
618 : : & +rhor_(ifft,3)*nhat(ifft,3) &
619 : 0 : & +rhor_(ifft,4)*nhat(ifft,4))/m_norm(ifft))
620 : : else
621 : : nhat_up(ifft)=half*(nhat(ifft,1) &
622 : 0 : & +sqrt(nhat(ifft,2)**2+nhat(ifft,3)**2+nhat(ifft,4)**2))
623 : : end if
624 : : end do
625 [ # # ]: 0 : rhocorval(:,1)=rhocorval(:,1)-nhat(:,1)
626 [ # # ]: 0 : rhocorval(:,2)=rhocorval(:,2)-nhat_up(:)
627 : : end if
628 : : end if
629 : :
630 : : ! rhonow will contain effective density (and gradients if GGA)
631 : : ! taunow will contain effective kinetic energy density (if MGGA)
632 : : ! lrhonow will contain the laplacian if we have a MGGA
633 : : ! d2rhonow will contain the 2nd derivatives if we have a MGGA and need the stress tensor
634 [ - + + - : 280265 : ABI_MALLOC(rhonow,(nfft,nspden_eff,ngrad*ngrad))
- + + - -
+ + - + -
+ - - + -
+ ]
635 [ - + + - : 224212 : ABI_MALLOC(lrhonow,(nfft,nspden_eff*uselaplacian))
+ + + - +
- + + - +
- + ]
636 [ - + + - : 280265 : ABI_MALLOC(taunow,(nfft,nspden_eff,usekden))
- + + - +
+ + - + -
+ + - + -
+ ]
637 [ + + + + ]: 56053 : if (compute_stress.and.uselaplacian>0) then
638 [ - + + - : 675 : ABI_MALLOC(d2rhonow,(nfft,nspden_eff,6*uselaplacian))
- + + - -
+ + - + -
+ - - + -
+ ]
639 : : end if
640 : :
641 : : ! ====================================================================
642 : : ! ====================================================================
643 : : ! Loop on unshifted or shifted grids
644 [ + + ]: 116309 : do ishift=0,xcdata%intxc
645 : :
646 : : ! Set up density on unshifted or shifted grid (will be in rhonow(:,:,1)),
647 : : ! as well as the gradient of the density, also on the unshifted
648 : : ! or shifted grid (will be in rhonow(:,:,2:4)), if needed.
649 [ + + ]: 60256 : if (uselaplacian==1) then
650 [ + - ]: 135 : if (compute_stress) then
651 : : call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_eff,&
652 [ - + ]: 135 : & qphon,rhocorval,rhonow,lrhonow=lrhonow,d2rhonow=d2rhonow)
653 : : else
654 : : call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_eff,&
655 [ # # ]: 0 : & qphon,rhocorval,rhonow,lrhonow=lrhonow)
656 : : end if
657 : : else
658 : : call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_eff,&
659 [ - + ]: 60121 : & qphon,rhocorval,rhonow)
660 : : end if
661 [ + + ]: 60256 : if (usekden==1) then
662 : : call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,1,nspden_eff,&
663 [ - + ]: 666 : & qphon,taucorval,taunow)
664 : : end if
665 : :
666 : : ! PAW+GGA: add "exact" gradients of compensation density
667 : : !if (test_nhat.and.usexcnhat==1.and.ishift==0) then
668 [ + + + + ]: 60256 : if (test_nhat.and.usexcnhat==1) then
669 [ + - ]: 730 : if (nspden==nspden_eff) then
670 [ + - + + : 57265608 : rhonow(:,1:nspden,1)=rhocorval(:,1:nspden)+nhat(:,1:nspden)
+ + + + +
+ ]
671 [ # # ]: 0 : else if (nspden==4) then
672 [ # # # # : 0 : rhonow(:,1,1)=rhocorval(:,1)+nhat(:,1)
# # ]
673 [ # # # # : 0 : rhonow(:,2,1)=rhocorval(:,2)+nhat_up(:)
# # ]
674 : : end if
675 [ + - + - : 730 : if (ngrad==2.and.nhatgrdim==1.and.nspden==nspden_eff) then
+ - ]
676 [ + + ]: 2920 : do ii=1,3
677 : 2190 : jj=ii+1
678 [ + + ]: 6982 : do ispden=1,nspden
679 [ + + ]: 85899507 : do ifft=1,nfft
680 : 85897317 : rhonow(ifft,ispden,jj)=rhonow(ifft,ispden,jj)+nhatgr(ifft,ispden,ii)
681 : : end do
682 : : end do
683 : : end do
684 : : end if
685 : : end if
686 : :
687 : : ! Deallocate temporary arrays
688 [ + + ]: 60256 : if (ishift==xcdata%intxc) then
689 [ + + + + : 56053 : if (n3xccc>0.or.test_nhat.or.nspden_eff/=nspden) then
+ + ]
690 [ - + ]: 31727 : ABI_FREE(rhocorval)
691 : : end if
692 [ + + + + : 56053 : if (usekden==1.and.(n3xccc>0.or.nspden_eff/=nspden)) then
- + ]
693 [ - + ]: 531 : ABI_FREE(taucorval)
694 : : end if
695 [ + + + - : 56053 : if (test_nhat.and.nspden/=nspden_eff.and.usexcnhat==1) then
- - ]
696 [ # # ]: 0 : ABI_FREE(nhat_up)
697 : : end if
698 : : end if
699 : :
700 : : ! In case of non-collinear magnetism, extract up and down density and gradients (if GGA)
701 [ + + + + ]: 60256 : if (nspden==4.and.nspden_eff==nspden) then
702 [ + - ]: 580 : if (ngrad==2) then
703 [ + + ]: 6565079 : do ifft=1,nfft
704 : 6564499 : gm_norm(1:3)=zero
705 [ + + ]: 6564499 : if(m_norm(ifft)>m_norm_min) then
706 : : ! if(m_norm(ifft)>rhonow(ifft,1,1)*tol10+tol14) then
707 [ + + ]: 20443468 : do jj=1,3 ! Compute here nabla(|m|)=(m.nabla(m))/|m| == (g|m| = m/|m| * gm)
708 [ + + ]: 66441271 : do ii=2,4
709 : 61330404 : gm_norm(jj)=gm_norm(jj)+rhonow(ifft,ii,1+jj)*rhonow(ifft,ii,1)
710 : : end do
711 : : end do
712 [ + + ]: 20443468 : gm_norm(1:3)=gm_norm(1:3)/m_norm(ifft)
713 : : end if
714 : 6564499 : rhonow(ifft,2,2)=half*(rhonow(ifft,1,2)+gm_norm(1))
715 : 6564499 : rhonow(ifft,2,3)=half*(rhonow(ifft,1,3)+gm_norm(2))
716 : 6565079 : rhonow(ifft,2,4)=half*(rhonow(ifft,1,4)+gm_norm(3))
717 : : end do
718 : : end if
719 [ + + ]: 6565079 : rhonow(:,2,1)=half*(rhonow(:,1,1)+m_norm(:))
720 [ - + - - ]: 580 : if (usekden==1) taunow(:,2,1)=half*(taunow(:,1,1)+m_norm(:))
721 : : end if
722 : : ! Make the density positive everywhere (but do not care about gradients)
723 : 60256 : call mkdenpos(iwarn,nfft,nspden_updn,1,rhonow(:,1:nspden_updn,1),xcdata%xc_denpos)
724 [ + + ]: 60256 : if (usekden==1) then
725 : 666 : call mkdenpos(iwarn,nfft,nspden_updn,1,taunow(:,1:nspden_updn,1),xcdata%xc_taupos)
726 : : end if
727 : :
728 : : ! write(std_out,*) 'rhonow',rhonow
729 : :
730 : : ! Uses a block formulation, in order to save simultaneously
731 : : ! CPU time and memory : xc routines
732 : : ! are called only once over mpts times, while the amount of allocated
733 : : ! space is kept at a low value, even if a lot of different
734 : : ! arrays are allocated, for use in different xc functionals.
735 : :
736 : 60256 : mpts=4000
737 [ + + ]: 60256 : if (usekden==1) mpts=nfft ! Why?
738 : :
739 [ - + - - : 60256 : do ifft=1,nfft,mpts
- + ]
740 : : ! npts=mpts
741 : : ! npts is the number of points to be treated in this bunch
742 : 355339 : npts=min(nfft-ifft+1,mpts)
743 : :
744 : : ! Allocation of mandatory arguments of drivexc
745 [ - + + - : 1066017 : ABI_MALLOC(exc_b,(npts))
+ - + - -
+ - + ]
746 [ - + + - : 710678 : ABI_MALLOC(tsxc_b,(npts))
- + - + ]
747 [ - + + - : 710678 : ABI_MALLOC(rho_b,(npts))
- + - + ]
748 [ - + + - : 1421356 : ABI_MALLOC(rho_b_updn,(npts,nspden_updn))
- + + - +
- + - - +
- + ]
749 [ - + + - : 1066017 : ABI_MALLOC(vxcrho_b_updn,(npts,nspden_updn))
- + + - -
+ - + ]
750 [ + + + + ]: 1548163597 : vxcrho_b_updn(:,:)=zero
751 : :
752 : : ! Allocation of optional arguments of drivexc
753 [ - + + - : 1421356 : ABI_MALLOC(grho2_b_updn,(npts,(2*nspden_updn-1)*usegradient))
+ + + - +
- + + - +
- + ]
754 [ - + + - : 1421356 : ABI_MALLOC(lrho_b_updn,(npts,nspden_updn*uselaplacian))
+ + + - +
- + + - +
- + ]
755 [ - + + - : 1421356 : ABI_MALLOC(tau_b_updn,(npts,nspden_updn*usekden))
+ + + - +
- + + - +
- + ]
756 [ - + + - : 1421356 : ABI_MALLOC(vxcgrho_b_updn,(npts,nvxcgrho))
+ + + - +
- + + - +
- + ]
757 [ - + + - : 1421356 : ABI_MALLOC(vxclrho_b_updn,(npts,nvxclrho))
+ + + - +
- + + - +
- + ]
758 [ - + + - : 1421356 : ABI_MALLOC(vxctau_b_updn,(npts,nvxctau))
+ + + - +
- + + - +
- + ]
759 [ - + + - : 1421356 : ABI_MALLOC(dvxc_b,(npts,ndvxc))
+ + + - +
- + + - +
- + ]
760 [ - + + - : 1421356 : ABI_MALLOC(d2vxc_b,(npts,nd2vxc))
+ + + - +
- + + - +
- + ]
761 [ + + + + : 1264794950 : if (nvxcgrho>0) vxcgrho_b_updn(:,:)=zero
+ + ]
762 [ + + + + : 7018412 : if (nvxclrho>0) vxclrho_b_updn(:,:)=zero
+ + ]
763 [ + + + + : 38141798 : if (nvxctau>0) vxctau_b_updn(:,:)=zero
+ + ]
764 : :
765 [ + + ]: 1343060029 : do ipts=ifft,ifft+npts-1
766 : : ! indx=ipts-ifft+1 varies from 1 to npts
767 : 1342704690 : indx=ipts-ifft+1
768 : 1342704690 : rho_b(indx)=rhonow(ipts,1,1)
769 [ + + ]: 1343060029 : if(nspden_updn==1)then
770 : 1138010766 : rho_b_updn(indx,1)=rhonow(ipts,1,1)*half
771 [ + + ]: 1138010766 : if (usegradient==1) grho2_b_updn(indx,1)=quarter*(rhonow(ipts,1,2)**2 &
772 : 355805484 : & +rhonow(ipts,1,3)**2+rhonow(ipts,1,4)**2)
773 [ + + ]: 1138010766 : if (usekden==1) tau_b_updn(indx,1)=taunow(ipts,1,1)*half
774 [ + + ]: 1138010766 : if (uselaplacian==1) lrho_b_updn(indx,1)=lrhonow(ipts,1)*half
775 : : else
776 : 204693924 : rho_b_updn(indx,1)=rhonow(ipts,2,1)
777 : 204693924 : rho_b_updn(indx,2)=rhonow(ipts,1,1)-rhonow(ipts,2,1)
778 [ + + ]: 204693924 : if(usegradient==1)then
779 : : grho2_b_updn(indx,1)=rhonow(ipts,2,2)**2+ &
780 : : & rhonow(ipts,2,3)**2+ &
781 : 67171506 : & rhonow(ipts,2,4)**2
782 : : grho2_b_updn(indx,2)=(rhonow(ipts,1,2)-rhonow(ipts,2,2))**2 + &
783 : : & (rhonow(ipts,1,3)-rhonow(ipts,2,3))**2 + &
784 : 67171506 : & (rhonow(ipts,1,4)-rhonow(ipts,2,4))**2
785 : : grho2_b_updn(indx,3)=rhonow(ipts,1,2)**2+ &
786 : : & rhonow(ipts,1,3)**2+ &
787 : 67171506 : & rhonow(ipts,1,4)**2
788 : : end if
789 [ + + ]: 204693924 : if (usekden==1) then
790 : 1322496 : tau_b_updn(indx,1)=taunow(ipts,2,1)
791 : 1322496 : tau_b_updn(indx,2)=taunow(ipts,1,1)-taunow(ipts,2,1)
792 : : end if
793 [ + + ]: 204693924 : if (uselaplacian==1) then
794 : 419904 : lrho_b_updn(indx,1)=lrhonow(ipts,2)
795 : 419904 : lrho_b_updn(indx,2)=lrhonow(ipts,1)-lrhonow(ipts,2)
796 : : end if
797 : : end if
798 : : end do
799 : : ! In case of a hybrid functional, if one needs to compute the auxiliary GGA Kxc,
800 : : ! a separate call to drivexc is first needed to compute Kxc using such auxiliary GGA,
801 : : ! before calling again drivexc using the correct functional for Exc and Vxc.
802 : :
803 [ + + + - ]: 355339 : if(xcdata%usefock==1 .and. auxc_ixc/=0)then
804 [ - + ]: 5851 : if (auxc_ixc<0) then
805 : 0 : call libxc_functionals_init(auxc_ixc,nspden,xc_functionals=xc_funcs_auxc)
806 : : end if
807 : : call drivexc(auxc_ixc,order,npts,nspden_updn,usegradient,0,0,&
808 : : & rho_b_updn,exc_b,tsxc_b,vxcrho_b_updn,nvxcgrho,0,0,ndvxc,nd2vxc,xcdata%tphysel, &
809 : : & grho2_updn=grho2_b_updn,vxcgrho=vxcgrho_b_updn,dvxc=dvxc_b, &
810 : 5851 : & hyb_mixing=xcdata%hyb_mixing,xc_funcs=xc_funcs_auxc)
811 : : ! Transfer the xc kernel
812 [ + + + - ]: 5851 : if (nkxc_eff==1.and.ndvxc==15) then
813 [ + + ]: 7434806 : kxc(ifft:ifft+npts-1,1)=half*(dvxc_b(1:npts,1)+dvxc_b(1:npts,9)+dvxc_b(1:npts,10))
814 [ + + + - ]: 3917 : else if (nkxc_eff==3.and.ndvxc==15) then
815 [ + + ]: 7984594 : kxc(ifft:ifft+npts-1,1)=dvxc_b(1:npts,1)+dvxc_b(1:npts,9)
816 [ + + ]: 7984594 : kxc(ifft:ifft+npts-1,2)=dvxc_b(1:npts,10)
817 [ + + ]: 7984594 : kxc(ifft:ifft+npts-1,3)=dvxc_b(1:npts,2)+dvxc_b(1:npts,11)
818 : : end if
819 [ - + ]: 5851 : if (auxc_ixc<0) then
820 : 0 : call libxc_functionals_end(xc_functionals=xc_funcs_auxc)
821 : : end if
822 : : end if
823 [ + + ]: 355339 : if (present(xc_funcs)) then
824 : : call libxc_functionals_get_hybridparams(hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,&
825 : 1024 : & hyb_range=hyb_range,xc_functionals=xc_funcs)
826 : : else
827 : : call libxc_functionals_get_hybridparams(hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,&
828 : 354315 : & hyb_range=hyb_range)
829 : : end if
830 : :
831 : : ! Call to main XC driver
832 [ + + ]: 355339 : if (present(xc_funcs)) then
833 : : call drivexc(ixc,order,npts,nspden_updn,&
834 : : & usegradient,uselaplacian,usekden,&
835 : : & rho_b_updn,exc_b,tsxc_b,vxcrho_b_updn,&
836 : : & nvxcgrho,nvxclrho,nvxctau,ndvxc,nd2vxc,xcdata%tphysel, &
837 : : & grho2_updn=grho2_b_updn,vxcgrho=vxcgrho_b_updn,&
838 : : & lrho_updn=lrho_b_updn,vxclrho=vxclrho_b_updn,&
839 : : & tau_updn=tau_b_updn,vxctau=vxctau_b_updn,&
840 : : & dvxc=dvxc_b,d2vxc=d2vxc_b,&
841 : : & hyb_mixing=xcdata%hyb_mixing,&
842 : 1024 : & xc_funcs=xc_funcs)
843 : : else
844 : : call drivexc(ixc,order,npts,nspden_updn,&
845 : : & usegradient,uselaplacian,usekden,&
846 : : & rho_b_updn,exc_b,tsxc_b,vxcrho_b_updn,&
847 : : & nvxcgrho,nvxclrho,nvxctau,ndvxc,nd2vxc,xcdata%tphysel, &
848 : : & grho2_updn=grho2_b_updn,vxcgrho=vxcgrho_b_updn,&
849 : : & lrho_updn=lrho_b_updn,vxclrho=vxclrho_b_updn,&
850 : : & tau_updn=tau_b_updn,vxctau=vxctau_b_updn,&
851 : : & dvxc=dvxc_b,d2vxc=d2vxc_b,&
852 : 354315 : & hyb_mixing=xcdata%hyb_mixing)
853 : : end if
854 : :
855 : : ! If fake meta-GGA, has to remove the core contribution
856 : : ! when electronic effective mass has been modified
857 [ + + + - : 355339 : if (n3xccc>0.and.(ixc==31.or.ixc==34.or.ixc==35)) then
+ - ]
858 [ # # ]: 0 : if (ixc==31.or.ixc==35) then
859 : 0 : coeff=one-(one/1.01_dp)
860 [ # # ]: 0 : if (nspden_updn==1) then
861 : 0 : coeff=half*coeff
862 [ # # ]: 0 : do ipts=1,npts
863 : : exc_b(ipts)=exc_b(ipts)-coeff*xcctau3d(ifft+ipts-1) &
864 : 0 : & /rho_b_updn(ipts,1)
865 : : end do
866 : : else
867 [ # # ]: 0 : do ipts=1,npts
868 : : exc_b(ipts)=exc_b(ipts)-coeff*xcctau3d(ifft+ipts-1) &
869 : 0 : & /(rho_b_updn(ipts,1)+rho_b_updn(ipts,2))
870 : : end do
871 : : end if
872 : : else
873 : 0 : message = 'MetaGGA ixc=34 is not yet allowed with a core kinetic energy density!'
874 : 0 : ABI_ERROR(message)
875 : : end if
876 : : end if
877 : :
878 : : ! Gradient Weiszacker correction to a Thomas-Fermi functional
879 [ + + ]: 355339 : if (my_add_tfw) then
880 [ + + + + ]: 2187756 : vxcgrho_b_updn(:,:)=zero
881 : : call xctfw(xcdata%tphysel,exc_b,tsxc_b,rho_b_updn,vxcrho_b_updn,npts,nspden_updn, &
882 : 189 : & vxcgrho_b_updn,nvxcgrho,grho2_b_updn)
883 : : end if
884 : :
885 : : ! Accumulate bigexc, bigsxc, strsxc and store vxc (and eventually kxc)
886 : 355339 : dstrsxc=zero
887 [ + + ]: 1343060029 : do ipts=ifft,ifft+npts-1
888 : 1342704690 : indx=ipts-ifft+1
889 : 1342704690 : rhoexc=rhoexc+rho_b(indx)*exc_b(indx) ! Will be normalized with respect to the volume later to get bigexc.
890 : 1342704690 : rhotsxc=rhotsxc+rho_b(indx)*tsxc_b(indx) ! Will be normalized with respect to the volume later to get bigsxc.
891 : 1342704690 : depsxc(ipts,1)=vxcrho_b_updn(indx,1)
892 [ + + ]: 1342704690 : if (nspden_updn==2) depsxc(ipts,2)=vxcrho_b_updn(indx,2)
893 [ + + ]: 1342704690 : if (compute_stress) then
894 : : ! Note for GGA/mGGA: the diagonal stress contribution is not complete because
895 : : ! this is not the complete Vxc. Only Int[n.dFxc/dn] is computed here. The
896 : : ! other terms are computed later
897 : : ! (i.e. -Int[n.grad(dFxc/grad(n))] and Int[n.Lapl(dFxc/Lapl(n))])
898 [ + + ]: 1324859266 : if(nspden_updn==1)then
899 : 1122177015 : strdiag=rho_b(indx)*(exc_b(indx)-vxcrho_b_updn(indx,1))
900 [ + + ]: 1122177015 : if (usekden==1) strdiag=strdiag-two*tau_b_updn(indx,1)*vxctau_b_updn(indx,1)
901 [ + - ]: 202682251 : else if(nspden_updn==2)then
902 : : strdiag=rho_b(indx)*exc_b(indx) &
903 : : & - rho_b_updn(indx,1)*vxcrho_b_updn(indx,1)&
904 : 202682251 : & - (rho_b(indx)-rho_b_updn(indx,1))*vxcrho_b_updn(indx,2)
905 [ + + ]: 202682251 : if (usekden==1) then
906 : : strdiag=strdiag - tau_b_updn(indx,1)*vxctau_b_updn(indx, 1) &
907 : 1322496 : & - tau_b_updn(indx,2)*vxctau_b_updn(indx,2)
908 : : end if
909 : : end if
910 : 1324859266 : dstrsxc=dstrsxc+strdiag
911 : : end if
912 : :
913 : : ! For GGAs, additional terms appear
914 : : ! (the LB functional does not lead to additional terms)
915 [ + + + - ]: 1342704690 : if(ngrad==2 .and. ixc/=13)then
916 : :
917 : : ! Treat explicitely spin up, spin down and total spin for spin-polarized
918 : : ! Will exit when ispden=1 is finished if non-spin-polarized
919 [ + + ]: 980296992 : do ispden=1,3
920 : :
921 [ + + + + ]: 913125486 : if(nspden_updn==1 .and. ispden>=2)exit
922 : :
923 : : ! If the norm of the gradient vanishes, then the different terms vanishes,
924 : : ! but the inverse of the gradient diverges, so skip the update.
925 [ + + ]: 557320002 : if(grho2_b_updn(indx,ispden) < 1.0d-24) then
926 : 4694730 : depsxc(ipts,ispden+nspden_updn)=zero
927 : 4694730 : cycle
928 : : end if
929 : :
930 : : ! Compute the derivative of n.e_xc wrt the
931 : : ! spin up, spin down, or total density. In the non-spin-polarized
932 : : ! case take the coefficient that will be multiplied by the
933 : : ! gradient of the total density
934 [ + + ]: 552625272 : if(nspden_updn==1)then
935 : : ! ! Definition of vxcgrho_b_updn changed in v3.3
936 [ + + ]: 355693870 : if (nvxcgrho == 3) then
937 : 351924288 : coeff=half*vxcgrho_b_updn(indx,1) + vxcgrho_b_updn(indx,3)
938 : : else
939 : 3769582 : coeff=half*vxcgrho_b_updn(indx,1)
940 : : end if
941 [ + - ]: 196931402 : else if(nspden_updn==2)then
942 [ + + ]: 196931402 : if (nvxcgrho == 3) then
943 : 193835624 : coeff=vxcgrho_b_updn(indx,ispden)
944 [ + + ]: 3095778 : else if (ispden /= 3) then
945 : 2063834 : coeff=vxcgrho_b_updn(indx,ispden)
946 : : else if (ispden == 3) then
947 : : coeff=zero
948 : : end if
949 : : end if
950 : 552625272 : depsxc(ipts,ispden+nspden_updn)=coeff
951 : :
952 : : ! In case of ixc 31 (mGGA functional fake 1),
953 : : ! skip the stress tensor to follow a LDA scheme (see doc/theory/MGGA/report_MGGA.pdf)
954 [ + + ]: 552625272 : if(ixc==31) cycle
955 : :
956 : : ! Compute the GGA contribution to the stress tensor, from the part of Vxc
957 : : ! coming from the derivative wrt the gradient (Eq (24) of PRB 50, 4327 (1994))
958 [ + + ]: 972616534 : if (compute_stress) then
959 [ + + + - : 1584550304 : if (nspden_updn==1.and.ispden==1) grho(1:3)=rhonow(ipts,1,2:4)
+ + ]
960 [ + + + + : 740084954 : if (nspden_updn==2.and.ispden==1) grho(1:3)=rhonow(ipts,2,2:4)
+ + ]
961 [ + + + + : 726433793 : if (nspden_updn==2.and.ispden==2) grho(1:3)=rhonow(ipts,1,2:4)-rhonow(ipts,2,2:4)
+ + ]
962 [ + + + + : 740082509 : if (nspden_updn==2.and.ispden==3) grho(1:3)=rhonow(ipts,1,2:4)
+ + ]
963 : 541593080 : s1=-grho(1)*grho(1)*coeff
964 : 541593080 : s2=-grho(2)*grho(2)*coeff
965 : 541593080 : s3=-grho(3)*grho(3)*coeff
966 : : !Diagonal part: +Int[n*Grad.dot.(dFxc/dgrad(n))] = -Int[grad(n).dot.(dFxc/dgrad(n))]
967 : : ! dFxc/dgrad(n) is 1/|grad(n)|.dFxc/d[grad(n)|*grad(n) = coeff * grho(:)
968 : 541593080 : dstrsxc=dstrsxc+s1+s2+s3
969 : : !Non-diagonal part: -Int[dn/dr_alpha.dFxc/grad_beta(n)]
970 : 541593080 : strsxc1_tot=strsxc1_tot+s1
971 : 541593080 : strsxc2_tot=strsxc2_tot+s2
972 : 541593080 : strsxc3_tot=strsxc3_tot+s3
973 : 541593080 : strsxc4_tot=strsxc4_tot-grho(3)*grho(2)*coeff
974 : 541593080 : strsxc5_tot=strsxc5_tot-grho(3)*grho(1)*coeff
975 : 541593080 : strsxc6_tot=strsxc6_tot-grho(2)*grho(1)*coeff
976 : : end if
977 : :
978 : : end do
979 : : end if
980 : :
981 : : ! For meta-GGAs, add the laplacian term (vxclrho)
982 : : ! and/or kinetic energy density term (vxctau)
983 [ + + + - ]: 1342704690 : if (usekden==1.and.with_vxctau) then
984 [ + + ]: 36463272 : if (nspden_updn==1)then
985 : 35140776 : vxctau(ipts,1,1) = vxctau_b_updn(indx,1)
986 [ + - ]: 1322496 : else if (nspden_updn==2)then
987 : 1322496 : vxctau(ipts,1,1) = vxctau_b_updn(indx,1)
988 : 1322496 : vxctau(ipts,2,1) = vxctau_b_updn(indx,2)
989 : : end if
990 : : end if
991 [ + + ]: 1343060029 : if (uselaplacian==1) then
992 [ + + ]: 6242728 : if (nspden_updn==1)then
993 : 5822824 : depsxc(ipts,3) = vxclrho_b_updn(indx,1)
994 [ + - ]: 419904 : else if (nspden_updn==2)then
995 : 419904 : depsxc(ipts,6) = vxclrho_b_updn(indx,1)
996 : 419904 : depsxc(ipts,7) = vxclrho_b_updn(indx,2)
997 : : end if
998 : : ! Compute the contribution to the stress tensor
999 [ + - ]: 6242728 : if (compute_stress) then
1000 [ + + ]: 12905360 : do ispden=1,2
1001 [ + + + + : 47422400 : if (nspden_updn==1.and.ispden==1) d2rho(1:6)=d2rhonow(ipts,1,1:6)
+ + ]
1002 [ + + + + ]: 12485456 : if (nspden_updn==1.and.ispden==2) exit
1003 [ + + + + : 9182056 : if (nspden_updn==2.and.ispden==1) d2rho(1:6)=d2rhonow(ipts,2,1:6)
+ + ]
1004 [ + + + + : 9182056 : if (nspden_updn==2.and.ispden==2) d2rho(1:6)=d2rhonow(ipts,1,1:6)-d2rhonow(ipts,2,1:6)
+ + ]
1005 : 6662632 : coeff=vxclrho_b_updn(indx,ispden)
1006 : : !Diagonal part: -Int[n.Lapl(dFxc/Lapl(n))] = -Int[Lapl(n).dFxc/grad(n)]
1007 : 6662632 : dstrsxc=dstrsxc-(d2rho(1)+d2rho(2)+d2rho(3))*coeff
1008 : : !Non-diagonal part: -2*Int[d2n/dr_alpha.dr_beta.dFxc/Lapl(n)]
1009 : 6662632 : strsxc1_tot=strsxc1_tot-two*d2rho(1)*coeff
1010 : 6662632 : strsxc2_tot=strsxc2_tot-two*d2rho(2)*coeff
1011 : 6662632 : strsxc3_tot=strsxc3_tot-two*d2rho(3)*coeff
1012 : 6662632 : strsxc4_tot=strsxc4_tot-two*d2rho(4)*coeff
1013 : 6662632 : strsxc5_tot=strsxc5_tot-two*d2rho(5)*coeff
1014 : 12905360 : strsxc6_tot=strsxc6_tot-two*d2rho(6)*coeff
1015 : : end do
1016 : : end if
1017 : : end if
1018 : :
1019 : : end do
1020 : :
1021 : : ! Additional electron-positron correlation terms
1022 [ + + ]: 355339 : if (ipositron==2) then
1023 : : ! Compute electron-positron XC energy per unit volume, potentials and derivatives
1024 [ - + ]: 10233 : ngr=0;if (ngrad_apn==2) ngr=npts
1025 [ - + + - : 30699 : ABI_MALLOC(fxc_apn,(npts))
+ - + - -
+ - + ]
1026 [ - + + - : 20466 : ABI_MALLOC(vxc_b_apn,(npts))
- + - + ]
1027 [ + - - - : 30699 : ABI_MALLOC(vxcgr_apn,(ngr))
+ - - + -
+ - + ]
1028 [ - + + - : 20466 : ABI_MALLOC(vxc_ep,(npts))
- + - + ]
1029 [ - + + - : 40932 : ABI_MALLOC(rhonow_apn,(npts,nspden_apn,1))
- + + - +
- + - - +
- + ]
1030 [ + - - - : 20466 : ABI_MALLOC(grho2_apn,(ngr))
- + - + ]
1031 [ + + ]: 40271769 : rhonow_apn(1:npts,1,1)=electronpositron%rhor_ep(ifft:ifft+npts-1,1)
1032 [ + + + + ]: 39218769 : if (usexcnhat==0) rhonow_apn(1:npts,1,1)=rhonow_apn(1:npts,1,1)-electronpositron%nhat_ep(ifft:ifft+npts-1,1)
1033 [ + + ]: 10233 : if (.not.electronpositron%posdensity0_limit) then
1034 : 49 : call mkdenpos(iwarnp,npts,nspden_apn,1,rhonow_apn(:,1,1),xcdata%xc_denpos)
1035 : : end if
1036 [ - + - - ]: 10233 : if (ngrad_apn==2.and.usegradient==1) then
1037 [ # # # # ]: 0 : if (nspden_apn==1) grho2_apn(:)=four*grho2_b_updn(:,1)
1038 [ # # # # ]: 0 : if (nspden_apn==2) grho2_apn(:)=grho2_b_updn(:,3)
1039 : : end if
1040 [ + + ]: 10233 : if (ndvxc==0) then
1041 : : call xcpositron(fxc_apn,grho2_apn,electronpositron%ixcpositron,ngr,npts,&
1042 : : & electronpositron%posdensity0_limit,rho_b,&
1043 : 2930 : & rhonow_apn(:,1,1),vxc_b_apn,vxcgr_apn,vxc_ep)
1044 : : else
1045 [ - + + - : 21909 : ABI_MALLOC(dvxc_apn,(npts))
+ - + - -
+ - + -
+ ]
1046 : : call xcpositron(fxc_apn,grho2_apn,electronpositron%ixcpositron,ngr,npts,&
1047 : : & electronpositron%posdensity0_limit,rho_b,&
1048 : 7303 : & rhonow_apn(:,1,1),vxc_b_apn,vxcgr_apn,vxc_ep,dvxce=dvxc_apn)
1049 : : end if
1050 : : ! Accumulate electron-positron XC energies
1051 : 10233 : s1=zero
1052 [ + + ]: 40271769 : do ipts=1,npts
1053 : 40271769 : s1=s1+fxc_apn(ipts)
1054 : : end do
1055 : 10233 : electronpositron%e_xc=electronpositron%e_xc+s1*ucvol/dble(nfftot)
1056 : : ! Add electron-positron dVxc_el/dRho_el to electron-electron one
1057 [ - + - - ]: 10233 : if (ndvxc==1) dvxc_b(:,1)=dvxc_b(:,1)+dvxc_apn(:)
1058 [ + + ]: 10233 : if (ndvxc==3) then
1059 [ + + ]: 28508224 : dvxc_b(:,1)=dvxc_b(:,1)+four*dvxc_apn(:)
1060 [ + + ]: 28508224 : dvxc_b(:,2)=dvxc_b(:,2)+four*dvxc_apn(:)
1061 [ + + ]: 28508224 : dvxc_b(:,3)=dvxc_b(:,3)+four*dvxc_apn(:)
1062 : : end if
1063 [ + + ]: 10233 : if (ndvxc==15) then
1064 [ + + ]: 243063 : dvxc_b(:, 9)=dvxc_b(:, 9)+four*dvxc_apn(:)
1065 [ + + ]: 243063 : dvxc_b(:,10)=dvxc_b(:,10)+four*dvxc_apn(:)
1066 [ + + ]: 243063 : dvxc_b(:,11)=dvxc_b(:,11)+four*dvxc_apn(:)
1067 : : end if
1068 : : ! Modify stresses - Compute factors for GGA
1069 [ + + ]: 40271769 : do ipts=ifft,ifft+npts-1
1070 : 40261536 : indx=ipts-ifft+1
1071 : 40261536 : depsxc_apn(ipts,1)=vxc_b_apn(indx)
1072 [ - + ]: 40261536 : if (ngrad_apn==2) depsxc_apn(ipts,2)=vxcgr_apn(indx)
1073 [ + - ]: 40271769 : if (compute_stress) then
1074 : 40261536 : dstrsxc=dstrsxc+fxc_apn(indx)-rho_b(indx)*vxc_b_apn(indx)
1075 [ - + ]: 40261536 : if (ngrad_apn==2) then
1076 [ # # # # ]: 0 : if (nspden_updn==1) grho(1:3)=rhonow(ipts,1,2:4)
1077 [ # # # # : 0 : if (ispden==1 .and. nspden_updn==2) grho(1:3)=rhonow(ipts,2,2:4)
# # ]
1078 [ # # # # : 0 : if (ispden==2 .and. nspden_updn==2) grho(1:3)=rhonow(ipts,1,2:4)-rhonow(ipts,2,2:4)
# # ]
1079 [ # # # # : 0 : if (ispden==3 .and. nspden_updn==2) grho(1:3)=rhonow(ipts,1,2:4)
# # ]
1080 : 0 : s1=-grho(1)*grho(1)*vxcgr_apn(indx)
1081 : 0 : s2=-grho(2)*grho(2)*vxcgr_apn(indx)
1082 : 0 : s3=-grho(3)*grho(3)*vxcgr_apn(indx)
1083 : 0 : dstrsxc=dstrsxc+s1+s2+s3
1084 : 0 : strsxc1_tot=strsxc1_tot+s1
1085 : 0 : strsxc2_tot=strsxc2_tot+s2
1086 : 0 : strsxc3_tot=strsxc3_tot+s3
1087 : 0 : strsxc4_tot=strsxc4_tot-grho(3)*grho(2)*vxcgr_apn(indx)
1088 : 0 : strsxc5_tot=strsxc5_tot-grho(3)*grho(1)*vxcgr_apn(indx)
1089 : 0 : strsxc6_tot=strsxc6_tot-grho(2)*grho(1)*vxcgr_apn(indx)
1090 : : end if ! GGA
1091 : : end if
1092 : : end do ! ipts
1093 : : ! Deallocations
1094 [ - + ]: 10233 : ABI_FREE(fxc_apn)
1095 [ - + ]: 10233 : ABI_FREE(vxc_b_apn)
1096 [ - + ]: 10233 : ABI_FREE(vxcgr_apn)
1097 [ - + ]: 10233 : ABI_FREE(vxc_ep)
1098 [ - + ]: 10233 : ABI_FREE(rhonow_apn)
1099 [ - + ]: 10233 : ABI_FREE(grho2_apn)
1100 [ + + ]: 10233 : if (ndvxc>0) then
1101 [ - + ]: 7303 : ABI_FREE(dvxc_apn)
1102 : : end if
1103 : : end if
1104 : :
1105 : : ! Transfer the xc kernel (if this must be done, and has not yet been done)
1106 [ + + + + : 355339 : if (nkxc_eff>0.and.ndvxc>0 .and. (xcdata%usefock==0 .or. auxc_ixc==0)) then
+ + - + ]
1107 [ + + + + ]: 129881 : if (nkxc_eff==1.and.ndvxc==15) then
1108 [ + + ]: 95546690 : kxc(ifft:ifft+npts-1,1)=half*(dvxc_b(1:npts,1)+dvxc_b(1:npts,9)+dvxc_b(1:npts,10))
1109 [ + + + + ]: 110422 : else if (nkxc_eff==3.and.ndvxc==15) then
1110 [ + + ]: 30410879 : kxc(ifft:ifft+npts-1,1)=dvxc_b(1:npts,1)+dvxc_b(1:npts,9)
1111 [ + + ]: 30410879 : kxc(ifft:ifft+npts-1,2)=dvxc_b(1:npts,10)
1112 [ + + ]: 30410879 : kxc(ifft:ifft+npts-1,3)=dvxc_b(1:npts,2)+dvxc_b(1:npts,11)
1113 [ + + + + ]: 102281 : else if (nkxc_eff==7.and.ndvxc==8) then
1114 [ + + ]: 9603 : kxc(ifft:ifft+npts-1,1)=half*dvxc_b(1:npts,1)
1115 [ + + ]: 9603 : kxc(ifft:ifft+npts-1,2)=half*dvxc_b(1:npts,3)
1116 [ + + ]: 9603 : kxc(ifft:ifft+npts-1,3)=quarter*dvxc_b(1:npts,5)
1117 [ + + ]: 9603 : kxc(ifft:ifft+npts-1,4)=eighth*dvxc_b(1:npts,7)
1118 [ + - ]: 1706 : else if (nkxc_eff==7.and.ndvxc==15) then
1119 [ + + ]: 6644962 : kxc(ifft:ifft+npts-1,1)=half*(dvxc_b(1:npts,1)+dvxc_b(1:npts,9)+dvxc_b(1:npts,10))
1120 [ + + ]: 6644962 : kxc(ifft:ifft+npts-1,2)=half*dvxc_b(1:npts,3)+dvxc_b(1:npts,12)
1121 [ + + ]: 6644962 : kxc(ifft:ifft+npts-1,3)=quarter*dvxc_b(1:npts,5)+dvxc_b(1:npts,13)
1122 [ + + ]: 6644962 : kxc(ifft:ifft+npts-1,4)=eighth*dvxc_b(1:npts,7)+dvxc_b(1:npts,15)
1123 [ + + + - ]: 100572 : else if (nkxc_eff==19.and.ndvxc==15) then
1124 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,1)=dvxc_b(1:npts,1)+dvxc_b(1:npts,9)
1125 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,2)=dvxc_b(1:npts,10)
1126 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,3)=dvxc_b(1:npts,2)+dvxc_b(1:npts,11)
1127 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,4)=dvxc_b(1:npts,3)
1128 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,5)=dvxc_b(1:npts,4)
1129 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,6)=dvxc_b(1:npts,5)
1130 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,7)=dvxc_b(1:npts,6)
1131 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,8)=dvxc_b(1:npts,7)
1132 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,9)=dvxc_b(1:npts,8)
1133 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,10)=dvxc_b(1:npts,12)
1134 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,11)=dvxc_b(1:npts,13)
1135 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,12)=dvxc_b(1:npts,14)
1136 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,13)=dvxc_b(1:npts,15)
1137 : : else ! All other cases
1138 [ + + + + ]: 492353314 : kxc(ifft:ifft+npts-1,1:nkxc_eff)=zero
1139 [ + + + + ]: 492161266 : kxc(ifft:ifft+npts-1,1:min(nkxc_eff,ndvxc))=dvxc_b(1:npts,1:min(nkxc_eff,ndvxc))
1140 : : end if
1141 [ + + ]: 129881 : if (nkxc_eff==7) then
1142 [ + + ]: 6654565 : kxc(ifft:ifft+npts-1,5)=rhonow(ifft:ifft+npts-1,1,2)
1143 [ + + ]: 6654565 : kxc(ifft:ifft+npts-1,6)=rhonow(ifft:ifft+npts-1,1,3)
1144 [ + + ]: 6654565 : kxc(ifft:ifft+npts-1,7)=rhonow(ifft:ifft+npts-1,1,4)
1145 [ + + ]: 128172 : else if (nkxc_eff==19) then
1146 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,14)=rhonow(ifft:ifft+npts-1,1,2)
1147 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,15)=rhonow(ifft:ifft+npts-1,2,2)
1148 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,16)=rhonow(ifft:ifft+npts-1,1,3)
1149 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,17)=rhonow(ifft:ifft+npts-1,2,3)
1150 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,18)=rhonow(ifft:ifft+npts-1,1,4)
1151 [ + + ]: 244412 : kxc(ifft:ifft+npts-1,19)=rhonow(ifft:ifft+npts-1,2,4)
1152 : : end if
1153 : : end if
1154 : :
1155 : : ! Transfer the XC 3rd-derivative
1156 [ + + + - : 355339 : if (abs(option)==3.and.order==3.and.nd2vxc>0) then
+ - ]
1157 [ + + + + ]: 278813 : k3xc(ifft:ifft+npts-1,1:nd2vxc)=d2vxc_b(1:npts,1:nd2vxc)
1158 : : end if
1159 : :
1160 : : ! Add the diagonal part to the xc stress
1161 [ + + ]: 355339 : if (compute_stress) then
1162 : 350335 : strsxc1_tot=strsxc1_tot+dstrsxc
1163 : 350335 : strsxc2_tot=strsxc2_tot+dstrsxc
1164 : 350335 : strsxc3_tot=strsxc3_tot+dstrsxc
1165 : : end if
1166 : :
1167 : : ! Accumulate integral of |Grad_rho|/Rho (to be used for TB09 XC)
1168 [ + + + + ]: 355339 : if (present(grho1_over_rho1).and.ixc<0) then
1169 [ - + ]: 17109 : if (present(xc_funcs)) then
1170 : 0 : test_tb09=libxc_functionals_is_tb09(xc_functionals=xc_funcs)
1171 : : else
1172 : 17109 : test_tb09=libxc_functionals_is_tb09()
1173 : : end if
1174 [ + + ]: 17109 : if (test_tb09) then
1175 [ - + ]: 10 : factor=merge(two,one,nspden_updn==1)
1176 [ - + ]: 10 : jj=merge(1,3,nspden_updn==1)
1177 [ + + ]: 218394 : do ipts=ifft,ifft+npts-1
1178 : 218384 : indx=ipts-ifft+1
1179 [ + - ]: 218394 : if (abs(rho_b(indx))>tol10) then
1180 : 218384 : grho1_over_rho1=grho1_over_rho1+factor*sqrt(grho2_b_updn(indx,jj))/rho_b(indx)
1181 : : end if
1182 : : end do
1183 : : end if
1184 : : end if
1185 : :
1186 [ - + ]: 355339 : ABI_FREE(exc_b)
1187 [ - + ]: 355339 : ABI_FREE(tsxc_b)
1188 [ - + ]: 355339 : ABI_FREE(rho_b)
1189 [ - + ]: 355339 : ABI_FREE(rho_b_updn)
1190 [ - + ]: 355339 : ABI_FREE(grho2_b_updn)
1191 [ - + ]: 355339 : ABI_FREE(vxcrho_b_updn)
1192 [ - + ]: 355339 : ABI_FREE(dvxc_b)
1193 [ - + ]: 355339 : ABI_FREE(d2vxc_b)
1194 [ - + ]: 355339 : ABI_FREE(vxcgrho_b_updn)
1195 [ - + ]: 355339 : ABI_FREE(vxclrho_b_updn)
1196 [ - + ]: 355339 : ABI_FREE(lrho_b_updn)
1197 [ - + ]: 355339 : ABI_FREE(tau_b_updn)
1198 [ - + + + ]: 355339 : ABI_FREE(vxctau_b_updn)
1199 : :
1200 : : ! End of the loop on blocks of data
1201 : : end do
1202 : :
1203 [ + + ]: 60256 : if (compute_stress) then
1204 : 59108 : strsxc(1)=strsxc1_tot
1205 : 59108 : strsxc(2)=strsxc2_tot
1206 : 59108 : strsxc(3)=strsxc3_tot
1207 : 59108 : strsxc(4)=strsxc4_tot
1208 : 59108 : strsxc(5)=strsxc5_tot
1209 : 59108 : strsxc(6)=strsxc6_tot
1210 : : end if
1211 : :
1212 : : ! If GGA, multiply the gradient of the density by the proper
1213 : : ! local partial derivatives of the XC functional
1214 : 60256 : rhonow_ptr => rhonow
1215 [ + + ]: 60256 : if (ipositron==2) then
1216 [ - + + - : 1440 : ABI_MALLOC(rhonow_ptr,(nfft,nspden_eff,ngrad*ngrad))
- + + - -
+ + - + -
+ - - + -
+ ]
1217 [ + - + + : 86842980 : rhonow_ptr=rhonow
+ + + + +
+ + + +
+ ]
1218 : : end if
1219 [ + + + - ]: 60256 : if(ngrad==2 .and. ixc/=13)then
1220 : 13952 : call xcmult(depsxc,nfft,ngrad,nspden_eff,nspgrad,rhonow_ptr)
1221 : : end if
1222 : :
1223 : : ! Compute contribution from this grid to vxc, and ADD to existing vxc
1224 [ + + ]: 60256 : if (nspden/=4) then
1225 [ + + ]: 57673 : if(with_vxctau)then
1226 : : call xcpot(cplex,gprimd,ishift,uselaplacian,mpi_enreg,nfft,ngfft,ngrad,nspden_eff,nspgrad,&
1227 [ + - - + : 666 : & qphon,depsxc=depsxc,rhonow=rhonow_ptr,vxc=vxc,vxctau=vxctau)
- - - - -
- - - + -
- + - - -
- - - ]
1228 : : else
1229 : : call xcpot(cplex,gprimd,ishift,uselaplacian,mpi_enreg,nfft,ngfft,ngrad,nspden_eff,nspgrad,&
1230 : 57007 : & qphon,depsxc=depsxc,rhonow=rhonow_ptr,vxc=vxc)
1231 : : end if
1232 : :
1233 : : else
1234 : :
1235 : : ! If non-collinear magnetism, restore potential in proper axis before adding it
1236 [ - + + - : 7749 : ABI_MALLOC(vxcrho_b_updn,(nfft,4))
+ - + - -
+ - + ]
1237 [ + + + + ]: 128054939 : vxcrho_b_updn=zero
1238 : : call xcpot(cplex,gprimd,ishift,uselaplacian,mpi_enreg,nfft,ngfft,ngrad,nspden_eff,nspgrad,&
1239 : 2583 : & qphon,depsxc=depsxc,rhonow=rhonow_ptr,vxc=vxcrho_b_updn)
1240 [ + + ]: 32013089 : do ifft=1,nfft
1241 : 32010506 : dvdn=half*(vxcrho_b_updn(ifft,1)+vxcrho_b_updn(ifft,2))
1242 [ + + ]: 32013089 : if(m_norm(ifft)>m_norm_min) then
1243 : : ! if(m_norm(ifft)>rhor_(ifft,1)*tol10+tol14) then
1244 : 27233195 : dvdz=half*(vxcrho_b_updn(ifft,1)-vxcrho_b_updn(ifft,2))/m_norm(ifft)
1245 : 27233195 : vxc(ifft,1)=vxc(ifft,1)+dvdn+rhor_(ifft,4)*dvdz
1246 : 27233195 : vxc(ifft,2)=vxc(ifft,2)+dvdn-rhor_(ifft,4)*dvdz
1247 : 27233195 : vxc(ifft,3)=vxc(ifft,3)+rhor_(ifft,2)*dvdz
1248 : 27233195 : vxc(ifft,4)=vxc(ifft,4)-rhor_(ifft,3)*dvdz
1249 : : else
1250 [ + + ]: 14331933 : vxc(ifft,1:2)=vxc(ifft,1:2)+dvdn
1251 : : end if
1252 : : end do
1253 [ - + ]: 2583 : ABI_FREE(vxcrho_b_updn)
1254 : : end if
1255 [ + + ]: 60256 : if (ipositron==2) then
1256 [ - + ]: 288 : ABI_FREE(rhonow_ptr)
1257 : : end if
1258 : 60256 : nullify(rhonow_ptr)
1259 : :
1260 : : ! Add electron-positron XC potential to electron-electron one
1261 : : ! Eventually compute GGA contribution
1262 : 56053 : if (ipositron==2) then
1263 [ - + + - : 1440 : ABI_MALLOC(rhonow_apn,(nfft,nspden_apn,ngrad_apn**2))
- + + - -
+ + - + -
+ - - + -
+ ]
1264 [ + + + + ]: 40262112 : rhonow_apn(1:nfft,1,1:ngrad_apn**2)=rhonow(1:nfft,1,1:ngrad_apn**2)
1265 [ - + ]: 288 : if (ngrad_apn==2) then
1266 : 0 : call xcmult(depsxc_apn,nfft,ngrad_apn,nspden_apn,ngrad_apn,rhonow_apn)
1267 : : end if
1268 [ - + + - : 1152 : ABI_MALLOC(vxc_apn,(nfft,nspden_apn))
- + + - +
- + - - +
- + ]
1269 [ + + + + ]: 40262112 : vxc_apn=zero
1270 : : call xcpot(cplex,gprimd,ishift,0,mpi_enreg,nfft,ngfft,ngrad_apn,&
1271 : 288 : & nspden_apn,ngrad_apn,qphon,depsxc=depsxc_apn,rhonow=rhonow_apn,vxc=vxc_apn)
1272 [ + + ]: 40261824 : vxc(:,1)=vxc(:,1)+vxc_apn(:,1)
1273 [ - + - - ]: 288 : if (nspden_updn==2) vxc(:,2)=vxc(:,2)+vxc_apn(:,1)
1274 : 288 : s1=zero
1275 [ + + ]: 40261824 : do ipts=1,nfft
1276 : 40261824 : s1=s1+vxc_apn(ipts,1)*rhonow(ipts,1,1)
1277 : : end do
1278 : 288 : electronpositron%e_xcdc=electronpositron%e_xcdc+s1*ucvol/dble(nfftot)
1279 [ - + ]: 288 : ABI_FREE(rhonow_apn)
1280 [ - + ]: 288 : ABI_FREE(vxc_apn)
1281 [ - + ]: 288 : ABI_FREE(depsxc_apn)
1282 : : end if
1283 : :
1284 : : ! End loop on unshifted or shifted grids
1285 : : end do
1286 : :
1287 : : ! Calculate van der Waals correction when requested
1288 : : #if defined DEV_YP_VDWXC
1289 : : if ( (xcdata%vdw_xc > 0) .and. (xcdata%vdw_xc < 3) .and. (xc_vdw_status()) ) then
1290 : : strsxc_vdw(:,:)=zero
1291 : : call xc_vdw_aggregate(ucvol,gprimd,nfft,nspden_updn,ngrad*ngrad, &
1292 : : & ngfft(1),ngfft(2),ngfft(3),rhonow, &
1293 : : & deltae_vdw,exc_vdw,decdrho_vdw,decdgrho_vdw,strsxc_vdw)
1294 : : end if
1295 : : #else
1296 [ - + ]: 56053 : if ( (xcdata%vdw_xc > 0) .and. (xcdata%vdw_xc < 3) ) then
1297 : : write(message,'(3a)')&
1298 : 0 : & 'vdW-DF functionals are not fully operational yet.',ch10,&
1299 : 0 : & 'Action : modify vdw_xc'
1300 : 0 : ABI_ERROR(message)
1301 : : end if
1302 : : #endif
1303 : : ! Normalize bigexc, bigsxc, strsxc and vxc
1304 : 56053 : divshft=one/dble(xcdata%intxc+1)
1305 : 56053 : bigexc=rhoexc*ucvol/dble(nfftot)*divshft
1306 : 56053 : bigsxc=rhotsxc*ucvol/dble(nfftot)*divshft/xcdata%tphysel
1307 [ + + + + ]: 1559239941 : vxc=vxc*divshft
1308 [ + + + + ]: 385801 : if (compute_stress) strsxc(:)=strsxc(:)/dble(nfftot)*divshft
1309 [ + + + + : 151204553 : if (with_vxctau) vxctau=vxctau*divshft
+ + + + ]
1310 [ + + ]: 56053 : if (present(grho1_over_rho1)) grho1_over_rho1=grho1_over_rho1*ucvol/dble(nfftot)*divshft
1311 : :
1312 : : ! Reduction in case of FFT distribution
1313 [ + + ]: 56053 : if (nproc_fft>1)then
1314 : 2920 : call timab(48,1,tsec)
1315 : 2920 : call xmpi_sum(bigexc,comm_fft,ierr)
1316 : 2920 : call xmpi_sum(bigsxc,comm_fft,ierr)
1317 [ + - ]: 2920 : if (compute_stress) then
1318 : 2920 : call xmpi_sum(strsxc,comm_fft,ierr)
1319 : : end if
1320 [ + + ]: 2920 : if (present(grho1_over_rho1)) then
1321 : 444 : call xmpi_sum(grho1_over_rho1,comm_fft ,ierr)
1322 : : end if
1323 [ - + ]: 2920 : if (ipositron==2) then
1324 : 0 : s1=electronpositron%e_xc;s2=electronpositron%e_xcdc
1325 : 0 : call xmpi_sum(s1,comm_fft ,ierr)
1326 : 0 : call xmpi_sum(s2,comm_fft ,ierr)
1327 : 0 : electronpositron%e_xc=s1;electronpositron%e_xcdc=s2
1328 : : end if
1329 : 2920 : call timab(48,2,tsec)
1330 : : end if
1331 : :
1332 : : ! Compute vxcavg
1333 : 56053 : call mean_fftr(vxc,vxcmean,nfft,nfftot,min(nspden,2),mpi_comm_sphgrid=comm_fft)
1334 [ + + ]: 56053 : if(nspden==1)then
1335 : 47554 : vxcavg=vxcmean(1)
1336 : : else
1337 : 8499 : vxcavg=half*(vxcmean(1)+vxcmean(2))
1338 : : end if
1339 : :
1340 [ - + ]: 56053 : ABI_FREE(depsxc)
1341 [ - + ]: 56053 : ABI_FREE(rhonow)
1342 [ - + ]: 56053 : ABI_FREE(lrhonow)
1343 [ - + ]: 56053 : ABI_FREE(taunow)
1344 [ + + + + ]: 56053 : if (compute_stress.and.uselaplacian>0) then
1345 [ - + ]: 135 : ABI_FREE(d2rhonow)
1346 : : end if
1347 [ + + + + ]: 56053 : if (need_nhat.or.non_magnetic_xc) then
1348 [ - + ]: 11163 : ABI_FREE(rhor_)
1349 : : end if
1350 [ + + - + ]: 56053 : if ((usekden==1).and.(non_magnetic_xc)) then
1351 [ # # ]: 0 : ABI_FREE(taur_)
1352 : : end if
1353 [ + + ]: 56053 : if (allocated(m_norm)) then
1354 : 2583 : ABI_FREE(m_norm)
1355 : : end if
1356 : :
1357 : : end if
1358 : :
1359 : : !Treat separately the Fermi-Amaldi correction.
1360 [ + + ]: 56830 : if (ixc==20 .or. ixc==21 .or. ixc==22) then
1361 [ + - ]: 44 : if(present(vhartr))then
1362 : :
1363 : : ! Fermi-Amaldi correction : minus Hartree divided by the
1364 : : ! number of electrons per unit cell. This is not size consistent, but
1365 : : ! interesting for isolated systems with a few electrons.
1366 : : ! nelect=ucvol*rhog(1,1)
1367 : 44 : factor=-one/xcdata%nelect
1368 [ + + ]: 20663556 : vxc(:,1)=factor*vhartr(:)
1369 [ + + + + ]: 17915948 : if(nspden>=2) vxc(:,2)=factor*vhartr(:)
1370 : :
1371 : : ! Compute corresponding xc energy and stress as well as vxcavg
1372 : 44 : call dotprod_vn(1,rhor,bigexc,doti,nfft,nfftot,1,1,vxc,ucvol,mpi_comm_sphgrid=comm_fft)
1373 : 44 : bigexc=half*bigexc
1374 [ + - + + ]: 176 : if (compute_stress) strsxc(1:3)=-bigexc/ucvol
1375 : :
1376 : : ! Compute average of vxc (one component only).
1377 : 44 : call mean_fftr(vxc,vxcmean,nfft,nfftot,1,mpi_comm_sphgrid=comm_fft)
1378 : 44 : vxcavg = vxcmean(1)
1379 : : ! For ixc=20, the local exchange-correlation kernel is zero, but the Hartree
1380 : : ! kernel will be modified in tddft. No other use of kxc should be made with ixc==20
1381 [ + + + + : 118848 : if(nkxc/=0 .and. ixc==20) kxc(:,:)=zero
+ + + + ]
1382 : : ! For ixc=21 or 22, the LDA (ixc=1) kernel has been computed previously.
1383 : :
1384 : : else
1385 : :
1386 : 0 : ABI_BUG('When ixc=20,21 or 22, vhartr needs to be present in the call to rhotoxc !')
1387 : :
1388 : : end if
1389 : :
1390 : : end if
1391 : :
1392 : : !Add van der Waals terms
1393 : : #if defined DEV_YP_VDWXC
1394 : : if ( (xcdata%vdw_xc > 0) .and. (xcdata%vdw_xc < 10) .and. (xc_vdw_status()) ) then
1395 : : bigexc = bigexc + exc_vdw + deltae_vdw
1396 : : do ispden=1,nspden
1397 : : vxc(:,ispden) = vxc(:,ispden) + decdrho_vdw(:,ispden)
1398 : : end do
1399 : : if (compute_stress) then
1400 : : strsxc(1) = strsxc(1) + strsxc_vdw(1,1)
1401 : : strsxc(2) = strsxc(2) + strsxc_vdw(2,2)
1402 : : strsxc(3) = strsxc(3) + strsxc_vdw(3,3)
1403 : : strsxc(4) = strsxc(4) + strsxc_vdw(3,2)
1404 : : strsxc(5) = strsxc(5) + strsxc_vdw(3,1)
1405 : : strsxc(6) = strsxc(6) + strsxc_vdw(2,1)
1406 : : end if
1407 : : end if
1408 : : #endif
1409 [ + + ]: 56830 : if ( present(exc_vdw_out) ) exc_vdw_out = exc_vdw
1410 : :
1411 : : !In case we have an entropy associated with XC contribution
1412 : : !(e.g. using finite-temperature exchange-correlation functionals),
1413 : : !we retrieve exchange-correlation internal energy bigexc using entropy bigsxc
1414 [ + + ]: 56830 : if(abs(bigsxc)>tiny(zero)) bigexc=bigexc+xcdata%tphysel*bigsxc
1415 : :
1416 : 56830 : call timab(81,2,tsec)
1417 : :
1418 : : DBG_EXIT("COLL")
1419 : :
1420 [ + + + + : 168981 : end subroutine rhotoxc
- + - + -
+ - + ]
1421 : : !!***
1422 : :
1423 : : end module m_rhotoxc
1424 : : !!***
|
