src/gromacs/gmxana/gmx_mindist.cpp

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