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