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