      subroutine assign_bas
      use basins_module
      use grid_module
      use orbitals_module
      implicit none
      integer i,j,k,l,m,irw,ix,iy,iz,npyz,neighb(6)
      integer nborder,iborder,lsave,iter,icircle
      integer ixa,iya,iza,indy
      integer*2 count
      integer low_limit
      integer npx,npy,npz,n_neighb(6)
      integer nchange,iacc
      integer new_attr
      integer i1,i2,i3
      integer l1,count2
      real t(2),elapsed,spent,thresh,oldf,spentt,expec
      double precision elf_min,elfd,step1,eps
      double precision xp,yp,zp
      double precision xold,yold,zold
c     double precision origin(3), edge(3),step(3),dv
      double precision dv
      double precision, dimension (:), allocatable :: xnew,ynew,znew
      double precision valp,delt,rmax
      double precision root(3)
      double precision pi,dist,rand
      character*1 molecule(40)
      character*40 filebas, fileelf, dumname
      logical func_type
      common /files/ molecule,fileelf,filebas
      common/ftype/func_type
      common/valence/thresh,low_limit
      common/accuracy/iacc
      equivalence (molecule(1),dumname)
      pi=dacos(-1.0d0)
      nsaddle=0
      allocate(xnew(maxbasin),ynew(maxbasin),znew(maxbasin))
      do i=1,maxbasin
         do j=1,maxbasin
            separatrices(i,j)=0
         enddo
      enddo
      count=0
      count2=0
      new_attr=0
      eps=1.0d-5
      origin(1)=xmin
      origin(2)=ymin
      origin(3)=zmin
      edge(1)=xmax-origin(1)
      edge(2)=ymax-origin(2)
      edge(3)=zmax-origin(3)
      npx=np(1)-1
      npy=np(2)-1
      npz=np(3)-1
      if(npx.eq.0) then
         step(1)=0.0d0
      else
         step(1)=edge(1)/npx
      endif
      if(npy.eq.0) then
         step(2)=0.0d0
      else
         step(2)=edge(2)/npy
      endif
      if(npz.eq.0) then
         step(3)=0.0d0
      else
         step(3)=edge(3)/npz
      endif
      dv=step(1)*step(2)*step(3)
      npts=np(1)*np(2)*np(3)
      npyz=np(2)*np(3)
*
*     assign attractor grid points
*
      do i=1,npts
         attract_code(i)=0
      enddo
      do i=1,attract_label
         call fill_bas(i)
      enddo
      nchange=0
      i=indx(1)
      if(attract_code(i).eq.0) then
         ix=(i-1)/npyz
         iy=(i-1-ix*npyz)/np(3)
         iz=i-1-ix*npyz-iy*np(3)
         xp=origin(1)+ix*step(1)
         yp=origin(2)+iy*step(2)
         zp=origin(3)+iz*step(3)
         irw=0
         do k=1,attract_label
            if(basin_type(k,1).ge.0) then
               valp=(x(k)-xp)**2+(y(k)-yp)**2+(z(k)-zp)**2
               valp=dsqrt(valp)
               if(valp.le.core_radius(k)) then
                 irw=k
               endif
            endif
         enddo
         if(irw.eq.0) then
            step1=0.02d0
            call attract_search(xp,yp,zp,root,elfd,eps)
            do k=1,attract_label
               if(basin_type(k,1).ge.0) then
                  valp=(x(k)-xp)**2+(y(k)-yp)**2+(z(k)-zp)**2
                  valp=dsqrt(valp)
                  if(valp.le.core_radius(k)) then
                    irw=k
                  endif
               endif
            enddo
         endif
         if(irw.eq.0) then
            delt=root(1)
            do k=2,3
               if(root(k).gt.delt) delt=root(k)
            enddo
            if(delt.lt.0.0d0) then
               if((xp.ge.xmin).and.(xp.le.xmax).and.
     .            (yp.ge.ymin).and.(yp.le.ymax).and.
     .            (zp.ge.zmin).and.(zp.le.zmax)) then
                  ixa=nint((xp-origin(1))/step(1))
                  iya=nint((yp-origin(2))/step(2))
                  iza=nint((zp-origin(3))/step(3))
                  indy=1+iza+np(3)*(iya+np(2)*ixa)
                  irw=attract_code(indy)
                  if(irw.eq.0) then
                     do k=1,attract_label
                        valp=(x(k)-xp)**2+(y(k)-yp)**2+(z(k)-zp)**2
                        valp=dsqrt(valp)
                        if(valp.le.1.0d-1) irw=k
                     enddo
                     if(irw.eq.0) then
                        new_attr=new_attr+1
                        xnew(new_attr)=xp
                        ynew(new_attr)=yp
                        znew(new_attr)=zp
                        irw=attract_label+new_attr
                        if(elfd.lt.dble(thresh)) then
                           new_attr=new_attr-1
                           irw=0
                        endif
                     endif
                  endif
               endif
            endif
         endif
         attract_code(i)=irw
      endif
      call xtime(t)
      elapsed=t(1)+t(2)-spent
      spent=t(1)+t(2)
      write(*,'(2x,''elapsed time in init step  '',f9.2,
     .             '' sec.'')') spent 
      nborder=0
      write(*,'(2x,''considered grid points '',i12)') low_limit
      l1=low_limit/10
      spentt=0.0
      j=1
      do j=2,low_limit
        if(mod(j,l1).eq.0) then
           call xtime(t)
           elapsed=t(1)+t(2)-spent
           spent=t(1)+t(2)
           spentt=spentt+elapsed
           expec=low_limit*spentt/j
           write(*,'(2x,i3,''% done'',f12.2,'' s. spent '',f12.2,
     .             '' s. expected '')') 10*j/l1,spentt,expec
        endif
        i=indx(j)
        ix=(i-1)/npyz
        iy=(i-1-ix*npyz)/np(3)
        iz=i-1-ix*npyz-iy*np(3)
        call neighbor(i,neighb)
        call neighcls(neighb,n_neighb,iborder)
*
*       iborder=0 new attractor
*
        if(iborder.eq.0) then
           if(attract_code(i).eq.0) then
              xp=origin(1)+ix*step(1)
              yp=origin(2)+iy*step(2)
              zp=origin(3)+iz*step(3)
              irw=0
              do k=1,attract_label
                 valp=(x(k)-xp)**2+(y(k)-yp)**2+(z(k)-zp)**2
                 valp=dsqrt(valp)
                 if(basin_type(k,1).eq.0) then
                    if(valp.le.core_radius(k)) then
                      irw=k
                    endif
                 else
                    delt=val_elf(k)-dble(elf(i))
                    if(dabs(delt).lt.5.0d-2) then
                       if(basin_type(k,1).eq.3) then
                          if((valp.ge.(0.8*core_radius(k))).and.
     .                   (valp.lt.(1.3*core_radius(k)))) then
                             irw=k
                          endif
                       else
                          if(valp.lt.0.6d0) then
                             irw=k
                          endif
                       endif
                    endif
                 endif
              enddo
              if(irw.eq.0) then
                step1=0.02d0
                call attract_search(xp,yp,zp,root,elfd,eps)
                delt=root(1)
                do k=2,3
                   if(root(k).gt.delt) then
                      delt=root(k)
                   endif
                enddo
                if(delt.gt.0.0d0) then 
                   xp=xp+0.1d0*(0.5d0-rand)*step(1)
                   yp=yp+0.1d0*(0.5d0-rand)*step(2)
                   zp=zp+0.1d0*(0.5d0-rand)*step(3)
                   call attract_search(xp,yp,zp,root,elfd,eps)
                endif
                do k=1,attract_label
                   icircle=iabs(basin_type(k,1))
                   valp=(x(k)-xp)**2+(y(k)-yp)**2+(z(k)-zp)**2
                   valp=dsqrt(valp)
                   if(icircle.ne.3) then
                      if(valp.le.core_radius(k)) then
                        irw=k
                      endif
                   else if (irw.eq.0) then
                      if((valp.ge.(0.8*core_radius(k))).and.
     .                   (valp.lt.(1.3*core_radius(k)))) then
                         irw=k
                      endif
                   endif
                enddo
              endif
              if(irw.eq.0) then
                 delt=root(1)
                 do k=2,3
                    if(root(k).gt.delt) delt=root(k)
                 enddo
                 if(delt.lt.-1.0d-5) then
                    if((xp.ge.xmin).and.(xp.le.xmax).and.
     .                 (yp.ge.ymin).and.(yp.le.ymax).and.
     .                 (zp.ge.zmin).and.(zp.le.zmax)) then
                        ixa=nint((xp-origin(1))/step(1))
                        iya=nint((yp-origin(2))/step(2))
                        iza=nint((zp-origin(3))/step(3))
                        indy=1+iza+np(3)*(iya+np(2)*ixa)
                        irw=attract_code(indy)
                        if(irw.eq.0) then
                          do k=1,attract_label
                             valp=(x(k)-xp)**2+(y(k)-yp)**2+(z(k)-zp)**2
                             valp=dsqrt(valp)
                             if(valp.le.1.0d-1) irw=k
                          enddo
                          if((irw.eq.0).and.(new_attr.gt.0)) then
                             do k=1,new_attr
                                valp=(xnew(k)-xp)**2+(ynew(k)-yp)**2+
     .                            (znew(k)-zp)**2
                                valp=dsqrt(valp)
                                if(valp.le.1.0d-1) irw=k+attract_label
                             enddo
                          endif
                          if(irw.eq.0) then
                             new_attr=new_attr+1
                             xnew(new_attr)=xp
                             ynew(new_attr)=yp
                             znew(new_attr)=zp
                             irw=attract_label+new_attr
                             if(elfd.lt.1.0d-1) then
                                new_attr=new_attr-1
                                irw=0
                             endif
                          endif
                       endif
                    else
                       new_attr=new_attr+1
                       xnew(new_attr)=xp
                       ynew(new_attr)=yp
                       znew(new_attr)=zp
                       irw=attract_label+new_attr
                       if(elfd.lt.1.0d-1) then
                          new_attr=new_attr-1
                          irw=0
                       endif
                    endif
                 else if(delt.lt.1.0d-4) then
                    do k=1,attract_label
                        elf_min=dabs(elfd-val_elf(k))
                        if(elf_min.le.1.0d-3) then
                           m=basin_type(k,3)
                           iter=basin_type(k,4)
                           step1=(xp-x(m))**2+(yp-y(m))**2+(zp-z(m))**2-
     .                           (x(k)-x(m))**2-(y(k)-y(m))**2-
     .                           (z(k)-z(m))**2
                           if(m.ne.iter) then
                               step1=step1+(xp-x(iter))**2+
     .                             (yp-y(iter))**2+(zp-z(iter))**2-
     .                             (x(k)-x(iter))**2-(y(k)-y(iter))**2-
     .                             (z(k)-z(iter))**2
                           endif
                           if(dabs(step1).lt.1.0d-4) irw=k
                        endif
                     enddo
                 endif
              endif
              attract_code(i)=irw
           endif
        endif
*
*       iborder=1 regular wandering point
*
        if(iborder.eq.1) then
           attract_code(i)=n_neighb(1)
        endif
*
*       iborder>1 separatrix point
*
        if(iborder.gt.1) then
           oldf=elf(i)
           nborder=nborder+1
           i1=n_neighb(1)
           i2=n_neighb(2)
           if(iborder.eq.3) then
              i3=n_neighb(3)
           else
              i3=0
           endif
           if(separatrices(i1,i2).eq.0) then
              nsaddle=nsaddle+1
              if(i1.lt.i2) then
                 separatrices(i1,i2)=i
                 separatrices(i2,i1)=nsaddle
              else
                 separatrices(i1,i2)=nsaddle
                 separatrices(i2,i1)=i
              endif
           endif
           lsave=0
           xp=origin(1)+ix*step(1)
           yp=origin(2)+iy*step(2)
           zp=origin(3)+iz*step(3)
           xold=xp
           yold=yp
           zold=zp
           valp=dble(oldf)
c          call ascent(xp,yp,zp,origin,step,valp,iter,irw)
           call ascent(xp,yp,zp,valp,iter,irw)
           if(irw.eq.0) then
              xp=xold
              yp=yold
              zp=zold
              call attract_search(xp,yp,zp,root,valp,eps)
              do l=1,attract_label
                 dist=(xp-x(l))**2+(yp-y(l))**2+(zp-z(l))**2
                 dist=dsqrt(dist)
                 if(dist.le.0.25d0) then
                    if(l.le.ncore) then
                       rmax=(x(l)-xold)**2+(y(l)-yold)**2+
     .                      (z(l)-zold)**2
                       rmax=dsqrt(rmax)
                       if(rmax.gt.(1.5d0*core_radius(l))) then
                          irw=0
                       else
                          irw=l 
                       endif
                    else
                       irw=l
                    endif
                 endif
              enddo
              count=count+1
           endif
          lsave=irw
          attract_code(i)=lsave
        endif
      enddo
      call xtime(t)
      elapsed=t(1)+t(2)-spent
      spent=t(1)+t(2)
      write(*,'(/)')
      write(*,'(2x,''elapsed time in search step  '',f9.2,
     .             '' sec.'')') spent
      write(*,'(2x,''number of attractors:    '',i10)') attract_label+
     .         new_attr
      write(*,'(2x,''border points:           '',2i10)') nborder,count
      write(*,'(/)')
      if(new_attr.gt.0) then
         do i=1,new_attr
            x(attract_label+i)=xnew(i)
            y(attract_label+i)=ynew(i)
            z(attract_label+i)=znew(i)
         enddo
      endif
      attract_label=attract_label+new_attr
      if(attract_label.gt.3000) then
         write(*,*) 'number of basins exceeds 3000',
     .              'the wave function may be wrong'
         stop
      endif
      do j=1,low_limit
        i=indx(j)
        if(attract_code(i).eq.0) then
           ix=(i-1)/npyz
           iy=(i-1-ix*npyz)/np(3)
           iz=i-1-ix*npyz-iy*np(3)
           call neighbor(i,neighb)
           call neighcls(neighb,n_neighb,iborder)
           if(iborder.ne.0) then
              attract_code(i)=n_neighb(1)
           endif
        endif
      enddo
      deallocate(xnew,ynew,znew)
      return
      end
