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:      >                               x_parameter,c_parameter,
       >                               xce,fn,fdn,fdtau)
        implicit none
  
        integer n2ft3d
        real*8 rho_in(n2ft3d,2),agr_in(n2ft3d,3),tau_in(n2ft3d,2)      !*** input  ****
        real*8 x_parameter,c_parameter
        real*8 xce(n2ft3d),fn(n2ft3d,2),fdn(n2ft3d,3),fdtau(n2ft3d,2)  !*** output ****
  
  ********* local declarations******
        integer i
        real*8 rhouu_in(n2ft3d,2),euu0(n2ft3d),fuu0n(n2ft3d,2)
        real*8 rhodd_in(n2ft3d,2),edd0(n2ft3d),fdd0n(n2ft3d,2)
        real*8 eud0(n2ft3d),fud0n(n2ft3d,2)
        real*8 nup,agrup,tauup
        real*8 ndn,agrdn,taudn
        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,GGup,dGup_dxup,dGup_dzup
        real*8 gammadn,xrdn,zrdn,GGdn,dGdn_dxdn,dGdn_dzdn
        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,fuu0_nupc,euuc
        real*8 fuu_nupc,fuu_dnupc,fuu_dtauupc
        real*8 edd0c,fdd0_ndnc,eddc
        real*8 fdd_ndnc,fdd_dndnc,fdd_dtaudnc
        real*8 eud0c,fud0_nupc,fud0_ndnc,eudc
        real*8 fud_nupc,fud_dnupc,fud_dtauupc
        real*8 fud_ndnc,fud_dndnc,fud_dtaudnc
        real*8 eupc,fnupc,fdnupc,fdtauupc
        real*8 ednc,fndnc,fdndnc,fdtaudnc
        real*8 ex,ec,n
  ********* constants ***************
        real*8 thrd,frthrd,fvthrd,etthrd
        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-10)
        parameter (ETA            =      1.0d-20)
  
  ********* VS98 constants ******
        real*8 cf,alpha,aa,bb,cc,dd,ee,ff
        parameter (cf = 9.115599720d0)         !***cf = 3.0d0/5.0d0*(6.0d0*pi*pi)**(2.0d0/3.0d0))***
  
        do i=1,n2ft3d
           rhouu_in(i,1) = rho_in(i,1)
           rhouu_in(i,2) = 0.0d0
           rhodd_in(i,1) = 0.0d0 
           rhodd_in(i,2) = rho_in(i,2)
        end do
        call gen_PW92_BW_unrestricted(n2ft3d,rhouu_in,euu0,fuu0n) 
        call gen_PW92_BW_unrestricted(n2ft3d,rhodd_in,edd0,fdd0n) 
        call gen_PW92_BW_unrestricted(n2ft3d,rho_in,eud0,fud0n) 
       
        do i=1,n2ft3d
  
          nup        = rho_in(i,1) + ETA
          agrup      = agr_in(i,1)
          tauup      = tau_in(i,1)
          euu0c      = euu0(i)
          fuu0_nupc  = fuu0n(i,1)
  
          ndn        = rho_in(i,2) + ETA
          agrdn      = agr_in(i,2)
          taudn      = tau_in(i,2)
          edd0c      = edd0(i)
          fdd0_ndnc  = fdd0n(i,2)
  
          eud0c      = eud0(i)
          fud0_nupc  = fud0n(i,1)
          fud0_ndnc  = fud0n(i,2)
  
          n          = nup + ndn
       
          agrup2  = agrup*agrup
          inv_nup = 1.0d0/nup
          nup_13  = nup**thrd
          nup_43  = nup_13*nup
          nup_53  = nup_43*nup_13
          nup_83  = nup_53*nup
  
          xup         = agrup2/nup_83
          dxup_dnup   = -etthrd*xup*inv_nup
          dxup_dagrup = 2.0d0*agrup/nup_83
          zup         = tauup/nup_53 - cf
          dzup_dnup   = -fvthrd*tauup/nup_83
          dzup_dtauup = 1.0d0/nup_53
          
          agrdn2  = agrdn*agrdn
          inv_ndn = 1.0d0/ndn
          ndn_13  = ndn**thrd
          ndn_43  = ndn_13*ndn
          ndn_53  = ndn_43*ndn_13
          ndn_83  = ndn_53*ndn
  
          xdn         = agrdn2/ndn_83
          dxdn_dndn   = -etthrd*xdn*inv_ndn
          dxdn_dagrdn = 2.0d0*agrdn/ndn_83
          zdn         = taudn/ndn_53 - cf
          dzdn_dndn   = -fvthrd*taudn/ndn_83
          dzdn_dtaudn = 1.0d0/ndn_53
  ********* VS98 Exchange *********
          aa    = -9.800683d-1
          bb    = -3.556788d-3
          cc    =  6.250326d-3
  
negative.code:         ee    = -1.269420d-3
          ff    =  1.296118d-3
          alpha =  0.00304966d0
  *       ***** alpha-beta *****
          rs     =  (0.75d0/(pi*n))**thrd
          drs_dn = -rs/(3.0d0*n)
  
          zeta       = (nup - ndn)/n
          dzeta_dnup = ( 1.0d0 - zeta)/n
          dzeta_dndn = (-1.0d0 - zeta)/n
  
          call gen_PW91_c_rz(tol,rs,zeta,eud0c,deud0c_drs,deud0c_dzeta)
          eud1c = eud0c - (euu0c*nup + edd0c*ndn)/n
          fud1c = n*eud1c
  
          xud     = xup + xdn
          zud     = zup + zdn
  
          gammaud = 1.0d0 + alpha*(xud + zud)
          xrud    = xud/gammaud
          zrud    = zud/gammaud
  