cellXC 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 - Aug.1997
 Parallel version written by J.Gale. June 1999.
 Argument dVxcdD added by J.Junquera. November 2000.
 Adapted for multiple functionals in the same run by J.Gale 2005
 Van der Waals functional added by J.M.Soler, Jan.2008, as explained in
   G.Roman-Perez and J.M.Soler, PRL 103, 096102 (2009)
 ************************* INPUT ***********************************
 integer  irel        : Relativistic exchange? (0=>no, 1=>yes)
 real(dp) cell(3,3)   : Unit cell vectors cell(ixyz,ivector)
 integer  nMesh(3)    : Total mesh divisions of each cell vector
 integer  lb1,lb2,lb3 : Lower bounds of arrays dens, Vxc, dVxcdD
 integer  ub1,ub2,ub3 : Upper bounds of arrays dens, Vxc, dVxcdD
 integer  nSpin       : nSpin=1 => unpolarized; nSpin=2 => polarized;
                        nSpin>2 => non-collinear polarization
 real(grid_p) dens(lb1:ub1,lb2:ub2,lb3:ub3,nSpin) : Total (nSpin=1) or 
                        spin (nSpin=2) electron density at mesh points
                        For non-collinear polarization, the density
                        matrix is given by: dens(1)=D11, dens(2)=D22,
                        dens(3)=Real(D12), dens(4)=Im(D12)
 ************************* OUTPUT **********************************
 real(dp) Ex             : Total exchange energy per unit cell
 real(dp) Ec             : Total correlation energy per unit cell
 real(dp) Dx             : IntegralOf( rho * (eps_x - v_x) ) in unit cell
 real(dp) Dc             : IntegralOf( rho * (eps_c - v_c) ) in unit cell
 real(dp) stress(3,3)    : xc contribution to the stress tensor, in unit
                           cell, assuming constant density (not charge),
                           i.e. r->r' => rho'(r') = rho(r)
                           For plane-wave and grid (finite diff) basis
                           sets, density rescaling gives an extra term
                           (not included) (Dx+Dc-Ex-Ec)/cell_volume for
                           the diagonal elements of stress. For other
                           basis sets, the extra term is, in general:
                           IntegralOf(v_xc * d_rho/d_strain) / cell_vol
 real(grid_p) Vxc(lb1:ub1,lb2:ub2,lb3:ub3,nSpin) : (Spin) xc potential
 ************************* OPTIONAL INPUT *************************
 logical keep_input_distribution :  Avoid internal load-balancer and lock
                                    the passed parallel distribution
 ************************* OPTIONAL OUTPUT *************************
 real(grid_p) dVxcdD(lb1:ub1,lb2:ub2,lb3:ub3,nSpin*nSpin) : Derivatives
                           of xc potential respect to charge density
                           Available only for LDA (libxc + built-in CA/PZ)
                           and collinear-spin
 *** OPTIONAL OUTPUT (must set keep_input_distribution=.true.) *****
 real(grid_p) epsxc(lb1:ub1,lb2:ub2,lb3:ub3) : XC energy density per electron
 real(grid_p) dexcdGD(lb1:ub1,lb2:ub2,lb3:ub3,3,nspin) 
                                 d(n*epsxc)/dGradD  (Notation: exc = n*epsxc)
 ************************ UNITS ************************************
 Distances in atomic units (Bohr).
 Densities in atomic units (electrons/Bohr**3)
 Energy unit depending of parameter EUnit below (currently Rydberg)
 Stress in EUnit/Bohr**3

 Arrays dens, Vxc, and dVxcdD may be alternatively 
 allocated and initialized with indexes beginning in 0 or 1, 
 or even use a single-index array, e.g.:
   real(grid_p),allocatable :: dens(:,:), Vxc(:,:)
   myMesh(1:3) = myBox(2,:) - myBox(1,:) + 1
   myPoints = myMesh(1)*myMesh(2)*myMesh(3)
   allocate( dens(myPoints,nSpin), Vxc(myPoints,nSpin) )
   iPoint = 0
   do i3 = myBox(1,3),myBox(2,3)
   do i2 = myBox(1,2),myBox(2,2)
   do i1 = myBox(1,1),myBox(2,1)
     iPoint = iPoint + 1
     do iSpin = 1,nSpin
       dens(iPoint,iSpin) = (spin)density at point (i1,i2,i3)
     end do
   end do
   end do
   end do
 but the call to cellXC must still be as given above, i.e. the
 arguments lb1,ub1,... must be the lower and upper bounds of the
 mesh box stored by each processor (not the allocated array bounds).
 However, the arrays size MUST be (ub1-lb1+1)*(ub2-lb2+1)*(ub3-lb3+1)
 The processor mesh boxes must not overlap, and they must cover the 
 entire unit cell mesh. This is NOT checked inside cellXC.

 ********* BEHAVIOUR ***********************************************
 - Stops and prints a warning if functl is not one of LDA, GGA, or VDW
 - The output values of Ex, Ec, Dx, Dc, and stress, are integrals over
   the whole unit cell, not over the mesh box of the local processor
 - Since the exchange and correlation part is usually a small fraction
   of a typical electronic structure calculation, this routine has
   been coded with emphasis on simplicity and functionality, not in
   efficiency.