program x
c----
c---- make a deshuffling index file based on a distance criterion
c---- input:
c---- 2 .gro files. First the atoms in the desired order, then the shuffled
c---- ones
c----

      implicit none

      integer i,j,nat,nr,iargc,last

      real x1(5000000),x2(50000000),dis,dx,dy,dz,boxx,boxy,boxz
      
      character*80 sdum

      logical found

      if (iargc().ne.2) then
         write(0,'(a)') 'Usage: deshuf sort.gro shuffled.gro > sort.ndx'
         write(0,'(a)')
         write(0,'(a,a)') 'writes an index file to deshuffle the atom', 
     &        ' order in shuffled.gro back to the one in sort.gro' 
         stop
      endif

      call getarg(1,sdum)
      open(1,file=sdum,status='old')
      call getarg(2,sdum)
      open(2,file=sdum,status='old')
      
      read(1,*)
      read(1,*) nat

      do i=1,nat
         read(1,'(20x,3f8.3)') x1(3*i-2),x1(3*i-1),x1(3*i)
      enddo
      read(1,*) boxx,boxy,boxz

      read(2,*)
      read(2,*) nat

      do i=1,nat
         read(2,'(20x,3f8.3)') x2(3*i-2),x2(3*i-1),x2(3*i)
      enddo

      write(6,'(a)') '[deshuf]'
      nr=0
      last=0
      do i=1,nat
         found=.false.
         do j=last+1,nat
            dx=x1(3*i-2)-x2(3*j-2)
            dy=x1(3*i-1)-x2(3*j-1)
            dz=x1(3*i)-x2(3*j)
            dis=dx**2+dy**2+dz**2
            if (dis.lt.0.005) then
               write(6,*) j
               nr=nr+1
               last=j
               found=.true.
               goto 10
            endif   
         enddo
         do j=1,max(1,last-1)
            dx=x1(3*i-2)-x2(3*j-2)
            dy=x1(3*i-1)-x2(3*j-1)
            dz=x1(3*i)-x2(3*j)
            dis=dx**2+dy**2+dz**2
            if (dis.lt.0.005) then
               write(6,*) j
               nr=nr+1
               last=j
               found=.true.
               goto 10
            endif   
         enddo
         if (.not.found) then
C-----
C----- check if atom is PBC-ed, rect box only
C-----
            do j=last+1,nat
               dx=x1(3*i-2)-x2(3*j-2)
               dy=x1(3*i-1)-x2(3*j-1)
               dz=x1(3*i)-x2(3*j)
               if (dx.gt.boxx*0.5) dx=dx-boxx
               if (dx.lt.boxx*-0.5) dx=dx+boxx
               if (dy.gt.boxy*0.5) dy=dy-boxy
               if (dy.lt.boxy*-0.5) dy=dy+boxy
               if (dz.gt.boxz*0.5) dz=dz-boxz
               if (dz.lt.boxz*-0.5) dz=dz+boxz
               dis=dx**2+dy**2+dz**2
               if (dis.lt.0.005) then
                  write(6,*) j
                  nr=nr+1
                  last=i
                  found=.true.
                  goto 10
               endif
            enddo
            do j=1,max(1,last-1)
               dx=x1(3*i-2)-x2(3*j-2)
               dy=x1(3*i-1)-x2(3*j-1)
               dz=x1(3*i)-x2(3*j)
               if (dx.gt.boxx*0.5) dx=dx-boxx
               if (dx.lt.boxx*-0.5) dx=dx+boxx
               if (dy.gt.boxy*0.5) dy=dy-boxy
               if (dy.lt.boxy*-0.5) dy=dy+boxy
               if (dz.gt.boxz*0.5) dz=dz-boxz
               if (dz.lt.boxz*-0.5) dz=dz+boxz
               dis=dx**2+dy**2+dz**2
               if (dis.lt.0.005) then
                  write(6,*) j
                  nr=nr+1
                  last=i
                  found=.true.
                  goto 10
               endif
            enddo
         endif
 10      continue
         if (.not.found) then
            last=0
            write(0,*) i,' not found!'
         endif
      enddo

      if (nr.ne.nat) then 
         write (0,*) 'WARNING: not all atoms could be assigned'
         write (0,*) 'only',nr,' out of',nat
      endif      

      end