anchor.diff:  *  to n, |grad n|, tau.
   
   *
  -*   Entry - n2ft3d   : number of grid points
  -*           rho_in(*) :  density (nup+ndn)
  -*           agr_in(*): |grad rho_in|
  -*           tau_in(*): tau
  -*           x_parameter: scale parameter for exchange
  -*           c_parameter: scale parameter for correlation
  +*     Entry - n2ft3d   : number of grid points
  +*             rho_in(*) :  density (nup+ndn)
  +*             agr_in(*): |grad rho_in|
  +*             tau_in(*): tau
  +*             x_parameter: scale parameter for exchange
  +*             c_parameter: scale parameter for correlation
   *
   *     Exit  - xce(n2ft3d) : VS98 exchange correlation energy density
   *             fn(n2ft3d)  : d(n*xce)/dn
  
positive.code: *  to nup, ndn, |grad nup|, |grad ndn|, tauup, taudn.
  
  *
  *   Entry - n2ft3d   : number of grid points
  *           rho_in(*,2) :  density (nup and ndn)
  *           agr_in(*,3): |grad rho_in| (nup, ndn and n)
  *           tau_in(*,2): tau (nup and ndn)
  *           x_parameter: scale parameter for exchange
  *           c_parameter: scale parameter for correlation
  *
  *     Exit  - xce(n2ft3d) : VS98 exchange correlation energy density
  *             fn(n2ft3d,2)  : d(n*xce)/dnup, d(n*xce)/dndn
  
negative.code:      >                               x_parameter,c_parameter,
       >                               xce,fn,fdn,fdtau)
        implicit none
  *     ***** input *****
        integer n2ft3d
        real*8 rho_in(n2ft3d,2),agr_in(n2ft3d,3),tau_in(n2ft3d,2)      
        real*8 x_parameter,c_parameter
  *     ***** output *****
        real*8 xce(n2ft3d),fn(n2ft3d,2),fdn(n2ft3d,3),fdtau(n2ft3d,2)  
  *     ***** local declarations *****
        integer i
        real*8 nup,agrup,tauup
        real*8 ndn,agrdn,taudn
        real*8 n,agr,tau
        real*8 inv_n,n_13,n_43,n_53,n_83,agr2
        real*8 x,dx_dn,dx_dagr,z,dz_dn,dz_dtau
        real*8 gamma,xr,zr,GG,dG_dx,dG_dz
        real*8 inv_nup,nup_13,nup_43,nup_53,nup_83,agrup2
        real*8 inv_ndn,ndn_13,ndn_43,ndn_53,ndn_83,agrdn2
        real*8 xup,dxup_dnup,dxup_dagrup,zup,dzup_dnup,dzup_dtauup
        real*8 xdn,dxdn_dndn,dxdn_dagrdn,zdn,dzdn_dndn,dzdn_dtaudn
        real*8 xud,zud
        real*8 gammaup,xrup,zrup
        real*8 gammadn,xrdn,zrdn
        real*8 gammaud,xrud,zrud,GGud,dGud_dxud,dGud_dzud
        real*8 dGud_dnup,dGud_dagrup,dGud_dtauup
        real*8 dGud_dndn,dGud_dagrdn,dGud_dtaudn
        real*8 eupx,fnupx,fdnupx,fdtauupx
        real*8 ednx,fndnx,fdndnx,fdtaudnx
        real*8 euu0c,deuu0c_dnup
        real*8 euuc,dfuuc_dnup,dfuuc_dagrup,dfuuc_dtauup
        real*8 edd0c,dedd0c_dndn
        real*8 eddc,dfddc_dndn,dfddc_dagrdn,dfddc_dtaudn
        real*8 rs,drs_dn,zeta,dzeta_dnup,dzeta_dndn
        real*8 eud0c,deud0c_drs,deud0c_dzeta
        real*8 eudc,eud1c,fud1c,dfud1c_dnup,dfud1c_dndn
        real*8 dfudc_dnup,dfudc_dagrup,dfudc_dtauup
        real*8 dfudc_dndn,dfudc_dagrdn,dfudc_dtaudn
        real*8 fnupc,fdnupc,fdtauupc
        real*8 fndnc,fdndnc,fdtaudnc
        real*8 ex,ec
  ********* constants ***************
        real*8 pi,thrd,frthrd,fvthrd,etthrd
        parameter (pi     =  3.14159265358979311599d0)
        parameter (thrd   = 1.0d0/3.0d0)
        parameter (frthrd = 4.0d0/3.0d0)
        parameter (fvthrd = 5.0d0/3.0d0)
        parameter (etthrd = 8.0d0/3.0d0)
  ********* density cutoff parameters******
        real*8 tol,ETA
        parameter (tol = 1.0d-18)
        parameter (ETA            =      1.0d-20)
  
  ********* VS98 constants ******
        real*8 cf,aa,bb,cc,dd,ee,ff,alpha
        parameter (cf = 9.115599720d0)        
  c     cf = 3.0d0/5.0d0*(6.0d0*pi*pi)**(2.0d0/3.0d0))
  
        do i=1,n2ft3d
          nup        = 0.50d0*rho_in(i,1) + ETA
          agrup      = 0.50d0*agr_in(i,1) + ETA
          tauup      = 0.50d0*tau_in(i,1) + ETA
          ndn        = 0.50d0*rho_in(i,2) + ETA
          agrdn      = 0.50d0*agr_in(i,2) + ETA
          taudn      = 0.50d0*tau_in(i,2) + ETA
  
  *       ***** VS98 Exchange *****
          aa    = -9.800683d-1
          bb    = -3.556788d-3
          cc    =  6.250326d-3
  