anchor.diff:  ******************************************************************
   *     *******calc. becke exchange energy density, fnx, fdnx*******
   *****************************************************************
  -       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 
  +          sd    = 1.0d0/(n_thrd*n) 
  +          chi   = two_thrd*agr*sd 
  +          if ((n.gt.tolrho)  .and.
  +     >        (agr.gt.minagr).and.
  +     >        (chi.lt.1.0d0)) then
  +             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 
  +          else
  +             xe = -n_thrd*(lda_c)/two_thrd
  +             fdnx = -(n_thrd/two_thrd)*dble(4.0d0/3.0d0)*(lda_c)
  +             fdagrx = 0.0d0
  +          end if
          end if
   
   
  
positive.code: 
  *      **************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 
  *      ************END UP*****************
  
  
  *      *************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
  
  
  *      ***********END DOWN****************
  
         xe = (xeup*nup + xedn*ndn)/n
  
  *      *******end excange part************
         end if
  
negative.code: 
  !$OMP DO
        do i=1,n2ft3d
         n      = rho_in(i) + ETA
         agr    = agr_in(i)
         n_thrd = n**thrd
  
         if (rho_in(i).lt.DNS_CUT) then
           xe     = 0.0d0
  