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