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