anchor.diff:  
         return
         end
  -
  - 
  +ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  +      subroutine LSDT2(a,a1,b1,b2,b3,b4,srs,ec,ec_rs)
  +      real*8 a,a1,b1,b2,b3,b4,srs,ec,ec_rs
  +      real*8 q0,q1,q1p,qd,ql
  +      q0= -2.0d0*a*(1.0d0+srs*srs)
  +      q1= 2.0d0*a*srs*(b1+srs*(b2+srs*(b3+srs*b4)))
  +      q1p= a*((b1/srs)+2.0d0*b2+srs*(3.0d0*b3+srs*4.0d0*b4))
  +      qd=1.0d0/(q1*q1+q1)
  +      ql= -dlog(qd*q1*q1)
  +      ec= q0*ql
  +      ec_rs= -2.0d0*a*a1*ql-q0*q1p*qd
  +      return
  +      end
  +      
  +      
  + 
  +
  + 
  
positive.code:          n     = 2.0d0*ndn
           agr   = 2.0d0*agrdn
  
           n_onethird = (3.0d0*n/pi)**onethird
           ex_lda     = -0.75d0*n_onethird
  
           kf = (3.0d0*pi*pi*n)**onethird
           s  = agr/(2.0d0*kf*n)
           P0 = 1.0d0 + (MU/KAPPA)*s*s
  
           F   = (1.0d0 + KAPPA - KAPPA/P0)
  
negative.code: 
        pi = 4.0d0*datan(1.0d0)
        rs_scale = (0.75d0/pi)**onethird
        fdnx_const = -3.0d0/(16.0d0*pi)
  
        
        do i=1,n2ft3d
           nup     = dn_in(i,1)
           agrup   = agr_in(i,1)
  
   
           ndn     = dn_in(i,2)
           agrdn   = agr_in(i,2)
  
           
  