TABLE OF CONTENTS


ABINIT/m_rhotoxc [ Modules ]

[ Top ] [ Modules ]

NAME

  m_rhotox

FUNCTION

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_rhotoxc
27 
28  use defs_basis
29  use defs_abitypes
30  use m_xmpi
31  use m_abicore
32  use m_errors
33  use m_cgtools
34  use m_xcdata
35  use m_xc_vdw
36  use libxc_functionals
37 
38  use m_time,             only : timab
39  use m_geometry,         only : metric
40  use m_electronpositron, only : electronpositron_type,electronpositron_calctype
41  use m_xcpositron,       only : xcpositron
42  use m_drivexc,          only : size_dvxc, drivexc_main, xcmult, mkdenpos
43  use m_xclda,            only : xctfw
44  use m_xctk,             only : xcden, xcpot
45 
46  implicit none
47 
48  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= true if usepawu==4
  n3xccc=dimension of the xccc3d array (0 or nfft or cplx*nfft).
  option=0 or 1 for xc only (exc, vxc, strsxc),
         2 for xc and kxc (no paramagnetic part if xcdata%nspden=1)
        10 for xc  and kxc with only partial derivatives wrt density part (d2Exc/drho^2)
        12 for xc and kxc with only partial derivatives wrt density part (d2Exc/drho^2)
              and, in the case of hybrid functionals, substitution of the hybrid functional
              by the related auxiliary GGA functional for the computation of the xc kernel (not for other quantities)
         3 for xc, kxc and k3xc
        -2 for xc and kxc (with paramagnetic part if xcdata%nspden=1)
  paral_kgb=Flag related to the kpoint-band-fft parallelism
  rhor(nfft,xcdata%nspden)=electron density in real space in electrons/bohr**3
   (total in first half and spin-up in second half if xcdata%nspden=2)
   (total in first comp. and magnetization in comp. 2 to 4 if xcdata%nspden=4)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
  [vhartr(nfft)=Hartree potential (only needed for Fermi-Amaldi functional)]
  xcdata <type(xcdata_type)>=storage for different input variables and derived parameters needed to compute the XC functional
  xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)

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

OUTPUT

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

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

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

SIDE EFFECTS

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

NOTES

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

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

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

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

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

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

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

PARENTS

      calc_vhxc_me,energy,m_kxc,nonlinear,nres2vres,odamix,prcref,prcref_PMA
      respfn,rhotov,scfcv,setvtr,xchybrid_ncpp_cc

CHILDREN

      dotprod_vn,drivexc_main,libxc_functionals_end
      libxc_functionals_get_hybridparams,libxc_functionals_init,mean_fftr
      metric,mkdenpos,size_dvxc,timab,xc_vdw_aggregate,xcden,xcmult
      xcpositron,xcpot,xctfw,xmpi_sum

SOURCE

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