subroutine laplac(r,gradup,graddn,i,imin,imax,done,
     j     lapl,laplup,lapldn)
c  calculates laplace(density), i.e. (del^2)(density)
c  spherical coordinates => (del^2)phi = 1/r^2*d/dr[r^2*d(phi)/dr]
c  input r              : radial coordinate (array)
c  input gradup         : grad(density), spin up (array)
c  input graddn         : grad(density), spin down (array)
c  input i              : counter, i.e. current position = r(i)
c  input imin           : min value of i
c  input imax           : max value of i
c  in/output done       : true if r^2*grad(density) is splined
c  output lapl          : laplace(density)
c  output laplup,lapldn : laplace(density), spin up and down
      implicit logical (a-z)
      logical done
      double precision r,gradup,graddn,lapl,laplup,lapldn,lapl2,temp
      double precision rj2,r2,deriv
      integer i,imin,imax,j
      common/ggalap/temp,lapl2
      dimension r(561),gradup(561),graddn(561)
      dimension lapl2(2,561),temp(2,561)
      if (.not.done) then
         do 10 j = imin,imax
            rj2 = r(j)**2
            temp(1,j) = rj2*gradup(j)
            temp(2,j) = rj2*graddn(j)
 10      continue
         call spline(r,imin+1,imax,temp,lapl2)
         done = .true.
      endif 
      if (i.eq.imin) then
         call laplac2(r,gradup,graddn,i+1,imin,imax,done,lapl,laplup,
     j        lapldn)
      else
         r2 = r(i)**2
         laplup = (deriv(r,temp,lapl2,1,i,imin,imax))/r2
         lapldn = (deriv(r,temp,lapl2,-1,i,imin,imax))/r2
         lapl = laplup + lapldn
      endif
      return
      end