src/gromacs/gmxana/gmx_traj.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 #include <string>
  46
  47 #include "gromacs/commandline/pargs.h"
  48 #include "gromacs/commandline/viewit.h"
  49 #include "gromacs/fileio/confio.h"
  50 #include "gromacs/fileio/trxio.h"
  51 #include "gromacs/fileio/xvgr.h"
  52 #include "gromacs/gmxana/gmx_ana.h"
  53 #include "gromacs/linearalgebra/nrjac.h"
  54 #include "gromacs/math/functions.h"
  55 #include "gromacs/math/units.h"
  56 #include "gromacs/math/vec.h"
  57 #include "gromacs/mdtypes/md_enums.h"
  58 #include "gromacs/pbcutil/pbc.h"
  59 #include "gromacs/pbcutil/rmpbc.h"
  60 #include "gromacs/topology/index.h"
  61 #include "gromacs/topology/topology.h"
  62 #include "gromacs/trajectory/trajectoryframe.h"
  63 #include "gromacs/utility/arraysize.h"
  64 #include "gromacs/utility/cstringutil.h"
  65 #include "gromacs/utility/fatalerror.h"
  66 #include "gromacs/utility/futil.h"
  67 #include "gromacs/utility/smalloc.h"
  68 #include "gromacs/utility/stringutil.h"
  69
  70 static void
  71 low_print_data(FILE* fp, real time, rvec x[], int n, const int* index, const gmx_bool bDim[], const char* sffmt)
  72 {
  73     int i, ii, d;
  74
  75     fprintf(fp, " %g", time);
  76     for (i = 0; i < n; i++)
  77     {
  78         if (index != nullptr)
  79         {
  80             ii = index[i];
  81         }
  82         else
  83         {
  84             ii = i;
  85         }
  86         for (d = 0; d < DIM; d++)
  87         {
  88             if (bDim[d])
  89             {
  90                 fprintf(fp, sffmt, x[ii][d]);
  91             }
  92         }
  93         if (bDim[DIM])
  94         {
  95             fprintf(fp, sffmt, norm(x[ii]));
  96         }
  97     }
  98     fprintf(fp, "\n");
  99 }
 100
 101 static void average_data(rvec x[], rvec xav[], const real* mass, int ngrps, const int isize[], int** index)
 102 {
 103     int    g, i, ind, d;
 104     real   m;
 105     rvec   tmp;
 106     double sum[DIM], mtot;
 107
 108     for (g = 0; g < ngrps; g++)
 109     {
 110         clear_dvec(sum);
 111         clear_rvec(xav[g]);
 112         mtot = 0;
 113         for (i = 0; i < isize[g]; i++)
 114         {
 115             ind = index[g][i];
 116             if (mass != nullptr)
 117             {
 118                 m = mass[ind];
 119                 svmul(m, x[ind], tmp);
 120                 for (d = 0; d < DIM; d++)
 121                 {
 122                     sum[d] += tmp[d];
 123                 }
 124                 mtot += m;
 125             }
 126             else
 127             {
 128                 for (d = 0; d < DIM; d++)
 129                 {
 130                     sum[d] += x[ind][d];
 131                 }
 132             }
 133         }
 134         if (mass != nullptr)
 135         {
 136             for (d = 0; d < DIM; d++)
 137             {
 138                 xav[g][d] = sum[d] / mtot;
 139             }
 140         }
 141         else
 142         {
 143             /* mass=NULL, so these are forces: sum only (do not average) */
 144             for (d = 0; d < DIM; d++)
 145             {
 146                 xav[g][d] = sum[d];
 147             }
 148         }
 149     }
 150 }
 151
 152 static void print_data(FILE*       fp,
 153                        real        time,
 154                        rvec        x[],
 155                        real*       mass,
 156                        gmx_bool    bCom,
 157                        int         ngrps,
 158                        int         isize[],
 159                        int**       index,
 160                        gmx_bool    bDim[],
 161                        const char* sffmt)
 162 {
 163     static rvec* xav = nullptr;
 164
 165     if (bCom)
 166     {
 167         if (xav == nullptr)
 168         {
 169             snew(xav, ngrps);
 170         }
 171         average_data(x, xav, mass, ngrps, isize, index);
 172         low_print_data(fp, time, xav, ngrps, nullptr, bDim, sffmt);
 173     }
 174     else
 175     {
 176         low_print_data(fp, time, x, isize[0], index[0], bDim, sffmt);
 177     }
 178 }
 179
 180 static void write_trx_x(t_trxstatus*      status,
 181                         const t_trxframe* fr,
 182                         real*             mass,
 183                         gmx_bool          bCom,
 184                         int               ngrps,
 185                         int               isize[],
 186                         int**             index)
 187 {
 188     static rvec*    xav   = nullptr;
 189     static t_atoms* atoms = nullptr;
 190     t_trxframe      fr_av;
 191     int             i;
 192
 193     if (bCom)
 194     {
 195         if (xav == nullptr)
 196         {
 197             snew(xav, ngrps);
 198             snew(atoms, 1);
 199             *atoms = *fr->atoms;
 200             snew(atoms->atom, ngrps);
 201             atoms->nr = ngrps;
 202             /* Note that atom and residue names will be the ones
 203              * of the first atom in each group.
 204              */
 205             for (i = 0; i < ngrps; i++)
 206             {
 207                 atoms->atom[i]     = fr->atoms->atom[index[i][0]];
 208                 atoms->atomname[i] = fr->atoms->atomname[index[i][0]];
 209             }
 210         }
 211         average_data(fr->x, xav, mass, ngrps, isize, index);
 212         fr_av        = *fr;
 213         fr_av.natoms = ngrps;
 214         fr_av.atoms  = atoms;
 215         fr_av.x      = xav;
 216         write_trxframe(status, &fr_av, nullptr);
 217     }
 218     else
 219     {
 220         write_trxframe_indexed(status, fr, isize[0], index[0], nullptr);
 221     }
 222 }
 223
 224 static void make_legend(FILE*                   fp,
 225                         int                     ngrps,
 226                         int                     isize,
 227                         int                     index[],
 228                         char**                  name,
 229                         gmx_bool                bCom,
 230                         gmx_bool                bMol,
 231                         const gmx_bool          bDim[],
 232                         const gmx_output_env_t* oenv)
 233 {
 234     char**      leg;
 235     const char* dimtxt[] = { " X", " Y", " Z", "" };
 236     int         n, i, j, d;
 237
 238     if (bCom)
 239     {
 240         n = ngrps;
 241     }
 242     else
 243     {
 244         n = isize;
 245     }
 246
 247     snew(leg, 4 * n);
 248     j = 0;
 249     for (i = 0; i < n; i++)
 250     {
 251         for (d = 0; d <= DIM; d++)
 252         {
 253             if (bDim[d])
 254             {
 255                 snew(leg[j], STRLEN);
 256                 if (bMol)
 257                 {
 258                     sprintf(leg[j], "mol %d%s", index[i] + 1, dimtxt[d]);
 259                 }
 260                 else if (bCom)
 261                 {
 262                     sprintf(leg[j], "%s%s", name[i], dimtxt[d]);
 263                 }
 264                 else
 265                 {
 266                     sprintf(leg[j], "atom %d%s", index[i] + 1, dimtxt[d]);
 267                 }
 268                 j++;
 269             }
 270         }
 271     }
 272     xvgr_legend(fp, j, leg, oenv);
 273
 274     for (i = 0; i < j; i++)
 275     {
 276         sfree(leg[i]);
 277     }
 278     sfree(leg);
 279 }
 280
 281 static real ekrot(rvec x[], rvec v[], const real mass[], int isize, const int index[])
 282 {
 283     real   TCM[5][5], L[5][5];
 284     double tm, m0, lxx, lxy, lxz, lyy, lyz, lzz, ekrot;
 285     rvec   a0, ocm;
 286     dvec   dx, b0;
 287     dvec   xcm, vcm, acm;
 288     int    i, j, m, n;
 289
 290     clear_dvec(xcm);
 291     clear_dvec(vcm);
 292     clear_dvec(acm);
 293     tm = 0.0;
 294     for (i = 0; i < isize; i++)
 295     {
 296         j  = index[i];
 297         m0 = mass[j];
 298         tm += m0;
 299         cprod(x[j], v[j], a0);
 300         for (m = 0; (m < DIM); m++)
 301         {
 302             xcm[m] += m0 * x[j][m]; /* c.o.m. position */
 303             vcm[m] += m0 * v[j][m]; /* c.o.m. velocity */
 304             acm[m] += m0 * a0[m];   /* rotational velocity around c.o.m. */
 305         }
 306     }
 307     dcprod(xcm, vcm, b0);
 308     for (m = 0; (m < DIM); m++)
 309     {
 310         xcm[m] /= tm;
 311         vcm[m] /= tm;
 312         acm[m] -= b0[m] / tm;
 313     }
 314
 315     lxx = lxy = lxz = lyy = lyz = lzz = 0.0;
 316     for (i = 0; i < isize; i++)
 317     {
 318         j  = index[i];
 319         m0 = mass[j];
 320         for (m = 0; m < DIM; m++)
 321         {
 322             dx[m] = x[j][m] - xcm[m];
 323         }
 324         lxx += dx[XX] * dx[XX] * m0;
 325         lxy += dx[XX] * dx[YY] * m0;
 326         lxz += dx[XX] * dx[ZZ] * m0;
 327         lyy += dx[YY] * dx[YY] * m0;
 328         lyz += dx[YY] * dx[ZZ] * m0;
 329         lzz += dx[ZZ] * dx[ZZ] * m0;
 330     }
 331
 332     L[XX][XX] = lyy + lzz;
 333     L[YY][XX] = -lxy;
 334     L[ZZ][XX] = -lxz;
 335     L[XX][YY] = -lxy;
 336     L[YY][YY] = lxx + lzz;
 337     L[ZZ][YY] = -lyz;
 338     L[XX][ZZ] = -lxz;
 339     L[YY][ZZ] = -lyz;
 340     L[ZZ][ZZ] = lxx + lyy;
 341
 342     m_inv_gen(&L[0][0], DIM, &TCM[0][0]);
 343
 344     /* Compute omega (hoeksnelheid) */
 345     clear_rvec(ocm);
 346     ekrot = 0;
 347     for (m = 0; m < DIM; m++)
 348     {
 349         for (n = 0; n < DIM; n++)
 350         {
 351             ocm[m] += TCM[m][n] * acm[n];
 352         }
 353         ekrot += 0.5 * ocm[m] * acm[m];
 354     }
 355
 356     return ekrot;
 357 }
 358
 359 static real ektrans(rvec v[], const real mass[], int isize, const int index[])
 360 {
 361     dvec   mvcom;
 362     double mtot;
 363     int    i, j, d;
 364
 365     clear_dvec(mvcom);
 366     mtot = 0;
 367     for (i = 0; i < isize; i++)
 368     {
 369         j = index[i];
 370         for (d = 0; d < DIM; d++)
 371         {
 372             mvcom[d] += mass[j] * v[j][d];
 373         }
 374         mtot += mass[j];
 375     }
 376
 377     return dnorm2(mvcom) / (mtot * 2);
 378 }
 379
 380 static real temp(rvec v[], const real mass[], int isize, const int index[])
 381 {
 382     double ekin2;
 383     int    i, j;
 384
 385     ekin2 = 0;
 386     for (i = 0; i < isize; i++)
 387     {
 388         j = index[i];
 389         ekin2 += mass[j] * norm2(v[j]);
 390     }
 391
 392     return ekin2 / (3 * isize * BOLTZ);
 393 }
 394
 395 static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[])
 396 {
 397     rvec hbox;
 398     int  d, i, m;
 399
 400     for (d = 0; d < DIM; d++)
 401     {
 402         hbox[d] = 0.5 * box[d][d];
 403     }
 404     for (i = 0; i < natoms; i++)
 405     {
 406         for (m = DIM - 1; m >= 0; m--)
 407         {
 408             while (x[i][m] - xp[i][m] <= -hbox[m])
 409             {
 410                 for (d = 0; d <= m; d++)
 411                 {
 412                     x[i][d] += box[m][d];
 413                 }
 414             }
 415             while (x[i][m] - xp[i][m] > hbox[m])
 416             {
 417                 for (d = 0; d <= m; d++)
 418                 {
 419                     x[i][d] -= box[m][d];
 420                 }
 421             }
 422         }
 423     }
 424 }
 425
 426 static void write_pdb_bfac(const char*             fname,
 427                            const char*             xname,
 428                            const char*             title,
 429                            t_atoms*                atoms,
 430                            PbcType                 pbcType,
 431                            matrix                  box,
 432                            int                     isize,
 433                            int*                    index,
 434                            int                     nfr_x,
 435                            rvec*                   x,
 436                            int                     nfr_v,
 437                            rvec*                   sum,
 438                            const gmx_bool          bDim[],
 439                            real                    scale_factor,
 440                            const gmx_output_env_t* oenv)
 441 {
 442     FILE* fp;
 443     real  max, len2, scale;
 444     int   maxi;
 445     int   i, m, onedim;
 446
 447     if ((nfr_x == 0) || (nfr_v == 0))
 448     {
 449         fprintf(stderr, "No frames found for %s, will not write %s\n", title, fname);
 450     }
 451     else
 452     {
 453         fprintf(stderr, "Used %d frames for %s\n", nfr_x, "coordinates");
 454         fprintf(stderr, "Used %d frames for %s\n", nfr_v, title);
 455         onedim = -1;
 456         if (!bDim[DIM])
 457         {
 458             m = 0;
 459             for (i = 0; i < DIM; i++)
 460             {
 461                 if (bDim[i])
 462                 {
 463                     onedim = i;
 464                     m++;
 465                 }
 466             }
 467             if (m != 1)
 468             {
 469                 onedim = -1;
 470             }
 471         }
 472         scale = 1.0 / nfr_v;
 473         for (i = 0; i < isize; i++)
 474         {
 475             svmul(scale, sum[index[i]], sum[index[i]]);
 476         }
 477
 478         fp = xvgropen(xname, title, "Atom", "Spatial component", oenv);
 479         for (i = 0; i < isize; i++)
 480         {
 481             fprintf(fp,
 482                     "%-5d  %10.3f  %10.3f  %10.3f\n",
 483                     1 + i,
 484                     sum[index[i]][XX],
 485                     sum[index[i]][YY],
 486                     sum[index[i]][ZZ]);
 487         }
 488         xvgrclose(fp);
 489         max  = 0;
 490         maxi = 0;
 491         for (i = 0; i < isize; i++)
 492         {
 493             len2 = 0;
 494             for (m = 0; m < DIM; m++)
 495             {
 496                 if (bDim[m] || bDim[DIM])
 497                 {
 498                     len2 += gmx::square(sum[index[i]][m]);
 499                 }
 500             }
 501             if (len2 > max)
 502             {
 503                 max  = len2;
 504                 maxi = index[i];
 505             }
 506         }
 507         if (scale_factor != 0)
 508         {
 509             scale = scale_factor;
 510         }
 511         else
 512         {
 513             if (max == 0)
 514             {
 515                 scale = 1;
 516             }
 517             else
 518             {
 519                 scale = 10.0 / std::sqrt(max);
 520             }
 521         }
 522
 523         printf("Maximum %s is %g on atom %d %s, res. %s %d\n",
 524                title,
 525                std::sqrt(max),
 526                maxi + 1,
 527                *(atoms->atomname[maxi]),
 528                *(atoms->resinfo[atoms->atom[maxi].resind].name),
 529                atoms->resinfo[atoms->atom[maxi].resind].nr);
 530
 531         if (atoms->pdbinfo == nullptr)
 532         {
 533             snew(atoms->pdbinfo, atoms->nr);
 534         }
 535         atoms->havePdbInfo = TRUE;
 536
 537         if (onedim == -1)
 538         {
 539             for (i = 0; i < isize; i++)
 540             {
 541                 len2 = 0;
 542                 for (m = 0; m < DIM; m++)
 543                 {
 544                     if (bDim[m] || bDim[DIM])
 545                     {
 546                         len2 += gmx::square(sum[index[i]][m]);
 547                     }
 548                 }
 549                 atoms->pdbinfo[index[i]].bfac = std::sqrt(len2) * scale;
 550             }
 551         }
 552         else
 553         {
 554             for (i = 0; i < isize; i++)
 555             {
 556                 atoms->pdbinfo[index[i]].bfac = sum[index[i]][onedim] * scale;
 557             }
 558         }
 559         write_sto_conf_indexed(fname, title, atoms, x, nullptr, pbcType, box, isize, index);
 560     }
 561 }
 562
 563 static void update_histo(int gnx, const int index[], rvec v[], int* nhisto, int** histo, real binwidth)
 564 {
 565     int  i, m, in, nnn;
 566     real vn, vnmax;
 567
 568     if (*histo == nullptr)
 569     {
 570         vnmax = 0;
 571         for (i = 0; (i < gnx); i++)
 572         {
 573             vn    = norm(v[index[i]]);
 574             vnmax = std::max(vn, vnmax);
 575         }
 576         vnmax *= 2;
 577         *nhisto = static_cast<int>(1 + (vnmax / binwidth));
 578         snew(*histo, *nhisto);
 579     }
 580     for (i = 0; (i < gnx); i++)
 581     {
 582         vn = norm(v[index[i]]);
 583         in = static_cast<int>(vn / binwidth);
 584         if (in >= *nhisto)
 585         {
 586             nnn = in + 100;
 587             fprintf(stderr, "Extending histogram from %d to %d\n", *nhisto, nnn);
 588
 589             srenew(*histo, nnn);
 590             for (m = *nhisto; (m < nnn); m++)
 591             {
 592                 (*histo)[m] = 0;
 593             }
 594             *nhisto = nnn;
 595         }
 596         (*histo)[in]++;
 597     }
 598 }
 599
 600 static void print_histo(const char* fn, int nhisto, int histo[], real binwidth, const gmx_output_env_t* oenv)
 601 {
 602     FILE* fp;
 603     int   i;
 604
 605     fp = xvgropen(fn, "Velocity distribution", "V (nm/ps)", "arbitrary units", oenv);
 606     for (i = 0; (i < nhisto); i++)
 607     {
 608         fprintf(fp, "%10.3e  %10d\n", i * binwidth, histo[i]);
 609     }
 610     xvgrclose(fp);
 611 }
 612
 613 int gmx_traj(int argc, char* argv[])
 614 {
 615     const char* desc[] = {
 616         "[THISMODULE] plots coordinates, velocities, forces and/or the box.",
 617         "With [TT]-com[tt] the coordinates, velocities and forces are",
 618         "calculated for the center of mass of each group.",
 619         "When [TT]-mol[tt] is set, the numbers in the index file are",
 620         "interpreted as molecule numbers and the same procedure as with",
 621         "[TT]-com[tt] is used for each molecule.[PAR]",
 622         "Option [TT]-ot[tt] plots the temperature of each group,",
 623         "provided velocities are present in the trajectory file.",
 624         "No corrections are made for constrained degrees of freedom!",
 625         "This implies [TT]-com[tt].[PAR]",
 626         "Options [TT]-ekt[tt] and [TT]-ekr[tt] plot the translational and",
 627         "rotational kinetic energy of each group,",
 628         "provided velocities are present in the trajectory file.",
 629         "This implies [TT]-com[tt].[PAR]",
 630         "Options [TT]-cv[tt] and [TT]-cf[tt] write the average velocities",
 631         "and average forces as temperature factors to a [REF].pdb[ref] file with",
 632         "the average coordinates or the coordinates at [TT]-ctime[tt].",
 633         "The temperature factors are scaled such that the maximum is 10.",
 634         "The scaling can be changed with the option [TT]-scale[tt].",
 635         "To get the velocities or forces of one",
 636         "frame set both [TT]-b[tt] and [TT]-e[tt] to the time of",
 637         "desired frame. When averaging over frames you might need to use",
 638         "the [TT]-nojump[tt] option to obtain the correct average coordinates.",
 639         "If you select either of these option the average force and velocity",
 640         "for each atom are written to an [REF].xvg[ref] file as well",
 641         "(specified with [TT]-av[tt] or [TT]-af[tt]).[PAR]",
 642         "Option [TT]-vd[tt] computes a velocity distribution, i.e. the",
 643         "norm of the vector is plotted. In addition in the same graph",
 644         "the kinetic energy distribution is given.",
 645         "",
 646         "See [gmx-trajectory] for plotting similar data for selections."
 647     };
 648     static gmx_bool bMol = FALSE, bCom = FALSE, bPBC = TRUE, bNoJump = FALSE;
 649     static gmx_bool bX = TRUE, bY = TRUE, bZ = TRUE, bNorm = FALSE, bFP = FALSE;
 650     static int      ngroups = 1;
 651     static real     ctime = -1, scale = 0, binwidth = 1;
 652     t_pargs         pa[] = {
 653         { "-com", FALSE, etBOOL, { &bCom }, "Plot data for the com of each group" },
 654         { "-pbc", FALSE, etBOOL, { &bPBC }, "Make molecules whole for COM" },
 655         { "-mol",
 656           FALSE,
 657           etBOOL,
 658           { &bMol },
 659           "Index contains molecule numbers instead of atom numbers" },
 660         { "-nojump", FALSE, etBOOL, { &bNoJump }, "Remove jumps of atoms across the box" },
 661         { "-x", FALSE, etBOOL, { &bX }, "Plot X-component" },
 662         { "-y", FALSE, etBOOL, { &bY }, "Plot Y-component" },
 663         { "-z", FALSE, etBOOL, { &bZ }, "Plot Z-component" },
 664         { "-ng", FALSE, etINT, { &ngroups }, "Number of groups to consider" },
 665         { "-len", FALSE, etBOOL, { &bNorm }, "Plot vector length" },
 666         { "-fp", FALSE, etBOOL, { &bFP }, "Full precision output" },
 667         { "-bin", FALSE, etREAL, { &binwidth }, "Binwidth for velocity histogram (nm/ps)" },
 668         { "-ctime",
 669           FALSE,
 670           etREAL,
 671           { &ctime },
 672           "Use frame at this time for x in [TT]-cv[tt] and [TT]-cf[tt] instead of the average x" },
 673         { "-scale",
 674           FALSE,
 675           etREAL,
 676           { &scale },
 677           "Scale factor for [REF].pdb[ref] output, 0 is autoscale" }
 678     };
 679     FILE *       outx = nullptr, *outv = nullptr, *outf = nullptr, *outb = nullptr, *outt = nullptr;
 680     FILE *       outekt = nullptr, *outekr = nullptr;
 681     t_topology   top;
 682     PbcType      pbcType;
 683     real *       mass, time;
 684     const char*  indexfn;
 685     t_trxframe   fr;
 686     int          flags, nvhisto = 0, *vhisto = nullptr;
 687     rvec *       xtop, *xp = nullptr;
 688     rvec *       sumx = nullptr, *sumv = nullptr, *sumf = nullptr;
 689     matrix       topbox;
 690     t_trxstatus* status;
 691     t_trxstatus* status_out = nullptr;
 692     gmx_rmpbc_t  gpbc       = nullptr;
 693     int          i, j;
 694     int          nr_xfr, nr_vfr, nr_ffr;
 695     char**       grpname;
 696     int *        isize0, *isize;
 697     int **       index0, **index;
 698     int*         atndx;
 699     t_block*     mols;
 700     gmx_bool     bTop, bOX, bOXT, bOV, bOF, bOB, bOT, bEKT, bEKR, bCV, bCF;
 701     gmx_bool     bDim[4], bDum[4], bVD;
 702     char         sffmt[STRLEN];
 703     const char*  box_leg[6] = { "XX", "YY", "ZZ", "YX", "ZX", "ZY" };
 704     gmx_output_env_t* oenv;
 705
 706     t_filenm fnm[] = {
 707         { efTRX, "-f", nullptr, ffREAD },       { efTPS, nullptr, nullptr, ffREAD },
 708         { efNDX, nullptr, nullptr, ffOPTRD },   { efXVG, "-ox", "coord", ffOPTWR },
 709         { efTRX, "-oxt", "coord", ffOPTWR },    { efXVG, "-ov", "veloc", ffOPTWR },
 710         { efXVG, "-of", "force", ffOPTWR },     { efXVG, "-ob", "box", ffOPTWR },
 711         { efXVG, "-ot", "temp", ffOPTWR },      { efXVG, "-ekt", "ektrans", ffOPTWR },
 712         { efXVG, "-ekr", "ekrot", ffOPTWR },    { efXVG, "-vd", "veldist", ffOPTWR },
 713         { efPDB, "-cv", "veloc", ffOPTWR },     { efPDB, "-cf", "force", ffOPTWR },
 714         { efXVG, "-av", "all_veloc", ffOPTWR }, { efXVG, "-af", "all_force", ffOPTWR }
 715     };
 716 #define NFILE asize(fnm)
 717
 718     if (!parse_common_args(&argc,
 719                            argv,
 720                            PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW,
 721                            NFILE,
 722                            fnm,
 723                            asize(pa),
 724                            pa,
 725                            asize(desc),
 726                            desc,
 727                            0,
 728                            nullptr,
 729                            &oenv))
 730     {
 731         return 0;
 732     }
 733
 734     if (bMol)
 735     {
 736         fprintf(stderr,
 737                 "Interpreting indexfile entries as molecules.\n"
 738                 "Using center of mass.\n");
 739     }
 740
 741     bOX  = opt2bSet("-ox", NFILE, fnm);
 742     bOXT = opt2bSet("-oxt", NFILE, fnm);
 743     bOV  = opt2bSet("-ov", NFILE, fnm);
 744     bOF  = opt2bSet("-of", NFILE, fnm);
 745     bOB  = opt2bSet("-ob", NFILE, fnm);
 746     bOT  = opt2bSet("-ot", NFILE, fnm);
 747     bEKT = opt2bSet("-ekt", NFILE, fnm);
 748     bEKR = opt2bSet("-ekr", NFILE, fnm);
 749     bCV  = opt2bSet("-cv", NFILE, fnm) || opt2bSet("-av", NFILE, fnm);
 750     bCF  = opt2bSet("-cf", NFILE, fnm) || opt2bSet("-af", NFILE, fnm);
 751     bVD  = opt2bSet("-vd", NFILE, fnm) || opt2parg_bSet("-bin", asize(pa), pa);
 752     if (bMol || bOT || bEKT || bEKR)
 753     {
 754         bCom = TRUE;
 755     }
 756
 757     bDim[XX]  = bX;
 758     bDim[YY]  = bY;
 759     bDim[ZZ]  = bZ;
 760     bDim[DIM] = bNorm;
 761
 762     if (bFP)
 763     {
 764         sprintf(sffmt, "\t%s", gmx_real_fullprecision_pfmt);
 765     }
 766     else
 767     {
 768         sprintf(sffmt, "\t%%g");
 769     }
 770     std::string sffmt6 = gmx::formatString("%s%s%s%s%s%s", sffmt, sffmt, sffmt, sffmt, sffmt, sffmt);
 771
 772     bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
 773                          &top,
 774                          &pbcType,
 775                          &xtop,
 776                          nullptr,
 777                          topbox,
 778                          bCom && (bOX || bOXT || bOV || bOT || bEKT || bEKR));
 779     sfree(xtop);
 780     if ((bMol || bCV || bCF) && !bTop)
 781     {
 782         gmx_fatal(FARGS, "Need a run input file for option -mol, -cv or -cf");
 783     }
 784
 785     if (bMol)
 786     {
 787         indexfn = ftp2fn(efNDX, NFILE, fnm);
 788     }
 789     else
 790     {
 791         indexfn = ftp2fn_null(efNDX, NFILE, fnm);
 792     }
 793
 794     if (!(bCom && !bMol))
 795     {
 796         ngroups = 1;
 797     }
 798     snew(grpname, ngroups);
 799     snew(isize0, ngroups);
 800     snew(index0, ngroups);
 801     get_index(&(top.atoms), indexfn, ngroups, isize0, index0, grpname);
 802
 803     if (bMol)
 804     {
 805         mols    = &(top.mols);
 806         atndx   = mols->index;
 807         ngroups = isize0[0];
 808         snew(isize, ngroups);
 809         snew(index, ngroups);
 810         for (i = 0; i < ngroups; i++)
 811         {
 812             if (index0[0][i] < 0 || index0[0][i] >= mols->nr)
 813             {
 814                 gmx_fatal(FARGS, "Molecule index (%d) is out of range (%d-%d)", index0[0][i] + 1, 1, mols->nr);
 815             }
 816             isize[i] = atndx[index0[0][i] + 1] - atndx[index0[0][i]];
 817             snew(index[i], isize[i]);
 818             for (j = 0; j < isize[i]; j++)
 819             {
 820                 index[i][j] = atndx[index0[0][i]] + j;
 821             }
 822         }
 823     }
 824     else
 825     {
 826         isize = isize0;
 827         index = index0;
 828     }
 829     if (bCom)
 830     {
 831         snew(mass, top.atoms.nr);
 832         for (i = 0; i < top.atoms.nr; i++)
 833         {
 834             mass[i] = top.atoms.atom[i].m;
 835         }
 836     }
 837     else
 838     {
 839         mass = nullptr;
 840     }
 841
 842     flags = 0;
 843     std::string label(output_env_get_xvgr_tlabel(oenv));
 844     if (bOX)
 845     {
 846         flags = flags | TRX_READ_X;
 847         outx  = xvgropen(opt2fn("-ox", NFILE, fnm),
 848                         bCom ? "Center of mass" : "Coordinate",
 849                         label,
 850                         "Coordinate (nm)",
 851                         oenv);
 852         make_legend(outx, ngroups, isize0[0], index0[0], grpname, bCom, bMol, bDim, oenv);
 853     }
 854     if (bOXT)
 855     {
 856         flags      = flags | TRX_READ_X;
 857         status_out = open_trx(opt2fn("-oxt", NFILE, fnm), "w");
 858     }
 859     if (bOV)
 860     {
 861         flags = flags | TRX_READ_V;
 862         outv  = xvgropen(opt2fn("-ov", NFILE, fnm),
 863                         bCom ? "Center of mass velocity" : "Velocity",
 864                         label,
 865                         "Velocity (nm/ps)",
 866                         oenv);
 867         make_legend(outv, ngroups, isize0[0], index0[0], grpname, bCom, bMol, bDim, oenv);
 868     }
 869     if (bOF)
 870     {
 871         flags = flags | TRX_READ_F;
 872         outf  = xvgropen(
 873                 opt2fn("-of", NFILE, fnm), "Force", label, "Force (kJ mol\\S-1\\N nm\\S-1\\N)", oenv);
 874         make_legend(outf, ngroups, isize0[0], index0[0], grpname, bCom, bMol, bDim, oenv);
 875     }
 876     if (bOB)
 877     {
 878         outb = xvgropen(opt2fn("-ob", NFILE, fnm), "Box vector elements", label, "(nm)", oenv);
 879
 880         xvgr_legend(outb, 6, box_leg, oenv);
 881     }
 882     if (bOT)
 883     {
 884         bDum[XX]  = FALSE;
 885         bDum[YY]  = FALSE;
 886         bDum[ZZ]  = FALSE;
 887         bDum[DIM] = TRUE;
 888         flags     = flags | TRX_READ_V;
 889         outt      = xvgropen(opt2fn("-ot", NFILE, fnm), "Temperature", label, "(K)", oenv);
 890         make_legend(outt, ngroups, isize[0], index[0], grpname, bCom, bMol, bDum, oenv);
 891     }
 892     if (bEKT)
 893     {
 894         bDum[XX]  = FALSE;
 895         bDum[YY]  = FALSE;
 896         bDum[ZZ]  = FALSE;
 897         bDum[DIM] = TRUE;
 898         flags     = flags | TRX_READ_V;
 899         outekt    = xvgropen(opt2fn("-ekt", NFILE, fnm),
 900                           "Center of mass translation",
 901                           label,
 902                           "Energy (kJ mol\\S-1\\N)",
 903                           oenv);
 904         make_legend(outekt, ngroups, isize[0], index[0], grpname, bCom, bMol, bDum, oenv);
 905     }
 906     if (bEKR)
 907     {
 908         bDum[XX]  = FALSE;
 909         bDum[YY]  = FALSE;
 910         bDum[ZZ]  = FALSE;
 911         bDum[DIM] = TRUE;
 912         flags     = flags | TRX_READ_X | TRX_READ_V;
 913         outekr    = xvgropen(opt2fn("-ekr", NFILE, fnm),
 914                           "Center of mass rotation",
 915                           label,
 916                           "Energy (kJ mol\\S-1\\N)",
 917                           oenv);
 918         make_legend(outekr, ngroups, isize[0], index[0], grpname, bCom, bMol, bDum, oenv);
 919     }
 920     if (bVD)
 921     {
 922         flags = flags | TRX_READ_V;
 923     }
 924     if (bCV)
 925     {
 926         flags = flags | TRX_READ_X | TRX_READ_V;
 927     }
 928     if (bCF)
 929     {
 930         flags = flags | TRX_READ_X | TRX_READ_F;
 931     }
 932     if ((flags == 0) && !bOB)
 933     {
 934         fprintf(stderr, "Please select one or more output file options\n");
 935         exit(0);
 936     }
 937
 938     read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
 939
 940
 941     if ((bOV || bOF) && fn2ftp(ftp2fn(efTRX, NFILE, fnm)) == efXTC)
 942     {
 943         gmx_fatal(FARGS,
 944                   "Cannot extract velocities or forces since your input XTC file does not contain "
 945                   "them.");
 946     }
 947
 948     if (bCV || bCF)
 949     {
 950         snew(sumx, fr.natoms);
 951     }
 952     if (bCV)
 953     {
 954         snew(sumv, fr.natoms);
 955     }
 956     if (bCF)
 957     {
 958         snew(sumf, fr.natoms);
 959     }
 960     nr_xfr = 0;
 961     nr_vfr = 0;
 962     nr_ffr = 0;
 963
 964     if (bCom && bPBC)
 965     {
 966         gpbc = gmx_rmpbc_init(&top.idef, pbcType, fr.natoms);
 967     }
 968
 969     do
 970     {
 971         time = output_env_conv_time(oenv, fr.time);
 972
 973         if (fr.bX && bNoJump && fr.bBox)
 974         {
 975             if (xp)
 976             {
 977                 remove_jump(fr.box, fr.natoms, xp, fr.x);
 978             }
 979             else
 980             {
 981                 snew(xp, fr.natoms);
 982             }
 983             for (i = 0; i < fr.natoms; i++)
 984             {
 985                 copy_rvec(fr.x[i], xp[i]);
 986             }
 987         }
 988
 989         if (fr.bX && bCom && bPBC)
 990         {
 991             gmx_rmpbc_trxfr(gpbc, &fr);
 992         }
 993
 994         if (bVD && fr.bV)
 995         {
 996             update_histo(isize[0], index[0], fr.v, &nvhisto, &vhisto, binwidth);
 997         }
 998
 999         if (bOX && fr.bX)
1000         {
1001             print_data(outx, time, fr.x, mass, bCom, ngroups, isize, index, bDim, sffmt);
1002         }
1003         if (bOXT && fr.bX)
1004         {
1005             t_trxframe frout = fr;
1006             if (!frout.bAtoms)
1007             {
1008                 frout.atoms  = &top.atoms;
1009                 frout.bAtoms = TRUE;
1010             }
1011             frout.bV = FALSE;
1012             frout.bF = FALSE;
1013             write_trx_x(status_out, &frout, mass, bCom, ngroups, isize, index);
1014         }
1015         if (bOV && fr.bV)
1016         {
1017             print_data(outv, time, fr.v, mass, bCom, ngroups, isize, index, bDim, sffmt);
1018         }
1019         if (bOF && fr.bF)
1020         {
1021             print_data(outf, time, fr.f, nullptr, bCom, ngroups, isize, index, bDim, sffmt);
1022         }
1023         if (bOB && fr.bBox)
1024         {
1025             fprintf(outb, "\t%g", fr.time);
1026             fprintf(outb,
1027                     sffmt6.c_str(),
1028                     fr.box[XX][XX],
1029                     fr.box[YY][YY],
1030                     fr.box[ZZ][ZZ],
1031                     fr.box[YY][XX],
1032                     fr.box[ZZ][XX],
1033                     fr.box[ZZ][YY]);
1034             fprintf(outb, "\n");
1035         }
1036         if (bOT && fr.bV)
1037         {
1038             fprintf(outt, " %g", time);
1039             for (i = 0; i < ngroups; i++)
1040             {
1041                 fprintf(outt, sffmt, temp(fr.v, mass, isize[i], index[i]));
1042             }
1043             fprintf(outt, "\n");
1044         }
1045         if (bEKT && fr.bV)
1046         {
1047             fprintf(outekt, " %g", time);
1048             for (i = 0; i < ngroups; i++)
1049             {
1050                 fprintf(outekt, sffmt, ektrans(fr.v, mass, isize[i], index[i]));
1051             }
1052             fprintf(outekt, "\n");
1053         }
1054         if (bEKR && fr.bX && fr.bV)
1055         {
1056             fprintf(outekr, " %g", time);
1057             for (i = 0; i < ngroups; i++)
1058             {
1059                 fprintf(outekr, sffmt, ekrot(fr.x, fr.v, mass, isize[i], index[i]));
1060             }
1061             fprintf(outekr, "\n");
1062         }
1063         if ((bCV || bCF) && fr.bX
1064             && (ctime < 0 || (fr.time >= ctime * 0.999999 && fr.time <= ctime * 1.000001)))
1065         {
1066             for (i = 0; i < fr.natoms; i++)
1067             {
1068                 rvec_inc(sumx[i], fr.x[i]);
1069             }
1070             nr_xfr++;
1071         }
1072         if (bCV && fr.bV)
1073         {
1074             for (i = 0; i < fr.natoms; i++)
1075             {
1076                 rvec_inc(sumv[i], fr.v[i]);
1077             }
1078             nr_vfr++;
1079         }
1080         if (bCF && fr.bF)
1081         {
1082             for (i = 0; i < fr.natoms; i++)
1083             {
1084                 rvec_inc(sumf[i], fr.f[i]);
1085             }
1086             nr_ffr++;
1087         }
1088
1089     } while (read_next_frame(oenv, status, &fr));
1090
1091     if (gpbc != nullptr)
1092     {
1093         gmx_rmpbc_done(gpbc);
1094     }
1095
1096     /* clean up a bit */
1097     close_trx(status);
1098
1099     if (bOX)
1100     {
1101         xvgrclose(outx);
1102     }
1103     if (bOXT)
1104     {
1105         close_trx(status_out);
1106     }
1107     if (bOV)
1108     {
1109         xvgrclose(outv);
1110     }
1111     if (bOF)
1112     {
1113         xvgrclose(outf);
1114     }
1115     if (bOB)
1116     {
1117         xvgrclose(outb);
1118     }
1119     if (bOT)
1120     {
1121         xvgrclose(outt);
1122     }
1123     if (bEKT)
1124     {
1125         xvgrclose(outekt);
1126     }
1127     if (bEKR)
1128     {
1129         xvgrclose(outekr);
1130     }
1131
1132     if (bVD)
1133     {
1134         print_histo(opt2fn("-vd", NFILE, fnm), nvhisto, vhisto, binwidth, oenv);
1135     }
1136
1137     if (bCV || bCF)
1138     {
1139         if (nr_xfr > 1)
1140         {
1141             if (pbcType != PbcType::No && !bNoJump)
1142             {
1143                 fprintf(stderr,
1144                         "\nWARNING: More than one frame was used for option -cv or -cf\n"
1145                         "If atoms jump across the box you should use the -nojump or -ctime "
1146                         "option\n\n");
1147             }
1148             for (i = 0; i < isize[0]; i++)
1149             {
1150                 svmul(1.0 / nr_xfr, sumx[index[0][i]], sumx[index[0][i]]);
1151             }
1152         }
1153         else if (nr_xfr == 0)
1154         {
1155             fprintf(stderr, "\nWARNING: No coordinate frames found for option -cv or -cf\n\n");
1156         }
1157     }
1158     if (bCV)
1159     {
1160         write_pdb_bfac(opt2fn("-cv", NFILE, fnm),
1161                        opt2fn("-av", NFILE, fnm),
1162                        "average velocity",
1163                        &(top.atoms),
1164                        pbcType,
1165                        topbox,
1166                        isize[0],
1167                        index[0],
1168                        nr_xfr,
1169                        sumx,
1170                        nr_vfr,
1171                        sumv,
1172                        bDim,
1173                        scale,
1174                        oenv);
1175     }
1176     if (bCF)
1177     {
1178         write_pdb_bfac(opt2fn("-cf", NFILE, fnm),
1179                        opt2fn("-af", NFILE, fnm),
1180                        "average force",
1181                        &(top.atoms),
1182                        pbcType,
1183                        topbox,
1184                        isize[0],
1185                        index[0],
1186                        nr_xfr,
1187                        sumx,
1188                        nr_ffr,
1189                        sumf,
1190                        bDim,
1191                        scale,
1192                        oenv);
1193     }
1194
1195     /* view it */
1196     view_all(oenv, NFILE, fnm);
1197
1198     done_top(&top);
1199     // Free index and isize only if they are distinct from index0 and isize0
1200     if (bMol)
1201     {
1202         for (int i = 0; i < ngroups; i++)
1203         {
1204             sfree(index[i]);
1205         }
1206         sfree(index);
1207         sfree(isize);
1208     }
1209     for (int i = 0; i < ngroups; i++)
1210     {
1211         sfree(index0[i]);
1212         sfree(grpname[i]);
1213     }
1214     sfree(index0);
1215     sfree(isize0);
1216     sfree(grpname);
1217     done_frame(&fr);
1218     output_env_done(oenv);
1219
1220     return 0;
1221 }