src/tools/g_dist.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_g_dist_c = "$Id$";
  30
  31 #include <typedefs.h>
  32 #include "smalloc.h"
  33 #include "macros.h"
  34 #include "math.h"
  35 #include "xvgr.h"
  36 #include "copyrite.h"
  37 #include "statutil.h"
  38 #include "string2.h"
  39 #include "vec.h"
  40 #include "rdgroup.h"
  41 #include "pbc.h"
  42 #include "fatal.h"
  43 #include "futil.h"
  44 #include "gstat.h"
  45 #include "pbc.h"
  46
  47 int main(int argc,char *argv[])
  48 {
  49   static char *desc[] = {
  50     "g_dist can calculate the distance between the centers of mass of two",
  51     "groups of atoms as a function of time.[PAR]",
  52     "Or when [TT]-dist[tt] is set, print all the atoms in group 2 that are",
  53     "closer than a certain distance to the center of mass of group 1."
  54   };
  55
  56   t_topology *top=NULL;
  57   real t,cut2,dist2;
  58   rvec *x=NULL,*v=NULL,dx;
  59   matrix box;
  60   int status;
  61   int natoms;
  62
  63   int g,d,i,j,res,teller=0;
  64   atom_id aid;
  65
  66   int     ngrps;     /* the number of index groups */
  67   atom_id **index,max;   /* the index for the atom numbers */
  68   int     *isize;    /* the size of each group */
  69   char    **grpname; /* the name of each group */
  70   rvec    *com;
  71   real    *mass;
  72   FILE    *fp=NULL;
  73   bool    bCutoff;
  74
  75   static real cut=0;
  76
  77   static t_pargs pa[] = {
  78     { "-dist",      FALSE, etREAL, {&cut},
  79       "Print all atoms in group 2 closer than dist to the center of mass of group 1" },
  80   };
  81 #define NPA asize(pa)
  82
  83   t_filenm fnm[] = {
  84     { efTRX, "-f", NULL, ffREAD },
  85     { efTPX, NULL, NULL, ffREAD },
  86     { efNDX, NULL, NULL, ffOPTRD },
  87     { efXVG, NULL, "dist", ffOPTWR },
  88   };
  89 #define NFILE asize(fnm)
  90
  91
  92   CopyRight(stdout,argv[0]);
  93
  94   parse_common_args(&argc,argv,PCA_CAN_TIME,TRUE,
  95                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
  96
  97   bCutoff=opt2parg_bSet("-dist",NPA,pa);
  98   cut2=cut*cut;
  99
 100   top=read_top(ftp2fn(efTPX,NFILE,fnm));
 101
 102   /* read index files */
 103   ngrps = 2;
 104   snew(com,ngrps);
 105   snew(grpname,ngrps);
 106   snew(index,ngrps);
 107   snew(isize,ngrps);
 108   get_index(&top->atoms,ftp2fn(efNDX,NFILE,fnm),ngrps,isize,index,grpname);
 109
 110   /* calculate mass */
 111   max=0;
 112   snew(mass,ngrps);
 113   for(g=0;(g<ngrps);g++) {
 114     mass[g]=0;
 115     for(i=0;(i<isize[g]);i++) {
 116       if (index[g][i]>max)
 117         max=index[g][i];
 118       if (index[g][i] >= top->atoms.nr)
 119         fatal_error(0,"Atom number %d, item %d of group %d, is larger than number of atoms in the topolgy (%d)\n",index[g][i]+1,i+1,g+1,top->atoms.nr+1);
 120       mass[g]+=top->atoms.atom[index[g][i]].m;
 121     }
 122   }
 123
 124   natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 125
 126   if (max>=natoms)
 127     fatal_error(0,"Atom number %d in an index group is larger than number of atoms in the trajectory (%d)\n",(int)max+1,natoms);
 128
 129   if (!bCutoff)
 130     /* open output file */
 131     fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),
 132                   "Distance","Time (ps)","Distance (nm)");
 133   else
 134     ngrps=1;
 135
 136   do {
 137     /* initialisation for correct distance calculations */
 138     init_pbc(box,FALSE);
 139     /* make molecules whole again */
 140     rm_pbc(&top->idef,natoms,box,x,x);
 141
 142     /* calculate center of masses */
 143     for(g=0;(g<ngrps);g++) {
 144       for(d=0;(d<DIM);d++) {
 145         com[g][d]=0;
 146         for(i=0;(i<isize[g]);i++) {
 147           com[g][d] += x[index[g][i]][d] * top->atoms.atom[index[g][i]].m;
 148         }
 149         com[g][d] /= mass[g];
 150       }
 151     }
 152
 153     if (!bCutoff) {
 154       /* write to output */
 155       fprintf(fp,"%8.3f ",t);
 156       for(g=0;(g<ngrps/2);g++) {
 157         pbc_dx(com[2*g],com[2*g+1],dx);
 158         fprintf(fp,"%10.5f %10.5f %10.5f %10.5f",
 159                 norm(dx),dx[XX],dx[YY],dx[ZZ]);
 160       }
 161       fprintf(fp,"\n");
 162     } else {
 163       for(i=0;(i<isize[1]);i++) {
 164         j=index[1][i];
 165         pbc_dx(x[j],com[0],dx);
 166         dist2 = norm2(dx);
 167         if (dist2<cut2) {
 168           res=top->atoms.atom[j].resnr;
 169           fprintf(stdout,"\rt: %g  %d %s %d %s  %g (nm)\n",
 170                   t,res+1,*top->atoms.resname[res],
 171                   j+1,*top->atoms.atomname[j],sqrt(dist2));
 172         }
 173       }
 174     }
 175
 176     teller++;
 177   } while (read_next_x(status,&t,natoms,x,box));
 178
 179   if (!bCutoff)
 180     fclose(fp);
 181
 182   close_trj(status);
 183
 184   thanx(stdout);
 185   return 0;
 186 }