src/tools/gmx_dist.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <typedefs.h>
  39
  40 #include "smalloc.h"
  41 #include "macros.h"
  42 #include "math.h"
  43 #include "xvgr.h"
  44 #include "copyrite.h"
  45 #include "statutil.h"
  46 #include "string2.h"
  47 #include "vec.h"
  48 #include "index.h"
  49 #include "pbc.h"
  50 #include "gmx_fatal.h"
  51 #include "futil.h"
  52 #include "gstat.h"
  53 #include "gmx_ana.h"
  54
  55
  56 static void add_contact_time(int **ccount,int *ccount_nalloc,int t)
  57 {
  58   int i;
  59
  60   if (t+2 >= *ccount_nalloc) {
  61     srenew(*ccount,t+2);
  62     for(i=*ccount_nalloc; i<t+2; i++)
  63       (*ccount)[i] = 0;
  64     *ccount_nalloc = t+2;
  65   }
  66   (*ccount)[t]++;
  67 }
  68
  69 int gmx_dist(int argc,char *argv[])
  70 {
  71   const char *desc[] = {
  72     "[TT]g_dist[tt] can calculate the distance between the centers of mass of two",
  73     "groups of atoms as a function of time. The total distance and its",
  74     "x, y and z components are plotted.[PAR]",
  75     "Or when [TT]-dist[tt] is set, print all the atoms in group 2 that are",
  76     "closer than a certain distance to the center of mass of group 1.[PAR]",
  77     "With options [TT]-lt[tt] and [TT]-dist[tt] the number of contacts",
  78     "of all atoms in group 2 that are closer than a certain distance",
  79     "to the center of mass of group 1 are plotted as a function of the time",
  80     "that the contact was continously present.[PAR]",
  81     "Other programs that calculate distances are [TT]g_mindist[tt]",
  82     "and [TT]g_bond[tt]."
  83   };
  84
  85   t_topology *top=NULL;
  86   int  ePBC;
  87   real t,t0,cut2,dist2;
  88   rvec *x=NULL,*v=NULL,dx;
  89   matrix box;
  90   t_trxstatus *status;
  91   int natoms;
  92
  93   int g,d,i,j,res,teller=0;
  94   atom_id aid;
  95
  96   int     ngrps;     /* the number of index groups */
  97   atom_id **index,max;   /* the index for the atom numbers */
  98   int     *isize;    /* the size of each group */
  99   char    **grpname; /* the name of each group */
 100   rvec    *com;
 101   real    *mass;
 102   FILE    *fp=NULL,*fplt=NULL;
 103   gmx_bool    bCutoff,bPrintDist,bLifeTime;
 104   t_pbc   *pbc;
 105   int     *contact_time=NULL,*ccount=NULL,ccount_nalloc=0,sum;
 106   char    buf[STRLEN];
 107   output_env_t oenv;
 108   gmx_rmpbc_t  gpbc=NULL;
 109
 110   const char *leg[4] = { "|d|","d\\sx\\N","d\\sy\\N","d\\sz\\N" };
 111
 112   static real cut=0;
 113
 114   static t_pargs pa[] = {
 115     { "-dist",      FALSE, etREAL, {&cut},
 116       "Print all atoms in group 2 closer than dist to the center of mass of group 1" }
 117   };
 118 #define NPA asize(pa)
 119
 120   t_filenm fnm[] = {
 121     { efTRX, "-f", NULL, ffREAD },
 122     { efTPX, NULL, NULL, ffREAD },
 123     { efNDX, NULL, NULL, ffOPTRD },
 124     { efXVG, NULL, "dist", ffOPTWR },
 125     { efXVG, "-lt", "lifetime", ffOPTWR },
 126   };
 127 #define NFILE asize(fnm)
 128
 129
 130   CopyRight(stderr,argv[0]);
 131
 132   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
 133                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 134
 135   bCutoff = opt2parg_bSet("-dist",NPA,pa);
 136   cut2 = cut*cut;
 137   bLifeTime = opt2bSet("-lt",NFILE,fnm);
 138   bPrintDist = (bCutoff && !bLifeTime);
 139
 140   top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
 141
 142   /* read index files */
 143   ngrps = 2;
 144   snew(com,ngrps);
 145   snew(grpname,ngrps);
 146   snew(index,ngrps);
 147   snew(isize,ngrps);
 148   get_index(&top->atoms,ftp2fn(efNDX,NFILE,fnm),ngrps,isize,index,grpname);
 149
 150   /* calculate mass */
 151   max=0;
 152   snew(mass,ngrps);
 153   for(g=0;(g<ngrps);g++) {
 154     mass[g]=0;
 155     for(i=0;(i<isize[g]);i++) {
 156       if (index[g][i]>max)
 157         max=index[g][i];
 158       if (index[g][i] >= top->atoms.nr)
 159         gmx_fatal(FARGS,"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);
 160       mass[g]+=top->atoms.atom[index[g][i]].m;
 161     }
 162   }
 163
 164   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 165   t0 = t;
 166
 167   if (max>=natoms)
 168     gmx_fatal(FARGS,"Atom number %d in an index group is larger than number of atoms in the trajectory (%d)\n",(int)max+1,natoms);
 169
 170   if (!bCutoff) {
 171     /* open output file */
 172     fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),
 173                   "Distance","Time (ps)","Distance (nm)",oenv);
 174     xvgr_legend(fp,4,leg,oenv);
 175   } else {
 176     ngrps = 1;
 177     if (bLifeTime)
 178       snew(contact_time,isize[1]);
 179   }
 180   if (ePBC != epbcNONE)
 181     snew(pbc,1);
 182   else
 183     pbc = NULL;
 184
 185   gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 186   do {
 187     /* initialisation for correct distance calculations */
 188     if (pbc) {
 189       set_pbc(pbc,ePBC,box);
 190       /* make molecules whole again */
 191       gmx_rmpbc(gpbc,natoms,box,x);
 192     }
 193     /* calculate center of masses */
 194     for(g=0;(g<ngrps);g++) {
 195       if (isize[g] == 1) {
 196         copy_rvec(x[index[g][0]],com[g]);
 197       } else {
 198         for(d=0;(d<DIM);d++) {
 199           com[g][d]=0;
 200           for(i=0;(i<isize[g]);i++) {
 201             com[g][d] += x[index[g][i]][d] * top->atoms.atom[index[g][i]].m;
 202           }
 203           com[g][d] /= mass[g];
 204         }
 205       }
 206     }
 207
 208     if (!bCutoff) {
 209       /* write to output */
 210       fprintf(fp,"%12.7f ",t);
 211       for(g=0;(g<ngrps/2);g++) {
 212         if (pbc)
 213           pbc_dx(pbc,com[2*g],com[2*g+1],dx);
 214         else
 215           rvec_sub(com[2*g],com[2*g+1],dx);
 216
 217         fprintf(fp,"%12.7f %12.7f %12.7f %12.7f",
 218                 norm(dx),dx[XX],dx[YY],dx[ZZ]);
 219       }
 220       fprintf(fp,"\n");
 221     } else {
 222       for(i=0;(i<isize[1]);i++) {
 223         j=index[1][i];
 224         if (pbc)
 225           pbc_dx(pbc,x[j],com[0],dx);
 226         else
 227           rvec_sub(x[j],com[0],dx);
 228
 229         dist2 = norm2(dx);
 230         if (dist2<cut2) {
 231           if (bPrintDist) {
 232             res=top->atoms.atom[j].resind;
 233             fprintf(stdout,"\rt: %g  %d %s %d %s  %g (nm)\n",
 234                     t,top->atoms.resinfo[res].nr,*top->atoms.resinfo[res].name,
 235                     j+1,*top->atoms.atomname[j],sqrt(dist2));
 236           }
 237           if (bLifeTime)
 238             contact_time[i]++;
 239         } else {
 240           if (bLifeTime) {
 241             if (contact_time[i]) {
 242               add_contact_time(&ccount,&ccount_nalloc,contact_time[i]-1);
 243               contact_time[i] = 0;
 244             }
 245           }
 246         }
 247       }
 248     }
 249
 250     teller++;
 251   } while (read_next_x(oenv,status,&t,natoms,x,box));
 252   gmx_rmpbc_done(gpbc);
 253
 254   if (!bCutoff)
 255     ffclose(fp);
 256
 257   close_trj(status);
 258
 259   if (bCutoff && bLifeTime) {
 260     /* Add the contacts still present in the last frame */
 261     for(i=0; i<isize[1]; i++)
 262       if (contact_time[i])
 263         add_contact_time(&ccount,&ccount_nalloc,contact_time[i]-1);
 264
 265     sprintf(buf,"%s - %s within %g nm",
 266             grpname[0],grpname[1],cut);
 267     fp = xvgropen(opt2fn("-lt",NFILE,fnm),
 268                   buf,"Time (ps)","Number of contacts",oenv);
 269     for(i=0; i<min(ccount_nalloc,teller-1); i++) {
 270       /* Account for all subintervals of longer intervals */
 271       sum = 0;
 272       for(j=i; j<ccount_nalloc; j++)
 273         sum += (j-i+1)*ccount[j];
 274
 275       fprintf(fp,"%10.3f %10.3f\n",i*(t-t0)/(teller-1),sum/(double)(teller-i));
 276     }
 277     ffclose(fp);
 278   }
 279
 280   thanx(stderr);
 281   return 0;
 282 }