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     {
  62         srenew(*ccount, t+2);
  63         for (i = *ccount_nalloc; i < t+2; i++)
  64         {
  65             (*ccount)[i] = 0;
  66         }
  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     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     {
 162         mass[g] = 0;
 163         for (i = 0; (i < isize[g]); i++)
 164         {
 165             if (index[g][i] > natoms_groups)
 166             {
 167                 natoms_groups = index[g][i];
 168             }
 169             if (index[g][i] >= top->atoms.nr)
 170             {
 171                 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);
 172             }
 173             mass[g] += top->atoms.atom[index[g][i]].m;
 174         }
 175     }
 176     /* The number is one more than the highest index */
 177     natoms_groups++;
 178
 179     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 180     t0     = t;
 181
 182     if (natoms_groups > natoms)
 183     {
 184         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);
 185     }
 186
 187     if (!bCutoff)
 188     {
 189         /* open output file */
 190         fp = xvgropen(ftp2fn(efXVG, NFILE, fnm),
 191                       "Distance", "Time (ps)", "Distance (nm)", oenv);
 192         xvgr_legend(fp, 4, leg, oenv);
 193     }
 194     else
 195     {
 196         ngrps = 1;
 197         if (bLifeTime)
 198         {
 199             snew(contact_time, isize[1]);
 200         }
 201     }
 202     if (ePBC != epbcNONE)
 203     {
 204         snew(pbc, 1);
 205     }
 206     else
 207     {
 208         pbc = NULL;
 209     }
 210
 211     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms, box);
 212     do
 213     {
 214         /* initialisation for correct distance calculations */
 215         if (pbc)
 216         {
 217             set_pbc(pbc, ePBC, box);
 218             /* make molecules whole again */
 219             gmx_rmpbc(gpbc, natoms, box, x);
 220         }
 221         /* calculate center of masses */
 222         for (g = 0; (g < ngrps); g++)
 223         {
 224             if (isize[g] == 1)
 225             {
 226                 copy_rvec(x[index[g][0]], com[g]);
 227             }
 228             else
 229             {
 230                 for (d = 0; (d < DIM); d++)
 231                 {
 232                     com[g][d] = 0;
 233                     for (i = 0; (i < isize[g]); i++)
 234                     {
 235                         com[g][d] += x[index[g][i]][d] * top->atoms.atom[index[g][i]].m;
 236                     }
 237                     com[g][d] /= mass[g];
 238                 }
 239             }
 240         }
 241
 242         if (!bCutoff)
 243         {
 244             /* write to output */
 245             fprintf(fp, "%12.7f ", t);
 246             for (g = 0; (g < ngrps/2); g++)
 247             {
 248                 if (pbc && (!bIntra))
 249                 {
 250                     pbc_dx(pbc, com[2*g], com[2*g+1], dx);
 251                 }
 252                 else
 253                 {
 254                     rvec_sub(com[2*g], com[2*g+1], dx);
 255                 }
 256
 257                 fprintf(fp, "%12.7f %12.7f %12.7f %12.7f",
 258                         norm(dx), dx[XX], dx[YY], dx[ZZ]);
 259             }
 260             fprintf(fp, "\n");
 261         }
 262         else
 263         {
 264             for (i = 0; (i < isize[1]); i++)
 265             {
 266                 j = index[1][i];
 267                 if (pbc && (!bIntra))
 268                 {
 269                     pbc_dx(pbc, x[j], com[0], dx);
 270                 }
 271                 else
 272                 {
 273                     rvec_sub(x[j], com[0], dx);
 274                 }
 275
 276                 dist2 = norm2(dx);
 277                 if (dist2 < cut2)
 278                 {
 279                     if (bPrintDist)
 280                     {
 281                         res = top->atoms.atom[j].resind;
 282                         fprintf(stdout, "\rt: %g  %d %s %d %s  %g (nm)\n",
 283                                 t, top->atoms.resinfo[res].nr, *top->atoms.resinfo[res].name,
 284                                 j+1, *top->atoms.atomname[j], sqrt(dist2));
 285                     }
 286                     if (bLifeTime)
 287                     {
 288                         contact_time[i]++;
 289                     }
 290                 }
 291                 else
 292                 {
 293                     if (bLifeTime)
 294                     {
 295                         if (contact_time[i])
 296                         {
 297                             add_contact_time(&ccount, &ccount_nalloc, contact_time[i]-1);
 298                             contact_time[i] = 0;
 299                         }
 300                     }
 301                 }
 302             }
 303         }
 304
 305         teller++;
 306     }
 307     while (read_next_x(oenv, status, &t, natoms, x, box));
 308     gmx_rmpbc_done(gpbc);
 309
 310     if (!bCutoff)
 311     {
 312         ffclose(fp);
 313     }
 314
 315     close_trj(status);
 316
 317     if (bCutoff && bLifeTime)
 318     {
 319         /* Add the contacts still present in the last frame */
 320         for (i = 0; i < isize[1]; i++)
 321         {
 322             if (contact_time[i])
 323             {
 324                 add_contact_time(&ccount, &ccount_nalloc, contact_time[i]-1);
 325             }
 326         }
 327
 328         sprintf(buf, "%s - %s within %g nm",
 329                 grpname[0], grpname[1], cut);
 330         fp = xvgropen(opt2fn("-lt", NFILE, fnm),
 331                       buf, "Time (ps)", "Number of contacts", oenv);
 332         for (i = 0; i < min(ccount_nalloc, teller-1); i++)
 333         {
 334             /* Account for all subintervals of longer intervals */
 335             sum = 0;
 336             for (j = i; j < ccount_nalloc; j++)
 337             {
 338                 sum += (j-i+1)*ccount[j];
 339             }
 340
 341             fprintf(fp, "%10.3f %10.3f\n", i*(t-t0)/(teller-1), sum/(double)(teller-i));
 342         }
 343         ffclose(fp);
 344     }
 345
 346     thanx(stderr);
 347     return 0;
 348 }