src/gromacs/gmxana/gmx_mindist.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
  39 #include <math.h>
  40 #include <stdlib.h>
  41
  42 #include "sysstuff.h"
  43 #include <string.h>
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "vec.h"
  48 #include "xvgr.h"
  49 #include "pbc.h"
  50 #include "copyrite.h"
  51 #include "futil.h"
  52 #include "statutil.h"
  53 #include "index.h"
  54 #include "tpxio.h"
  55 #include "rmpbc.h"
  56 #include "xtcio.h"
  57 #include "gmx_ana.h"
  58
  59
  60 static void periodic_dist(matrix box, rvec x[], int n, atom_id index[],
  61                           real *rmin, real *rmax, int *min_ind)
  62 {
  63 #define NSHIFT 26
  64     int  sx, sy, sz, i, j, s;
  65     real sqr_box, r2min, r2max, r2;
  66     rvec shift[NSHIFT], d0, d;
  67
  68     sqr_box = sqr(min(norm(box[XX]), min(norm(box[YY]), norm(box[ZZ]))));
  69
  70     s = 0;
  71     for (sz = -1; sz <= 1; sz++)
  72     {
  73         for (sy = -1; sy <= 1; sy++)
  74         {
  75             for (sx = -1; sx <= 1; sx++)
  76             {
  77                 if (sx != 0 || sy != 0 || sz != 0)
  78                 {
  79                     for (i = 0; i < DIM; i++)
  80                     {
  81                         shift[s][i] = sx*box[XX][i]+sy*box[YY][i]+sz*box[ZZ][i];
  82                     }
  83                     s++;
  84                 }
  85             }
  86         }
  87     }
  88
  89     r2min = sqr_box;
  90     r2max = 0;
  91
  92     for (i = 0; i < n; i++)
  93     {
  94         for (j = i+1; j < n; j++)
  95         {
  96             rvec_sub(x[index[i]], x[index[j]], d0);
  97             r2 = norm2(d0);
  98             if (r2 > r2max)
  99             {
 100                 r2max = r2;
 101             }
 102             for (s = 0; s < NSHIFT; s++)
 103             {
 104                 rvec_add(d0, shift[s], d);
 105                 r2 = norm2(d);
 106                 if (r2 < r2min)
 107                 {
 108                     r2min      = r2;
 109                     min_ind[0] = i;
 110                     min_ind[1] = j;
 111                 }
 112             }
 113         }
 114     }
 115
 116     *rmin = sqrt(r2min);
 117     *rmax = sqrt(r2max);
 118 }
 119
 120 static void periodic_mindist_plot(const char *trxfn, const char *outfn,
 121                                   t_topology *top, int ePBC,
 122                                   int n, atom_id index[], gmx_bool bSplit,
 123                                   const output_env_t oenv)
 124 {
 125     FILE        *out;
 126     const char  *leg[5] = { "min per.", "max int.", "box1", "box2", "box3" };
 127     t_trxstatus *status;
 128     real         t;
 129     rvec        *x;
 130     matrix       box;
 131     int          natoms, ind_min[2] = {0, 0}, ind_mini = 0, ind_minj = 0;
 132     real         r, rmin, rmax, rmint, tmint;
 133     gmx_bool     bFirst;
 134     gmx_rmpbc_t  gpbc = NULL;
 135
 136     natoms = read_first_x(oenv, &status, trxfn, &t, &x, box);
 137
 138     check_index(NULL, n, index, NULL, natoms);
 139
 140     out = xvgropen(outfn, "Minimum distance to periodic image",
 141                    output_env_get_time_label(oenv), "Distance (nm)", oenv);
 142     if (output_env_get_print_xvgr_codes(oenv))
 143     {
 144         fprintf(out, "@ subtitle \"and maximum internal distance\"\n");
 145     }
 146     xvgr_legend(out, 5, leg, oenv);
 147
 148     rmint = box[XX][XX];
 149     tmint = 0;
 150
 151     if (NULL != top)
 152     {
 153         gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms, box);
 154     }
 155
 156     bFirst = TRUE;
 157     do
 158     {
 159         if (NULL != top)
 160         {
 161             gmx_rmpbc(gpbc, natoms, box, x);
 162         }
 163
 164         periodic_dist(box, x, n, index, &rmin, &rmax, ind_min);
 165         if (rmin < rmint)
 166         {
 167             rmint    = rmin;
 168             tmint    = t;
 169             ind_mini = ind_min[0];
 170             ind_minj = ind_min[1];
 171         }
 172         if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5)
 173         {
 174             fprintf(out, "&\n");
 175         }
 176         fprintf(out, "\t%g\t%6.3f %6.3f %6.3f %6.3f %6.3f\n",
 177                 output_env_conv_time(oenv, t), rmin, rmax, norm(box[0]), norm(box[1]), norm(box[2]));
 178         bFirst = FALSE;
 179     }
 180     while (read_next_x(oenv, status, &t, natoms, x, box));
 181
 182     if (NULL != top)
 183     {
 184         gmx_rmpbc_done(gpbc);
 185     }
 186
 187     ffclose(out);
 188
 189     fprintf(stdout,
 190             "\nThe shortest periodic distance is %g (nm) at time %g (%s),\n"
 191             "between atoms %d and %d\n",
 192             rmint, output_env_conv_time(oenv, tmint), output_env_get_time_unit(oenv),
 193             index[ind_mini]+1, index[ind_minj]+1);
 194 }
 195
 196 static void calc_dist(real rcut, gmx_bool bPBC, int ePBC, matrix box, rvec x[],
 197                       int nx1, int nx2, atom_id index1[], atom_id index2[],
 198                       gmx_bool bGroup,
 199                       real *rmin, real *rmax, int *nmin, int *nmax,
 200                       int *ixmin, int *jxmin, int *ixmax, int *jxmax)
 201 {
 202     int      i, j, i0 = 0, j1;
 203     int      ix, jx;
 204     atom_id *index3;
 205     rvec     dx;
 206     real     r2, rmin2, rmax2, rcut2;
 207     t_pbc    pbc;
 208     int      nmin_j, nmax_j;
 209
 210     *ixmin = -1;
 211     *jxmin = -1;
 212     *ixmax = -1;
 213     *jxmax = -1;
 214     *nmin  = 0;
 215     *nmax  = 0;
 216
 217     rcut2 = sqr(rcut);
 218
 219     /* Must init pbc every step because of pressure coupling */
 220     if (bPBC)
 221     {
 222         set_pbc(&pbc, ePBC, box);
 223     }
 224     if (index2)
 225     {
 226         i0     = 0;
 227         j1     = nx2;
 228         index3 = index2;
 229     }
 230     else
 231     {
 232         j1     = nx1;
 233         index3 = index1;
 234     }
 235
 236     rmin2 = 1e12;
 237     rmax2 = -1e12;
 238
 239     for (j = 0; (j < j1); j++)
 240     {
 241         jx = index3[j];
 242         if (index2 == NULL)
 243         {
 244             i0 = j + 1;
 245         }
 246         nmin_j = 0;
 247         nmax_j = 0;
 248         for (i = i0; (i < nx1); i++)
 249         {
 250             ix = index1[i];
 251             if (ix != jx)
 252             {
 253                 if (bPBC)
 254                 {
 255                     pbc_dx(&pbc, x[ix], x[jx], dx);
 256                 }
 257                 else
 258                 {
 259                     rvec_sub(x[ix], x[jx], dx);
 260                 }
 261                 r2 = iprod(dx, dx);
 262                 if (r2 < rmin2)
 263                 {
 264                     rmin2  = r2;
 265                     *ixmin = ix;
 266                     *jxmin = jx;
 267                 }
 268                 if (r2 > rmax2)
 269                 {
 270                     rmax2  = r2;
 271                     *ixmax = ix;
 272                     *jxmax = jx;
 273                 }
 274                 if (r2 <= rcut2)
 275                 {
 276                     nmin_j++;
 277                 }
 278                 else if (r2 > rcut2)
 279                 {
 280                     nmax_j++;
 281                 }
 282             }
 283         }
 284         if (bGroup)
 285         {
 286             if (nmin_j > 0)
 287             {
 288                 (*nmin)++;
 289             }
 290             if (nmax_j > 0)
 291             {
 292                 (*nmax)++;
 293             }
 294         }
 295         else
 296         {
 297             *nmin += nmin_j;
 298             *nmax += nmax_j;
 299         }
 300     }
 301     *rmin = sqrt(rmin2);
 302     *rmax = sqrt(rmax2);
 303 }
 304
 305 void dist_plot(const char *fn, const char *afile, const char *dfile,
 306                const char *nfile, const char *rfile, const char *xfile,
 307                real rcut, gmx_bool bMat, t_atoms *atoms,
 308                int ng, atom_id *index[], int gnx[], char *grpn[], gmx_bool bSplit,
 309                gmx_bool bMin, int nres, atom_id *residue, gmx_bool bPBC, int ePBC,
 310                gmx_bool bGroup, gmx_bool bEachResEachTime, gmx_bool bPrintResName,
 311                const output_env_t oenv)
 312 {
 313     FILE            *atm, *dist, *num;
 314     t_trxstatus     *trxout;
 315     char             buf[256];
 316     char           **leg;
 317     real             t, dmin, dmax, **mindres = NULL, **maxdres = NULL;
 318     int              nmin, nmax;
 319     t_trxstatus     *status;
 320     int              i = -1, j, k, natoms;
 321     int              min1, min2, max1, max2, min1r, min2r, max1r, max2r;
 322     atom_id          oindex[2];
 323     rvec            *x0;
 324     matrix           box;
 325     t_trxframe       frout;
 326     gmx_bool         bFirst;
 327     FILE            *respertime = NULL;
 328
 329     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 330     {
 331         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 332     }
 333
 334     sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
 335     dist = xvgropen(dfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
 336     sprintf(buf, "Number of Contacts %s %g nm", bMin ? "<" : ">", rcut);
 337     num    = nfile ? xvgropen(nfile, buf, output_env_get_time_label(oenv), "Number", oenv) : NULL;
 338     atm    = afile ? ffopen(afile, "w") : NULL;
 339     trxout = xfile ? open_trx(xfile, "w") : NULL;
 340
 341     if (bMat)
 342     {
 343         if (ng == 1)
 344         {
 345             snew(leg, 1);
 346             sprintf(buf, "Internal in %s", grpn[0]);
 347             leg[0] = strdup(buf);
 348             xvgr_legend(dist, 0, (const char**)leg, oenv);
 349             if (num)
 350             {
 351                 xvgr_legend(num, 0, (const char**)leg, oenv);
 352             }
 353         }
 354         else
 355         {
 356             snew(leg, (ng*(ng-1))/2);
 357             for (i = j = 0; (i < ng-1); i++)
 358             {
 359                 for (k = i+1; (k < ng); k++, j++)
 360                 {
 361                     sprintf(buf, "%s-%s", grpn[i], grpn[k]);
 362                     leg[j] = strdup(buf);
 363                 }
 364             }
 365             xvgr_legend(dist, j, (const char**)leg, oenv);
 366             if (num)
 367             {
 368                 xvgr_legend(num, j, (const char**)leg, oenv);
 369             }
 370         }
 371     }
 372     else
 373     {
 374         snew(leg, ng-1);
 375         for (i = 0; (i < ng-1); i++)
 376         {
 377             sprintf(buf, "%s-%s", grpn[0], grpn[i+1]);
 378             leg[i] = strdup(buf);
 379         }
 380         xvgr_legend(dist, ng-1, (const char**)leg, oenv);
 381         if (num)
 382         {
 383             xvgr_legend(num, ng-1, (const char**)leg, oenv);
 384         }
 385     }
 386
 387     if (bEachResEachTime)
 388     {
 389         sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
 390         respertime = xvgropen(rfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
 391         xvgr_legend(respertime, ng-1, (const char**)leg, oenv);
 392         if (bPrintResName)
 393         {
 394             fprintf(respertime, "# ");
 395         }
 396         for (j = 0; j < nres; j++)
 397         {
 398             fprintf(respertime, "%s%d ", *(atoms->resinfo[atoms->atom[index[0][residue[j]]].resind].name), atoms->atom[index[0][residue[j]]].resind);
 399         }
 400         fprintf(respertime, "\n");
 401
 402     }
 403
 404     j = 0;
 405     if (nres)
 406     {
 407         snew(mindres, ng-1);
 408         snew(maxdres, ng-1);
 409         for (i = 1; i < ng; i++)
 410         {
 411             snew(mindres[i-1], nres);
 412             snew(maxdres[i-1], nres);
 413             for (j = 0; j < nres; j++)
 414             {
 415                 mindres[i-1][j] = 1e6;
 416             }
 417             /* maxdres[*][*] is already 0 */
 418         }
 419     }
 420     bFirst = TRUE;
 421     do
 422     {
 423         if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5)
 424         {
 425             fprintf(dist, "&\n");
 426             if (num)
 427             {
 428                 fprintf(num, "&\n");
 429             }
 430             if (atm)
 431             {
 432                 fprintf(atm, "&\n");
 433             }
 434         }
 435         fprintf(dist, "%12e", output_env_conv_time(oenv, t));
 436         if (num)
 437         {
 438             fprintf(num, "%12e", output_env_conv_time(oenv, t));
 439         }
 440
 441         if (bMat)
 442         {
 443             if (ng == 1)
 444             {
 445                 calc_dist(rcut, bPBC, ePBC, box, x0, gnx[0], gnx[0], index[0], index[0], bGroup,
 446                           &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
 447                 fprintf(dist, "  %12e", bMin ? dmin : dmax);
 448                 if (num)
 449                 {
 450                     fprintf(num, "  %8d", bMin ? nmin : nmax);
 451                 }
 452             }
 453             else
 454             {
 455                 for (i = 0; (i < ng-1); i++)
 456                 {
 457                     for (k = i+1; (k < ng); k++)
 458                     {
 459                         calc_dist(rcut, bPBC, ePBC, box, x0, gnx[i], gnx[k], index[i], index[k],
 460                                   bGroup, &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
 461                         fprintf(dist, "  %12e", bMin ? dmin : dmax);
 462                         if (num)
 463                         {
 464                             fprintf(num, "  %8d", bMin ? nmin : nmax);
 465                         }
 466                     }
 467                 }
 468             }
 469         }
 470         else
 471         {
 472             for (i = 1; (i < ng); i++)
 473             {
 474                 calc_dist(rcut, bPBC, ePBC, box, x0, gnx[0], gnx[i], index[0], index[i], bGroup,
 475                           &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
 476                 fprintf(dist, "  %12e", bMin ? dmin : dmax);
 477                 if (num)
 478                 {
 479                     fprintf(num, "  %8d", bMin ? nmin : nmax);
 480                 }
 481                 if (nres)
 482                 {
 483                     for (j = 0; j < nres; j++)
 484                     {
 485                         calc_dist(rcut, bPBC, ePBC, box, x0, residue[j+1]-residue[j], gnx[i],
 486                                   &(index[0][residue[j]]), index[i], bGroup,
 487                                   &dmin, &dmax, &nmin, &nmax, &min1r, &min2r, &max1r, &max2r);
 488                         mindres[i-1][j] = min(mindres[i-1][j], dmin);
 489                         maxdres[i-1][j] = max(maxdres[i-1][j], dmax);
 490                     }
 491                 }
 492             }
 493         }
 494         fprintf(dist, "\n");
 495         if (num)
 496         {
 497             fprintf(num, "\n");
 498         }
 499         if ( (bMin ? min1 : max1) != -1)
 500         {
 501             if (atm)
 502             {
 503                 fprintf(atm, "%12e  %12d  %12d\n",
 504                         output_env_conv_time(oenv, t), 1+(bMin ? min1 : max1),
 505                         1+(bMin ? min2 : max2));
 506             }
 507         }
 508
 509         if (trxout)
 510         {
 511             oindex[0] = bMin ? min1 : max1;
 512             oindex[1] = bMin ? min2 : max2;
 513             write_trx(trxout, 2, oindex, atoms, i, t, box, x0, NULL, NULL);
 514         }
 515         bFirst = FALSE;
 516         /*dmin should be minimum distance for residue and group*/
 517         if (bEachResEachTime)
 518         {
 519             fprintf(respertime, "%12e", t);
 520             for (i = 1; i < ng; i++)
 521             {
 522                 for (j = 0; j < nres; j++)
 523                 {
 524                     fprintf(respertime, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
 525                     /*reset distances for next time point*/
 526                     mindres[i-1][j] = 1e6;
 527                     maxdres[i-1][j] = 0;
 528                 }
 529             }
 530             fprintf(respertime, "\n");
 531         }
 532     }
 533     while (read_next_x(oenv, status, &t, natoms, x0, box));
 534
 535     close_trj(status);
 536     ffclose(dist);
 537     if (num)
 538     {
 539         ffclose(num);
 540     }
 541     if (atm)
 542     {
 543         ffclose(atm);
 544     }
 545     if (trxout)
 546     {
 547         close_trx(trxout);
 548     }
 549
 550     if (nres && !bEachResEachTime)
 551     {
 552         FILE *res;
 553
 554         sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
 555         res = xvgropen(rfile, buf, "Residue (#)", "Distance (nm)", oenv);
 556         xvgr_legend(res, ng-1, (const char**)leg, oenv);
 557         for (j = 0; j < nres; j++)
 558         {
 559             fprintf(res, "%4d", j+1);
 560             for (i = 1; i < ng; i++)
 561             {
 562                 fprintf(res, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
 563             }
 564             fprintf(res, "\n");
 565         }
 566     }
 567
 568     sfree(x0);
 569 }
 570
 571 int find_residues(t_atoms *atoms, int n, atom_id index[], atom_id **resindex)
 572 {
 573     int  i;
 574     int  nres = 0, resnr, presnr;
 575     int *residx;
 576
 577     /* build index of first atom numbers for each residue */
 578     presnr = NOTSET;
 579     snew(residx, atoms->nres+1);
 580     for (i = 0; i < n; i++)
 581     {
 582         resnr = atoms->atom[index[i]].resind;
 583         if (resnr != presnr)
 584         {
 585             residx[nres] = i;
 586             nres++;
 587             presnr = resnr;
 588         }
 589     }
 590     if (debug)
 591     {
 592         printf("Found %d residues out of %d (%d/%d atoms)\n",
 593                nres, atoms->nres, atoms->nr, n);
 594     }
 595     srenew(residx, nres+1);
 596     /* mark end of last residue */
 597     residx[nres] = n;
 598     *resindex    = residx;
 599     return nres;
 600 }
 601
 602 void dump_res(FILE *out, int nres, atom_id *resindex, int n, atom_id index[])
 603 {
 604     int i, j;
 605
 606     for (i = 0; i < nres-1; i++)
 607     {
 608         fprintf(out, "Res %d (%d):", i, resindex[i+1]-resindex[i]);
 609         for (j = resindex[i]; j < resindex[i+1]; j++)
 610         {
 611             fprintf(out, " %d(%d)", j, index[j]);
 612         }
 613         fprintf(out, "\n");
 614     }
 615 }
 616
 617 int gmx_mindist(int argc, char *argv[])
 618 {
 619     const char     *desc[] = {
 620         "[TT]g_mindist[tt] computes the distance between one group and a number of",
 621         "other groups. Both the minimum distance",
 622         "(between any pair of atoms from the respective groups)",
 623         "and the number of contacts within a given",
 624         "distance are written to two separate output files.",
 625         "With the [TT]-group[tt] option a contact of an atom in another group",
 626         "with multiple atoms in the first group is counted as one contact",
 627         "instead of as multiple contacts.",
 628         "With [TT]-or[tt], minimum distances to each residue in the first",
 629         "group are determined and plotted as a function of residue number.[PAR]",
 630         "With option [TT]-pi[tt] the minimum distance of a group to its",
 631         "periodic image is plotted. This is useful for checking if a protein",
 632         "has seen its periodic image during a simulation. Only one shift in",
 633         "each direction is considered, giving a total of 26 shifts.",
 634         "It also plots the maximum distance within the group and the lengths",
 635         "of the three box vectors.[PAR]",
 636         "Other programs that calculate distances are [TT]g_dist[tt]",
 637         "and [TT]g_bond[tt]."
 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     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     trxfnm  = ftp2fn(efTRX, NFILE, fnm);
 702     ndxfnm  = ftp2fn_null(efNDX, NFILE, fnm);
 703     distfnm = opt2fn("-od", NFILE, fnm);
 704     numfnm  = opt2fn_null("-on", NFILE, fnm);
 705     atmfnm  = ftp2fn_null(efOUT, NFILE, fnm);
 706     oxfnm   = opt2fn_null("-ox", NFILE, fnm);
 707     resfnm  = opt2fn_null("-or", NFILE, fnm);
 708     if (bPI || resfnm != NULL)
 709     {
 710         /* We need a tps file */
 711         tpsfnm = ftp2fn(efTPS, NFILE, fnm);
 712     }
 713     else
 714     {
 715         tpsfnm = ftp2fn_null(efTPS, NFILE, fnm);
 716     }
 717
 718     if (!tpsfnm && !ndxfnm)
 719     {
 720         gmx_fatal(FARGS, "You have to specify either the index file or a tpr file");
 721     }
 722
 723     if (bPI)
 724     {
 725         ng = 1;
 726         fprintf(stderr, "Choose a group for distance calculation\n");
 727     }
 728     else if (!bMat)
 729     {
 730         ng++;
 731     }
 732
 733     snew(gnx, ng);
 734     snew(index, ng);
 735     snew(grpname, ng);
 736
 737     if (tpsfnm || resfnm || !ndxfnm)
 738     {
 739         snew(top, 1);
 740         bTop = read_tps_conf(tpsfnm, title, top, &ePBC, &x, NULL, box, FALSE);
 741         if (bPI && !bTop)
 742         {
 743             printf("\nWARNING: Without a run input file a trajectory with broken molecules will not give the correct periodic image distance\n\n");
 744         }
 745     }
 746     get_index(top ? &(top->atoms) : NULL, ndxfnm, ng, gnx, index, grpname);
 747
 748     if (bMat && (ng == 1))
 749     {
 750         ng = gnx[0];
 751         printf("Special case: making distance matrix between all atoms in group %s\n",
 752                grpname[0]);
 753         srenew(gnx, ng);
 754         srenew(index, ng);
 755         srenew(grpname, ng);
 756         for (i = 1; (i < ng); i++)
 757         {
 758             gnx[i]      = 1;
 759             grpname[i]  = grpname[0];
 760             snew(index[i], 1);
 761             index[i][0] = index[0][i];
 762         }
 763         gnx[0] = 1;
 764     }
 765
 766     if (resfnm)
 767     {
 768         nres = find_residues(top ? &(top->atoms) : NULL,
 769                              gnx[0], index[0], &residues);
 770         if (debug)
 771         {
 772             dump_res(debug, nres, residues, gnx[0], index[0]);
 773         }
 774     }
 775
 776     if (bPI)
 777     {
 778         periodic_mindist_plot(trxfnm, distfnm, top, ePBC, gnx[0], index[0], bSplit, oenv);
 779     }
 780     else
 781     {
 782         dist_plot(trxfnm, atmfnm, distfnm, numfnm, resfnm, oxfnm,
 783                   rcutoff, bMat, top ? &(top->atoms) : NULL,
 784                   ng, index, gnx, grpname, bSplit, !bMax, nres, residues, bPBC, ePBC,
 785                   bGroup, bEachResEachTime, bPrintResName, oenv);
 786     }
 787
 788     do_view(oenv, distfnm, "-nxy");
 789     if (!bPI)
 790     {
 791         do_view(oenv, numfnm, "-nxy");
 792     }
 793
 794     thanx(stderr);
 795
 796     return 0;
 797 }