anchor.diff:        parameter (ETA            =      1.0d-20)
   
         real*8 tolrho,minagr
  -      parameter (tolrho=2.0e-8,minagr=1.0e-10)
  +      parameter (tolrho=2.0e-11,minagr=1.0e-12)
   ****** Becke constants *****************************
         real*8 lda_c,beta
         parameter (beta = 0.0042d0)
  
positive.code: 
  *      *************DOWN******************
         chidn = agrdn/ndn**(4.0d0/3.0d0)
         chidn2 = chidn*chidn
         chidnSQ = dsqrt(1.0d0+chidn2)
  
         Kdn = 6.0d0*beta*dlog(chidn+chidnSQ)
         F1dn = chidn2/(1.0d0 + chidn*Kdn)
         xedn = -ndn**thrd*(lda_c + beta*F1dn)
         F2dn = (2.0d0 + chidn*Kdn
       &           - 6.0d0*beta*chidn2/chidnSQ)
       &           /(1.0d0+chidn*Kdn)**2.0d0
         fdnxdn = -ndn**(thrd)*(4.0d0/3.0d0)
       &             *(lda_c+beta*(F1dn-chidn2*F2dn))
         fdagrxdn = -beta*chidn*F2dn
         if ((fdnxdn-xedn).gt.0.0d0) then
            xedn     = -ndn**thrd*(lda_c)
            fdnxdn   = -ndn**(thrd)*(4.0d0/3.0d0)*(lda_c)
            fdagrxdn = 0.0d0
  
negative.code: 
  *      **************UP*******************
         chiup = agrup/nup**(4.0d0/3.0d0)
         if ((nup.gt.tolrho)  .and.
       >     (agrup.gt.minagr)) then
            chiup2 = chiup*chiup
            chiupSQ = dsqrt(1.0d0+chiup2)
  
            Kup = 6.0d0*beta*dlog(chiup+chiupSQ)
            F1up = chiup2/(1.0d0 + chiup*Kup)
            xeup = -nup**thrd*(lda_c + beta*F1up)
            F2up = (2.0d0 + chiup*Kup 
       &           - 6.0d0*beta*chiup2/chiupSQ)
       &           /(1.0d0+chiup*Kup)**2.0d0
            fdnxup = -nup**(thrd)*(4.0d0/3.0d0)
       &             *(lda_c+beta*(F1up-chiup2*F2up))
            fdagrxup = -beta*chiup*F2up 
           
            !if ((fdnxup-xeup).gt.0.0d0) then
            !   xeup     = -nup**thrd*(lda_c)
            !   fdnxup   = -nup**(thrd)*(4.0d0/3.0d0)*(lda_c)
            !   fdagrxup = 0.0d0
            !end if
         else
            xeup     = -nup**thrd*(lda_c)
            fdnxup   = -nup**(thrd)*(4.0d0/3.0d0)*(lda_c)
            fdagrxup = 0.0d0
  