      subroutine attractors(xp,yp,zp,i,j,vpoint,vcircle)
      use grid_module
      use orbitals_module
      use basins_module
      implicit none
      integer i,j,k,l,iacc
      integer vpoint,vcircle
      double precision elfd,rad
      double precision xp,yp,zp
      double precision scale
      double precision radi,radj,ah,rcircle
      double precision root(3),xold,yold,zold
      double precision root_max
      double precision dist,stepc,eps
      logical found1
      common/accuracy/iacc
      eps=1.0d-6
      stepc=0.17d0
      call attract_search(xp,yp,zp,root,elfd,eps)
      xs(attract_label+1)=xp
      ys(attract_label+1)=yp
      zs(attract_label+1)=zp
      xold=xp
      yold=yp
      zold=zp
      radi=(xp-x(i))**2+(yp-y(i))**2+(zp-z(i))**2
*
*      find the largest eigenvalue
*
      root_max=root(1)
      do l=2,3
         if(root(l).gt.root_max) then
            root_max=root(l)
         endif
      enddo   
      found1=.true.
      if((root_max.gt.-5.d-6).and.(root_max.lt.1.0d-3)) then
         rad=1.0D6
         l=0
         do k=1,natoms
            dist=(x(i)-xat(k))**2+(y(i)-yat(k))**2+(z(i)-zat(k))**2
            if((dist.le.rad).and.(dist.gt.1.0d0)) then
               rad=dist
               l=k
            endif
         enddo
         rad=dsqrt(rad)
         radj=(xp-xat(l))**2+(yp-yat(l))**2+(zp-zat(l))**2
         ah=0.5d0*(rad+(radi-radj)/rad)
         rcircle=dsqrt(radi-ah*ah)
         xp=x(i)+ah*(xat(l)-x(i))/rad
         yp=y(i)+ah*(yat(l)-y(i))/rad
         zp=z(i)+ah*(zat(l)-z(i))/rad
         found1=.false.
      endif
      k=0
      dist=1.0d6
      do l=1,attract_label
         if(basin_type(l,1).ge.0) then
            scale=(xold-x(l))**2+(yold-y(l))**2+(zold-z(l))**2
            scale=dsqrt(scale)
            if(basin_type(l,1).eq.3) then
               scale=dabs(scale-core_radius(l))
            endif
            if(scale.lt.dist) then
               dist=scale
            endif
         endif
      enddo
      if(dist.le.2.0d-02) then
         found1=.true.
         return
      endif
      do l=1,attract_label
         dist=(xp-x(l))**2+(yp-y(l))**2+(zp-z(l))**2
         dist=dsqrt(dist)
         if(basin_type(l,1).eq.0) then
            scale=(xold-x(l))**2+(yold-y(l))**2+(zold-z(l))**2
            scale=dsqrt(scale)
            if(vpoint.eq.2) then
               if(scale.le.core_radius(l)) then
                  k=-l-j
               endif
            endif
         endif
         if(basin_type(l,1).ne.vcircle) then
            if(found1) then
               if(dist.le.stepc) then
                  k=l
               else if(dist.gt.1.0d3) then
                  k=-l-j
               endif
               scale=core_radius(l)
               if(vpoint.lt.0) scale=0.1d0
               if(dist.le.scale) then
                  k=-l-j
               endif
            endif
            if((vcircle.gt.0).and.(basin_type(l,1).eq.0)) then
               if((dist.le.core_radius(l)).and.
     .           (rcircle.le.core_radius(l))) then
               k=-l-j
               endif
            endif
         else if (.not.found1) then
            if(dist.le.stepc) then
               k=l
            endif
            if(dist.le.core_radius(l)) then
               k=-l-j
            endif
         else 
            if(dist.gt.1.0d3) then
               k=-l-j
            else if(dabs(core_radius(l)-dist).le.stepc) then
               if(root_max.gt.-5.0d-4) then
                  k=l
               else
                 k=0
               endif
            endif
            if(dist.le.stepc) then
               k=l
            endif
         endif
         if(basin_type(l,1).eq.1) then
            if(dist.le.1.5d0) then
               if(dabs(elfd-val_elf(l)).le.1.0d-3) then
                  k=l
               endif
            endif
         endif
      enddo
      if((k.eq.0).and.(root_max.lt.1.0d-3)) then
         attract_label=attract_label+1
         if(attract_label.gt.3000) then
            write(*,*) 'the number of basins excedds 3000',
     .                 'the wave function may be wrong'
            stop
         endif
         k=attract_label
         val_elf(attract_label)=elfd
         aux(attract_label,1)=root(1)
         aux(attract_label,2)=root(2)
         aux(attract_label,3)=root(3)
         x(attract_label)=xp
         y(attract_label)=yp
         z(attract_label)=zp
         if(root_max.gt.-5.0d-4) then
            basin_type(attract_label,1)=vcircle
            basin_type(attract_label,2)=0
            basin_type(attract_label,3)=i
            if(natoms.eq.1) then
               core_radius(attract_label)=dsqrt(radi)
            else
               core_radius(attract_label)=dsqrt(radi-ah*ah)+0.01
            endif
         else
            basin_type(attract_label,1)=vpoint
            basin_type(attract_label,2)=0
            core_radius(attract_label)=0.2d0
         endif   
      endif
      return
      end
