anchor.diff:  *   This routine does small n (larger gradient) correction in which
   * the function is slowly switched to LDA when (fdnx-xe)>0.
   *
  -      subroutine Becke_smalln_correction(n,n_thrd,agr,beta,lda_c,
  +      subroutine Becke_smalln_correction(n,n_thrd,agr,fac,beta,lda_c,
        >                                   chi,chi2,chiSQ,K,F1,F2,
        >                                   xe,fdnx,fdagrx)
         implicit none
  -      real*8 n,n_thrd,agr,beta,lda_c,chi,chi2,chiSQ,K,F1,F2
  +      real*8 n,n_thrd,agr,fac,beta,lda_c,chi,chi2,chiSQ,K,F1,F2
         real*8 xe,fdnx,fdagrx
   
   *     **** constants ****
  -      real*8 pi,thrd,two_thrd
  +      real*8 pi,thrd,frthrd,two_thrd
         parameter (pi=3.14159265358979311599d0)
         parameter (thrd=1.0d0/3.0d0)
  +      parameter (frthrd=4.0d0/3.0d0)
         parameter (two_thrd=1.25992104989487319066d0)  ! two_thrd = 2**thrd
   
   *     **** local variables ****
  +      real*8 nf,nf_thrd,agrf
         real*8 x,s,dsdx,tmp1,tmp2,dchi,dF1,dF2
         real*8 dfdnxdagr,ddf0,dxdn,dxdagr,dsdn,dsdagr,xeold
   
  -      x = 0.85d0*(fdnx-xe)/(n_thrd)
  +      nf = n/fac
  +      nf_thrd = n_thrd/(fac**thrd)
  +      agrf = agr/fac
  +
  +      x = 0.85d0*(fdnx-xe)/(nf_thrd)
         s = dtanh(x)
         dsdx = 0.0d0
         if (x<100.0d0) then
  
positive.code:              fdagrxup = -beta*chiup*F2up 
           
               if ((fdnxup-xeup).gt.0.0d0) then
               call Becke_smalln_correction(nup,nup_thrd,agrup,beta,lda_c,
       >                               chiup,chiup2,chiupSQ,Kup,F1up,F2up,
       >                               xeup,fdnxup,fdagrxup)
               end if
  
negative.code:      &             *(lda_c+beta*(F1dn-chidn2*F2dn))
               fdagrxdn = -beta*chidn*F2dn
               if ((fdnxdn-xedn).gt.0.0d0) then
               call Becke_smalln_correction(ndn,ndn_thrd,agrdn,1.0d0,beta,
       >                               lda_c,
       >                               chidn,chidn2,chidnSQ,Kdn,F1dn,F2dn,
       >                               xedn,fdnxdn,fdagrxdn)
               end if 
            end if 
  *         ***********END DOWN****************
  
             xe = (xeup*nup + xedn*ndn)/n
  
         end if
  *      *******end excange part************
  