anchor.diff:  
   *      **************UP*******************
          chiup = agrup/nup**(4.0d0/3.0d0)
  -       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 ((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
  +          !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
  
positive.code: *****************************************************************
            sd    = 1.0d0/(n_thrd*n) 
            chi   = two_thrd*agr*sd 
            chi2  = chi*chi
            chiSQ = dsqrt(1.0d0+chi2)   
  
            K  = 6.0d0*beta*dlog(chi+chiSQ)
            F1 = chi2/(1.0d0+chi*K)
            xe = -n_thrd*(lda_c+beta*F1)/two_thrd
            F2 = (2.0d0 + chi*K-(chi2)*6.0d0*beta
       &             /chiSQ)
       &           /((1.0d0+chi*K)*(1.0d0+chi*K))
            fdnx = -(n_thrd/two_thrd)*dble(4.0d0/3.0d0)
       &               *(lda_c+beta*(F1-chi2*F2))
            fdagrx = -beta*chi*F2 
            if ((fdnx-xe).gt.0.0d0) then
               xe = -n_thrd*(lda_c)/two_thrd
               fdnx = -(n_thrd/two_thrd)*dble(4.0d0/3.0d0)*(lda_c)
               fdagrx = 0.0d0
  
negative.code:       parameter (DNS_CUT        =      1.0d-18)
  
        real*8 tolrho,minagr
        parameter (tolrho=2.0e-11,minagr=1.0e-12)
        !parameter (tolrho=2.0e-8,minagr=1.0e-10)
  
  *     ***LYP parameters****************
        real*8 a,b,c,d
  