atomXC Subroutine

 This version implements the Local (spin) Density Approximation and
 the Generalized-Gradient-Aproximation with the 'explicit mesh 
 functional' approach of White & Bird, PRB 50, 4954 (1994).
 Gradients are 'defined' by numerical derivatives, using 2*nn+1 mesh
   points, where nn is a parameter defined below
 Ref: L.C.Balbas et al, PRB 64, 165110 (2001)
 Written by J.M.Soler using algorithms developed by 
   L.C.Balbas, J.L.Martins and J.M.Soler, Dec.1996
 Van der Waals functional added by J.M.Soler, Jul.2008, as explained in
   G.Roman-Perez and J.M.Soler, PRL 103, 096102 (2009)
 ************************* INPUT ***********************************
 INTEGER irel         : Relativistic exchange? (0=>no, 1=>yes)
 INTEGER nr           : Number of radial mesh points
 INTEGER maxr         : Physical first dimension of Dens and Vxc
 REAL*8  rmesh(nr)    : Radial mesh points. Must be nr.le.maxr
 INTEGER nSpin        : nSpin=1 => unpolarized; nSpin=2 => polarized
 REAL*8  Dens(maxr,nSpin) : Total (nSpin=1) or spin (nSpin=2) electron
                            density at mesh points
 ************************* OUTPUT **********************************
 REAL*8  Ex              : Total exchange energy
 REAL*8  Ec              : Total correlation energy
 REAL*8  Dx              : IntegralOf( rho * (eps_x - v_x) )
 REAL*8  Dc              : IntegralOf( rho * (eps_c - v_c) )
 REAL*8  Vxc(maxr,nSpin) : (Spin) exch-corr potential
 ************************ UNITS ************************************
 Distances in atomic units (Bohr).
 Densities in atomic units (electrons/Bohr**3)
 Energy unit depending of parameter Eunit below (currently Rydberg)