anchor.diff:           fnxdn  = fourthird*(exdn - ex_lda*Fs*s)
            fdnxdn = fdnx_const*Fs
   
  -         n = nup + ndn
  -         ex = (exup*nup+ exdn*ndn)/n
  -         
  +         n = nup+ndn
  +
  +         ex = (exup*nup+ exdn*ndn)/ n
  +                  
   c        *******************************************************************
   c        ***** calculate polarized correlation energies and potentials ***** 
   c        *******************************************************************
            agr   = agr_in(i,3)
  -         n  = max(n,1.d-30)
  +
            zet = (nup-ndn)/n
  +c         if (zet.gt.0.0d0) zet = zet - ETA2
  +c         if (zet.lt.0.0d0) zet = zet + ETA2
  +c        if (dabs(dn_in(i,2)).gt.DNS_CUT) zet_nup =  2*ndn/n**2
  +c        if (dabs(dn_in(i,1)).gt.DNS_CUT) zet_ndn = -2*nup/n**2
  +c        if (dabs(dn_in(i,2)).gt.DNS_CUT) zet_nup =  2*ndn/n
  +c        zet_nup =  2*ndn/n
  +c        zet_ndn = -2*nup/n
            zet_nup = -(zet - 1.0d0)
            zet_ndn = -(zet + 1.0d0)
            zetpm_1_3 = (1.0d0+zet*alpha_zeta)**onethirdm
  
positive.code: *     **** Density cutoff parameter ****
        real*8 DNS_CUT,ETA,ETA2,alpha_zeta,alpha_zeta2
        parameter (DNS_CUT = 1.0d-20)
        parameter (ETA=1.0d-12)
        parameter (ETA2=1.0d-6)
        parameter (alpha_zeta=(1.0d0-ETA2))
        parameter (alpha_zeta2=(1.0d0-ETA2))
  
  c     ***** revPBE GGA exchange constants ******
        real*8 MU,KAPPA
        parameter (MU    = 0.2195149727645171d0)
  c      parameter (KAPPA = 0.8040000000000000d0)
        parameter (KAPPA = 1.2450000000000000d0)
   
  c     ****** PBE96 GGA correlation constants ******
  
negative.code:          n     = 2.0d0*nup
           agr   = 2.0d0*agrup
  
           n_onethird = (3.0d0*n/pi)**onethird
           ex_lda     = -0.75d0*n_onethird
  
           kf = (3.0d0*pi*pi*n)**onethird
           s  = agr/(2.0d0*kf*n)
           P0 = 1.0d0 + (MU/KAPPA)*s*s
  
           F   = (1.0d0 + KAPPA - KAPPA/P0)
  