anchor.diff:       &               *(lda_c+beta*(F1-chi2*F2))
             fdagrx = -beta*chi*F2 
             if ((fdnx-xe).gt.0.0d0) then
  -              call Becke_smalln_correction(n,n_thrd,agr,beta,lda_c,
  +              call Becke_smalln_correction(n,n_thrd,agr,2.0d0,beta,
  +     >                                   lda_c,
        >                                   chi,chi2,chiSQ,K,F1,F2,
        >                                   xe,fdnx,fdagrx)
             end if
  
positive.code:       parameter (c = 0.2533d0)
        parameter (d = 0.349d0)
  *     ***Becke88 parameters*************
        real*8 beta,Cf
        parameter (beta = 0.0042d0)
        parameter (Cf =2.87123400018819108225d0)     ! Cf = (3/10)*(3*pi*pi)**(2/3)
        
        real*8 thrd
        real*8 twthrd,frthrd,fvthrd,snthrd,etthrd
  
negative.code: *   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,fac,beta,lda_c,
       >                                   chi,chi2,chiSQ,K,F1,F2,
       >                                   xe,fdnx,fdagrx)
        implicit none
        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,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
  
        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
  