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     "[IT]x[it]-, [IT]y[it]-, and [IT]z[it]-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 continuously present. The [TT]-intra[tt] switch enables",
  81     "calculations of intramolecular distances avoiding distance calculation to its",
  82     "periodic images. For a proper function, the molecule in the input trajectory",
  83     "should be whole (e.g. by preprocessing with [TT]trjconv -pbc[tt]) or a matching",
  84     "topology should be provided. The [TT]-intra[tt] switch will only give",
  85     "meaningful results for intramolecular and not intermolecular distances.[PAR]",
  86     "Other programs that calculate distances are [TT]g_mindist[tt]",
  87     "and [TT]g_bond[tt]."
  88   };
  89
  90   t_topology *top=NULL;
  91   int  ePBC;
  92   real t,t0,cut2,dist2;
  93   rvec *x=NULL,*v=NULL,dx;
  94   matrix box;
  95   t_trxstatus *status;
  96   int natoms;
  97
  98   int g,d,i,j,res,teller=0;
  99   atom_id aid;
 100
 101   int     ngrps;     /* the number of index groups */
 102   atom_id **index, natoms_groups;   /* the index for the atom numbers */
 103   int     *isize;    /* the size of each group */
 104   char    **grpname; /* the name of each group */
 105   rvec    *com;
 106   real    *mass;
 107   FILE    *fp=NULL,*fplt=NULL;
 108   gmx_bool    bCutoff,bPrintDist,bLifeTime,bIntra=FALSE;
 109   t_pbc   *pbc;
 110   int     *contact_time=NULL,*ccount=NULL,ccount_nalloc=0,sum;
 111   char    buf[STRLEN];
 112   output_env_t oenv;
 113   gmx_rmpbc_t  gpbc=NULL;
 114
 115   const char *leg[4] = { "|d|","d\\sx\\N","d\\sy\\N","d\\sz\\N" };
 116
 117   static real cut=0;
 118
 119   t_pargs pa[] = {
 120     { "-intra",      FALSE, etBOOL, {&bIntra},
 121       "Calculate distances without considering periodic boundaries, e.g. intramolecular." },
 122     { "-dist",      FALSE, etREAL, {&cut},
 123       "Print all atoms in group 2 closer than dist to the center of mass of group 1" }
 124   };
 125 #define NPA asize(pa)
 126
 127   t_filenm fnm[] = {
 128     { efTRX, "-f", NULL, ffREAD },
 129     { efTPX, NULL, NULL, ffREAD },
 130     { efNDX, NULL, NULL, ffOPTRD },
 131     { efXVG, NULL, "dist", ffOPTWR },
 132     { efXVG, "-lt", "lifetime", ffOPTWR },
 133   };
 134 #define NFILE asize(fnm)
 135
 136
 137   CopyRight(stderr,argv[0]);
 138
 139   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
 140                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 141
 142   bCutoff = opt2parg_bSet("-dist",NPA,pa);
 143   cut2 = cut*cut;
 144   bLifeTime = opt2bSet("-lt",NFILE,fnm);
 145   bPrintDist = (bCutoff && !bLifeTime);
 146
 147   top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
 148
 149   /* read index files */
 150   ngrps = 2;
 151   snew(com,ngrps);
 152   snew(grpname,ngrps);
 153   snew(index,ngrps);
 154   snew(isize,ngrps);
 155   get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),ngrps,isize,index,grpname);
 156
 157   /* calculate mass */
 158   natoms_groups=0;
 159   snew(mass,ngrps);
 160   for(g=0;(g<ngrps);g++) {
 161     mass[g]=0;
 162     for(i=0;(i<isize[g]);i++) {
 163       if (index[g][i]>natoms_groups)
 164         natoms_groups=index[g][i];
 165       if (index[g][i] >= top->atoms.nr)
 166         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);
 167       mass[g]+=top->atoms.atom[index[g][i]].m;
 168     }
 169   }
 170   /* The number is one more than the highest index */
 171   natoms_groups++;
 172
 173   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 174   t0 = t;
 175
 176   if (natoms_groups>natoms)
 177     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);
 178
 179   if (!bCutoff) {
 180     /* open output file */
 181     fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),
 182                   "Distance","Time (ps)","Distance (nm)",oenv);
 183     xvgr_legend(fp,4,leg,oenv);
 184   } else {
 185     ngrps = 1;
 186     if (bLifeTime)
 187       snew(contact_time,isize[1]);
 188   }
 189   if (ePBC != epbcNONE)
 190     snew(pbc,1);
 191   else
 192     pbc = NULL;
 193
 194   gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 195   do {
 196     /* initialisation for correct distance calculations */
 197     if (pbc) {
 198       set_pbc(pbc,ePBC,box);
 199       /* make molecules whole again */
 200       gmx_rmpbc(gpbc,natoms,box,x);
 201     }
 202     /* calculate center of masses */
 203     for(g=0;(g<ngrps);g++) {
 204       if (isize[g] == 1) {
 205         copy_rvec(x[index[g][0]],com[g]);
 206       } else {
 207         for(d=0;(d<DIM);d++) {
 208           com[g][d]=0;
 209           for(i=0;(i<isize[g]);i++) {
 210             com[g][d] += x[index[g][i]][d] * top->atoms.atom[index[g][i]].m;
 211           }
 212           com[g][d] /= mass[g];
 213         }
 214       }
 215     }
 216
 217     if (!bCutoff) {
 218       /* write to output */
 219       fprintf(fp,"%12.7f ",t);
 220       for(g=0;(g<ngrps/2);g++) {
 221         if (pbc && (!bIntra))
 222           pbc_dx(pbc,com[2*g],com[2*g+1],dx);
 223         else
 224           rvec_sub(com[2*g],com[2*g+1],dx);
 225
 226         fprintf(fp,"%12.7f %12.7f %12.7f %12.7f",
 227                 norm(dx),dx[XX],dx[YY],dx[ZZ]);
 228       }
 229       fprintf(fp,"\n");
 230     } else {
 231       for(i=0;(i<isize[1]);i++) {
 232         j=index[1][i];
 233         if (pbc && (!bIntra))
 234           pbc_dx(pbc,x[j],com[0],dx);
 235         else
 236           rvec_sub(x[j],com[0],dx);
 237
 238         dist2 = norm2(dx);
 239         if (dist2<cut2) {
 240           if (bPrintDist) {
 241             res=top->atoms.atom[j].resind;
 242             fprintf(stdout,"\rt: %g  %d %s %d %s  %g (nm)\n",
 243                     t,top->atoms.resinfo[res].nr,*top->atoms.resinfo[res].name,
 244                     j+1,*top->atoms.atomname[j],sqrt(dist2));
 245           }
 246           if (bLifeTime)
 247             contact_time[i]++;
 248         } else {
 249           if (bLifeTime) {
 250             if (contact_time[i]) {
 251               add_contact_time(&ccount,&ccount_nalloc,contact_time[i]-1);
 252               contact_time[i] = 0;
 253             }
 254           }
 255         }
 256       }
 257     }
 258
 259     teller++;
 260   } while (read_next_x(oenv,status,&t,natoms,x,box));
 261   gmx_rmpbc_done(gpbc);
 262
 263   if (!bCutoff)
 264     ffclose(fp);
 265
 266   close_trj(status);
 267
 268   if (bCutoff && bLifeTime) {
 269     /* Add the contacts still present in the last frame */
 270     for(i=0; i<isize[1]; i++)
 271       if (contact_time[i])
 272         add_contact_time(&ccount,&ccount_nalloc,contact_time[i]-1);
 273
 274     sprintf(buf,"%s - %s within %g nm",
 275             grpname[0],grpname[1],cut);
 276     fp = xvgropen(opt2fn("-lt",NFILE,fnm),
 277                   buf,"Time (ps)","Number of contacts",oenv);
 278     for(i=0; i<min(ccount_nalloc,teller-1); i++) {
 279       /* Account for all subintervals of longer intervals */
 280       sum = 0;
 281       for(j=i; j<ccount_nalloc; j++)
 282         sum += (j-i+1)*ccount[j];
 283
 284       fprintf(fp,"%10.3f %10.3f\n",i*(t-t0)/(teller-1),sum/(double)(teller-i));
 285     }
 286     ffclose(fp);
 287   }
 288
 289   thanx(stderr);
 290   return 0;
 291 }