anchor.diff:           n     = 2.0d0*ndn
            agr   = 2.0d0*agrdn
   
  -         n_onethird = (3.0d0*n/pi)**onethird
  +         n13 = n**onethird
  +         n13 = max(n13,1.e-14)
  +         n_onethird = n13*(3.0d0/pi)**onethird
            ex_lda     = -0.75d0*n_onethird
   
  -         kf = (3.0d0*pi*pi*n)**onethird
  -         s  = agr/(2.0d0*kf*n)
  +         kf = (3.0d0*pi*pi)**onethird * n13
  +         s  = agr/(2.0d0*kf*n13*n13*n13)
            P0 = 1.0d0 + (MU/KAPPA)*s*s
   
            F   = (1.0d0 + KAPPA - KAPPA/P0)
  
positive.code: c     **** other constants ****
        real*8 onethird,fourthird,fivethird,onesixthm
        real*8 twothird,sevensixthm
        real*8 onethirdm
        parameter (onethird=1.0d0/3.0d0)
        parameter (onethirdm=-1.0d0/3.0d0)
        parameter (twothird=2.0d0/3.0d0)
  
negative.code: *     **** Density cutoff parameter ****
        real*8 DNS_CUT,ETA,ETA2,alpha_zeta,alpha_zeta2
        parameter (DNS_CUT = 1.0d-20)
        parameter (ETA=1.0d-12)
        parameter (ETA2=1.0d-6)
        parameter (alpha_zeta=(1.0d0-ETA2))
        parameter (alpha_zeta2=(1.0d0-ETA2))
  
  