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