      subroutine at_shell(xp,yp,zp,dx,dy,dz,value)
      implicit none
      double precision xp,yp,zp,value
      double precision rho,drho(3),d2rho(6)
      double precision elfd,delf(3),d2elf(6),stepc
      double precision dx,dy,dz
      integer level
      stepc=dx*dx+dy*dy+dz*dz
      stepc=dsqrt(stepc)
      level=0
      do while (stepc.ge.1.d-3)
         elfd=1.0d0
         value=0.0d0
         do while(elfd.ge.value)
            call elf_eval(rho,drho,d2rho,elfd,delf,d2elf,xp,yp,zp,level)
            if(elfd.ge.value) then
               value=elfd
               xp=xp+dx
               yp=yp+dy
               zp=zp+dz
            endif
         enddo
         xp=xp-2.0d0*dx
         yp=yp-2.0d0*dy
         zp=zp-2.0d0*dz
         dx=0.38197d0*dx
         dy=0.38197d0*dy
         dz=0.38197d0*dz
         stepc=0.38197d0*stepc
      enddo
      value=elfd
      return
      end
