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