function deriv(r,x,x2,spin,i,imin,imax)
c  deriv calculates the derivative of x with respect to r, at
c  position r(i).
c  input r        : radial coordinate (array)
c  input x        : function to be derivated  (array)
c  input x2       : second derivative of x (array)
c  input spin     : -1 if spin-down, 1 if spin-up, 0 otherwise
c  input i        : counter, i.e. current position = r(i)
c  input imin     : min value of i
c  input imax     : max value of i
      implicit logical (a-z)
      double precision deriv,r,x,x2
      double precision h,d,dforw,dback,dforup,dfordn,dbacup,dbacdn
      integer spin,i,imin,imax
      dimension r(561),x(2,561),x2(2,561)
      if (i.eq.imin) then
         h = (r(i+1)-r(i))/2.d0
         if (spin.eq.-1) then
            d = x(2,i)
            call splint(r,x,x2,-1,imin,imax,r(i)+h,dforw)
         elseif (spin.eq.1) then
            d = x(1,i)
            call splint(r,x,x2,1,imin,imax,r(i)+h,dforw)
         else
            d = (x(1,i) + x(2,i))
            call splint(r,x,x2,1,imin,imax,r(i)+h,dforup)
            call splint(r,x,x2,-1,imin,imax,r(i)+h,dfordn)
            dforw = dforup + dfordn
         endif   
         deriv = (dforw-d)/h
      elseif (i.eq.imax) then
         h = (r(i)-r(i-1))/2.d0
         if (spin.eq.-1) then
            call splint(r,x,x2,-1,imin,imax,r(i)-h,dback)
            d = x(2,i)
         elseif (spin.eq.1) then
            call splint(r,x,x2,1,imin,imax,r(i)-h,dback)
            d = x(1,i)
         else
            call splint(r,x,x2,1,imin,imax,r(i)-h,dbacup)
            call splint(r,x,x2,-1,imin,imax,r(i)-h,dbacdn)
            dback = dbacup + dbacdn
            d = x(1,i) + x(2,i)
         endif   
         deriv = (d-dback)/h
      else
         h = (r(i)-r(i-1))/2.d0
         if (spin.eq.-1) then
            call splint(r,x,x2,-1,imin,imax,r(i)+h,dforw)
            call splint(r,x,x2,-1,imin,imax,r(i)-h,dback)
         elseif (spin.eq.1) then
            call splint(r,x,x2,1,imin,imax,r(i)+h,dforw)
            call splint(r,x,x2,1,imin,imax,r(i)-h,dback)
         else
            call splint(r,x,x2,-1,imin,imax,r(i)+h,dfordn)
            call splint(r,x,x2,-1,imin,imax,r(i)-h,dbacdn)
            call splint(r,x,x2,1,imin,imax,r(i)+h,dforup)
            call splint(r,x,x2,1,imin,imax,r(i)-h,dbacup)
            dback = dbacup + dbacdn
            dforw = dforup + dfordn
         endif   
         deriv = (dforw-dback)/(2.d0*h)
      endif
      return
      end