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