src/gromacs/gmxana/gmx_mindist.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  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <math.h>
  42 #include <stdlib.h>
  43
  44 #include "sysstuff.h"
  45 #include <string.h>
  46 #include "typedefs.h"
  47 #include "smalloc.h"
  48 #include "macros.h"
  49 #include "vec.h"
  50 #include "xvgr.h"
  51 #include "pbc.h"
  52 #include "gromacs/fileio/futil.h"
  53 #include "gromacs/commandline/pargs.h"
  54 #include "index.h"
  55 #include "gromacs/fileio/tpxio.h"
  56 #include "gromacs/fileio/trxio.h"
  57 #include "rmpbc.h"
  58 #include "gmx_ana.h"
  59
  60
  61 static void periodic_dist(matrix box, rvec x[], int n, atom_id index[],
  62                           real *rmin, real *rmax, int *min_ind)
  63 {
  64 #define NSHIFT 26
  65     int  sx, sy, sz, i, j, s;
  66     real sqr_box, r2min, r2max, r2;
  67     rvec shift[NSHIFT], d0, d;
  68
  69     sqr_box = sqr(min(norm(box[XX]), min(norm(box[YY]), norm(box[ZZ]))));
  70
  71     s = 0;
  72     for (sz = -1; sz <= 1; sz++)
  73     {
  74         for (sy = -1; sy <= 1; sy++)
  75         {
  76             for (sx = -1; sx <= 1; sx++)
  77             {
  78                 if (sx != 0 || sy != 0 || sz != 0)
  79                 {
  80                     for (i = 0; i < DIM; i++)
  81                     {
  82                         shift[s][i] = sx*box[XX][i]+sy*box[YY][i]+sz*box[ZZ][i];
  83                     }
  84                     s++;
  85                 }
  86             }
  87         }
  88     }
  89
  90     r2min = sqr_box;
  91     r2max = 0;
  92
  93     for (i = 0; i < n; i++)
  94     {
  95         for (j = i+1; j < n; j++)
  96         {
  97             rvec_sub(x[index[i]], x[index[j]], d0);
  98             r2 = norm2(d0);
  99             if (r2 > r2max)
 100             {
 101                 r2max = r2;
 102             }
 103             for (s = 0; s < NSHIFT; s++)
 104             {
 105                 rvec_add(d0, shift[s], d);
 106                 r2 = norm2(d);
 107                 if (r2 < r2min)
 108                 {
 109                     r2min      = r2;
 110                     min_ind[0] = i;
 111                     min_ind[1] = j;
 112                 }
 113             }
 114         }
 115     }
 116
 117     *rmin = sqrt(r2min);
 118     *rmax = sqrt(r2max);
 119 }
 120
 121 static void periodic_mindist_plot(const char *trxfn, const char *outfn,
 122                                   t_topology *top, int ePBC,
 123                                   int n, atom_id index[], gmx_bool bSplit,
 124                                   const output_env_t oenv)
 125 {
 126     FILE        *out;
 127     const char  *leg[5] = { "min per.", "max int.", "box1", "box2", "box3" };
 128     t_trxstatus *status;
 129     real         t;
 130     rvec        *x;
 131     matrix       box;
 132     int          natoms, ind_min[2] = {0, 0}, ind_mini = 0, ind_minj = 0;
 133     real         r, rmin, rmax, rmint, tmint;
 134     gmx_bool     bFirst;
 135     gmx_rmpbc_t  gpbc = NULL;
 136
 137     natoms = read_first_x(oenv, &status, trxfn, &t, &x, box);
 138
 139     check_index(NULL, n, index, NULL, natoms);
 140
 141     out = xvgropen(outfn, "Minimum distance to periodic image",
 142                    output_env_get_time_label(oenv), "Distance (nm)", oenv);
 143     if (output_env_get_print_xvgr_codes(oenv))
 144     {
 145         fprintf(out, "@ subtitle \"and maximum internal distance\"\n");
 146     }
 147     xvgr_legend(out, 5, leg, oenv);
 148
 149     rmint = box[XX][XX];
 150     tmint = 0;
 151
 152     if (NULL != top)
 153     {
 154         gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
 155     }
 156
 157     bFirst = TRUE;
 158     do
 159     {
 160         if (NULL != top)
 161         {
 162             gmx_rmpbc(gpbc, natoms, box, x);
 163         }
 164
 165         periodic_dist(box, x, n, index, &rmin, &rmax, ind_min);
 166         if (rmin < rmint)
 167         {
 168             rmint    = rmin;
 169             tmint    = t;
 170             ind_mini = ind_min[0];
 171             ind_minj = ind_min[1];
 172         }
 173         if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5)
 174         {
 175             fprintf(out, "&\n");
 176         }
 177         fprintf(out, "\t%g\t%6.3f %6.3f %6.3f %6.3f %6.3f\n",
 178                 output_env_conv_time(oenv, t), rmin, rmax, norm(box[0]), norm(box[1]), norm(box[2]));
 179         bFirst = FALSE;
 180     }
 181     while (read_next_x(oenv, status, &t, x, box));
 182
 183     if (NULL != top)
 184     {
 185         gmx_rmpbc_done(gpbc);
 186     }
 187
 188     gmx_ffclose(out);
 189
 190     fprintf(stdout,
 191             "\nThe shortest periodic distance is %g (nm) at time %g (%s),\n"
 192             "between atoms %d and %d\n",
 193             rmint, output_env_conv_time(oenv, tmint), output_env_get_time_unit(oenv),
 194             index[ind_mini]+1, index[ind_minj]+1);
 195 }
 196
 197 static void calc_dist(real rcut, gmx_bool bPBC, int ePBC, matrix box, rvec x[],
 198                       int nx1, int nx2, atom_id index1[], atom_id index2[],
 199                       gmx_bool bGroup,
 200                       real *rmin, real *rmax, int *nmin, int *nmax,
 201                       int *ixmin, int *jxmin, int *ixmax, int *jxmax)
 202 {
 203     int      i, j, i0 = 0, j1;
 204     int      ix, jx;
 205     atom_id *index3;
 206     rvec     dx;
 207     real     r2, rmin2, rmax2, rcut2;
 208     t_pbc    pbc;
 209     int      nmin_j, nmax_j;
 210
 211     *ixmin = -1;
 212     *jxmin = -1;
 213     *ixmax = -1;
 214     *jxmax = -1;
 215     *nmin  = 0;
 216     *nmax  = 0;
 217
 218     rcut2 = sqr(rcut);
 219
 220     /* Must init pbc every step because of pressure coupling */
 221     if (bPBC)
 222     {
 223         set_pbc(&pbc, ePBC, box);
 224     }
 225     if (index2)
 226     {
 227         i0     = 0;
 228         j1     = nx2;
 229         index3 = index2;
 230     }
 231     else
 232     {
 233         j1     = nx1;
 234         index3 = index1;
 235     }
 236
 237     rmin2 = 1e12;
 238     rmax2 = -1e12;
 239
 240     for (j = 0; (j < j1); j++)
 241     {
 242         jx = index3[j];
 243         if (index2 == NULL)
 244         {
 245             i0 = j + 1;
 246         }
 247         nmin_j = 0;
 248         nmax_j = 0;
 249         for (i = i0; (i < nx1); i++)
 250         {
 251             ix = index1[i];
 252             if (ix != jx)
 253             {
 254                 if (bPBC)
 255                 {
 256                     pbc_dx(&pbc, x[ix], x[jx], dx);
 257                 }
 258                 else
 259                 {
 260                     rvec_sub(x[ix], x[jx], dx);
 261                 }
 262                 r2 = iprod(dx, dx);
 263                 if (r2 < rmin2)
 264                 {
 265                     rmin2  = r2;
 266                     *ixmin = ix;
 267                     *jxmin = jx;
 268                 }
 269                 if (r2 > rmax2)
 270                 {
 271                     rmax2  = r2;
 272                     *ixmax = ix;
 273                     *jxmax = jx;
 274                 }
 275                 if (r2 <= rcut2)
 276                 {
 277                     nmin_j++;
 278                 }
 279                 else if (r2 > rcut2)
 280                 {
 281                     nmax_j++;
 282                 }
 283             }
 284         }
 285         if (bGroup)
 286         {
 287             if (nmin_j > 0)
 288             {
 289                 (*nmin)++;
 290             }
 291             if (nmax_j > 0)
 292             {
 293                 (*nmax)++;
 294             }
 295         }
 296         else
 297         {
 298             *nmin += nmin_j;
 299             *nmax += nmax_j;
 300         }
 301     }
 302     *rmin = sqrt(rmin2);
 303     *rmax = sqrt(rmax2);
 304 }
 305
 306 void dist_plot(const char *fn, const char *afile, const char *dfile,
 307                const char *nfile, const char *rfile, const char *xfile,
 308                real rcut, gmx_bool bMat, t_atoms *atoms,
 309                int ng, atom_id *index[], int gnx[], char *grpn[], gmx_bool bSplit,
 310                gmx_bool bMin, int nres, atom_id *residue, gmx_bool bPBC, int ePBC,
 311                gmx_bool bGroup, gmx_bool bEachResEachTime, gmx_bool bPrintResName,
 312                const output_env_t oenv)
 313 {
 314     FILE            *atm, *dist, *num;
 315     t_trxstatus     *trxout;
 316     char             buf[256];
 317     char           **leg;
 318     real             t, dmin, dmax, **mindres = NULL, **maxdres = NULL;
 319     int              nmin, nmax;
 320     t_trxstatus     *status;
 321     int              i = -1, j, k, natoms;
 322     int              min1, min2, max1, max2, min1r, min2r, max1r, max2r;
 323     atom_id          oindex[2];
 324     rvec            *x0;
 325     matrix           box;
 326     t_trxframe       frout;
 327     gmx_bool         bFirst;
 328     FILE            *respertime = NULL;
 329
 330     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 331     {
 332         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 333     }
 334
 335     sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
 336     dist = xvgropen(dfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
 337     sprintf(buf, "Number of Contacts %s %g nm", bMin ? "<" : ">", rcut);
 338     num    = nfile ? xvgropen(nfile, buf, output_env_get_time_label(oenv), "Number", oenv) : NULL;
 339     atm    = afile ? gmx_ffopen(afile, "w") : NULL;
 340     trxout = xfile ? open_trx(xfile, "w") : NULL;
 341
 342     if (bMat)
 343     {
 344         if (ng == 1)
 345         {
 346             snew(leg, 1);
 347             sprintf(buf, "Internal in %s", grpn[0]);
 348             leg[0] = strdup(buf);
 349             xvgr_legend(dist, 0, (const char**)leg, oenv);
 350             if (num)
 351             {
 352                 xvgr_legend(num, 0, (const char**)leg, oenv);
 353             }
 354         }
 355         else
 356         {
 357             snew(leg, (ng*(ng-1))/2);
 358             for (i = j = 0; (i < ng-1); i++)
 359             {
 360                 for (k = i+1; (k < ng); k++, j++)
 361                 {
 362                     sprintf(buf, "%s-%s", grpn[i], grpn[k]);
 363                     leg[j] = strdup(buf);
 364                 }
 365             }
 366             xvgr_legend(dist, j, (const char**)leg, oenv);
 367             if (num)
 368             {
 369                 xvgr_legend(num, j, (const char**)leg, oenv);
 370             }
 371         }
 372     }
 373     else
 374     {
 375         snew(leg, ng-1);
 376         for (i = 0; (i < ng-1); i++)
 377         {
 378             sprintf(buf, "%s-%s", grpn[0], grpn[i+1]);
 379             leg[i] = strdup(buf);
 380         }
 381         xvgr_legend(dist, ng-1, (const char**)leg, oenv);
 382         if (num)
 383         {
 384             xvgr_legend(num, ng-1, (const char**)leg, oenv);
 385         }
 386     }
 387
 388     if (bEachResEachTime)
 389     {
 390         sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
 391         respertime = xvgropen(rfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
 392         xvgr_legend(respertime, ng-1, (const char**)leg, oenv);
 393         if (bPrintResName)
 394         {
 395             fprintf(respertime, "# ");
 396         }
 397         for (j = 0; j < nres; j++)
 398         {
 399             fprintf(respertime, "%s%d ", *(atoms->resinfo[atoms->atom[index[0][residue[j]]].resind].name), atoms->atom[index[0][residue[j]]].resind);
 400         }
 401         fprintf(respertime, "\n");
 402
 403     }
 404
 405     j = 0;
 406     if (nres)
 407     {
 408         snew(mindres, ng-1);
 409         snew(maxdres, ng-1);
 410         for (i = 1; i < ng; i++)
 411         {
 412             snew(mindres[i-1], nres);
 413             snew(maxdres[i-1], nres);
 414             for (j = 0; j < nres; j++)
 415             {
 416                 mindres[i-1][j] = 1e6;
 417             }
 418             /* maxdres[*][*] is already 0 */
 419         }
 420     }
 421     bFirst = TRUE;
 422     do
 423     {
 424         if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5)
 425         {
 426             fprintf(dist, "&\n");
 427             if (num)
 428             {
 429                 fprintf(num, "&\n");
 430             }
 431             if (atm)
 432             {
 433                 fprintf(atm, "&\n");
 434             }
 435         }
 436         fprintf(dist, "%12e", output_env_conv_time(oenv, t));
 437         if (num)
 438         {
 439             fprintf(num, "%12e", output_env_conv_time(oenv, t));
 440         }
 441
 442         if (bMat)
 443         {
 444             if (ng == 1)
 445             {
 446                 calc_dist(rcut, bPBC, ePBC, box, x0, gnx[0], gnx[0], index[0], index[0], bGroup,
 447                           &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
 448                 fprintf(dist, "  %12e", bMin ? dmin : dmax);
 449                 if (num)
 450                 {
 451                     fprintf(num, "  %8d", bMin ? nmin : nmax);
 452                 }
 453             }
 454             else
 455             {
 456                 for (i = 0; (i < ng-1); i++)
 457                 {
 458                     for (k = i+1; (k < ng); k++)
 459                     {
 460                         calc_dist(rcut, bPBC, ePBC, box, x0, gnx[i], gnx[k], index[i], index[k],
 461                                   bGroup, &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
 462                         fprintf(dist, "  %12e", bMin ? dmin : dmax);
 463                         if (num)
 464                         {
 465                             fprintf(num, "  %8d", bMin ? nmin : nmax);
 466                         }
 467                     }
 468                 }
 469             }
 470         }
 471         else
 472         {
 473             for (i = 1; (i < ng); i++)
 474             {
 475                 calc_dist(rcut, bPBC, ePBC, box, x0, gnx[0], gnx[i], index[0], index[i], bGroup,
 476                           &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
 477                 fprintf(dist, "  %12e", bMin ? dmin : dmax);
 478                 if (num)
 479                 {
 480                     fprintf(num, "  %8d", bMin ? nmin : nmax);
 481                 }
 482                 if (nres)
 483                 {
 484                     for (j = 0; j < nres; j++)
 485                     {
 486                         calc_dist(rcut, bPBC, ePBC, box, x0, residue[j+1]-residue[j], gnx[i],
 487                                   &(index[0][residue[j]]), index[i], bGroup,
 488                                   &dmin, &dmax, &nmin, &nmax, &min1r, &min2r, &max1r, &max2r);
 489                         mindres[i-1][j] = min(mindres[i-1][j], dmin);
 490                         maxdres[i-1][j] = max(maxdres[i-1][j], dmax);
 491                     }
 492                 }
 493             }
 494         }
 495         fprintf(dist, "\n");
 496         if (num)
 497         {
 498             fprintf(num, "\n");
 499         }
 500         if ( (bMin ? min1 : max1) != -1)
 501         {
 502             if (atm)
 503             {
 504                 fprintf(atm, "%12e  %12d  %12d\n",
 505                         output_env_conv_time(oenv, t), 1+(bMin ? min1 : max1),
 506                         1+(bMin ? min2 : max2));
 507             }
 508         }
 509
 510         if (trxout)
 511         {
 512             oindex[0] = bMin ? min1 : max1;
 513             oindex[1] = bMin ? min2 : max2;
 514             write_trx(trxout, 2, oindex, atoms, i, t, box, x0, NULL, NULL);
 515         }
 516         bFirst = FALSE;
 517         /*dmin should be minimum distance for residue and group*/
 518         if (bEachResEachTime)
 519         {
 520             fprintf(respertime, "%12e", t);
 521             for (i = 1; i < ng; i++)
 522             {
 523                 for (j = 0; j < nres; j++)
 524                 {
 525                     fprintf(respertime, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
 526                     /*reset distances for next time point*/
 527                     mindres[i-1][j] = 1e6;
 528                     maxdres[i-1][j] = 0;
 529                 }
 530             }
 531             fprintf(respertime, "\n");
 532         }
 533     }
 534     while (read_next_x(oenv, status, &t, x0, box));
 535
 536     close_trj(status);
 537     gmx_ffclose(dist);
 538     if (num)
 539     {
 540         gmx_ffclose(num);
 541     }
 542     if (atm)
 543     {
 544         gmx_ffclose(atm);
 545     }
 546     if (trxout)
 547     {
 548         close_trx(trxout);
 549     }
 550
 551     if (nres && !bEachResEachTime)
 552     {
 553         FILE *res;
 554
 555         sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
 556         res = xvgropen(rfile, buf, "Residue (#)", "Distance (nm)", oenv);
 557         xvgr_legend(res, ng-1, (const char**)leg, oenv);
 558         for (j = 0; j < nres; j++)
 559         {
 560             fprintf(res, "%4d", j+1);
 561             for (i = 1; i < ng; i++)
 562             {
 563                 fprintf(res, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
 564             }
 565             fprintf(res, "\n");
 566         }
 567     }
 568
 569     sfree(x0);
 570 }
 571
 572 int find_residues(t_atoms *atoms, int n, atom_id index[], atom_id **resindex)
 573 {
 574     int  i;
 575     int  nres = 0, resnr, presnr;
 576     int *residx;
 577
 578     /* build index of first atom numbers for each residue */
 579     presnr = NOTSET;
 580     snew(residx, atoms->nres+1);
 581     for (i = 0; i < n; i++)
 582     {
 583         resnr = atoms->atom[index[i]].resind;
 584         if (resnr != presnr)
 585         {
 586             residx[nres] = i;
 587             nres++;
 588             presnr = resnr;
 589         }
 590     }
 591     if (debug)
 592     {
 593         printf("Found %d residues out of %d (%d/%d atoms)\n",
 594                nres, atoms->nres, atoms->nr, n);
 595     }
 596     srenew(residx, nres+1);
 597     /* mark end of last residue */
 598     residx[nres] = n;
 599     *resindex    = residx;
 600     return nres;
 601 }
 602
 603 void dump_res(FILE *out, int nres, atom_id *resindex, atom_id index[])
 604 {
 605     int i, j;
 606
 607     for (i = 0; i < nres-1; i++)
 608     {
 609         fprintf(out, "Res %d (%d):", i, resindex[i+1]-resindex[i]);
 610         for (j = resindex[i]; j < resindex[i+1]; j++)
 611         {
 612             fprintf(out, " %d(%d)", j, index[j]);
 613         }
 614         fprintf(out, "\n");
 615     }
 616 }
 617
 618 int gmx_mindist(int argc, char *argv[])
 619 {
 620     const char     *desc[] = {
 621         "[THISMODULE] computes the distance between one group and a number of",
 622         "other groups. Both the minimum distance",
 623         "(between any pair of atoms from the respective groups)",
 624         "and the number of contacts within a given",
 625         "distance are written to two separate output files.",
 626         "With the [TT]-group[tt] option a contact of an atom in another group",
 627         "with multiple atoms in the first group is counted as one contact",
 628         "instead of as multiple contacts.",
 629         "With [TT]-or[tt], minimum distances to each residue in the first",
 630         "group are determined and plotted as a function of residue number.[PAR]",
 631         "With option [TT]-pi[tt] the minimum distance of a group to its",
 632         "periodic image is plotted. This is useful for checking if a protein",
 633         "has seen its periodic image during a simulation. Only one shift in",
 634         "each direction is considered, giving a total of 26 shifts.",
 635         "It also plots the maximum distance within the group and the lengths",
 636         "of the three box vectors.[PAR]",
 637         "Also [gmx-distance] calculates distances."
 638     };
 639     const char     *bugs[] = {
 640         "The [TT]-pi[tt] option is very slow."
 641     };
 642
 643     static gmx_bool bMat             = FALSE, bPI = FALSE, bSplit = FALSE, bMax = FALSE, bPBC = TRUE;
 644     static gmx_bool bGroup           = FALSE;
 645     static real     rcutoff          = 0.6;
 646     static int      ng               = 1;
 647     static gmx_bool bEachResEachTime = FALSE, bPrintResName = FALSE;
 648     t_pargs         pa[]             = {
 649         { "-matrix", FALSE, etBOOL, {&bMat},
 650           "Calculate half a matrix of group-group distances" },
 651         { "-max",    FALSE, etBOOL, {&bMax},
 652           "Calculate *maximum* distance instead of minimum" },
 653         { "-d",      FALSE, etREAL, {&rcutoff},
 654           "Distance for contacts" },
 655         { "-group",      FALSE, etBOOL, {&bGroup},
 656           "Count contacts with multiple atoms in the first group as one" },
 657         { "-pi",     FALSE, etBOOL, {&bPI},
 658           "Calculate minimum distance with periodic images" },
 659         { "-split",  FALSE, etBOOL, {&bSplit},
 660           "Split graph where time is zero" },
 661         { "-ng",       FALSE, etINT, {&ng},
 662           "Number of secondary groups to compute distance to a central group" },
 663         { "-pbc",    FALSE, etBOOL, {&bPBC},
 664           "Take periodic boundary conditions into account" },
 665         { "-respertime",  FALSE, etBOOL, {&bEachResEachTime},
 666           "When writing per-residue distances, write distance for each time point" },
 667         { "-printresname",  FALSE, etBOOL, {&bPrintResName},
 668           "Write residue names" }
 669     };
 670     output_env_t    oenv;
 671     t_topology     *top  = NULL;
 672     int             ePBC = -1;
 673     char            title[256];
 674     real            t;
 675     rvec           *x;
 676     matrix          box;
 677     gmx_bool        bTop = FALSE;
 678
 679     FILE           *atm;
 680     int             i, j, nres = 0;
 681     const char     *trxfnm, *tpsfnm, *ndxfnm, *distfnm, *numfnm, *atmfnm, *oxfnm, *resfnm;
 682     char          **grpname;
 683     int            *gnx;
 684     atom_id       **index, *residues = NULL;
 685     t_filenm        fnm[] = {
 686         { efTRX, "-f",  NULL,      ffREAD },
 687         { efTPS,  NULL, NULL,      ffOPTRD },
 688         { efNDX,  NULL, NULL,      ffOPTRD },
 689         { efXVG, "-od", "mindist",  ffWRITE },
 690         { efXVG, "-on", "numcont",  ffOPTWR },
 691         { efOUT, "-o", "atm-pair", ffOPTWR },
 692         { efTRO, "-ox", "mindist",  ffOPTWR },
 693         { efXVG, "-or", "mindistres", ffOPTWR }
 694     };
 695 #define NFILE asize(fnm)
 696
 697     if (!parse_common_args(&argc, argv,
 698                            PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
 699                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 700     {
 701         return 0;
 702     }
 703
 704     trxfnm  = ftp2fn(efTRX, NFILE, fnm);
 705     ndxfnm  = ftp2fn_null(efNDX, NFILE, fnm);
 706     distfnm = opt2fn("-od", NFILE, fnm);
 707     numfnm  = opt2fn_null("-on", NFILE, fnm);
 708     atmfnm  = ftp2fn_null(efOUT, NFILE, fnm);
 709     oxfnm   = opt2fn_null("-ox", NFILE, fnm);
 710     resfnm  = opt2fn_null("-or", NFILE, fnm);
 711     if (bPI || resfnm != NULL)
 712     {
 713         /* We need a tps file */
 714         tpsfnm = ftp2fn(efTPS, NFILE, fnm);
 715     }
 716     else
 717     {
 718         tpsfnm = ftp2fn_null(efTPS, NFILE, fnm);
 719     }
 720
 721     if (!tpsfnm && !ndxfnm)
 722     {
 723         gmx_fatal(FARGS, "You have to specify either the index file or a tpr file");
 724     }
 725
 726     if (bPI)
 727     {
 728         ng = 1;
 729         fprintf(stderr, "Choose a group for distance calculation\n");
 730     }
 731     else if (!bMat)
 732     {
 733         ng++;
 734     }
 735
 736     snew(gnx, ng);
 737     snew(index, ng);
 738     snew(grpname, ng);
 739
 740     if (tpsfnm || resfnm || !ndxfnm)
 741     {
 742         snew(top, 1);
 743         bTop = read_tps_conf(tpsfnm, title, top, &ePBC, &x, NULL, box, FALSE);
 744         if (bPI && !bTop)
 745         {
 746             printf("\nWARNING: Without a run input file a trajectory with broken molecules will not give the correct periodic image distance\n\n");
 747         }
 748     }
 749     get_index(top ? &(top->atoms) : NULL, ndxfnm, ng, gnx, index, grpname);
 750
 751     if (bMat && (ng == 1))
 752     {
 753         ng = gnx[0];
 754         printf("Special case: making distance matrix between all atoms in group %s\n",
 755                grpname[0]);
 756         srenew(gnx, ng);
 757         srenew(index, ng);
 758         srenew(grpname, ng);
 759         for (i = 1; (i < ng); i++)
 760         {
 761             gnx[i]      = 1;
 762             grpname[i]  = grpname[0];
 763             snew(index[i], 1);
 764             index[i][0] = index[0][i];
 765         }
 766         gnx[0] = 1;
 767     }
 768
 769     if (resfnm)
 770     {
 771         nres = find_residues(top ? &(top->atoms) : NULL,
 772                              gnx[0], index[0], &residues);
 773         if (debug)
 774         {
 775             dump_res(debug, nres, residues, index[0]);
 776         }
 777     }
 778
 779     if (bPI)
 780     {
 781         periodic_mindist_plot(trxfnm, distfnm, top, ePBC, gnx[0], index[0], bSplit, oenv);
 782     }
 783     else
 784     {
 785         dist_plot(trxfnm, atmfnm, distfnm, numfnm, resfnm, oxfnm,
 786                   rcutoff, bMat, top ? &(top->atoms) : NULL,
 787                   ng, index, gnx, grpname, bSplit, !bMax, nres, residues, bPBC, ePBC,
 788                   bGroup, bEachResEachTime, bPrintResName, oenv);
 789     }
 790
 791     do_view(oenv, distfnm, "-nxy");
 792     if (!bPI)
 793     {
 794         do_view(oenv, numfnm, "-nxy");
 795     }
 796
 797     return 0;
 798 }