[gmx-users] shuffle/sort

[Bert de Groot]

Hi,

to return to the recent shuffle/sort discussion, I just found back a small
(and stupid) program to generate a deshuffle/desort index file based on
two .gro files: one with the atoms in the desired (normal sorting) and the
second one the shuffled one. The two files should contain the same coordinates,
but in a different order. You can e.g. use the input structure file for your
md for the first (usually confout.gro from a minimisation or posres run),
and the shuffled structure at time t=0 from your traj.??? for the second.

can be useful if you lost your deshuf.ndx or want to follow that special
solvent/lipid/whatever molecule...

Bert

____________________________________________________________________________
Dr. Bert de Groot

Max Planck Institute for Biophysical Chemistry
Theoretical molecular biophysics group
Am Fassberg 11 
37077 Goettingen, Germany

tel: +49-551-2011306, fax: +49-551-2011089

email: bgroot at gwdg.de
http://www.mpibpc.gwdg.de/abteilungen/071/bgroot
____________________________________________________________________________
-------------- next part --------------
      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

      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
      do i=1,nat
         found=.false.
         do j=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
               found=.true.
               goto 10
            endif
         enddo
         if (.not.found) then
C-----
C----- check if atom is PBC-ed, rect box only
C-----
            do j=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
                  found=.true.
                  goto 10
               endif
            enddo
         endif
 10      continue
         if (.not.found) write(0,*) i,' not found!'
      enddo

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

      end