TABLE OF CONTENTS


ABINIT/m_rhotoxc [ Modules ]

[ Top ] [ Modules ]

NAME

  m_rhotox

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2022 ABINIT group (DCA, XG, GMR, MF, GZ, DRH, MT, SPr)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

15 #if defined HAVE_CONFIG_H
16 #include "config.h"
17 #endif
18 
19 #include "abi_common.h"
20 
21 module m_rhotoxc
22 
23  use defs_basis
24  use m_xmpi
25  use m_abicore
26  use m_errors
27  use m_cgtools
28  use m_xcdata
29  use m_xc_vdw
30  use libxc_functionals
31 
32  use defs_abitypes,      only : MPI_type
33  use m_time,             only : timab
34  use m_geometry,         only : metric
35  use m_electronpositron, only : electronpositron_type,electronpositron_calctype
36  use m_xcpositron,       only : xcpositron
37  use m_drivexc,          only : size_dvxc,drivexc,xcmult,mkdenpos
38  use m_xclda,            only : xctfw
39  use m_xctk,             only : xcden, xcpot
40 
41  implicit none
42 
43  private

ABINIT/rhotoxc [ Functions ]

[ Top ] [ Functions ]

NAME

 rhotoxc

FUNCTION

 Start from the density or spin-density, and
 compute xc correlation potential and energies.
 Eventually compute xc kernel (if option=-2, 2, 3, 10 or 12).
 Cannot be used with wavelets.

INPUTS

  mpi_enreg=information about MPI parallelization
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nhat(nfft,xcdata%nspden*nhatdim)= -PAW only- compensation density
  nhatdim= -PAW only- 0 if nhat array is not used ; 1 otherwise
  nhatgr(nfft,xcdata%nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density
  nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise
  nkxc=second dimension of the kxc array. If /=0,
   the exchange-correlation kernel must be computed.
  non_magnetic_xc= if true, handle density/potential as non-magnetic (even if it is)
  n3xccc=dimension of the xccc3d array (0 or nfft or cplx*nfft).
  option=0 or 1 for xc only (exc, vxc, strsxc),
         2 for xc and kxc (no paramagnetic part if xcdata%nspden=1)
        10 for xc  and kxc with only partial derivatives wrt density part (d2Exc/drho^2)
        12 for xc and kxc with only partial derivatives wrt density part (d2Exc/drho^2)
              and, in the case of hybrid functionals, substitution of the hybrid functional
              by the related auxiliary GGA functional for the computation of the xc kernel (not for other quantities)
         3 for xc, kxc and k3xc
        -2 for xc and kxc (with paramagnetic part if xcdata%nspden=1)
  rhor(nfft,xcdata%nspden)=electron density in real space in electrons/bohr**3
   (total in first half and spin-up in second half if xcdata%nspden=2)
   (total in first comp. and magnetization in comp. 2 to 4 if xcdata%nspden=4)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
  [vhartr(nfft)=Hartree potential (only needed for Fermi-Amaldi functional)]
  xcdata <type(xcdata_type)>=storage for different input variables and derived parameters needed to compute the XC functional
  xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)

  === optional inputs ===
  [add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy
  [taur(nfftf,xcdata%nspden*xcdata%usekden)]=array for kinetic energy density
  [xc_funcs(2)]= <type(libxc_functional_type)>, optional : libxc XC functionals. Must be coherent with xcdata.
  [xccctau3d(n3xccc)]=3D core electron kinetic energy density for XC core correction (bohr^-3)

OUTPUT

  enxc=returned exchange and correlation energy (hartree).
  strsxc(6)= contribution of xc to stress tensor (hartree/bohr^3),
   given in order (1,1), (2,2), (3,3), (3,2), (3,1), (2,1).
   (note: fxc is rho*exc in the following)
   Explicitely : strsxc(mu,nu) = (1/N) Sum(i=1,N)
    ( delta(mu,nu) * [  exc(i)rhotot(i)
               - depsxc_drho(up,i)*rhor(up,i)-depsxc_drho(dn,i)*rhor(dn,i)]
     - gradrho(up,mu)*gradrho(up,nu) * depsxc_dgradrho(up,i) / gradrho(up,i)
     - gradrho(dn,mu)*gradrho(dn,nu) * depsxc_dgradrho(dn,i) / gradrho(dn,i) )
  vxc(nfft,xcdata%nspden)=xc potential
    (spin up in first half and spin down in second half if xcdata%nspden=2)
    (v^11, v^22, Re[V^12], Im[V^12] if xcdata%nspden=4)
  vxcavg=<Vxc>=unit cell average of Vxc = (1/ucvol) Int [Vxc(r) d^3 r].

  === Only if abs(option)=2, -2, 3, 10, 12 (in case 12, for hybrids, substitution of the related GGA) ===
  kxc(nfft,nkxc)=exchange and correlation kernel (returned only if nkxc/=0)
    Content of Kxc array:
   ===== if LDA
    if xcdata%nspden==1: kxc(:,1)= d2Exc/drho2
              that is 1/2 ( d2Exc/drho_up drho_up + d2Exc/drho_up drho_dn )
                         kxc(:,2)= d2Exc/drho_up drho_dn
    if xcdata%nspden>=2: kxc(:,1)= d2Exc/drho_up drho_up
                         kxc(:,2)= d2Exc/drho_up drho_dn
                         kxc(:,3)= d2Exc/drho_dn drho_dn
   ===== if GGA or mGGA
    if xcdata%nspden==1:
       kxc(:,1)= d2Exc/drho2
       kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)|
       kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho
       kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| )
       kxc(:,5)= gradx(rho)
       kxc(:,6)= grady(rho)
       kxc(:,7)= gradz(rho)
    if xcdata%nspden>=2:
       kxc(:,1)= d2Exc/drho_up drho_up
       kxc(:,2)= d2Exc/drho_up drho_dn
       kxc(:,3)= d2Exc/drho_dn drho_dn
       kxc(:,4)= 1/|grad(rho_up)| dEx/d|grad(rho_up)|
       kxc(:,5)= 1/|grad(rho_dn)| dEx/d|grad(rho_dn)|
       kxc(:,6)= 1/|grad(rho_up)| d2Ex/d|grad(rho_up)| drho_up
       kxc(:,7)= 1/|grad(rho_dn)| d2Ex/d|grad(rho_dn)| drho_dn
       kxc(:,8)= 1/|grad(rho_up)| * d/d|grad(rho_up)| ( 1/|grad(rho_up)| dEx/d|grad(rho_up)| )
       kxc(:,9)= 1/|grad(rho_dn)| * d/d|grad(rho_dn)| ( 1/|grad(rho_dn)| dEx/d|grad(rho_dn)| )
       kxc(:,10)=1/|grad(rho)| dEc/d|grad(rho)|
       kxc(:,11)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_up
       kxc(:,12)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_dn
       kxc(:,13)=1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dEc/d|grad(rho)| )
       kxc(:,14)=gradx(rho_up)
       kxc(:,15)=gradx(rho_dn)
       kxc(:,16)=grady(rho_up)
       kxc(:,17)=grady(rho_dn)
       kxc(:,18)=gradz(rho_up)
       kxc(:,19)=gradz(rho_dn)
    Note about mGGA: 2nd derivatives involving Tau or Laplacian are not output

  === Only if abs(option)=3 ===
  [k3xc(nfft,nk3xc)]= -- optional -- third derivative of the XC energy functional of the density,
    at each point of the real space grid (only in the LDA or LSDA)
    Content of K3xc array:
    ===== if LDA
    if xcdata%nspden==1: k3xc(:,1)= d3Exc/drho3
    if xcdata%nspden>=2, k3xc(:,1)= d3Exc/drho_up drho_up drho_up
                  k3xc(:,2)= d3Exc/drho_up drho_up drho_dn
                  k3xc(:,3)= d3Exc/drho_up drho_dn drho_dn
                  k3xc(:,4)= d3Exc/drho_dn drho_dn drho_dn
 === Additional optional output ===
  [exc_vdw_out]= vdW-DF contribution to enxc (hartree)
  [vxctau(nfft,xcdata%nspden,4*xcdata%usekden)]=(only for meta-GGA)=
    vxctau(:,:,1): derivative of XC energy density with respect to kinetic energy density (depsxcdtau).
    vxctau(:,:,2:4): gradient of vxctau (gvxctau)

SIDE EFFECTS

  electronpositron <type(electronpositron_type)>= -- optional argument -- quantities for the electron-positron annihilation

NOTES

 Start from the density, and compute Hartree (if option>=1) and xc correlation potential and energies.
 Eventually compute xc kernel (if option=-2, 2, 3, 10 or 12 - in the latter case, substitution by the related GGA kernel).
 Allows a variety of exchange-correlation functionals
 according to ixc. Here is a list of allowed values.
                                                    subroutine name
   <0 means use of libxc
    0 means no xc applied (usually for testing)
 *LDA,LSD
    1 means new Teter (4/93) with spin-pol option        xcspol
    2 means Perdew-Zunger-Ceperley-Alder                 xcpzca
    3 means old Teter (4/91) fit to Ceperley-Alder data  xctetr
    4 means Wigner                                       xcwign
    5 means Hedin-Lundqvist                              xchelu
    6 means "X-alpha" xc                                 xcxalp
    7 mean Perdew-Wang 92 LSD fit to Ceperley-Alder data xcpbe
    8 mean Perdew-Wang 92 LSD , exchange-only            xcpbe
    9 mean Perdew-Wang 92 Ex+Ec_RPA  energy              xcpbe
   10 means RPA LSD energy (only the energy !!)          xcpbe
 *GGA
   11 means Perdew-Burke-Ernzerhof GGA functional        xcpbe
   12 means x-only Perdew-Burke-Ernzerhof GGA functional xcpbe
   13 means LDA (ixc==7), except that the xc potential
      is given within the van Leeuwen-Baerends GGA       xclb
   14 means revPBE GGA functional                        xcpbe
   15 means RPBE GGA functional                          xcpbe
   16 means HCTH GGA functional                          xchcth
   23 means WC GGA functional                            xcpbe
   24 means C09x GGA exchange functional                 xcpbe
 *Fermi-Amaldi
   20 means Fermi-Amaldi correction
   21 means Fermi-Amaldi correction with LDA(ixc=1) kernel
   22 means Fermi-Amaldi correction with hybrid BPG kernel
 *Hybrid GGA
   41 means PBE0-1/4                                     xcpbe
   42 means PBE0-1/3                                     xcpbe
 *Other
   50 means IIT xc                                       xciit

 NOTE: please update echo_xc_name.F90 if you add new functional (apart from libxc)

 Allow for improved xc quadrature (intxc=1) by using the usual FFT grid
 as well as another, shifted, grid, and combining both results.
 Spin-polarization is allowed only with ixc=0, 1, and GGAs until now.

 To make the variable names easier to understand, a rule notation is tentatively proposed here:
   rho ---> means density
   tau ---> means kinetic energy density
   exc ---> means exchange-correlation energy density per particle
   epsxc ---> means rho*exc == exchange-correlation energy density
   vxc ---> means exchange-correlation potential
   bigexc ---> means exchange-correlation energy E_xc (for the moment it is named "enxc")
   m_norm ---> means norm of magnetization

   g... --> means gradient of something (e.g. : grho --> means gradient of electron density)
   g...2 -> means square norm of gradient of something (e.g. : grho2 -> means square norm of gradient of electron density)
   l... --> means laplacian of something (e.g. : lrho --> means laplacian of electron density)
   d...d... --> means derivative of something with regards to something else.
   (d2...d...d...  ---> means second derivative of ... with regards to ... and to ...) etc...
   d... --> without the occurence of the second "d" means that this is an array of
            several derivative of the same quantity (e.g. : depsxc)

   ..._b ----> means a block of the quantity "..." (use in mpi loops which treat the data block by block)
   ..._updn -> means that spin up and spin down is available in that array
               as (..,1) and (..,2). (if xcdata%nspden >=2 of course).
   ..._apn --> in case of positrons are concerned.

   for more details about notations please see pdf in /doc/theory/MGGA/

SOURCE

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