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 112194 : subroutine rhotoxc(bigexc,bigsxc,kxc,mpi_enreg,nfft,ngfft, &
259 56097 : & nhat,nhatdim,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,option, &
260 112194 : & rhor,rprimd,usexcnhat,vxc,vxcavg,xccc3d,xcdata, &
261 110685 : & 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 56097 : real(dp),allocatable :: d2rhonow(:,:,:)
306 56097 : real(dp),allocatable :: d2vxc_b(:,:),depsxc(:,:),depsxc_apn(:,:),dvxc_apn(:),dvxc_b(:,:)
307 56097 : real(dp),allocatable :: exc_b(:),tsxc_b(:),fxc_apn(:),grho2_apn(:),grho2_b_updn(:,:)
308 56097 : real(dp),allocatable :: lrhonow(:,:),lrho_b_updn(:,:)
309 56097 : real(dp),allocatable :: m_norm(:),nhat_up(:),rho_b_updn(:,:),rho_b(:),rhonow_apn(:,:,:)
310 56097 : real(dp),allocatable :: tau_b_updn(:,:),vxc_apn(:,:),vxcgr_apn(:),vxcgrho_b_updn(:,:),vxcrho_b_updn(:,:)
311 56097 : real(dp),allocatable :: vxc_b_apn(:),vxc_ep(:),vxctau_b_updn(:,:),vxclrho_b_updn(:,:)
312 56097 : real(dp),allocatable,target :: rhonow(:,:,:),taunow(:,:,:)
313 56097 : real(dp),pointer :: rhocorval(:,:),rhor_(:,:),taucorval(:,:),taur_(:,:)
314 56097 : real(dp),contiguous, pointer :: rhonow_ptr(:,:,:)
315 112194 : real(dp) :: decdrho_vdw(nfft,xcdata%nspden),decdgrho_vdw(nfft,3,xcdata%nspden)
316 : real(dp) :: strsxc_vdw(3,3)
317 168291 : 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 56097 : call timab(81,1,tsec)
326 :
327 : !Optional arguments
328 56097 : compute_stress=present(strsxc)
329 56097 : my_add_tfw=.false.;if (present(add_tfw)) my_add_tfw=add_tfw
330 :
331 : !Useful scalars
332 56097 : nspden=xcdata%nspden
333 56097 : ixc=xcdata%ixc
334 56097 : auxc_ixc=xcdata%auxc_ixc
335 56097 : n3xctau=0
336 :
337 : !nspden_updn: 1 for non-polarized, 2 for polarized
338 56097 : 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 56097 : order=1
347 56097 : if(option==2.or.option==10.or.option==12)order=2
348 56097 : if(option==-2)order=-2
349 56097 : if(option==3)order=3
350 :
351 : !Sizes of local arrays
352 56097 : 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 55576 : & 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 56097 : 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 56097 : 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 56097 : nkxc_eff=nkxc;if (option==10.or.option==12) nkxc_eff=min(nkxc,3)
375 :
376 : !Check options
377 56097 : 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 56097 : 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 56097 : if(ixc<0) then
387 12499 : if (present(xc_funcs)) then
388 521 : ixc_from_lib=libxc_functionals_ixc(xc_functionals=xc_funcs)
389 : else
390 11978 : 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 12499 : 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 56097 : with_vxctau=(present(vxctau))
404 218178 : if (with_vxctau) with_vxctau=(size(vxctau)>0)
405 56097 : 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 56097 : 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 56097 : 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 56097 : 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 56097 : cplex=1
440 56097 : qphon(:)=zero
441 56097 : iwarn=0
442 56097 : nfftot=ngfft(1)*ngfft(2)*ngfft(3)
443 :
444 : !Initializations
445 56097 : bigexc=zero
446 56097 : bigsxc=zero
447 56097 : rhoexc=zero
448 56097 : rhotsxc=zero
449 1553019157 : vxc(:,:)=zero
450 56097 : vxcavg=zero
451 56097 : if (compute_stress) then
452 54982 : strsxc(:)=zero
453 54982 : strsxc1_tot=zero;strsxc2_tot=zero;strsxc3_tot=zero
454 54982 : strsxc4_tot=zero;strsxc5_tot=zero;strsxc6_tot=zero
455 : strsxc_vdw(:,:)=zero
456 : end if
457 151204597 : if (with_vxctau) vxctau(:,:,:)=zero
458 729381405 : if (nkxc/=0) kxc(:,:)=zero
459 370818 : if(abs(option)==3.and.nk3xc/=0) k3xc(:,:)=zero
460 56097 : ipositron=0
461 56097 : 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 56097 : deltae_vdw = zero
469 56097 : exc_vdw = zero
470 1553019157 : decdrho_vdw(:,:) = zero
471 4659015063 : decdgrho_vdw(:,:,:) = zero
472 56097 : if (present(grho1_over_rho1)) grho1_over_rho1=zero
473 :
474 56097 : 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 55362 : else if (ixc/=20) then
479 :
480 : ! Test: has a compensation density to be added/substracted (PAW) ?
481 55320 : need_nhat=(nhatdim==1.and.usexcnhat==0)
482 55320 : need_nhatgr=(nhatdim==1.and.nhatgrdim==1.and.ngrad==2.and.xcdata%intxc==0)
483 55320 : 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 55320 : nspgrad=nspden_updn*ngrad;if(nspden_updn==2.and.ngrad==2)nspgrad=5
503 55320 : if(uselaplacian==1) nspgrad=nspgrad+nspden_updn
504 221280 : ABI_MALLOC(depsxc,(nfft,nspgrad))
505 1934124459 : depsxc(:,:)=zero
506 :
507 : ! PAW: select the valence density (and magnetization) to use:
508 : ! link the correct density, according to usexcnhat option
509 55320 : if ((.not.need_nhat).and.(.not.non_magnetic_xc)) then
510 44890 : rhor_ => rhor
511 : else
512 41720 : ABI_MALLOC(rhor_,(nfft,nspden))
513 10430 : if (need_nhat) then
514 22288 : do ispden=1,nspden
515 373058406 : do ifft=1,nfft
516 373048263 : 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 10430 : 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 55320 : 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 55320 : 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 55320 : m_norm_min=tol8 ! EB: tol14 is still too small, tests are underway
559 55320 : 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 55320 : if (n3xccc>0.or.test_nhat.or.nspden_eff/=nspden) then
570 123976 : ABI_MALLOC(rhocorval,(nfft,nspden_eff))
571 30994 : if (nspden==nspden_eff) then
572 1627532391 : 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 55320 : 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 54789 : taucorval => taur_
591 : end if
592 :
593 : ! Add core electron density to effective density
594 55320 : if (n3xccc>0) then
595 621275912 : rhocorval(:,1)=rhocorval(:,1)+xccc3d(:)
596 26017 : if(nspden_eff==2) then
597 87540113 : rhocorval(:,2)=rhocorval(:,2)+half*xccc3d(:)
598 : end if
599 : end if
600 55320 : 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 55320 : 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 276600 : ABI_MALLOC(rhonow,(nfft,nspden_eff,ngrad*ngrad))
635 221280 : ABI_MALLOC(lrhonow,(nfft,nspden_eff*uselaplacian))
636 276600 : ABI_MALLOC(taunow,(nfft,nspden_eff,usekden))
637 55320 : 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 114843 : 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 59523 : 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 59388 : & qphon,rhocorval,rhonow)
660 : end if
661 59523 : 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 59523 : 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 59523 : if (ishift==xcdata%intxc) then
689 55320 : if (n3xccc>0.or.test_nhat.or.nspden_eff/=nspden) then
690 30994 : ABI_FREE(rhocorval)
691 : end if
692 55320 : if (usekden==1.and.(n3xccc>0.or.nspden_eff/=nspden)) then
693 531 : ABI_FREE(taucorval)
694 : end if
695 55320 : 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 59523 : 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 59523 : call mkdenpos(iwarn,nfft,nspden_updn,1,rhonow(:,1:nspden_updn,1),xcdata%xc_denpos)
724 59523 : 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 59523 : mpts=4000
737 59523 : if (usekden==1) mpts=nfft ! Why?
738 :
739 59523 : do ifft=1,nfft,mpts
740 : ! npts=mpts
741 : ! npts is the number of points to be treated in this bunch
742 338939 : npts=min(nfft-ifft+1,mpts)
743 :
744 : ! Allocation of mandatory arguments of drivexc
745 1016817 : ABI_MALLOC(exc_b,(npts))
746 677878 : ABI_MALLOC(tsxc_b,(npts))
747 677878 : ABI_MALLOC(rho_b,(npts))
748 1355756 : ABI_MALLOC(rho_b_updn,(npts,nspden_updn))
749 1016817 : ABI_MALLOC(vxcrho_b_updn,(npts,nspden_updn))
750 1483185645 : vxcrho_b_updn(:,:)=zero
751 :
752 : ! Allocation of optional arguments of drivexc
753 1355756 : ABI_MALLOC(grho2_b_updn,(npts,(2*nspden_updn-1)*usegradient))
754 1355756 : ABI_MALLOC(lrho_b_updn,(npts,nspden_updn*uselaplacian))
755 1355756 : ABI_MALLOC(tau_b_updn,(npts,nspden_updn*usekden))
756 1355756 : ABI_MALLOC(vxcgrho_b_updn,(npts,nvxcgrho))
757 1355756 : ABI_MALLOC(vxclrho_b_updn,(npts,nvxclrho))
758 1355756 : ABI_MALLOC(vxctau_b_updn,(npts,nvxctau))
759 1355756 : ABI_MALLOC(dvxc_b,(npts,ndvxc))
760 1355756 : ABI_MALLOC(d2vxc_b,(npts,nd2vxc))
761 1092574470 : if (nvxcgrho>0) vxcgrho_b_updn(:,:)=zero
762 7002012 : if (nvxclrho>0) vxclrho_b_updn(:,:)=zero
763 38125398 : if (nvxctau>0) vxctau_b_updn(:,:)=zero
764 :
765 1278098477 : do ipts=ifft,ifft+npts-1
766 : ! indx=ipts-ifft+1 varies from 1 to npts
767 1277759538 : indx=ipts-ifft+1
768 1277759538 : rho_b(indx)=rhonow(ipts,1,1)
769 1278098477 : if(nspden_updn==1)then
770 1073065614 : rho_b_updn(indx,1)=rhonow(ipts,1,1)*half
771 1073065614 : if (usegradient==1) grho2_b_updn(indx,1)=quarter*(rhonow(ipts,1,2)**2 &
772 298418604 : & +rhonow(ipts,1,3)**2+rhonow(ipts,1,4)**2)
773 1073065614 : if (usekden==1) tau_b_updn(indx,1)=taunow(ipts,1,1)*half
774 1073065614 : 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 338939 : 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 338939 : 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 337915 : & hyb_range=hyb_range)
829 : end if
830 :
831 : ! Call to main XC driver
832 338939 : 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 337915 : & 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 338939 : 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 338939 : 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 338939 : dstrsxc=zero
887 1278098477 : do ipts=ifft,ifft+npts-1
888 1277759538 : indx=ipts-ifft+1
889 1277759538 : rhoexc=rhoexc+rho_b(indx)*exc_b(indx) ! Will be normalized with respect to the volume later to get bigexc.
890 1277759538 : rhotsxc=rhotsxc+rho_b(indx)*tsxc_b(indx) ! Will be normalized with respect to the volume later to get bigsxc.
891 1277759538 : depsxc(ipts,1)=vxcrho_b_updn(indx,1)
892 1277759538 : if (nspden_updn==2) depsxc(ipts,2)=vxcrho_b_updn(indx,2)
893 1277759538 : 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 1263914114 : if(nspden_updn==1)then
899 1061231863 : strdiag=rho_b(indx)*(exc_b(indx)-vxcrho_b_updn(indx,1))
900 1061231863 : 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 1263914114 : dstrsxc=dstrsxc+strdiag
911 : end if
912 :
913 : ! For GGAs, additional terms appear
914 : ! (the LB functional does not lead to additional terms)
915 1277759538 : 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 865523232 : do ispden=1,3
920 :
921 798351726 : 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 499933122 : if(grho2_b_updn(indx,ispden) < 1.0d-24) then
926 4694640 : depsxc(ipts,ispden+nspden_updn)=zero
927 4694640 : 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 495238482 : if(nspden_updn==1)then
935 : ! ! Definition of vxcgrho_b_updn changed in v3.3
936 298307080 : if (nvxcgrho == 3) then
937 294537498 : 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 495238482 : 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 495238482 : 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 857842864 : if (compute_stress) then
959 1371003032 : if (nspden_updn==1.and.ispden==1) grho(1:3)=rhonow(ipts,1,2:4)
960 686698136 : if (nspden_updn==2.and.ispden==1) grho(1:3)=rhonow(ipts,2,2:4)
961 673046975 : if (nspden_updn==2.and.ispden==2) grho(1:3)=rhonow(ipts,1,2:4)-rhonow(ipts,2,2:4)
962 686695691 : if (nspden_updn==2.and.ispden==3) grho(1:3)=rhonow(ipts,1,2:4)
963 488206262 : s1=-grho(1)*grho(1)*coeff
964 488206262 : s2=-grho(2)*grho(2)*coeff
965 488206262 : 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 488206262 : dstrsxc=dstrsxc+s1+s2+s3
969 : !Non-diagonal part: -Int[dn/dr_alpha.dFxc/grad_beta(n)]
970 488206262 : strsxc1_tot=strsxc1_tot+s1
971 488206262 : strsxc2_tot=strsxc2_tot+s2
972 488206262 : strsxc3_tot=strsxc3_tot+s3
973 488206262 : strsxc4_tot=strsxc4_tot-grho(3)*grho(2)*coeff
974 488206262 : strsxc5_tot=strsxc5_tot-grho(3)*grho(1)*coeff
975 488206262 : 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 1277759538 : 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 1278098477 : 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 338939 : 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 338939 : if (nkxc_eff>0.and.ndvxc>0 .and. (xcdata%usefock==0 .or. auxc_ixc==0)) then
1107 127901 : if (nkxc_eff==1.and.ndvxc==15) then
1108 92158830 : kxc(ifft:ifft+npts-1,1)=half*(dvxc_b(1:npts,1)+dvxc_b(1:npts,9)+dvxc_b(1:npts,10))
1109 109422 : 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 101281 : 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 706 : else if (nkxc_eff==7.and.ndvxc==15) then
1119 2643962 : kxc(ifft:ifft+npts-1,1)=half*(dvxc_b(1:npts,1)+dvxc_b(1:npts,9)+dvxc_b(1:npts,10))
1120 2643962 : kxc(ifft:ifft+npts-1,2)=half*dvxc_b(1:npts,3)+dvxc_b(1:npts,12)
1121 2643962 : kxc(ifft:ifft+npts-1,3)=quarter*dvxc_b(1:npts,5)+dvxc_b(1:npts,13)
1122 2643962 : 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 127901 : if (nkxc_eff==7) then
1142 2653565 : kxc(ifft:ifft+npts-1,5)=rhonow(ifft:ifft+npts-1,1,2)
1143 2653565 : kxc(ifft:ifft+npts-1,6)=rhonow(ifft:ifft+npts-1,1,3)
1144 2653565 : kxc(ifft:ifft+npts-1,7)=rhonow(ifft:ifft+npts-1,1,4)
1145 127192 : 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 338939 : 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 338939 : if (compute_stress) then
1162 334935 : strsxc1_tot=strsxc1_tot+dstrsxc
1163 334935 : strsxc2_tot=strsxc2_tot+dstrsxc
1164 334935 : strsxc3_tot=strsxc3_tot+dstrsxc
1165 : end if
1166 :
1167 : ! Accumulate integral of |Grad_rho|/Rho (to be used for TB09 XC)
1168 338939 : if (present(grho1_over_rho1).and.ixc<0) then
1169 14621 : if (present(xc_funcs)) then
1170 0 : test_tb09=libxc_functionals_is_tb09(xc_functionals=xc_funcs)
1171 : else
1172 14621 : test_tb09=libxc_functionals_is_tb09()
1173 : end if
1174 14621 : 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 338939 : ABI_FREE(exc_b)
1187 338939 : ABI_FREE(tsxc_b)
1188 338939 : ABI_FREE(rho_b)
1189 338939 : ABI_FREE(rho_b_updn)
1190 338939 : ABI_FREE(grho2_b_updn)
1191 338939 : ABI_FREE(vxcrho_b_updn)
1192 338939 : ABI_FREE(dvxc_b)
1193 338939 : ABI_FREE(d2vxc_b)
1194 338939 : ABI_FREE(vxcgrho_b_updn)
1195 338939 : ABI_FREE(vxclrho_b_updn)
1196 338939 : ABI_FREE(lrho_b_updn)
1197 338939 : ABI_FREE(tau_b_updn)
1198 338939 : ABI_FREE(vxctau_b_updn)
1199 :
1200 : ! End of the loop on blocks of data
1201 : end do
1202 :
1203 59523 : if (compute_stress) then
1204 58379 : strsxc(1)=strsxc1_tot
1205 58379 : strsxc(2)=strsxc2_tot
1206 58379 : strsxc(3)=strsxc3_tot
1207 58379 : strsxc(4)=strsxc4_tot
1208 58379 : strsxc(5)=strsxc5_tot
1209 58379 : 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 59523 : rhonow_ptr => rhonow
1215 59523 : if (ipositron==2) then
1216 1440 : ABI_MALLOC(rhonow_ptr,(nfft,nspden_eff,ngrad*ngrad))
1217 86842980 : rhonow_ptr=rhonow
1218 : end if
1219 59523 : if(ngrad==2 .and. ixc/=13)then
1220 13603 : 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 59523 : if (nspden/=4) then
1225 56940 : 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 56274 : & 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 59523 : if (ipositron==2) then
1256 288 : ABI_FREE(rhonow_ptr)
1257 : end if
1258 59523 : nullify(rhonow_ptr)
1259 :
1260 : ! Add electron-positron XC potential to electron-electron one
1261 : ! Eventually compute GGA contribution
1262 55320 : 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 55320 : 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 55320 : divshft=one/dble(xcdata%intxc+1)
1305 55320 : bigexc=rhoexc*ucvol/dble(nfftot)*divshft
1306 55320 : bigsxc=rhotsxc*ucvol/dble(nfftot)*divshft/xcdata%tphysel
1307 1494293323 : vxc=vxc*divshft
1308 380694 : if (compute_stress) strsxc(:)=strsxc(:)/dble(nfftot)*divshft
1309 151203820 : if (with_vxctau) vxctau=vxctau*divshft
1310 55320 : if (present(grho1_over_rho1)) grho1_over_rho1=grho1_over_rho1*ucvol/dble(nfftot)*divshft
1311 :
1312 : ! Reduction in case of FFT distribution
1313 55320 : if (nproc_fft>1)then
1314 2820 : call timab(48,1,tsec)
1315 2820 : call xmpi_sum(bigexc,comm_fft,ierr)
1316 2820 : call xmpi_sum(bigsxc,comm_fft,ierr)
1317 2820 : if (compute_stress) then
1318 2820 : call xmpi_sum(strsxc,comm_fft,ierr)
1319 : end if
1320 2820 : if (present(grho1_over_rho1)) then
1321 440 : call xmpi_sum(grho1_over_rho1,comm_fft ,ierr)
1322 : end if
1323 2820 : 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 2820 : call timab(48,2,tsec)
1330 : end if
1331 :
1332 : ! Compute vxcavg
1333 55320 : call mean_fftr(vxc,vxcmean,nfft,nfftot,min(nspden,2),mpi_comm_sphgrid=comm_fft)
1334 55320 : if(nspden==1)then
1335 46821 : vxcavg=vxcmean(1)
1336 : else
1337 8499 : vxcavg=half*(vxcmean(1)+vxcmean(2))
1338 : end if
1339 :
1340 55320 : ABI_FREE(depsxc)
1341 55320 : ABI_FREE(rhonow)
1342 55320 : ABI_FREE(lrhonow)
1343 55320 : ABI_FREE(taunow)
1344 55320 : if (compute_stress.and.uselaplacian>0) then
1345 135 : ABI_FREE(d2rhonow)
1346 : end if
1347 55320 : if (need_nhat.or.non_magnetic_xc) then
1348 10430 : ABI_FREE(rhor_)
1349 : end if
1350 55320 : if ((usekden==1).and.(non_magnetic_xc)) then
1351 0 : ABI_FREE(taur_)
1352 : end if
1353 55320 : 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 56097 : 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 56097 : 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 56097 : if(abs(bigsxc)>tiny(zero)) bigexc=bigexc+xcdata%tphysel*bigsxc
1415 :
1416 56097 : call timab(81,2,tsec)
1417 :
1418 : DBG_EXIT("COLL")
1419 :
1420 166782 : end subroutine rhotoxc
1421 : !!***
1422 :
1423 : end module m_rhotoxc
1424 : !!***
|
