      use grid_module
      use confs_module
      use orbitals_module
      use eval_module
      use basins_module
      implicit none
      integer i,niter,nshell(0:104),ier,j
      integer atom_list(0:104), list_at, level
      integer ns,k,l,m
      integer ns2,ix
      integer saddles,n_neighb
      integer vpoint,vcircle,n_fixed,iacc,munu
      integer nchange,nchange1
      real t(2),elapsed,spent
      integer, dimension (:,:), allocatable :: dummy,index_basin
      integer, dimension (:), allocatable :: val_label
      double precision, dimension (:), allocatable :: volume
      double precision xp,yp,zp,stepc,xx,yy,zz
      double precision rad,radius,elfold,ds(0:5)
      double precision root(3),dv
      double precision rho,drho(3),d2rho(6),elfd,gelf(3),d2elf(6)
      double precision stepold
      double precision dist,xv(26),yv(26),zv(26)
      double precision eps
      double precision dx,dy,dz,scale
      double precision radi,radj,radc,ah,rv(12),pi,rand
      character*4, dimension (:,:), allocatable :: basin_name
      character*4 N_At
      character*40 filebas, fileelf,filein,filegam,molecule
      character*40 name1
      character*3 func_name
      character*2 Bq
      character*1 ok
      logical func_type, do_attrib,rho_and_elf
      common/ftype/func_type
      common /files/ molecule,filein,filegam,fileelf,filebas,rho_and_elf
      common/accuracy/iacc
      external func,func2,mfunc
      data nshell/3*0,8*1,8*2,18*3,18*4,32*5,18*6/
      data atom_list/105*0/
      data N_At /'N_At'/
      data Bq /'Bq'/
      open(unit=99,file='temp.bas',status='unknown')
      open(unit=45,file='elf09_bas.gjf',status='unknown')
      write(*,*) ' function: elf/rho '
      read(*,'(a3)') func_name
      write(99,'(a3)') func_name
      eps=1.0d-4
      ah=rand
      pi=dacos(-1.0d0)
      if(func_name.eq.'elf') then
        func_type=.true.
      else
        func_type=.false.
      endif
      do_attrib=.true.
      saddles=0
      niter=0
      ds(0)=0.275d0
      ds(1)=0.1d0
      ds(2)=0.05d0
      ds(3)=0.025d0
      ds(4)=0.012d0
      ds(5)=0.006d0
      call xtime(t)
      call read_wfn
      allocate(indij(nmol*(nmol+1)/2))
      allocate(g(ngto),gx(ngto),gy(ngto),gz(ngto),extest(ngto))
      allocate(fxx(ngto),fxy(ngto),fxz(ngto),fyy(ngto),fyz(ngto))
      allocate(fzz(ngto),fxxx(ngto),fxxy(ngto),fxxz(ngto),fxyy(ngto))
      allocate(fxyz(ngto),fxzz(ngto),fyyy(ngto),fyyz(ngto),fyzz(ngto))
      allocate(fzzz(ngto))
      allocate(phi(nmol),dphi(nmol,3),d2phi(nmol,6),d3phi(nmol,10))
      allocate(rhorb(nmol*(nmol+1)/2),drhorb(nmol*(nmol+1)/2,3))
      allocate(d2rhorb(nmol*(nmol+1)/2,6),dotorb(nmol*(nmol+1)/2,3))
      allocate(ddotorb(nmol*(nmol+1)/2,9),d2dotorb(nmol*(nmol+1)/2,18))
      maxbasin=3000
      allocate(x(maxbasin),y(maxbasin),z(maxbasin))
      allocate(xs(maxbasin),ys(maxbasin),zs(maxbasin))
      allocate(val_elf(maxbasin),core_radius(maxbasin),volume(maxbasin))
      allocate(name_basin(maxbasin),basin_type(maxbasin,4))
      allocate(dummy(maxbasin,12),index_basin(maxbasin,5))
      allocate(aux(maxbasin,3))
      allocate(separatrices(maxbasin,maxbasin),val_label(maxbasin))
      allocate(basin_name(maxbasin,12))
      call read_elf
      open(unit=1,file=filebas,form='unformatted',status='unknown')
      do i=1,maxbasin
         do j=1,12
            basin_name(i,j)='    '
         enddo
      enddo
      write(99,'(a40)') filein
      munu=0
*      write(*,*) ' accuracy: 0 very high, 1 high, 2 medium '
*      read(*,*) iacc
        iacc = 1
      write(99,*) iacc
      if(func_type) then
*      write(*,*) ' find external core shell attractors?'
*      read(*,'(a)') ok
       ok = 'n'
      write(99,'(a)') ok
      write(*,'(/," Core basins")')
      write(*,'(/,"atom",t22,"attractor position",t45,"ELF",t55,
     .            "eigenvalues",t74,"radius",/)')
***
*     find core basins basin_type(i,1)=0
***
      attract_label=0
      do i=1,natoms
         ns=nshell(nat(i))
         ns2=ns
         if(atom_list(nat(i)).eq.0) then
         xp=xat(i)
         yp=yat(i)
         zp=zat(i)
         level=0
         call elf_eval(rho,drho,d2rho,elfd,gelf,d2elf,xp,yp,zp,level)
         if(ns.eq.0) then
            if(elfd.le.0.3d0) then
            attract_label=attract_label+1
            if(attract_label.gt.3000) then
               write(*,*) 'number of basins exceeds 3000, ',
     .                    'wave function probably wrong'
               stop
            endif
            basin_type(attract_label,1)=0
            basin_type(attract_label,2)=0
            aux(attract_label,1)=0.0d0
            aux(attract_label,2)=0.0d0
            aux(attract_label,3)=0.0d0
            x(attract_label)=xp
            y(attract_label)=yp
            z(attract_label)=zp
            xs(attract_label)=xp
            ys(attract_label)=yp
            zs(attract_label)=zp
            basin_name(attract_label,1)=atom_name(i)
            val_elf(attract_label)=elfd
            write(*,'(a4,12x,7f8.3)') atom_name(i),xp,yp,zp,elfd,root
            endif
         endif
         if(ns.ge.1) then
            radius=1.0d2
            n_neighb=0
            do k=1,natoms
               if(k.ne.i) then
                  dist=(xat(i)-xat(k))**2+(yat(i)-yat(k))**2+
     .                 (zat(i)-zat(k))**2
                  if(dist.le.radius) then
                     radius=dist
                     n_neighb=k
                  endif
               endif
            enddo
            if(n_neighb.eq.0) then
               radius=1.0d2
            endif
            scale=dsqrt(radius)
            root(1)=0.0d0
            root(2)=0.0d0
            root(3)=0.0d0
            if(elfd.ge.0.3d0) then
               call attract_search(xp,yp,zp,root,elfd,eps)
            else
               write(*,'(/,a4," pseudopotential type: 1 large core",
     .                   /,27x,"2 small core")') atom_name(i)
               read(*,*) ns2
               write(99,'(i4)') ns2
            endif
            attract_label=attract_label+1
            if(attract_label.gt.3000) then
               write(*,*) 'number of basins exceeds 3000, ',
     .                    'wave function probably wrong'
               stop
            endif
            atom_list(nat(i))=attract_label
            x(attract_label)=xp
            y(attract_label)=yp
            z(attract_label)=zp
            xs(attract_label)=xp
            ys(attract_label)=yp
            zs(attract_label)=zp
            basin_type(attract_label,1)=0
            basin_type(attract_label,2)=0
            basin_name(attract_label,1)=atom_name(i)
            write(*,'(a4,12x,7f8.3)') atom_name(i),xp,yp,zp,elfd,root
            aux(attract_label,1)=root(1)
            aux(attract_label,2)=root(2)
            aux(attract_label,3)=root(3)
            core_radius(attract_label)=0.2d0
            val_elf(attract_label)=elfd
            rad=0.0d0
            if(ns2.ne.0) then
               do j=1,ns2
                  elfold=elfd
                  stepc=ds(ns2-j)/scale
                  stepold=stepc
                  if(n_neighb.ne.0) then
                     dx=stepc*(xat(n_neighb)-xat(i))
                     dy=stepc*(yat(n_neighb)-yat(i))
                     dz=stepc*(zat(n_neighb)-zat(i))
                     xp=xat(i)+rad*(xat(n_neighb)-xat(i))/scale+dx
                     yp=yat(i)+rad*(yat(n_neighb)-yat(i))/scale+dy
                     zp=zat(i)+rad*(zat(n_neighb)-zat(i))/scale+dz
                  else
                     xp=xat(i)+5.0d0*stepc
                     yp=yat(i)
                     zp=zat(i)
                     dx=stepold
                  endif
                  call sep_core(xp,yp,zp,dx,dy,dz,elfold)
                  radius=(x(attract_label)-xp)**2+
     .                   (y(attract_label)-yp)**2+
     .                   (z(attract_label)-zp)**2
                  radius=dsqrt(radius)
                  write(*,'(16x,4f8.3,24x,f6.3)') xp,yp,zp,elfold,radius
                  if(ok.ne.'y') then
                     core_radius(attract_label)=radius
                  endif
                  if(j.ne.ns2) then
                     core_radius(attract_label)=radius
                  endif
                  if(j.lt.ns2) then
                     rad=0.0d0
                     stepc=ds(ns2-j)/scale
                     do while((rad-core_radius(attract_label)).le.0.01)
                        stepold=stepc
                        if(n_neighb.ne.0) then
                           dx=stepc*(xat(n_neighb)-xat(i))
                           dy=stepc*(yat(n_neighb)-yat(i))
                           dz=stepc*(zat(n_neighb)-zat(i))
                           xx=xp+dx
                           yy=yp+dy
                           zz=zp+dz
                        else
                           xx=xx+5.d0*stepold
                        endif
                        call at_shell(xx,yy,zz,dx,dy,dz,elfold)
                        eps=2.0d0
                        call attract_search(xx,yy,zz,root,elfd,eps)
                        eps=1.0d-4
                        radius=(x(attract_label)-xx)**2+
     .                         (y(attract_label)-yy)**2+
     .                         (z(attract_label)-zz)**2
                        rad=dsqrt(radius)
                        if(rad.lt.core_radius(attract_label)) then
                          stepc=1.404d0*stepold
                        else
                          write(*,'(16x,7f8.3,f6.3)') xx,yy,zz,elfd,
     .                                                root,rad
                          xp=xx
                          yp=yy
                          zp=zz
                        endif
                     enddo
                  endif
               enddo
***
*       find external shell core basins of heavy atoms (Z>21)
*       basin_type(i,1)=-1
***
               if((nat(i).ge.19).and.(ok.eq.'y')) then
                  call shifts(rad,xv,yv,zv)
                  vpoint=-1
                  vcircle=-3
                  do l=1,26
                     xp=xat(i)+xv(l)+0.05d0*(1.0d0-rand)
                     yp=yat(i)+yv(l)+0.05d0*(1.0d0-rand)
                     zp=zat(i)+zv(l)+0.05d0*(1.0d0-rand)
                     call attractors(xp,yp,zp,i,j,vpoint,vcircle)
                     basin_name(attract_label,1)=atom_name(i)
                  enddo
               else
                  core_radius(attract_label)=radius
               endif
            else
            endif
         endif
         else
            list_at=atom_list(nat(i))
            attract_label=attract_label+1
            if(attract_label.gt.3000) then
               write(*,*) 'number of basins exceeds 3000, ',
     .                    'wave function probably wrong'
               stop
            endif
            xp=xat(i)
            yp=yat(i)
            zp=zat(i)
            x(attract_label)=xp
            y(attract_label)=yp
            z(attract_label)=zp
            xs(attract_label)=xp
            ys(attract_label)=yp
            zs(attract_label)=zp
            basin_type(attract_label,1)=0
            basin_type(attract_label,2)=0
            basin_name(attract_label,1)=atom_name(i)
            aux(attract_label,1)=aux(list_at,1)
            aux(attract_label,2)=aux(list_at,2)
            aux(attract_label,3)=aux(list_at,3)
            core_radius(attract_label)=core_radius(list_at)
            val_elf(attract_label)=val_elf(list_at)
            root(1)=aux(attract_label,1)
            root(2)=aux(attract_label,2)
            root(3)=aux(attract_label,3)
            elfd=val_elf(attract_label)
            write(*,'(a4,12x,7f8.3,f6.3)') atom_name(i),xp,yp,zp,elfd,
     .                                  root,core_radius(attract_label)
         endif
      enddo
      ncore=attract_label
***
*     find protonated basins basin_type(i,1)=1
***
      do i=1,natoms
         ns=nshell(nat(i))
         if(ns.eq.0) then
            xp=xat(i)
            yp=yat(i)
            zp=zat(i)
***
*           check for pseudopotential
***
            level=0
            call elf_eval(rho,drho,d2rho,elfd,gelf,d2elf,xp,yp,zp,level)
            root(1)=0.0d0
            root(2)=0.0d0
            root(3)=0.0d0
            if(elfd.ge.0.3d0) then
               attract_label=attract_label+1
               if(attract_label.gt.3000) then
                  write(*,*) 'number of basins exceeds 3000, ',
     .                       'wave function probably wrong'
                  stop
               endif
               call attract_search(xp,yp,zp,root,elfd,eps)
               radius=dsqrt((xp-xat(i))**2+(yp-yat(i))**2+
     .                      (zp-zat(i))**2)
               basin_type(attract_label,1)=1
               basin_type(attract_label,2)=1
               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
               xs(attract_label)=xp
               ys(attract_label)=yp
               zs(attract_label)=zp
               basin_name(attract_label,1)=atom_name(i)
               core_radius(attract_label)=radius
               val_elf(attract_label)=elfd
            endif
         endif
      enddo
      n_fixed=attract_label
***
*     find valence basins basin_type(i,1)=2 if not protonated
*                         basin_type(i,1)=3 if circular attractor
*                         basin_type(i,1)=4 if circular attractor
***
      vpoint=2
      vcircle=3
      write(*,'(/)')
*      write(*,'(" search mode: 0 automatic, >0 number of attractors",
*     .          " in input, <0 fast")')
*      read(*,*) ier
       ier = 0
      if(ier.gt.0) then
         do_attrib=.false.
      endif
      if(ier.eq.0) then
         do i=1,ncore
            if(basin_type(i,1).eq.0) then
               rad=1.85d0*core_radius(i)
               call shifts(rad,xv,yv,zv)
               do j=1,26
                  xp=x(i)+xv(j)+0.05d0*(1.0d0-rand)
                  yp=y(i)+yv(j)+0.05d0*(1.0d0-rand)
                  zp=z(i)+zv(j)+0.05d0*(1.0d0-rand)
                  call attractors(xp,yp,zp,i,j,vpoint,vcircle)
               enddo
            endif
         enddo
         do i=1,ncore-1
            if(basin_type(i,1).eq.0) then
               do j=i+1,ncore
                  if(basin_type(j,1).eq.0) then
                     rad=4.0d0*(core_radius(i)+core_radius(j))
                     radc=(x(i)-x(j))**2+(y(i)-y(j))**2+(z(i)-z(j))**2
                     if(dsqrt(radc).lt.rad) then
                        xp=x(i)+0.3d0*(x(j)-x(i))
                        yp=y(i)+0.3d0*(y(j)-y(i))
                        zp=z(i)+0.3d0*(z(j)-z(i))
                        call attractors(xp,yp,zp,i,j,vpoint,vcircle)
                        xp=x(i)+0.7d0*(x(j)-x(i))
                        yp=y(i)+0.7d0*(y(j)-y(i))
                        zp=z(i)+0.7d0*(z(j)-z(i))
                        call attractors(xp,yp,zp,i,j,vpoint,vcircle)
                     endif
                  endif
               enddo
            endif
         enddo
      endif
      do j=1,attract_label
         xp=x(j)
         yp=y(j)
         zp=z(j)
         level=1
         call elf_eval(rho,drho,d2rho,elfold,gelf,d2elf,
     .                          xp,yp,zp,level)
      enddo
      if(ier.lt.0) then
         call valence_attr
      else
         do j=1,ier
            attract_label=attract_label+1
            if(attract_label.gt.3000) then
               write(*,*) 'number of basins exceeds 3000, ',
     .                    'wave function probably wrong'
               stop
            endif
            write(*,*) ' input coordinates of attractor ',j
            read(*,*) x(attract_label),y(attract_label),z(attract_label)
            write(*,*) ' input type (0 nuclear, 1 point, 3 circle or ',
     .                 '4 sphere and atom label defining axis direction'
            read(*,*) basin_type(attract_label,1),
     .                basin_type(attract_label,3)
            write(*,*) ' synaptic order '
            read(*,*) basin_type(attract_label,2)
            write(*,*) ' input radius'
            read(*,*) core_radius(attract_label)
            write(*,*) ' input basin name 12 characters'
            read(*,'(3a4)') (basin_name(attract_label,k),k=1,3)
            do_attrib=.false.
            xp=x(attract_label)
            yp=y(attract_label)
            zp=z(attract_label)
            if(basin_type(attract_label,1).eq.4) then
               xp=xp+core_radius(attract_label)
            endif
            level=0
            call elf_eval(rho,drho,d2rho,elfold,gelf,d2elf,
     .                          xp,yp,zp,level)
            val_elf(attract_label)=elfold
         enddo
      endif
      call xtime(t)
      elapsed=t(1)+t(2)-spent
      spent=t(1)+t(2)
      write(*,'(/,2x,''elapsed time in attractor search'',f9.2,
     .               '' sec.'')') elapsed
      write(*,'(/)')
      write(*,*) attract_label,' attractors found'
*      write(*,*) 'assign grid points? '
*      read(*,'(a)') ok
        ok ='y'
      if(ok.eq.'y') then
         call assign_bas
      endif
      call xtime(t)
      spent=t(1)+t(2)
      eps=1.0d-3
      do i=ncore+1,attract_label
         k=0
         if(basin_type(i,1).eq.1) k=1
         if(do_attrib.or.basin_type(i,1).eq.1) then
         if(natoms.eq.1) then
            basin_name(i,1)=basin_name(1,1)
         else
            do j=1,ncore
               if(basin_type(j,1).eq.0) then
                  rad=(xs(i)-x(j))**2+(ys(i)-y(j))**2+(zs(i)-z(j))**2
                  rad=1.2d0*core_radius(j)/dsqrt(rad)
                  xp=x(j)+rad*(xs(i)-x(j))+0.005d0*(1.0d0-rand)
                  yp=y(j)+rad*(ys(i)-y(j))+0.005d0*(1.0d0-rand)
                  zp=z(j)+rad*(zs(i)-z(j))+0.005d0*(1.0d0-rand)
                  xx=xp
                  yy=yp
                  zz=zp
                  level=0
                  call elf_eval(rho,drho,d2rho,elfold,gelf,d2elf,
     .                          xp,yp,zp,level)
                  eps=1.0d-4
                  call attract_search(xp,yp,zp,root,elfd,eps)
                  if(elfd.gt.elfold) then
                     dist=(xp-x(i))**2+(yp-y(i))**2+(zp-z(i))**2
                     dist=dsqrt(dist)
                     stepc=0.1d0
                     if(dist.lt.stepc) then
                        k=k+1
                        basin_type(i,2)=basin_type(i,2)+1
                        basin_name(i,k)=basin_name(j,1)
                        dummy(i,k)=j
                     endif
                     if(basin_type(i,1).eq.3) then
                        if(dabs(dist-core_radius(i)).le.0.2d0) then
                           k=k+1
                           basin_type(i,2)=basin_type(i,2)+1
                           basin_name(i,k)=basin_name(j,1)
                           dummy(i,k)=j
                        endif
                     endif
                  endif
               endif
            enddo
         endif
      endif
      enddo
      call xtime(t)
      elapsed=t(1)+t(2)-spent
      spent=t(1)+t(2)
      write(*,'(/,2x,''elapsed time in attractor assignment'',f9.2,
     .               '' sec.'',/)') elapsed
      do i=1,ncore
         name1='C('//basin_name(i,1)//')'
         call short_name(name1,name_basin(i))
      enddo
      do i=ncore+1,attract_label
         j=basin_type(i,1)
         l=basin_type(i,2)
         if(j.le.2) then
            if(l.lt.1) then
               basin_name(i,1)='Asyn'
               name1='V(Asyn)'
            endif
            if(l.eq.1) then
               name1='V('//basin_name(i,1)//')'
            endif
            if(l.eq.2) then
               name1='V('//basin_name(i,1)//','//basin_name(i,2)
     .                      //')'
            endif
            if(l.eq.3) then
               name1='V('//basin_name(i,1)//','//basin_name(i,2)
     .                   //','//basin_name(i,3)//')'
            endif
            call short_name(name1,name_basin(i))
         else if(j.ge.3) then
            k=basin_type(i,1)
            if(l.le.1) then
               name1='V('//basin_name(i,1)//')'
            endif
            if(l.eq.2) then
               name1='V('//basin_name(i,1)//','//basin_name(i,2)
     .                      //')'
            endif
            if(l.eq.3) then
               name1='V('//basin_name(i,1)//','//basin_name(i,2)
     .                   //','//basin_name(i,3)//')'
            endif
            call short_name(name1,name_basin(i))
         endif
      enddo
      write(*,'(/)')
      npts=(np(1)-1)*(np(2)-1)*(np(3)-1)
      dv=0.0d0
      if(npts.gt.0) then
         dv=(xmax-xmin)*(ymax-ymin)*(zmax-zmin)/npts
      endif
      npts=np(1)*np(2)*np(3)
      do i=1,attract_label
         index_basin(i,4)=0
         index_basin(i,5)=0
      enddo
      nchange=0
      nchange1=0
      do i=1,npts
         j=attract_code(i)
         if(j.gt.0) then
            index_basin(j,5)=index_basin(j,5)+1
            if(elf(i).ge.0.02) then
               index_basin(j,4)=index_basin(j,4)+1
            endif
         else
            nchange=nchange+1
         endif
      enddo
      do i=1,attract_label
         j=basin_type(i,1)
         if(j.le.0) then
            val_label(i)=5
         else if (j.eq.1) then
            val_label(i)=2
         else if (basin_type(i,2).eq.0) then
            val_label(i)=4
         else if (basin_type(i,2).eq.1) then
            val_label(i)=1
         else
            val_label(i)=3
         endif
      enddo
    		write(45,*) '# '
      		write(45,*) '   '
      		write(45,*) 'ELF basins. Edit atom types for core/valence'
      		write(45,*) '   '
      		write(45,*) '0  1'
      do i=1,attract_label
         volume(i)=dv*index_basin(i,4)
         write(*,'(i4,4f8.3,2x,a16,2i7,f8.2,i3)') i,x(i),y(i),z(i),
     .            val_elf(i),name_basin(i),index_basin(i,5),
     .            index_basin(i,4),volume(i),val_label(i)
      write(45,'(a2,1x,a16,3f15.4)') Bq, 
     .            name_basin(i),x(i),y(i),z(i)
      enddo
      write(45,*) '   '
      write (*,*) ' File  elf09_bas.gjf written for Gaussian display'
      write(99,'(i4)') attract_label-n_fixed
      do i=n_fixed+1,attract_label
         write(99,'(3f20.14)') x(i),y(i),z(i)
         write(99,'(2i4)') basin_type(i,1),basin_type(i,3)
         write(99,'(i4)') basin_type(i,2)
         write(99,'(f20.14)') core_radius(i)
         write(99,'(3a4)') (basin_name(i,k),k=1,3)
      enddo
      write(*,'(/)')
      write(*,'(2x,''unassigned grid points'',2i10)') nchange,nchange1 
      if(ok.eq.'y') then
      write(1) np(1),np(2),np(3)
      write(1) xmin,xmax,ymin,ymax,zmin,zmax
      write(1) natoms,attract_label
      do i=1,natoms
         write(1) nat(i),atom_name(i),xat(i),yat(i),zat(i)
      enddo
      do i=1,attract_label
         write(1) val_label(i),name_basin(i),volume(i),x(i),y(i),z(i)
      enddo
      write(1)(attract_code(i),i=1,npts)
      close(unit=1)
      endif
***
      write(*,'(/," distances from nuclei (A)",/)')
      do i=ncore+1,attract_label
         j=basin_type(i,2)
         m=1
         if(basin_type(i,1).eq.1) then
            m=2
         endif
         do k=m,j
            l=dummy(i,k)
            dist=dsqrt((x(i)-x(l))**2+(y(i)-y(l))**2+(z(i)-z(l))**2)
            rv(k)=0.52917715d0*dist
         enddo
         write(*,'(a16,6(4x,a4))') name_basin(i),
     .                              (basin_name(i,k),k=m,j)
         write(*,'(16x,6f8.3)') (rv(k),k=m,j)
      enddo
      if(attract_label.gt.(ncore+1)) then
         write(*,'(/," angles around core attractors ",/)')
         do i=1,ncore
            if(basin_type(i,1).eq.0) then
               do j=ncore+1,attract_label-1
                  ix=0
                  do l=1,basin_type(j,2)
                     if(dummy(j,l).eq.i) ix=j
                  enddo
                  if(ix.ne.0) then
                     radi=dsqrt((x(ix)-x(i))**2+(y(ix)-y(i))**2+
     .                         (z(ix)-z(i))**2)
                     do l=j+1,attract_label
                        do m=1,basin_type(l,2)
                           if(dummy(l,m).eq.i) then
                              radj=dsqrt((x(l)-x(i))**2+(y(l)-y(i))**2+
     .                            (z(l)-z(i))**2)
                              radc=(x(l)-x(ix))**2+(y(l)-y(ix))**2
     .                             +(z(l)-z(ix))**2
                              ah=(radi*radi+radj*radj-radc)/(radi*radj)
                              rad=acos(0.49999999999999d0*ah)*180.0d0/pi
                              write(*,'(3a16,f12.3)') name_basin(ix),
     .                          name_basin(i),name_basin(l),rad
                           endif
                        enddo
                     enddo
                  endif
               enddo
            endif
         enddo
      endif
      write(*,'(/," angles around polysynaptic attractors ",/)')
      do i=ncore+1,attract_label
         if(basin_type(i,1).ge.2) then
            niter=basin_type(i,2)
            do j=1,niter-1
               l=dummy(i,j)
               radi=dsqrt((x(l)-x(i))**2+(y(l)-y(i))**2+(z(l)-z(i))**2)
               do k=j+1,niter
                  m=dummy(i,k)
                  radj=dsqrt((x(m)-x(i))**2+(y(m)-y(i))**2+
     .                        (z(m)-z(i))**2)
                  radc=(x(m)-x(l))**2+(y(m)-y(l))**2+(z(m)-z(l))**2
                  ah=(radi*radi+radj*radj-radc)/(radi*radj)
                  rad=acos(0.49999999999999d0*ah)*180.0d0/pi
                  write(*,'(a4,a16,a4,f12.3)') basin_name(i,j),
     .                   name_basin(i),basin_name(i,k),rad
               enddo
            enddo
         endif
      enddo
      else
         attract_label=natoms
         do i=1,natoms
            x(i)=xat(i)
            y(i)=yat(i)
            z(i)=zat(i)
            xp=x(i)
            yp=y(i)
            zp=z(i)
            core_radius(i)=0.75d0
            level=0
            call elf_eval(rho,drho,d2rho,elfd,gelf,d2elf,xp,yp,zp,level)
            dist=drho(1)**2+drho(2)**2+drho(3)**2
            dist=dist/rho
            if(elfd.le.0.1d0) then
               write(*,'(/," pseudopotential: input core radius for",
     .                   /,2x,a4)') atom_name(i)
               core_radius(i)=1.4d0
            else if((dist.gt.1.0d-06).and.(rho.lt.1.0d0)) then
               xp=x(i)
               yp=y(i)
               zp=z(i)
               call attract_search(xp,yp,zp,root,elfd,eps)
               x(i)=xp
               y(i)=yp
               z(i)=zp
            endif
            basin_name(i,1)=atom_name(i)
            basin_type(i,1)=0
            name_basin(i)='Atom('//basin_name(i,1)//')'
            write(*,'(a16,7f8.3,2i3)') name_basin(i),x(i),y(i),z(i)
         enddo
      write(*,'(/)')
      write(*,*) 'assign grid points? '
      read(*,'(a)') ok
      if(ok.eq.'y') then
         call assign_bas
      endif
      if(attract_label.gt.natoms) then
         do i=natoms+1,attract_label
            name_basin(i)='Atom('//N_At//')'
         enddo
      endif
      call xtime(t)
      spent=t(1)+t(2)
      write(*,'(/)')
      npts=(np(1)-1)*(np(2)-1)*(np(3)-1)
      dv=(xmax-xmin)*(ymax-ymin)*(zmax-zmin)/npts
      npts=np(1)*np(2)*np(3)
      do i=1,attract_label
         index_basin(i,4)=0
         index_basin(i,5)=0
      enddo
      nchange=0
      nchange1=0
      do i=1,npts
         j=attract_code(i)
         if(j.gt.0) then
            index_basin(j,5)=index_basin(j,5)+1
            if(elf(i).ge.1.0d-10) then
               index_basin(j,4)=index_basin(j,4)+1
            endif
         else
            nchange=nchange+1
         endif
      enddo
      do i=1,attract_label
         j=basin_type(i,1)
         if(j.le.0) then
            val_label(i)=5
         else if (j.eq.1) then
            val_label(i)=2
         else if (basin_type(i,2).eq.1) then
            val_label(i)=1
         else
            val_label(i)=3
         endif
      enddo
      do i=1,attract_label
         volume(i)=dv*index_basin(i,4)
         write(*,'(i4,4f8.3,2x,a16,2i7,f8.2,i3)') i,x(i),y(i),z(i),
     .            val_elf(i),name_basin(i),index_basin(i,5),
     .            index_basin(i,4),volume(i),val_label(i)
      enddo
      write(*,'(/)')
      write(*,'(2x,''unassigned grid points'',2i10)') nchange,nchange1 
      if(ok.eq.'y') then
*     open(unit=1,file=filebas,form='unformatted',status='unknown')
      write(1) np(1),np(2),np(3)
      write(1) xmin,xmax,ymin,ymax,zmin,zmax
      write(1) natoms,attract_label
      do i=1,natoms
         write(1) nat(i),atom_name(i),xat(i),yat(i),zat(i)
      enddo
      do i=1,attract_label
         write(1) val_label(i),name_basin(i),volume(i),x(i),y(i),z(i)
      enddo
      write(1)(attract_code(i),i=1,npts)
      j=0
      do i=1,npts
         if(attract_code(i).gt.0) then
            j=j+1
         endif
      enddo
      close(unit=1)
      endif
      endif
      write(99,'("y")')
      close(unit=99)
      j=0
      do i=1,npts
         if(attract_code(i).gt.0) then
            j=j+1
         endif
      enddo
      end
