TABLE OF CONTENTS
ABINIT/m_rhotoxc [ Modules ]
NAME
m_rhotox
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, MF, GZ, DRH, MT, SPr) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
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
ABINIT/rhotoxc [ Functions ]
NAME
rhotoxc
FUNCTION
Start from the density or spin-density, and compute xc correlation potential and energies. Eventually compute xc kernel (if option=-2, 2, 3, 10 or 12). Cannot be used with wavelets.
INPUTS
mpi_enreg=information about MPI parallelization nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nhat(nfft,xcdata%nspden*nhatdim)= -PAW only- compensation density nhatdim= -PAW only- 0 if nhat array is not used ; 1 otherwise nhatgr(nfft,xcdata%nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise nkxc=second dimension of the kxc array. If /=0, the exchange-correlation kernel must be computed. non_magnetic_xc= if true, handle density/potential as non-magnetic (even if it is) n3xccc=dimension of the xccc3d array (0 or nfft or cplx*nfft). option=0 or 1 for xc only (exc, vxc, strsxc), 2 for xc and kxc (no paramagnetic part if xcdata%nspden=1) 10 for xc and kxc with only partial derivatives wrt density part (d2Exc/drho^2) 12 for xc and kxc with only partial derivatives wrt density part (d2Exc/drho^2) and, in the case of hybrid functionals, substitution of the hybrid functional by the related auxiliary GGA functional for the computation of the xc kernel (not for other quantities) 3 for xc, kxc and k3xc -2 for xc and kxc (with paramagnetic part if xcdata%nspden=1) rhor(nfft,xcdata%nspden)=electron density in real space in electrons/bohr**3 (total in first half and spin-up in second half if xcdata%nspden=2) (total in first comp. and magnetization in comp. 2 to 4 if xcdata%nspden=4) rprimd(3,3)=dimensional primitive translations in real space (bohr) usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc [vhartr(nfft)=Hartree potential (only needed for Fermi-Amaldi functional)] xcdata <type(xcdata_type)>=storage for different input variables and derived parameters needed to compute the XC functional xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3) === optional inputs === [add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy [taur(nfftf,xcdata%nspden*xcdata%usekden)]=array for kinetic energy density [xc_funcs(2)]= <type(libxc_functional_type)>, optional : libxc XC functionals. Must be coherent with xcdata. [xccctau3d(n3xccc)]=3D core electron kinetic energy density for XC core correction (bohr^-3)
OUTPUT
enxc=returned exchange and correlation energy (hartree). strsxc(6)= contribution of xc to stress tensor (hartree/bohr^3), given in order (1,1), (2,2), (3,3), (3,2), (3,1), (2,1). (note: fxc is rho*exc in the following) Explicitely : strsxc(mu,nu) = (1/N) Sum(i=1,N) ( delta(mu,nu) * [ exc(i)rhotot(i) - depsxc_drho(up,i)*rhor(up,i)-depsxc_drho(dn,i)*rhor(dn,i)] - gradrho(up,mu)*gradrho(up,nu) * depsxc_dgradrho(up,i) / gradrho(up,i) - gradrho(dn,mu)*gradrho(dn,nu) * depsxc_dgradrho(dn,i) / gradrho(dn,i) ) vxc(nfft,xcdata%nspden)=xc potential (spin up in first half and spin down in second half if xcdata%nspden=2) (v^11, v^22, Re[V^12], Im[V^12] if xcdata%nspden=4) vxcavg=<Vxc>=unit cell average of Vxc = (1/ucvol) Int [Vxc(r) d^3 r]. === Only if abs(option)=2, -2, 3, 10, 12 (in case 12, for hybrids, substitution of the related GGA) === kxc(nfft,nkxc)=exchange and correlation kernel (returned only if nkxc/=0) Content of Kxc array: ===== if LDA if xcdata%nspden==1: kxc(:,1)= d2Exc/drho2 that is 1/2 ( d2Exc/drho_up drho_up + d2Exc/drho_up drho_dn ) kxc(:,2)= d2Exc/drho_up drho_dn if xcdata%nspden>=2: kxc(:,1)= d2Exc/drho_up drho_up kxc(:,2)= d2Exc/drho_up drho_dn kxc(:,3)= d2Exc/drho_dn drho_dn ===== if GGA or mGGA if xcdata%nspden==1: kxc(:,1)= d2Exc/drho2 kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)| kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| ) kxc(:,5)= gradx(rho) kxc(:,6)= grady(rho) kxc(:,7)= gradz(rho) if xcdata%nspden>=2: kxc(:,1)= d2Exc/drho_up drho_up kxc(:,2)= d2Exc/drho_up drho_dn kxc(:,3)= d2Exc/drho_dn drho_dn kxc(:,4)= 1/|grad(rho_up)| dEx/d|grad(rho_up)| kxc(:,5)= 1/|grad(rho_dn)| dEx/d|grad(rho_dn)| kxc(:,6)= 1/|grad(rho_up)| d2Ex/d|grad(rho_up)| drho_up kxc(:,7)= 1/|grad(rho_dn)| d2Ex/d|grad(rho_dn)| drho_dn kxc(:,8)= 1/|grad(rho_up)| * d/d|grad(rho_up)| ( 1/|grad(rho_up)| dEx/d|grad(rho_up)| ) kxc(:,9)= 1/|grad(rho_dn)| * d/d|grad(rho_dn)| ( 1/|grad(rho_dn)| dEx/d|grad(rho_dn)| ) kxc(:,10)=1/|grad(rho)| dEc/d|grad(rho)| kxc(:,11)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_up kxc(:,12)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_dn kxc(:,13)=1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dEc/d|grad(rho)| ) kxc(:,14)=gradx(rho_up) kxc(:,15)=gradx(rho_dn) kxc(:,16)=grady(rho_up) kxc(:,17)=grady(rho_dn) kxc(:,18)=gradz(rho_up) kxc(:,19)=gradz(rho_dn) Note about mGGA: 2nd derivatives involving Tau or Laplacian are not output === Only if abs(option)=3 === [k3xc(nfft,nk3xc)]= -- optional -- third derivative of the XC energy functional of the density, at each point of the real space grid (only in the LDA or LSDA) Content of K3xc array: ===== if LDA if xcdata%nspden==1: k3xc(:,1)= d3Exc/drho3 if xcdata%nspden>=2, k3xc(:,1)= d3Exc/drho_up drho_up drho_up k3xc(:,2)= d3Exc/drho_up drho_up drho_dn k3xc(:,3)= d3Exc/drho_up drho_dn drho_dn k3xc(:,4)= d3Exc/drho_dn drho_dn drho_dn === Additional optional output === [exc_vdw_out]= vdW-DF contribution to enxc (hartree) [vxctau(nfft,xcdata%nspden,4*xcdata%usekden)]=(only for meta-GGA)= vxctau(:,:,1): derivative of XC energy density with respect to kinetic energy density (depsxcdtau). vxctau(:,:,2:4): gradient of vxctau (gvxctau) === For the TB09 XC functional (modified Becke-Johnson) === [grho1_over_rho1]=Integral of |Grad(rho^1)|/rho^1 over the augmentation region Used to compute the c parameter of the TB09 XC functional
SIDE EFFECTS
electronpositron <type(electronpositron_type)>= -- optional argument -- quantities for the electron-positron annihilation
NOTES
Start from the density, and compute Hartree (if option>=1) and xc correlation potential and energies. Eventually compute xc kernel (if option=-2, 2, 3, 10 or 12 - in the latter case, substitution by the related GGA kernel). Allows a variety of exchange-correlation functionals according to ixc. Here is a list of allowed values. subroutine name <0 means use of libxc 0 means no xc applied (usually for testing) *LDA,LSD 1 means new Teter (4/93) with spin-pol option xcspol 2 means Perdew-Zunger-Ceperley-Alder xcpzca 3 means old Teter (4/91) fit to Ceperley-Alder data xctetr 4 means Wigner xcwign 5 means Hedin-Lundqvist xchelu 6 means "X-alpha" xc xcxalp 7 mean Perdew-Wang 92 LSD fit to Ceperley-Alder data xcpbe 8 mean Perdew-Wang 92 LSD , exchange-only xcpbe 9 mean Perdew-Wang 92 Ex+Ec_RPA energy xcpbe 10 means RPA LSD energy (only the energy !!) xcpbe *GGA 11 means Perdew-Burke-Ernzerhof GGA functional xcpbe 12 means x-only Perdew-Burke-Ernzerhof GGA functional xcpbe 13 means LDA (ixc==7), except that the xc potential is given within the van Leeuwen-Baerends GGA xclb 14 means revPBE GGA functional xcpbe 15 means RPBE GGA functional xcpbe 16 means HCTH GGA functional xchcth 23 means WC GGA functional xcpbe 24 means C09x GGA exchange functional xcpbe *Fermi-Amaldi 20 means Fermi-Amaldi correction 21 means Fermi-Amaldi correction with LDA(ixc=1) kernel 22 means Fermi-Amaldi correction with hybrid BPG kernel *Hybrid GGA 41 means PBE0-1/4 xcpbe 42 means PBE0-1/3 xcpbe *Other 50 means IIT xc xciit NOTE: please update echo_xc_name.F90 if you add new functional (apart from libxc) Allow for improved xc quadrature (intxc=1) by using the usual FFT grid as well as another, shifted, grid, and combining both results. Spin-polarization is allowed only with ixc=0, 1, and GGAs until now. To make the variable names easier to understand, a rule notation is tentatively proposed here: rho ---> means density tau ---> means kinetic energy density exc ---> means exchange-correlation energy density per particle epsxc ---> means rho*exc == exchange-correlation energy density vxc ---> means exchange-correlation potential bigexc ---> means exchange-correlation energy E_xc (for the moment it is named "enxc") m_norm ---> means norm of magnetization g... --> means gradient of something (e.g. : grho --> means gradient of electron density) g...2 -> means square norm of gradient of something (e.g. : grho2 -> means square norm of gradient of electron density) l... --> means laplacian of something (e.g. : lrho --> means laplacian of electron density) d...d... --> means derivative of something with regards to something else. (d2...d...d... ---> means second derivative of ... with regards to ... and to ...) etc... d... --> without the occurence of the second "d" means that this is an array of several derivative of the same quantity (e.g. : depsxc) ..._b ----> means a block of the quantity "..." (use in mpi loops which treat the data block by block) ..._updn -> means that spin up and spin down is available in that array as (..,1) and (..,2). (if xcdata%nspden >=2 of course). ..._apn --> in case of positrons are concerned. for more details about notations please see pdf in /doc/theory/MGGA/
SOURCE
247 subroutine rhotoxc(enxc,kxc,mpi_enreg,nfft,ngfft, & 248 & nhat,nhatdim,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,option, & 249 & rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata, & 250 & add_tfw,exc_vdw_out,grho1_over_rho1,electronpositron,k3xc,taur,vhartr,vxctau,xc_funcs,xcctau3d) ! optional arguments 251 252 !Arguments ------------------------------------ 253 !scalars 254 integer,intent(in) :: nk3xc,n3xccc,nfft,nhatdim,nhatgrdim,nkxc,option 255 integer,intent(in) :: usexcnhat 256 logical,intent(in) :: non_magnetic_xc 257 logical,intent(in),optional :: add_tfw 258 real(dp),intent(out) :: enxc,vxcavg 259 real(dp),intent(out),optional :: exc_vdw_out,grho1_over_rho1 260 type(MPI_type),intent(in) :: mpi_enreg 261 type(electronpositron_type),pointer,optional :: electronpositron 262 type(xcdata_type), intent(in) :: xcdata 263 !arrays 264 integer,intent(in) :: ngfft(18) 265 real(dp),intent(in) :: nhat(nfft,xcdata%nspden*nhatdim) 266 real(dp),intent(in) :: nhatgr(nfft,xcdata%nspden,3*nhatgrdim) 267 real(dp),intent(in),target :: rhor(nfft,xcdata%nspden) 268 real(dp),intent(in) :: rprimd(3,3),xccc3d(n3xccc) 269 real(dp),intent(in),optional :: xcctau3d(:) 270 real(dp),intent(out) :: kxc(nfft,nkxc),strsxc(6),vxc(nfft,xcdata%nspden) 271 real(dp),intent(in),optional :: vhartr(nfft) 272 real(dp),intent(in),target,optional :: taur(:,:) 273 real(dp),intent(out),optional :: k3xc(1:nfft,1:nk3xc),vxctau(:,:,:) 274 type(libxc_functional_type),intent(inout),optional :: xc_funcs(2) 275 276 !Local variables------------------------------- 277 !scalars 278 integer :: auxc_ixc,cplex,ierr,ifft,ii,ixc,ixc_from_lib,indx,ipositron,ipts,ishift,ispden,iwarn,iwarnp 279 integer :: jj,mpts,ndvxc,nd2vxc,nfftot,ngr,ngrad,ngrad_apn,nkxc_eff,npts 280 integer :: nspden,nspden_apn,nspden_eff,nspden_updn,nspgrad,nvxcgrho,nvxclrho,nvxctau 281 integer :: n3xctau,order,usefxc,nproc_fft,comm_fft,usegradient,usekden,uselaplacian 282 logical :: my_add_tfw 283 real(dp),parameter :: mot=-one/3.0_dp 284 real(dp) :: coeff,divshft,doti,dstrsxc,dvdn,dvdz,epsxc,exc_str,factor,m_norm_min,s1,s2,s3 285 real(dp) :: strdiag,strsxc1_tot,strsxc2_tot,strsxc3_tot,strsxc4_tot 286 real(dp) :: strsxc5_tot,strsxc6_tot,ucvol 287 real(dp) :: deltae_vdw,exc_vdw 288 logical :: test_nhat,test_tb09,need_nhat,need_nhatgr,with_vxctau 289 character(len=500) :: message 290 real(dp) :: hyb_mixing, hyb_mixing_sr, hyb_range 291 !arrays 292 real(dp) :: gm_norm(3),grho(3),gmet(3,3),gprimd(3,3),qphon(3),rmet(3,3) 293 real(dp) :: tsec(2),vxcmean(4) 294 real(dp),allocatable :: d2vxc_b(:,:),depsxc(:,:),depsxc_apn(:,:),dvxc_apn(:),dvxc_b(:,:) 295 real(dp),allocatable :: exc_b(:),fxc_b(:),fxc_apn(:),grho2_apn(:),grho2_b_updn(:,:),lrhonow(:,:),lrho_b_updn(:,:) 296 real(dp),allocatable :: m_norm(:),nhat_up(:),rho_b_updn(:,:),rho_b(:),rhonow_apn(:,:,:) 297 real(dp),allocatable :: tau_b_updn(:,:),vxc_apn(:,:),vxcgr_apn(:),vxcgrho_b_updn(:,:),vxcrho_b_updn(:,:) 298 real(dp),allocatable :: vxc_b_apn(:),vxc_ep(:),vxctau_b_updn(:,:),vxclrho_b_updn(:,:) 299 real(dp),allocatable,target :: rhonow(:,:,:),taunow(:,:,:) 300 real(dp),pointer :: rhocorval(:,:),rhor_(:,:),taucorval(:,:),taur_(:,:) 301 real(dp),ABI_CONTIGUOUS pointer :: rhonow_ptr(:,:,:) 302 real(dp) :: decdrho_vdw(nfft,xcdata%nspden),decdgrho_vdw(nfft,3,xcdata%nspden) 303 real(dp) :: strsxc_vdw(3,3) 304 type(libxc_functional_type) :: xc_funcs_auxc(2) 305 306 ! ************************************************************************* 307 308 ! Note: the following cases seem to never be tested (should be fixed) 309 ! - ipositron==2 and ngrad_apn==2 310 ! - usewvl/=0 311 ! - test_nhat and usexcnhat==1 and nspden==4 312 313 call timab(81,1,tsec) 314 315 !Optional arguments 316 my_add_tfw=.false.;if (present(add_tfw)) my_add_tfw=add_tfw 317 318 !Useful scalars 319 nspden=xcdata%nspden 320 ixc=xcdata%ixc 321 auxc_ixc=xcdata%auxc_ixc 322 n3xctau=0 323 324 !nspden_updn: 1 for non-polarized, 2 for polarized 325 nspden_updn=min(nspden,2) 326 327 !The variable order indicates to which derivative of the energy 328 !the computation must be done. Computing exc and vxc needs order=1 . 329 !Meaningful values are 1, 2, 3. Lower than 1 is the same as 1, and larger 330 !than 3 is the same as 3. 331 !order=1 or 2 supported for all LSD and GGA ixc 332 !order=3 supported only for ixc=3 and ixc=7 333 order=1 334 if(option==2.or.option==10.or.option==12)order=2 335 if(option==-2)order=-2 336 if(option==3)order=3 337 338 !Sizes of local arrays 339 if (present(xc_funcs)) then 340 call size_dvxc(ixc,order,nspden_updn,& 341 & usegradient=usegradient,uselaplacian=uselaplacian,usekden=usekden,& 342 & nvxcgrho=nvxcgrho,nvxclrho=nvxclrho,nvxctau=nvxctau,& 343 & ndvxc=ndvxc,nd2vxc=nd2vxc,add_tfw=my_add_tfw,xc_funcs=xc_funcs) 344 else 345 call size_dvxc(ixc,order,nspden_updn,& 346 & usegradient=usegradient,uselaplacian=uselaplacian,usekden=usekden,& 347 & nvxcgrho=nvxcgrho,nvxclrho=nvxclrho,nvxctau=nvxctau,& 348 & ndvxc=ndvxc,nd2vxc=nd2vxc,add_tfw=my_add_tfw) 349 end if 350 351 !ngrad=1 is for LDAs or LSDs, ngrad=2 is for GGAs/mGGAs 352 ngrad=1;if(xcdata%xclevel==2.or.usegradient==1) ngrad=2 353 354 !nspden_eff: effective value of nspden used to compute gradients of density: 355 ! 1 for non-polarized system, 356 ! 2 for collinear polarized system or LDA (can be reduced to a collinear system) 357 ! 4 for non-collinear polarized system and GGA 358 nspden_eff=nspden_updn;if (nspden==4.and.ngrad==2) nspden_eff=4 359 360 !Number of kcxc components depends on option (force LDA type if option==10 or 12) 361 nkxc_eff=nkxc;if (option==10.or.option==12) nkxc_eff=min(nkxc,3) 362 363 !Check options 364 if(option==3.and.nd2vxc==0.and.ixc/=0)then 365 write(message, '(3a,i0)' )& 366 & 'Third-order xc kernel can only be computed for ixc = 0, 3, 7 or 8,',ch10,& 367 & 'while it is found to be ',ixc 368 ABI_ERROR(message) 369 end if 370 if(nspden==4.and.xcdata%xclevel==2.and.(abs(option)==2))then 371 ABI_BUG('When nspden==4 and GGA, the absolute value of option cannot be 2 !') 372 end if 373 if(ixc<0) then 374 if (present(xc_funcs)) then 375 ixc_from_lib=libxc_functionals_ixc(xc_functionals=xc_funcs) 376 else 377 ixc_from_lib=libxc_functionals_ixc() 378 end if 379 ! Check consistency between ixc passed in input and the one used to initialize the library. 380 if (ixc /= ixc_from_lib) then 381 write(message, '(a,i0,2a,i0,2a)')& 382 & 'The value of ixc specified in input, ixc = ',ixc,ch10,& 383 & 'differs from the one used to initialize the functional ',ixc_from_lib,ch10,& 384 & 'Action: reinitialize the global structure funcs, see NOTES in m_libxc_functionals' 385 ABI_BUG(message) 386 end if 387 end if 388 389 !Handling of mGGA functionals 390 with_vxctau=(present(vxctau)) 391 if (with_vxctau) with_vxctau=(size(vxctau)>0) 392 if (usekden==1) then 393 if (.not.present(taur)) then 394 message=' For mGGA functionals, kinetic energy density is needed. Set input variable usekden to 1.' 395 message=trim(message)//' Also use NC pseudopotentials without non-linear XC core correction.' 396 ABI_BUG(message) 397 else if (size(taur)/=nfft*nspden) then 398 ABI_BUG('Invalid size for taur!') 399 end if 400 if (present(xcctau3d)) then 401 n3xctau=size(xcctau3d) 402 if (n3xctau/=0.and.n3xctau/=nfft) then 403 ABI_BUG('Invalid size for xccctau3d!') 404 end if 405 end if 406 if (with_vxctau) then 407 if (size(vxctau)/=nfft*nspden*4) then 408 ABI_BUG('Invalid size for vxctau!') 409 end if 410 end if 411 end if 412 if((usekden==1.or.uselaplacian==1).and.nspden==4)then 413 !mGGA en NC-magnetism: how do we rotate tau kinetic energy density? 414 message=' At present, meta-GGA (usekden=1 or uselaplacian=1) is not compatible with non-collinear magnetism (nspden=4).' 415 ABI_ERROR(message) 416 end if 417 418 !MPI FFT communicator 419 comm_fft = mpi_enreg%comm_fft; nproc_fft = mpi_enreg%nproc_fft 420 421 !Compute different geometric tensor, as well as ucvol, from rprimd 422 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 423 424 !In this routine, hartre, xcden and xcpot are called for real 425 !densities and potentials, corresponding to zero wavevector 426 cplex=1 427 qphon(:)=zero 428 iwarn=0 429 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 430 usefxc=0;if (ixc==50) usefxc=1 431 432 !Initializations 433 enxc=zero 434 epsxc=zero 435 vxc(:,:)=zero 436 vxcavg=zero 437 strsxc(:)=zero 438 strsxc1_tot=zero;strsxc2_tot=zero;strsxc3_tot=zero 439 strsxc4_tot=zero;strsxc5_tot=zero;strsxc6_tot=zero 440 if (with_vxctau) vxctau(:,:,:)=zero 441 if (nkxc/=0) kxc(:,:)=zero 442 if(abs(option)==3.and.nk3xc/=0) k3xc(:,:)=zero 443 ipositron=0 444 if (present(electronpositron)) then 445 ipositron=electronpositron_calctype(electronpositron) 446 if (ipositron==2) then 447 electronpositron%e_xc =zero 448 electronpositron%e_xcdc=zero 449 end if 450 end if 451 deltae_vdw = zero 452 exc_vdw = zero 453 decdrho_vdw(:,:) = zero 454 decdgrho_vdw(:,:,:) = zero 455 strsxc_vdw(:,:) = zero 456 if (present(grho1_over_rho1)) grho1_over_rho1=zero 457 458 if ((xcdata%xclevel==0.or.ixc==0).and.(.not.my_add_tfw)) then 459 ! No xc at all is applied (usually for testing) 460 ABI_WARNING('Note that no xc is applied (ixc=0).') 461 462 else if (ixc/=20) then 463 464 ! Test: has a compensation density to be added/substracted (PAW) ? 465 need_nhat=(nhatdim==1.and.usexcnhat==0) 466 need_nhatgr=(nhatdim==1.and.nhatgrdim==1.and.ngrad==2.and.xcdata%intxc==0) 467 test_nhat=(need_nhat.or.need_nhatgr) 468 469 ! The different components of depsxc will be 470 ! for nspden=1, depsxc(:,1)=d(rho.exc)/d(rho) == (depsxcdrho) == (vxcrho) 471 ! and if ngrad=2, depsxc(:,2)=1/2*1/|grad rho_up|*d(rho.exc)/d(|grad rho_up|) 472 ! +1/|grad rho|*d(rho.exc)/d(|grad rho|) 473 ! == (1/2 * 1/|grho_up| * depsxcd|grho_up|) + 1/|grho| * depsxcd|grho| 474 ! (vxcgrho=1/|grho| * depsxcd|grho|) 475 ! (do not forget : |grad rho| /= |grad rho_up| + |grad rho_down| 476 ! and if use_laplacian, depsxc(:,3)=d(rho.exc)/d(lapl rho) == (depsxcdlrho) == (vxclrho) 477 ! 478 ! for nspden>=2, depsxc(:,1)=d(rho.exc)/d(rho_up) == (depsxcdrho_up) == (vxcrho_up) 479 ! depsxc(:,2)=d(rho.exc)/d(rho_down) == (depsxcdrho_dn) == (vxcrho_dn) 480 ! 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) 481 ! depsxc(:,4)=1/|grad rho_down|*d(rho.exc)/d(|grad rho_down|) == (1/|grho_dn| * depsxcd|grho_dn|) == (vxcgrho_dn) 482 ! depsxc(:,5)=1/|grad rho|*d(rho.exc)/d(|grad rho|) == (1/|grho| * depsxcd|grho|) == (vxcgrho) 483 ! and if use_laplacian, depsxc(:,6)=d(rho.exc)/d(lapl rho_up) == (depsxcdlrho_up) == (vxclrho_up) 484 ! depsxc(:,7)=d(rho.exc)/d(lapl rho_dn) == (depsxcdlrho_dn) == (vxclrho_dn) 485 ! Note: if nspden=4, rho_up=(rho+|m|)/2, rho_down=(rho-|m|)/2 486 nspgrad=nspden_updn*ngrad;if(nspden_updn==2.and.ngrad==2)nspgrad=5 487 if(uselaplacian==1) nspgrad=nspgrad+nspden_updn 488 ABI_MALLOC(depsxc,(nfft,nspgrad)) 489 depsxc(:,:)=zero 490 491 ! PAW: select the valence density (and magnetization) to use: 492 ! link the correct density, according to usexcnhat option 493 if ((.not.need_nhat).and.(.not.non_magnetic_xc)) then 494 rhor_ => rhor 495 else 496 ABI_MALLOC(rhor_,(nfft,nspden)) 497 if (need_nhat) then 498 do ispden=1,nspden 499 do ifft=1,nfft 500 rhor_(ifft,ispden)=rhor(ifft,ispden)-nhat(ifft,ispden) 501 end do 502 end do 503 else 504 do ispden=1,nspden 505 do ifft=1,nfft 506 rhor_(ifft,ispden)=rhor(ifft,ispden) 507 end do 508 end do 509 end if 510 if(non_magnetic_xc) then 511 if(nspden==2) rhor_(:,2)=rhor_(:,1)*half 512 if(nspden==4) rhor_(:,2:4)=zero 513 endif 514 end if 515 if (usekden==1) then 516 if(non_magnetic_xc) then 517 ABI_MALLOC(taur_,(nfft,nspden)) 518 if(nspden==2) taur_(:,2)=taur_(:,1)*half 519 if(nspden==4) taur_(:,2:4)=zero 520 else 521 taur_ => taur 522 end if 523 end if 524 525 ! Some initializations for the electron-positron correlation 526 if (ipositron==2) then 527 nspden_apn=1;ngrad_apn=1;iwarnp=1 528 if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad_apn=2 529 if (ngrad_apn==2.and.xcdata%xclevel<2) then 530 message = 'GGA for the positron can only be performed with GGA pseudopotentials for the electron !' 531 ABI_ERROR(message) 532 end if 533 if (ngrad_apn>1.and.option/=0.and.option/=1.and.option/=10.and.option/=12) then 534 message = 'You cannot compute full GGA XC kernel for electrons-positron systems !' 535 ABI_ERROR(message) 536 end if 537 ABI_MALLOC(depsxc_apn,(nfft,ngrad_apn)) 538 end if 539 540 ! Non-collinear magnetism: store norm of magnetization 541 ! m_norm_min= EPSILON(0.0_dp)**2 ! EB: TOO SMALL!!! 542 m_norm_min=tol8 ! EB: tol14 is still too small, tests are underway 543 if (nspden==4) then 544 ABI_MALLOC(m_norm,(nfft)) 545 m_norm(:)=sqrt(rhor_(:,2)**2+rhor_(:,3)**2+rhor_(:,4)**2) 546 end if 547 548 ! rhocorval will contain effective density used to compute gradients: 549 ! - with core density (if NLCC) 550 ! - without compensation density (if PAW under certain conditions) 551 ! - in (up+dn,up) or (n,mx,my,mz) format according to collinearity 552 ! of polarization and use of gradients (GGA) 553 if (n3xccc>0.or.test_nhat.or.nspden_eff/=nspden) then 554 ABI_MALLOC(rhocorval,(nfft,nspden_eff)) 555 if (nspden==nspden_eff) then 556 rhocorval(:,1:nspden)=rhor_(:,1:nspden) 557 else if (nspden==4) then 558 rhocorval(:,1)=rhor_(:,1) 559 rhocorval(:,2)=half*(rhor_(:,1)+m_norm(:)) 560 else 561 rhocorval=zero 562 end if 563 else 564 rhocorval => rhor_ 565 end if 566 if (usekden==1.and.(n3xctau>0.or.nspden_eff/=nspden)) then 567 ABI_MALLOC(taucorval,(nfft,nspden_eff)) 568 if (nspden==nspden_eff) then 569 taucorval(:,1:nspden)=taur_(:,1:nspden) 570 else 571 taucorval=zero 572 end if 573 else 574 taucorval => taur_ 575 end if 576 577 ! Add core electron density to effective density 578 if (n3xccc>0) then 579 rhocorval(:,1)=rhocorval(:,1)+xccc3d(:) 580 if(nspden_eff==2) then 581 rhocorval(:,2)=rhocorval(:,2)+half*xccc3d(:) 582 end if 583 end if 584 if (n3xctau>0) then 585 taucorval(:,1)=taucorval(:,1)+xcctau3d(:) 586 if(nspden_eff==2) then 587 taucorval(:,2)=taucorval(:,2)+half*xcctau3d(:) 588 end if 589 end if 590 591 ! If PAW, substract compensation density from effective density: 592 ! - if GGA, because nhat gradients are computed separately 593 if (test_nhat.and.usexcnhat==1) then 594 if (nspden==nspden_eff) then 595 rhocorval(:,1:nspden)=rhocorval(:,1:nspden)-nhat(:,1:nspden) 596 else if (nspden==4) then 597 ABI_MALLOC(nhat_up,(nfft)) 598 do ifft=1,nfft 599 if (m_norm(ifft)>m_norm_min) then 600 nhat_up(ifft)=half*(nhat(ifft,1) & 601 & +(rhor_(ifft,2)*nhat(ifft,2) & 602 & +rhor_(ifft,3)*nhat(ifft,3) & 603 & +rhor_(ifft,4)*nhat(ifft,4))/m_norm(ifft)) 604 else 605 nhat_up(ifft)=half*(nhat(ifft,1) & 606 & +sqrt(nhat(ifft,2)**2+nhat(ifft,3)**2+nhat(ifft,4)**2)) 607 end if 608 end do 609 rhocorval(:,1)=rhocorval(:,1)-nhat(:,1) 610 rhocorval(:,2)=rhocorval(:,2)-nhat_up(:) 611 end if 612 end if 613 614 ! rhonow will contain effective density (and gradients if GGA) 615 ! taunow will contain effective kinetic energy density (if MetaGGA) 616 ! lrhonow will contain the laplacian if we have a MGGA 617 ABI_MALLOC(rhonow,(nfft,nspden_eff,ngrad*ngrad)) 618 ABI_MALLOC(lrhonow,(nfft,nspden_eff*uselaplacian)) 619 ABI_MALLOC(taunow,(nfft,nspden_eff,usekden)) 620 621 ! ==================================================================== 622 ! ==================================================================== 623 ! Loop on unshifted or shifted grids 624 do ishift=0,xcdata%intxc 625 626 ! Set up density on unshifted or shifted grid (will be in rhonow(:,:,1)), 627 ! as well as the gradient of the density, also on the unshifted 628 ! or shifted grid (will be in rhonow(:,:,2:4)), if needed. 629 if (uselaplacian==1) then 630 call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_eff,& 631 & qphon,rhocorval,rhonow,lrhonow=lrhonow) 632 else 633 call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_eff,& 634 & qphon,rhocorval,rhonow) 635 end if 636 if (usekden==1) then 637 call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,1,nspden_eff,& 638 & qphon,taucorval,taunow) 639 end if 640 641 ! PAW+GGA: add "exact" gradients of compensation density 642 !if (test_nhat.and.usexcnhat==1.and.ishift==0) then 643 if (test_nhat.and.usexcnhat==1) then 644 if (nspden==nspden_eff) then 645 rhonow(:,1:nspden,1)=rhocorval(:,1:nspden)+nhat(:,1:nspden) 646 else if (nspden==4) then 647 rhonow(:,1,1)=rhocorval(:,1)+nhat(:,1) 648 rhonow(:,2,1)=rhocorval(:,2)+nhat_up(:) 649 end if 650 if (ngrad==2.and.nhatgrdim==1.and.nspden==nspden_eff) then 651 do ii=1,3 652 jj=ii+1 653 do ispden=1,nspden 654 do ifft=1,nfft 655 rhonow(ifft,ispden,jj)=rhonow(ifft,ispden,jj)+nhatgr(ifft,ispden,ii) 656 end do 657 end do 658 end do 659 end if 660 end if 661 662 ! Deallocate temporary arrays 663 if (ishift==xcdata%intxc) then 664 if (n3xccc>0.or.test_nhat.or.nspden_eff/=nspden) then 665 ABI_FREE(rhocorval) 666 end if 667 if (usekden==1.and.(n3xccc>0.or.nspden_eff/=nspden)) then 668 ABI_FREE(taucorval) 669 end if 670 if (test_nhat.and.nspden/=nspden_eff.and.usexcnhat==1) then 671 ABI_FREE(nhat_up) 672 end if 673 end if 674 675 ! In case of non-collinear magnetism, extract up and down density and gradients (if GGA) 676 if (nspden==4.and.nspden_eff==nspden) then 677 if (ngrad==2) then 678 do ifft=1,nfft 679 gm_norm(1:3)=zero 680 if(m_norm(ifft)>m_norm_min) then 681 ! if(m_norm(ifft)>rhonow(ifft,1,1)*tol10+tol14) then 682 do jj=1,3 ! Compute here nabla(|m|)=(m.nabla(m))/|m| == (g|m| = m/|m| * gm) 683 do ii=2,4 684 gm_norm(jj)=gm_norm(jj)+rhonow(ifft,ii,1+jj)*rhonow(ifft,ii,1) 685 end do 686 end do 687 gm_norm(1:3)=gm_norm(1:3)/m_norm(ifft) 688 end if 689 rhonow(ifft,2,2)=half*(rhonow(ifft,1,2)+gm_norm(1)) 690 rhonow(ifft,2,3)=half*(rhonow(ifft,1,3)+gm_norm(2)) 691 rhonow(ifft,2,4)=half*(rhonow(ifft,1,4)+gm_norm(3)) 692 end do 693 end if 694 rhonow(:,2,1)=half*(rhonow(:,1,1)+m_norm(:)) 695 if (usekden==1) taunow(:,2,1)=half*(taunow(:,1,1)+m_norm(:)) 696 end if 697 ! Make the density positive everywhere (but do not care about gradients) 698 call mkdenpos(iwarn,nfft,nspden_updn,1,rhonow(:,1:nspden_updn,1),xcdata%xc_denpos) 699 if (usekden==1) then 700 call mkdenpos(iwarn,nfft,nspden_updn,1,taunow(:,1:nspden_updn,1),xcdata%xc_taupos) 701 end if 702 703 ! write(std_out,*) 'rhonow',rhonow 704 705 ! Uses a block formulation, in order to save simultaneously 706 ! CPU time and memory : xc routines 707 ! are called only once over mpts times, while the amount of allocated 708 ! space is kept at a low value, even if a lot of different 709 ! arrays are allocated, for use in different xc functionals. 710 711 mpts=4000 712 if (usekden==1) mpts=nfft ! Why? 713 714 do ifft=1,nfft,mpts 715 ! npts=mpts 716 ! npts is the number of points to be treated in this bunch 717 npts=min(nfft-ifft+1,mpts) 718 719 ! Allocation of mandatory arguments of drivexc 720 ABI_MALLOC(exc_b,(npts)) 721 ABI_MALLOC(rho_b,(npts)) 722 ABI_MALLOC(rho_b_updn,(npts,nspden_updn)) 723 ABI_MALLOC(vxcrho_b_updn,(npts,nspden_updn)) 724 vxcrho_b_updn(:,:)=zero 725 726 ! Allocation of optional arguments of drivexc 727 ABI_MALLOC(grho2_b_updn,(npts,(2*nspden_updn-1)*usegradient)) 728 ABI_MALLOC(lrho_b_updn,(npts,nspden_updn*uselaplacian)) 729 ABI_MALLOC(tau_b_updn,(npts,nspden_updn*usekden)) 730 ABI_MALLOC(vxcgrho_b_updn,(npts,nvxcgrho)) 731 ABI_MALLOC(vxclrho_b_updn,(npts,nvxclrho)) 732 ABI_MALLOC(vxctau_b_updn,(npts,nvxctau)) 733 ABI_MALLOC(dvxc_b,(npts,ndvxc)) 734 ABI_MALLOC(d2vxc_b,(npts,nd2vxc)) 735 ABI_MALLOC(fxc_b,(npts*usefxc)) 736 if (nvxcgrho>0) vxcgrho_b_updn(:,:)=zero 737 if (nvxclrho>0) vxclrho_b_updn(:,:)=zero 738 if (nvxctau>0) vxctau_b_updn(:,:)=zero 739 740 do ipts=ifft,ifft+npts-1 741 ! indx=ipts-ifft+1 varies from 1 to npts 742 indx=ipts-ifft+1 743 rho_b(indx)=rhonow(ipts,1,1) 744 if(nspden_updn==1)then 745 rho_b_updn(indx,1)=rhonow(ipts,1,1)*half 746 if (usegradient==1) grho2_b_updn(indx,1)=quarter*(rhonow(ipts,1,2)**2 & 747 & +rhonow(ipts,1,3)**2+rhonow(ipts,1,4)**2) 748 if (usekden==1) tau_b_updn(indx,1)=taunow(ipts,1,1)*half 749 if (uselaplacian==1) lrho_b_updn(indx,1)=lrhonow(ipts,1)*half 750 else 751 rho_b_updn(indx,1)=rhonow(ipts,2,1) 752 rho_b_updn(indx,2)=rhonow(ipts,1,1)-rhonow(ipts,2,1) 753 if(usegradient==1)then 754 grho2_b_updn(indx,1)=rhonow(ipts,2,2)**2+ & 755 & rhonow(ipts,2,3)**2+ & 756 & rhonow(ipts,2,4)**2 757 grho2_b_updn(indx,2)=(rhonow(ipts,1,2)-rhonow(ipts,2,2))**2 + & 758 & (rhonow(ipts,1,3)-rhonow(ipts,2,3))**2 + & 759 & (rhonow(ipts,1,4)-rhonow(ipts,2,4))**2 760 grho2_b_updn(indx,3)=rhonow(ipts,1,2)**2+ & 761 & rhonow(ipts,1,3)**2+ & 762 & rhonow(ipts,1,4)**2 763 end if 764 if (usekden==1) then 765 tau_b_updn(indx,1)=taunow(ipts,2,1) 766 tau_b_updn(indx,2)=taunow(ipts,1,1)-taunow(ipts,2,1) 767 end if 768 if (uselaplacian==1) then 769 lrho_b_updn(indx,1)=lrhonow(ipts,2) 770 lrho_b_updn(indx,2)=lrhonow(ipts,1)-lrhonow(ipts,2) 771 end if 772 end if 773 end do 774 ! In case of a hybrid functional, if one needs to compute the auxiliary GGA Kxc, 775 ! a separate call to drivexc is first needed to compute Kxc using such auxiliary GGA, 776 ! before calling again drivexc using the correct functional for Exc and Vxc. 777 778 if(xcdata%usefock==1 .and. auxc_ixc/=0)then 779 if (auxc_ixc<0) then 780 call libxc_functionals_init(auxc_ixc,nspden,xc_functionals=xc_funcs_auxc) 781 end if 782 call drivexc(auxc_ixc,order,npts,nspden_updn,usegradient,0,0,& 783 & rho_b_updn,exc_b,vxcrho_b_updn,nvxcgrho,0,0,ndvxc,nd2vxc, & 784 & grho2_updn=grho2_b_updn,vxcgrho=vxcgrho_b_updn,dvxc=dvxc_b, & 785 & fxcT=fxc_b,hyb_mixing=xcdata%hyb_mixing,el_temp=xcdata%tphysel,& 786 & xc_funcs=xc_funcs_auxc) 787 ! Transfer the xc kernel 788 if (nkxc_eff==1.and.ndvxc==15) then 789 kxc(ifft:ifft+npts-1,1)=half*(dvxc_b(1:npts,1)+dvxc_b(1:npts,9)+dvxc_b(1:npts,10)) 790 else if (nkxc_eff==3.and.ndvxc==15) then 791 kxc(ifft:ifft+npts-1,1)=dvxc_b(1:npts,1)+dvxc_b(1:npts,9) 792 kxc(ifft:ifft+npts-1,2)=dvxc_b(1:npts,10) 793 kxc(ifft:ifft+npts-1,3)=dvxc_b(1:npts,2)+dvxc_b(1:npts,11) 794 end if 795 if (auxc_ixc<0) then 796 call libxc_functionals_end(xc_functionals=xc_funcs_auxc) 797 end if 798 end if 799 if (present(xc_funcs)) then 800 call libxc_functionals_get_hybridparams(hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,& 801 & hyb_range=hyb_range,xc_functionals=xc_funcs) 802 else 803 call libxc_functionals_get_hybridparams(hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,& 804 & hyb_range=hyb_range) 805 end if 806 807 ! Call to main XC driver 808 if (present(xc_funcs)) then 809 call drivexc(ixc,order,npts,nspden_updn,& 810 & usegradient,uselaplacian,usekden,& 811 & rho_b_updn,exc_b,vxcrho_b_updn,& 812 & nvxcgrho,nvxclrho,nvxctau,ndvxc,nd2vxc, & 813 & grho2_updn=grho2_b_updn,vxcgrho=vxcgrho_b_updn,& 814 & lrho_updn=lrho_b_updn,vxclrho=vxclrho_b_updn,& 815 & tau_updn=tau_b_updn,vxctau=vxctau_b_updn,& 816 & dvxc=dvxc_b,d2vxc=d2vxc_b,el_temp=xcdata%tphysel,fxcT=fxc_b,& 817 & hyb_mixing=xcdata%hyb_mixing,& 818 & xc_funcs=xc_funcs) 819 else 820 call drivexc(ixc,order,npts,nspden_updn,& 821 & usegradient,uselaplacian,usekden,& 822 & rho_b_updn,exc_b,vxcrho_b_updn,& 823 & nvxcgrho,nvxclrho,nvxctau,ndvxc,nd2vxc, & 824 & grho2_updn=grho2_b_updn,vxcgrho=vxcgrho_b_updn,& 825 & lrho_updn=lrho_b_updn,vxclrho=vxclrho_b_updn,& 826 & tau_updn=tau_b_updn,vxctau=vxctau_b_updn,& 827 & dvxc=dvxc_b,d2vxc=d2vxc_b,el_temp=xcdata%tphysel,fxcT=fxc_b,& 828 & hyb_mixing=xcdata%hyb_mixing) 829 end if 830 831 ! If fake meta-GGA, has to remove the core contribution 832 ! when electronic effective mass has been modified 833 if (n3xccc>0.and.(ixc==31.or.ixc==34.or.ixc==35)) then 834 if (ixc==31.or.ixc==35) then 835 coeff=one-(one/1.01_dp) 836 if (nspden_updn==1) then 837 coeff=half*coeff 838 do ipts=1,npts 839 exc_b(ipts)=exc_b(ipts)-coeff*xcctau3d(ifft+ipts-1) & 840 & /rho_b_updn(ipts,1) 841 end do 842 else 843 do ipts=1,npts 844 exc_b(ipts)=exc_b(ipts)-coeff*xcctau3d(ifft+ipts-1) & 845 & /(rho_b_updn(ipts,1)+rho_b_updn(ipts,2)) 846 end do 847 end if 848 else 849 message = 'MetaGGA ixc=34 is not yet allowed with a core kinetic energy density!' 850 ABI_ERROR(message) 851 end if 852 end if 853 854 ! Gradient Weiszacker correction to a Thomas-Fermi functional 855 if (my_add_tfw) then 856 vxcgrho_b_updn(:,:)=zero 857 call xctfw(xcdata%tphysel,exc_b,fxc_b,usefxc,rho_b_updn,vxcrho_b_updn,npts,nspden_updn, & 858 & vxcgrho_b_updn,nvxcgrho,grho2_b_updn) 859 end if 860 861 ! Accumulate enxc, strsxc and store vxc (and eventually kxc) 862 dstrsxc=zero 863 do ipts=ifft,ifft+npts-1 864 indx=ipts-ifft+1 865 epsxc=epsxc+rho_b(indx)*exc_b(indx) !will be normalized with respect to the volume later to get enxc ("bigexc"). 866 depsxc(ipts,1)=vxcrho_b_updn(indx,1) 867 exc_str=exc_b(indx);if(usefxc==1) exc_str=fxc_b(indx) 868 if(nspden_updn==1)then 869 strdiag=rho_b(indx)*(exc_str-vxcrho_b_updn(indx,1)) 870 else if(nspden_updn==2)then 871 depsxc(ipts,2)=vxcrho_b_updn(indx,2) 872 ! Note : this is not the complete Vxc in the GGA case 873 strdiag=rho_b(indx)*exc_str & 874 & -rho_b_updn(indx,1)*vxcrho_b_updn(indx,1)& 875 & -(rho_b(indx)-rho_b_updn(indx,1))*vxcrho_b_updn(indx,2) 876 end if 877 dstrsxc=dstrsxc+strdiag 878 879 ! For GGAs, additional terms appear 880 ! (the LB functional does not lead to additional terms) 881 if(ngrad==2 .and. ixc/=13)then 882 883 ! Treat explicitely spin up, spin down and total spin for spin-polarized 884 ! Will exit when ispden=1 is finished if non-spin-polarized 885 do ispden=1,3 886 887 if(nspden_updn==1 .and. ispden>=2)exit 888 889 ! If the norm of the gradient vanishes, then the different terms vanishes, 890 ! but the inverse of the gradient diverges, so skip the update. 891 if(grho2_b_updn(indx,ispden) < 1.0d-24) then 892 depsxc(ipts,ispden+nspden_updn)=zero 893 cycle 894 end if 895 896 ! Compute the derivative of n.e_xc wrt the 897 ! spin up, spin down, or total density. In the non-spin-polarized 898 ! case take the coefficient that will be multiplied by the 899 ! gradient of the total density 900 if(nspden_updn==1)then 901 ! ! Definition of vxcgrho_b_updn changed in v3.3 902 if (nvxcgrho == 3) then 903 coeff=half*vxcgrho_b_updn(indx,1) + vxcgrho_b_updn(indx,3) 904 else 905 coeff=half*vxcgrho_b_updn(indx,1) 906 end if 907 else if(nspden_updn==2)then 908 if (nvxcgrho == 3) then 909 coeff=vxcgrho_b_updn(indx,ispden) 910 else if (ispden /= 3) then 911 coeff=vxcgrho_b_updn(indx,ispden) 912 else if (ispden == 3) then 913 coeff=zero 914 end if 915 end if 916 depsxc(ipts,ispden+nspden_updn)=coeff 917 918 ! Store the gradient of up, down or total density, depending on ispden and nspden, at point ipts 919 if(nspden_updn==1)then 920 grho(1:3)=rhonow(ipts,1,2:4) 921 else if(ispden==1 .and. nspden_updn==2)then 922 grho(1:3)=rhonow(ipts,2,2:4) 923 else if(ispden==2 .and. nspden_updn==2)then 924 grho(1:3)=rhonow(ipts,1,2:4)-rhonow(ipts,2,2:4) 925 else if(ispden==3 .and. nspden_updn==2)then 926 grho(1:3)=rhonow(ipts,1,2:4) 927 end if 928 929 ! In case of ixc 31 (mGGA functional fake 1), 930 ! skip the stress tensor to follow a LDA scheme (see doc/theory/MGGA/report_MGGA.pdf) 931 if(ixc==31) cycle 932 933 ! Compute the contribution to the stress tensor 934 s1=-grho(1)*grho(1)*coeff 935 s2=-grho(2)*grho(2)*coeff 936 s3=-grho(3)*grho(3)*coeff 937 ! The contribution of the next line comes from the part of Vxc 938 ! obtained from the derivative wrt the gradient 939 dstrsxc=dstrsxc+s1+s2+s3 940 strsxc1_tot=strsxc1_tot+s1 941 strsxc2_tot=strsxc2_tot+s2 942 strsxc3_tot=strsxc3_tot+s3 943 strsxc4_tot=strsxc4_tot-grho(3)*grho(2)*coeff 944 strsxc5_tot=strsxc5_tot-grho(3)*grho(1)*coeff 945 strsxc6_tot=strsxc6_tot-grho(2)*grho(1)*coeff 946 947 end do 948 end if 949 950 ! For meta-GGAs, add the laplacian term (vxclrho) and/or kinetic energy density term (vxctau) 951 if (usekden==1.and.with_vxctau) then 952 if (nspden_updn==1)then 953 vxctau(ipts,1,1) = vxctau_b_updn(indx,1) 954 else if (nspden_updn==2)then 955 vxctau(ipts,1,1) = vxctau_b_updn(indx,1) 956 vxctau(ipts,2,1) = vxctau_b_updn(indx,2) 957 end if 958 end if 959 if (uselaplacian==1) then 960 if (nspden_updn==1)then 961 depsxc(ipts,3) = vxclrho_b_updn(indx,1) 962 else if (nspden_updn==2)then 963 depsxc(ipts,6) = vxclrho_b_updn(indx,1) 964 depsxc(ipts,7) = vxclrho_b_updn(indx,2) 965 end if 966 end if 967 968 end do 969 970 ! Additional electron-positron correlation terms 971 if (ipositron==2) then 972 ! Compute electron-positron XC energy per unit volume, potentials and derivatives 973 ngr=0;if (ngrad_apn==2) ngr=npts 974 ABI_MALLOC(fxc_apn,(npts)) 975 ABI_MALLOC(vxc_b_apn,(npts)) 976 ABI_MALLOC(vxcgr_apn,(ngr)) 977 ABI_MALLOC(vxc_ep,(npts)) 978 ABI_MALLOC(rhonow_apn,(npts,nspden_apn,1)) 979 ABI_MALLOC(grho2_apn,(ngr)) 980 rhonow_apn(1:npts,1,1)=electronpositron%rhor_ep(ifft:ifft+npts-1,1) 981 if (usexcnhat==0) rhonow_apn(1:npts,1,1)=rhonow_apn(1:npts,1,1)-electronpositron%nhat_ep(ifft:ifft+npts-1,1) 982 if (.not.electronpositron%posdensity0_limit) then 983 call mkdenpos(iwarnp,npts,nspden_apn,1,rhonow_apn(:,1,1),xcdata%xc_denpos) 984 end if 985 if (ngrad_apn==2.and.usegradient==1) then 986 if (nspden_apn==1) grho2_apn(:)=four*grho2_b_updn(:,1) 987 if (nspden_apn==2) grho2_apn(:)=grho2_b_updn(:,3) 988 end if 989 if (ndvxc==0) then 990 call xcpositron(fxc_apn,grho2_apn,electronpositron%ixcpositron,ngr,npts,& 991 & electronpositron%posdensity0_limit,rho_b,& 992 & rhonow_apn(:,1,1),vxc_b_apn,vxcgr_apn,vxc_ep) 993 else 994 ABI_MALLOC(dvxc_apn,(npts)) 995 call xcpositron(fxc_apn,grho2_apn,electronpositron%ixcpositron,ngr,npts,& 996 & electronpositron%posdensity0_limit,rho_b,& 997 & rhonow_apn(:,1,1),vxc_b_apn,vxcgr_apn,vxc_ep,dvxce=dvxc_apn) 998 end if 999 ABI_FREE(vxc_ep) 1000 ABI_FREE(rhonow_apn) 1001 ABI_FREE(grho2_apn) 1002 ! Accumulate electron-positron XC energies 1003 s1=zero 1004 do ipts=1,npts 1005 s1=s1+fxc_apn(ipts) 1006 end do 1007 electronpositron%e_xc=electronpositron%e_xc+s1*ucvol/dble(nfftot) 1008 ! Add electron-positron dVxc_el/dRho_el to electron-electron one 1009 if (ndvxc==1) dvxc_b(:,1)=dvxc_b(:,1)+dvxc_apn(:) 1010 if (ndvxc==3) then 1011 dvxc_b(:,1)=dvxc_b(:,1)+four*dvxc_apn(:) 1012 dvxc_b(:,2)=dvxc_b(:,2)+four*dvxc_apn(:) 1013 dvxc_b(:,3)=dvxc_b(:,3)+four*dvxc_apn(:) 1014 end if 1015 if (ndvxc==15) then 1016 dvxc_b(:, 9)=dvxc_b(:, 9)+four*dvxc_apn(:) 1017 dvxc_b(:,10)=dvxc_b(:,10)+four*dvxc_apn(:) 1018 dvxc_b(:,11)=dvxc_b(:,11)+four*dvxc_apn(:) 1019 end if 1020 ! Modify stresses - Compute factors for GGA 1021 do ipts=ifft,ifft+npts-1 1022 indx=ipts-ifft+1 1023 depsxc_apn(ipts,1)=vxc_b_apn(indx) 1024 dstrsxc=dstrsxc+fxc_apn(indx)-rho_b(indx)*vxc_b_apn(indx) 1025 if (ngrad_apn==2) then 1026 depsxc_apn(ipts,2)=vxcgr_apn(indx) 1027 s1=-grho(1)*grho(1)*vxcgr_apn(indx) 1028 s2=-grho(2)*grho(2)*vxcgr_apn(indx) 1029 s3=-grho(3)*grho(3)*vxcgr_apn(indx) 1030 dstrsxc=dstrsxc+s1+s2+s3 1031 strsxc1_tot=strsxc1_tot+s1 1032 strsxc2_tot=strsxc2_tot+s2 1033 strsxc3_tot=strsxc3_tot+s3 1034 strsxc4_tot=strsxc4_tot-grho(3)*grho(2)*vxcgr_apn(indx) 1035 strsxc5_tot=strsxc5_tot-grho(3)*grho(1)*vxcgr_apn(indx) 1036 strsxc6_tot=strsxc6_tot-grho(2)*grho(1)*vxcgr_apn(indx) 1037 end if ! GGA 1038 end do ! ipts 1039 ! Deallocations 1040 ABI_FREE(fxc_apn) 1041 ABI_FREE(vxc_b_apn) 1042 ABI_FREE(vxcgr_apn) 1043 if (ndvxc>0) then 1044 ABI_FREE(dvxc_apn) 1045 end if 1046 end if 1047 1048 ! Transfer the xc kernel (if this must be done, and has not yet been done) 1049 if (nkxc_eff>0.and.ndvxc>0 .and. (xcdata%usefock==0 .or. auxc_ixc==0)) then 1050 if (nkxc_eff==1.and.ndvxc==15) then 1051 kxc(ifft:ifft+npts-1,1)=half*(dvxc_b(1:npts,1)+dvxc_b(1:npts,9)+dvxc_b(1:npts,10)) 1052 else if (nkxc_eff==3.and.ndvxc==15) then 1053 kxc(ifft:ifft+npts-1,1)=dvxc_b(1:npts,1)+dvxc_b(1:npts,9) 1054 kxc(ifft:ifft+npts-1,2)=dvxc_b(1:npts,10) 1055 kxc(ifft:ifft+npts-1,3)=dvxc_b(1:npts,2)+dvxc_b(1:npts,11) 1056 else if (nkxc_eff==7.and.ndvxc==8) then 1057 kxc(ifft:ifft+npts-1,1)=half*dvxc_b(1:npts,1) 1058 kxc(ifft:ifft+npts-1,2)=half*dvxc_b(1:npts,3) 1059 kxc(ifft:ifft+npts-1,3)=quarter*dvxc_b(1:npts,5) 1060 kxc(ifft:ifft+npts-1,4)=eighth*dvxc_b(1:npts,7) 1061 else if (nkxc_eff==7.and.ndvxc==15) then 1062 kxc(ifft:ifft+npts-1,1)=half*(dvxc_b(1:npts,1)+dvxc_b(1:npts,9)+dvxc_b(1:npts,10)) 1063 kxc(ifft:ifft+npts-1,2)=half*dvxc_b(1:npts,3)+dvxc_b(1:npts,12) 1064 kxc(ifft:ifft+npts-1,3)=quarter*dvxc_b(1:npts,5)+dvxc_b(1:npts,13) 1065 kxc(ifft:ifft+npts-1,4)=eighth*dvxc_b(1:npts,7)+dvxc_b(1:npts,15) 1066 else if (nkxc_eff==19.and.ndvxc==15) then 1067 kxc(ifft:ifft+npts-1,1)=dvxc_b(1:npts,1)+dvxc_b(1:npts,9) 1068 kxc(ifft:ifft+npts-1,2)=dvxc_b(1:npts,10) 1069 kxc(ifft:ifft+npts-1,3)=dvxc_b(1:npts,2)+dvxc_b(1:npts,11) 1070 kxc(ifft:ifft+npts-1,4)=dvxc_b(1:npts,3) 1071 kxc(ifft:ifft+npts-1,5)=dvxc_b(1:npts,4) 1072 kxc(ifft:ifft+npts-1,6)=dvxc_b(1:npts,5) 1073 kxc(ifft:ifft+npts-1,7)=dvxc_b(1:npts,6) 1074 kxc(ifft:ifft+npts-1,8)=dvxc_b(1:npts,7) 1075 kxc(ifft:ifft+npts-1,9)=dvxc_b(1:npts,8) 1076 kxc(ifft:ifft+npts-1,10)=dvxc_b(1:npts,12) 1077 kxc(ifft:ifft+npts-1,11)=dvxc_b(1:npts,13) 1078 kxc(ifft:ifft+npts-1,12)=dvxc_b(1:npts,14) 1079 kxc(ifft:ifft+npts-1,13)=dvxc_b(1:npts,15) 1080 else ! All other cases 1081 kxc(ifft:ifft+npts-1,1:nkxc_eff)=zero 1082 kxc(ifft:ifft+npts-1,1:min(nkxc_eff,ndvxc))=dvxc_b(1:npts,1:min(nkxc_eff,ndvxc)) 1083 end if 1084 if (nkxc_eff==7) then 1085 kxc(ifft:ifft+npts-1,5)=rhonow(ifft:ifft+npts-1,1,2) 1086 kxc(ifft:ifft+npts-1,6)=rhonow(ifft:ifft+npts-1,1,3) 1087 kxc(ifft:ifft+npts-1,7)=rhonow(ifft:ifft+npts-1,1,4) 1088 else if (nkxc_eff==19) then 1089 kxc(ifft:ifft+npts-1,14)=rhonow(ifft:ifft+npts-1,1,2) 1090 kxc(ifft:ifft+npts-1,15)=rhonow(ifft:ifft+npts-1,2,2) 1091 kxc(ifft:ifft+npts-1,16)=rhonow(ifft:ifft+npts-1,1,3) 1092 kxc(ifft:ifft+npts-1,17)=rhonow(ifft:ifft+npts-1,2,3) 1093 kxc(ifft:ifft+npts-1,18)=rhonow(ifft:ifft+npts-1,1,4) 1094 kxc(ifft:ifft+npts-1,19)=rhonow(ifft:ifft+npts-1,2,4) 1095 end if 1096 end if 1097 1098 ! Transfer the XC 3rd-derivative 1099 if (abs(option)==3.and.order==3.and.nd2vxc>0) then 1100 k3xc(ifft:ifft+npts-1,1:nd2vxc)=d2vxc_b(1:npts,1:nd2vxc) 1101 end if 1102 1103 ! Add the diagonal part to the xc stress 1104 strsxc1_tot=strsxc1_tot+dstrsxc 1105 strsxc2_tot=strsxc2_tot+dstrsxc 1106 strsxc3_tot=strsxc3_tot+dstrsxc 1107 1108 ! Accumulate integral of |Grad_rho|/Rho (to be used for TB09 XC) 1109 if (present(grho1_over_rho1).and.ixc<0) then 1110 if (present(xc_funcs)) then 1111 test_tb09=libxc_functionals_is_tb09(xc_functionals=xc_funcs) 1112 else 1113 test_tb09=libxc_functionals_is_tb09() 1114 end if 1115 if (test_tb09) then 1116 factor=merge(two,one,nspden_updn==1) 1117 jj=merge(1,3,nspden_updn==1) 1118 do ipts=ifft,ifft+npts-1 1119 indx=ipts-ifft+1 1120 if (abs(rho_b(indx))>tol10) then 1121 grho1_over_rho1=grho1_over_rho1+factor*sqrt(grho2_b_updn(indx,jj))/rho_b(indx) 1122 end if 1123 end do 1124 end if 1125 end if 1126 1127 ABI_FREE(exc_b) 1128 ABI_FREE(rho_b) 1129 ABI_FREE(rho_b_updn) 1130 ABI_FREE(grho2_b_updn) 1131 ABI_FREE(vxcrho_b_updn) 1132 ABI_FREE(dvxc_b) 1133 ABI_FREE(d2vxc_b) 1134 ABI_FREE(vxcgrho_b_updn) 1135 ABI_FREE(fxc_b) 1136 ABI_FREE(vxclrho_b_updn) 1137 ABI_FREE(lrho_b_updn) 1138 ABI_FREE(tau_b_updn) 1139 ABI_FREE(vxctau_b_updn) 1140 1141 ! End of the loop on blocks of data 1142 end do 1143 1144 strsxc(1)=strsxc1_tot 1145 strsxc(2)=strsxc2_tot 1146 strsxc(3)=strsxc3_tot 1147 strsxc(4)=strsxc4_tot 1148 strsxc(5)=strsxc5_tot 1149 strsxc(6)=strsxc6_tot 1150 1151 ! If GGA, multiply the gradient of the density by the proper 1152 ! local partial derivatives of the XC functional 1153 rhonow_ptr => rhonow 1154 if (ipositron==2) then 1155 ABI_MALLOC(rhonow_ptr,(nfft,nspden_eff,ngrad*ngrad)) 1156 rhonow_ptr=rhonow 1157 end if 1158 if(ngrad==2 .and. ixc/=13)then 1159 call xcmult(depsxc,nfft,ngrad,nspden_eff,nspgrad,rhonow_ptr) 1160 end if 1161 1162 ! Compute contribution from this grid to vxc, and ADD to existing vxc 1163 if (nspden/=4) then 1164 if(with_vxctau)then 1165 call xcpot(cplex,gprimd,ishift,uselaplacian,mpi_enreg,nfft,ngfft,ngrad,nspden_eff,nspgrad,& 1166 & qphon,depsxc=depsxc,rhonow=rhonow_ptr,vxc=vxc,vxctau=vxctau) 1167 else 1168 call xcpot(cplex,gprimd,ishift,uselaplacian,mpi_enreg,nfft,ngfft,ngrad,nspden_eff,nspgrad,& 1169 & qphon,depsxc=depsxc,rhonow=rhonow_ptr,vxc=vxc) 1170 end if 1171 1172 else 1173 1174 ! If non-collinear magnetism, restore potential in proper axis before adding it 1175 ABI_MALLOC(vxcrho_b_updn,(nfft,4)) 1176 vxcrho_b_updn=zero 1177 call xcpot(cplex,gprimd,ishift,uselaplacian,mpi_enreg,nfft,ngfft,ngrad,nspden_eff,nspgrad,& 1178 & qphon,depsxc=depsxc,rhonow=rhonow_ptr,vxc=vxcrho_b_updn) 1179 do ifft=1,nfft 1180 dvdn=half*(vxcrho_b_updn(ifft,1)+vxcrho_b_updn(ifft,2)) 1181 if(m_norm(ifft)>m_norm_min) then 1182 ! if(m_norm(ifft)>rhor_(ifft,1)*tol10+tol14) then 1183 dvdz=half*(vxcrho_b_updn(ifft,1)-vxcrho_b_updn(ifft,2))/m_norm(ifft) 1184 vxc(ifft,1)=vxc(ifft,1)+dvdn+rhor_(ifft,4)*dvdz 1185 vxc(ifft,2)=vxc(ifft,2)+dvdn-rhor_(ifft,4)*dvdz 1186 vxc(ifft,3)=vxc(ifft,3)+rhor_(ifft,2)*dvdz 1187 vxc(ifft,4)=vxc(ifft,4)-rhor_(ifft,3)*dvdz 1188 else 1189 vxc(ifft,1:2)=vxc(ifft,1:2)+dvdn 1190 end if 1191 end do 1192 ABI_FREE(vxcrho_b_updn) 1193 end if 1194 if (ipositron==2) then 1195 ABI_FREE(rhonow_ptr) 1196 end if 1197 nullify(rhonow_ptr) 1198 1199 ! Add electron-positron XC potential to electron-electron one 1200 ! Eventually compute GGA contribution 1201 if (ipositron==2) then 1202 ABI_MALLOC(rhonow_apn,(nfft,nspden_apn,ngrad_apn**2)) 1203 rhonow_apn(1:nfft,1,1:ngrad_apn**2)=rhonow(1:nfft,1,1:ngrad_apn**2) 1204 if (ngrad_apn==2) then 1205 call xcmult(depsxc_apn,nfft,ngrad_apn,nspden_apn,ngrad_apn,rhonow_apn) 1206 end if 1207 ABI_MALLOC(vxc_apn,(nfft,nspden_apn)) 1208 vxc_apn=zero 1209 call xcpot(cplex,gprimd,ishift,0,mpi_enreg,nfft,ngfft,ngrad_apn,& 1210 & nspden_apn,ngrad_apn,qphon,depsxc=depsxc_apn,rhonow=rhonow_apn,vxc=vxc_apn) 1211 vxc(:,1)=vxc(:,1)+vxc_apn(:,1) 1212 if (nspden_updn==2) vxc(:,2)=vxc(:,2)+vxc_apn(:,1) 1213 s1=zero 1214 do ipts=1,nfft 1215 s1=s1+vxc_apn(ipts,1)*rhonow(ipts,1,1) 1216 end do 1217 electronpositron%e_xcdc=electronpositron%e_xcdc+s1*ucvol/dble(nfftot) 1218 ABI_FREE(rhonow_apn) 1219 ABI_FREE(vxc_apn) 1220 ABI_FREE(depsxc_apn) 1221 end if 1222 1223 ! End loop on unshifted or shifted grids 1224 end do 1225 1226 ! Calculate van der Waals correction when requested 1227 #if defined DEV_YP_VDWXC 1228 if ( (xcdata%vdw_xc > 0) .and. (xcdata%vdw_xc < 3) .and. (xc_vdw_status()) ) then 1229 call xc_vdw_aggregate(ucvol,gprimd,nfft,nspden_updn,ngrad*ngrad, & 1230 & ngfft(1),ngfft(2),ngfft(3),rhonow, & 1231 & deltae_vdw,exc_vdw,decdrho_vdw,decdgrho_vdw,strsxc_vdw) 1232 end if 1233 #else 1234 if ( (xcdata%vdw_xc > 0) .and. (xcdata%vdw_xc < 3) ) then 1235 write(message,'(3a)')& 1236 & 'vdW-DF functionals are not fully operational yet.',ch10,& 1237 & 'Action : modify vdw_xc' 1238 ABI_ERROR(message) 1239 end if 1240 #endif 1241 ! Normalize enxc, strsxc and vxc 1242 divshft=one/dble(xcdata%intxc+1) 1243 strsxc(:)=strsxc(:)/dble(nfftot)*divshft 1244 enxc=epsxc*ucvol/dble(nfftot)*divshft 1245 vxc=vxc*divshft 1246 if (with_vxctau) vxctau=vxctau*divshft 1247 if (present(grho1_over_rho1)) grho1_over_rho1=grho1_over_rho1*ucvol/dble(nfftot)*divshft 1248 1249 ! Reduction in case of FFT distribution 1250 if (nproc_fft>1)then 1251 call timab(48,1,tsec) 1252 call xmpi_sum(strsxc,comm_fft ,ierr) 1253 call xmpi_sum(enxc ,comm_fft ,ierr) 1254 if (present(grho1_over_rho1)) then 1255 call xmpi_sum(grho1_over_rho1,comm_fft ,ierr) 1256 end if 1257 if (ipositron==2) then 1258 s1=electronpositron%e_xc;s2=electronpositron%e_xcdc 1259 call xmpi_sum(s1,comm_fft ,ierr) 1260 call xmpi_sum(s2,comm_fft ,ierr) 1261 electronpositron%e_xc=s1;electronpositron%e_xcdc=s2 1262 end if 1263 call timab(48,2,tsec) 1264 end if 1265 1266 ! Compute vxcavg 1267 call mean_fftr(vxc,vxcmean,nfft,nfftot,min(nspden,2),mpi_comm_sphgrid=comm_fft) 1268 if(nspden==1)then 1269 vxcavg=vxcmean(1) 1270 else 1271 vxcavg=half*(vxcmean(1)+vxcmean(2)) 1272 end if 1273 1274 ABI_FREE(depsxc) 1275 ABI_FREE(rhonow) 1276 ABI_FREE(lrhonow) 1277 ABI_FREE(taunow) 1278 if (need_nhat.or.non_magnetic_xc) then 1279 ABI_FREE(rhor_) 1280 end if 1281 if ((usekden==1).and.(non_magnetic_xc)) then 1282 ABI_FREE(taur_) 1283 end if 1284 if (allocated(m_norm)) then 1285 ABI_FREE(m_norm) 1286 end if 1287 1288 end if 1289 1290 !Treat separately the Fermi-Amaldi correction. 1291 if (ixc==20 .or. ixc==21 .or. ixc==22) then 1292 if(present(vhartr))then 1293 1294 ! Fermi-Amaldi correction : minus Hartree divided by the 1295 ! number of electrons per unit cell. This is not size consistent, but 1296 ! interesting for isolated systems with a few electrons. 1297 ! nelect=ucvol*rhog(1,1) 1298 factor=-one/xcdata%nelect 1299 vxc(:,1)=factor*vhartr(:) 1300 if(nspden>=2) vxc(:,2)=factor*vhartr(:) 1301 1302 ! Compute corresponding xc energy and stress as well as vxcavg 1303 call dotprod_vn(1,rhor,enxc,doti,nfft,nfftot,1,1,vxc,ucvol,mpi_comm_sphgrid=comm_fft) 1304 enxc=half*enxc 1305 strsxc(1:3)=-enxc/ucvol 1306 1307 ! Compute average of vxc (one component only). 1308 call mean_fftr(vxc,vxcmean,nfft,nfftot,1,mpi_comm_sphgrid=comm_fft) 1309 vxcavg = vxcmean(1) 1310 ! For ixc=20, the local exchange-correlation kernel is zero, but the Hartree 1311 ! kernel will be modified in tddft. No other use of kxc should be made with ixc==20 1312 if(nkxc/=0 .and. ixc==20) kxc(:,:)=zero 1313 ! For ixc=21 or 22, the LDA (ixc=1) kernel has been computed previously. 1314 1315 else 1316 1317 ABI_BUG('When ixc=20,21 or 22, vhartr needs to be present in the call to rhotoxc !') 1318 1319 end if 1320 1321 end if 1322 1323 !Add van der Waals terms 1324 #if defined DEV_YP_VDWXC 1325 if ( (xcdata%vdw_xc > 0) .and. (xcdata%vdw_xc < 10) .and. (xc_vdw_status()) ) then 1326 enxc = enxc + exc_vdw + deltae_vdw 1327 do ispden=1,nspden 1328 vxc(:,ispden) = vxc(:,ispden) + decdrho_vdw(:,ispden) 1329 end do 1330 strsxc(1) = strsxc(1) + strsxc_vdw(1,1) 1331 strsxc(2) = strsxc(2) + strsxc_vdw(2,2) 1332 strsxc(3) = strsxc(3) + strsxc_vdw(3,3) 1333 strsxc(4) = strsxc(4) + strsxc_vdw(3,2) 1334 strsxc(5) = strsxc(5) + strsxc_vdw(3,1) 1335 strsxc(6) = strsxc(6) + strsxc_vdw(2,1) 1336 end if 1337 #endif 1338 if ( present(exc_vdw_out) ) exc_vdw_out = exc_vdw 1339 1340 call timab(81,2,tsec) 1341 1342 DBG_EXIT("COLL") 1343 1344 end subroutine rhotoxc