src/tools/gmx_dist.c

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