implicit real*8(a-h,o-z)
      parameter (nat=1000)
      parameter (ntypd = 4)
      parameter (nshelld=30)
      parameter (natmd=30)
      dimension num(5)
      dimension a(3,3),natom(3),neigh(ntypd,natmd,ntypd,natmd)
      dimension r(ntypd,natmd,ntypd,natmd,nat),rr(nshelld)
      dimension no(ntypd,natmd,ntypd,natmd,nat),mo(nshelld)
      dimension itype(ntypd,natmd,ntypd,natmd,nat)
      dimension t(ntypd,natmd,3)
      open(8,file='CONTCAR')
      open(22,file='neig-con.dat')

      write(6,*) 'enter the type of atom'
      read(5,*) ntype
      if (ntype .gt. ntypd) stop ' ntype too large '
      read(8,600) temp
 600  format(20a4)
      read (8,*) aa
 
c      *** read lattice constant from POSCAR**
 
      do i=1,3
      read (8,*) (a(i,j),j=1,3)
 500  format (3f12.8)
      enddo
      do i=1,3
       do j=1,3
       a(i,j)=aa*a(i,j)
       enddo
      enddo
    
      read(8,*) (natom(i),i=1,ntype)
      do i = 1,ntype
      if (natom(i) .gt. natmd) stop ' natom too large '
      enddo
      write(6,*) ' selective dynamic ? (1/yes 0/no) '

      read(5,*) isel
      if (isel .eq. 1) read(8,600) temp
      read(8,600) temp
      do i=1,ntype
       do j=1,natom(i)
       read(8,*) t(i,j,1),t(i,j,2),t(i,j,3)
       enddo
      enddo

c     *** read lattice vector from POSCAR***

      do 200 i = 1,ntype
       do 200 j = 1,natom(i)
       x = t(i,j,1)*a(1,1)+t(i,j,2)*a(2,1)+t(i,j,3)*a(3,1)
       y = t(i,j,1)*a(1,2)+t(i,j,2)*a(2,2)+t(i,j,3)*a(3,2)
       z = t(i,j,1)*a(1,3)+t(i,j,2)*a(2,3)+t(i,j,3)*a(3,3)
      
        do 100 ip = 1,ntype
         do 100 jp = 1,natom(ip)
         nn = 0
          do 50 n1=-5,5
           do 50 n2=-5,5
            do 50 n3=-5,5
            dtij1 = t(i,j,1) - t(ip,jp,1)
            dtij2 = t(i,j,2) - t(ip,jp,2)
            dtij3 = t(i,j,3) - t(ip,jp,3)
            xp = (n1+dtij1)*a(1,1) + (n2+dtij2)*a(2,1) +
     &                                  (n3+dtij3)*a(3,1)
            yp = (n1+dtij1)*a(1,2) + (n2+dtij2)*a(2,2) +
     &                                  (n3+dtij3)*a(3,2)
            zp = (n1+dtij1)*a(1,3) + (n2+dtij2)*a(2,3) +
     &                                  (n3+dtij3)*a(3,3)
            dr = sqrt(xp**2+yp**2+zp**2)
1234	 format(4i3,f10.3,3i3,3f10.3)
             if ((dr .gt. 12.0) .or. (dr .eq. 0.0)) go to 50
             nn = nn+ 1
             if (nn .gt. nat) stop ' nn too large '
             r(i,j,ip,jp,nn) = dr
             itype(i,j,ip,jp,nn) = ip
  50  continue

      neigh(i,j,ip,jp) = nn
 100  continue
 200  continue

      do 1000 i = 1,ntype
       do 1000 j = 1,natom(i)
        do 1000 ip = 1,ntype
         do 1000 jp = 1,natom(ip)
           do 60 m = 1,neigh(i,j,ip,jp)-1
            do 60 n = m+1,neigh(i,j,ip,jp)
             if (r(i,j,ip,jp,n) .lt. r(i,j,ip,jp,m)) then
             dump = r(i,j,ip,jp,n)
             r(i,j,ip,jp,n) = r(i,j,ip,jp,m)
             r(i,j,ip,jp,m) = dump
             idump = itype(i,j,ip,jp,n)
             itype(i,j,ip,jp,n) = itype(i,j,ip,jp,m)
             itype(i,j,ip,jp,m) = idump
             endif
  60  continue        
 1000 continue

      do 2000 i = 1,ntype 
       do 2000 j = 1,natom(i) 
        do 2000 ip = 1,ntype
         do 2000 jp = 1,natom(ip)
      
         nshell = 1
          do n= 2,neigh(i,j,ip,jp)
           if (abs(r(i,j,ip,jp,n)-r(i,j,ip,jp,n-1)) .gt. 0.0001) then
           nshell = nshell + 1
            if(nshell .gt. nshelld) stop ' nshell too large '
            rr(nshell) = r(i,j,ip,jp,n)
            no(i,j,ip,jp,n) = n
            mo(nshell) = no(i,j,ip,jp,n)
           endif
          enddo
          rr(1) = r(i,j,ip,jp,(mo(2)-1))
          mo(1) = mo(2)-1 
          mo(nshell+1) = neigh(i,j,ip,jp)
         
      num(1) = mo(1)
      do m = 2,nshell-1 
      num(m) = mo(m+1)-mo(m)
      num(nshell) = mo(nshell+1)-mo(nshell)+1
 111  format(4i4,5(f8.3,i4))
      enddo
      write(22,111) i,j,ip,jp,(rr(ii),num(ii),ii=1,5)
2000  continue

      stop
      end