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