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, "%-5d  %10.3f  %10.3f  %10.3f\n", 1 + i, sum[index[i]][XX],
 482                     sum[index[i]][YY], sum[index[i]][ZZ]);
 483         }
 484         xvgrclose(fp);
 485         max  = 0;
 486         maxi = 0;
 487         for (i = 0; i < isize; i++)
 488         {
 489             len2 = 0;
 490             for (m = 0; m < DIM; m++)
 491             {
 492                 if (bDim[m] || bDim[DIM])
 493                 {
 494                     len2 += gmx::square(sum[index[i]][m]);
 495                 }
 496             }
 497             if (len2 > max)
 498             {
 499                 max  = len2;
 500                 maxi = index[i];
 501             }
 502         }
 503         if (scale_factor != 0)
 504         {
 505             scale = scale_factor;
 506         }
 507         else
 508         {
 509             if (max == 0)
 510             {
 511                 scale = 1;
 512             }
 513             else
 514             {
 515                 scale = 10.0 / std::sqrt(max);
 516             }
 517         }
 518
 519         printf("Maximum %s is %g on atom %d %s, res. %s %d\n", title, std::sqrt(max), maxi + 1,
 520                *(atoms->atomname[maxi]), *(atoms->resinfo[atoms->atom[maxi].resind].name),
 521                atoms->resinfo[atoms->atom[maxi].resind].nr);
 522
 523         if (atoms->pdbinfo == nullptr)
 524         {
 525             snew(atoms->pdbinfo, atoms->nr);
 526         }
 527         atoms->havePdbInfo = TRUE;
 528
 529         if (onedim == -1)
 530         {
 531             for (i = 0; i < isize; i++)
 532             {
 533                 len2 = 0;
 534                 for (m = 0; m < DIM; m++)
 535                 {
 536                     if (bDim[m] || bDim[DIM])
 537                     {
 538                         len2 += gmx::square(sum[index[i]][m]);
 539                     }
 540                 }
 541                 atoms->pdbinfo[index[i]].bfac = std::sqrt(len2) * scale;
 542             }
 543         }
 544         else
 545         {
 546             for (i = 0; i < isize; i++)
 547             {
 548                 atoms->pdbinfo[index[i]].bfac = sum[index[i]][onedim] * scale;
 549             }
 550         }
 551         write_sto_conf_indexed(fname, title, atoms, x, nullptr, pbcType, box, isize, index);
 552     }
 553 }
 554
 555 static void update_histo(int gnx, const int index[], rvec v[], int* nhisto, int** histo, real binwidth)
 556 {
 557     int  i, m, in, nnn;
 558     real vn, vnmax;
 559
 560     if (*histo == nullptr)
 561     {
 562         vnmax = 0;
 563         for (i = 0; (i < gnx); i++)
 564         {
 565             vn    = norm(v[index[i]]);
 566             vnmax = std::max(vn, vnmax);
 567         }
 568         vnmax *= 2;
 569         *nhisto = static_cast<int>(1 + (vnmax / binwidth));
 570         snew(*histo, *nhisto);
 571     }
 572     for (i = 0; (i < gnx); i++)
 573     {
 574         vn = norm(v[index[i]]);
 575         in = static_cast<int>(vn / binwidth);
 576         if (in >= *nhisto)
 577         {
 578             nnn = in + 100;
 579             fprintf(stderr, "Extending histogram from %d to %d\n", *nhisto, nnn);
 580
 581             srenew(*histo, nnn);
 582             for (m = *nhisto; (m < nnn); m++)
 583             {
 584                 (*histo)[m] = 0;
 585             }
 586             *nhisto = nnn;
 587         }
 588         (*histo)[in]++;
 589     }
 590 }
 591
 592 static void print_histo(const char* fn, int nhisto, int histo[], real binwidth, const gmx_output_env_t* oenv)
 593 {
 594     FILE* fp;
 595     int   i;
 596
 597     fp = xvgropen(fn, "Velocity distribution", "V (nm/ps)", "arbitrary units", oenv);
 598     for (i = 0; (i < nhisto); i++)
 599     {
 600         fprintf(fp, "%10.3e  %10d\n", i * binwidth, histo[i]);
 601     }
 602     xvgrclose(fp);
 603 }
 604
 605 int gmx_traj(int argc, char* argv[])
 606 {
 607     const char* desc[] = {
 608         "[THISMODULE] plots coordinates, velocities, forces and/or the box.",
 609         "With [TT]-com[tt] the coordinates, velocities and forces are",
 610         "calculated for the center of mass of each group.",
 611         "When [TT]-mol[tt] is set, the numbers in the index file are",
 612         "interpreted as molecule numbers and the same procedure as with",
 613         "[TT]-com[tt] is used for each molecule.[PAR]",
 614         "Option [TT]-ot[tt] plots the temperature of each group,",
 615         "provided velocities are present in the trajectory file.",
 616         "No corrections are made for constrained degrees of freedom!",
 617         "This implies [TT]-com[tt].[PAR]",
 618         "Options [TT]-ekt[tt] and [TT]-ekr[tt] plot the translational and",
 619         "rotational kinetic energy of each group,",
 620         "provided velocities are present in the trajectory file.",
 621         "This implies [TT]-com[tt].[PAR]",
 622         "Options [TT]-cv[tt] and [TT]-cf[tt] write the average velocities",
 623         "and average forces as temperature factors to a [REF].pdb[ref] file with",
 624         "the average coordinates or the coordinates at [TT]-ctime[tt].",
 625         "The temperature factors are scaled such that the maximum is 10.",
 626         "The scaling can be changed with the option [TT]-scale[tt].",
 627         "To get the velocities or forces of one",
 628         "frame set both [TT]-b[tt] and [TT]-e[tt] to the time of",
 629         "desired frame. When averaging over frames you might need to use",
 630         "the [TT]-nojump[tt] option to obtain the correct average coordinates.",
 631         "If you select either of these option the average force and velocity",
 632         "for each atom are written to an [REF].xvg[ref] file as well",
 633         "(specified with [TT]-av[tt] or [TT]-af[tt]).[PAR]",
 634         "Option [TT]-vd[tt] computes a velocity distribution, i.e. the",
 635         "norm of the vector is plotted. In addition in the same graph",
 636         "the kinetic energy distribution is given.",
 637         "",
 638         "See [gmx-trajectory] for plotting similar data for selections."
 639     };
 640     static gmx_bool bMol = FALSE, bCom = FALSE, bPBC = TRUE, bNoJump = FALSE;
 641     static gmx_bool bX = TRUE, bY = TRUE, bZ = TRUE, bNorm = FALSE, bFP = FALSE;
 642     static int      ngroups = 1;
 643     static real     ctime = -1, scale = 0, binwidth = 1;
 644     t_pargs         pa[] = {
 645         { "-com", FALSE, etBOOL, { &bCom }, "Plot data for the com of each group" },
 646         { "-pbc", FALSE, etBOOL, { &bPBC }, "Make molecules whole for COM" },
 647         { "-mol",
 648           FALSE,
 649           etBOOL,
 650           { &bMol },
 651           "Index contains molecule numbers instead of atom numbers" },
 652         { "-nojump", FALSE, etBOOL, { &bNoJump }, "Remove jumps of atoms across the box" },
 653         { "-x", FALSE, etBOOL, { &bX }, "Plot X-component" },
 654         { "-y", FALSE, etBOOL, { &bY }, "Plot Y-component" },
 655         { "-z", FALSE, etBOOL, { &bZ }, "Plot Z-component" },
 656         { "-ng", FALSE, etINT, { &ngroups }, "Number of groups to consider" },
 657         { "-len", FALSE, etBOOL, { &bNorm }, "Plot vector length" },
 658         { "-fp", FALSE, etBOOL, { &bFP }, "Full precision output" },
 659         { "-bin", FALSE, etREAL, { &binwidth }, "Binwidth for velocity histogram (nm/ps)" },
 660         { "-ctime",
 661           FALSE,
 662           etREAL,
 663           { &ctime },
 664           "Use frame at this time for x in [TT]-cv[tt] and [TT]-cf[tt] instead of the average x" },
 665         { "-scale",
 666           FALSE,
 667           etREAL,
 668           { &scale },
 669           "Scale factor for [REF].pdb[ref] output, 0 is autoscale" }
 670     };
 671     FILE *       outx = nullptr, *outv = nullptr, *outf = nullptr, *outb = nullptr, *outt = nullptr;
 672     FILE *       outekt = nullptr, *outekr = nullptr;
 673     t_topology   top;
 674     PbcType      pbcType;
 675     real *       mass, time;
 676     const char*  indexfn;
 677     t_trxframe   fr;
 678     int          flags, nvhisto = 0, *vhisto = nullptr;
 679     rvec *       xtop, *xp = nullptr;
 680     rvec *       sumx = nullptr, *sumv = nullptr, *sumf = nullptr;
 681     matrix       topbox;
 682     t_trxstatus* status;
 683     t_trxstatus* status_out = nullptr;
 684     gmx_rmpbc_t  gpbc       = nullptr;
 685     int          i, j;
 686     int          nr_xfr, nr_vfr, nr_ffr;
 687     char**       grpname;
 688     int *        isize0, *isize;
 689     int **       index0, **index;
 690     int*         atndx;
 691     t_block*     mols;
 692     gmx_bool     bTop, bOX, bOXT, bOV, bOF, bOB, bOT, bEKT, bEKR, bCV, bCF;
 693     gmx_bool     bDim[4], bDum[4], bVD;
 694     char         sffmt[STRLEN];
 695     const char*  box_leg[6] = { "XX", "YY", "ZZ", "YX", "ZX", "ZY" };
 696     gmx_output_env_t* oenv;
 697
 698     t_filenm fnm[] = {
 699         { efTRX, "-f", nullptr, ffREAD },       { efTPS, nullptr, nullptr, ffREAD },
 700         { efNDX, nullptr, nullptr, ffOPTRD },   { efXVG, "-ox", "coord", ffOPTWR },
 701         { efTRX, "-oxt", "coord", ffOPTWR },    { efXVG, "-ov", "veloc", ffOPTWR },
 702         { efXVG, "-of", "force", ffOPTWR },     { efXVG, "-ob", "box", ffOPTWR },
 703         { efXVG, "-ot", "temp", ffOPTWR },      { efXVG, "-ekt", "ektrans", ffOPTWR },
 704         { efXVG, "-ekr", "ekrot", ffOPTWR },    { efXVG, "-vd", "veldist", ffOPTWR },
 705         { efPDB, "-cv", "veloc", ffOPTWR },     { efPDB, "-cf", "force", ffOPTWR },
 706         { efXVG, "-av", "all_veloc", ffOPTWR }, { efXVG, "-af", "all_force", ffOPTWR }
 707     };
 708 #define NFILE asize(fnm)
 709
 710     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW, NFILE, fnm,
 711                            asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 712     {
 713         return 0;
 714     }
 715
 716     if (bMol)
 717     {
 718         fprintf(stderr,
 719                 "Interpreting indexfile entries as molecules.\n"
 720                 "Using center of mass.\n");
 721     }
 722
 723     bOX  = opt2bSet("-ox", NFILE, fnm);
 724     bOXT = opt2bSet("-oxt", NFILE, fnm);
 725     bOV  = opt2bSet("-ov", NFILE, fnm);
 726     bOF  = opt2bSet("-of", NFILE, fnm);
 727     bOB  = opt2bSet("-ob", NFILE, fnm);
 728     bOT  = opt2bSet("-ot", NFILE, fnm);
 729     bEKT = opt2bSet("-ekt", NFILE, fnm);
 730     bEKR = opt2bSet("-ekr", NFILE, fnm);
 731     bCV  = opt2bSet("-cv", NFILE, fnm) || opt2bSet("-av", NFILE, fnm);
 732     bCF  = opt2bSet("-cf", NFILE, fnm) || opt2bSet("-af", NFILE, fnm);
 733     bVD  = opt2bSet("-vd", NFILE, fnm) || opt2parg_bSet("-bin", asize(pa), pa);
 734     if (bMol || bOT || bEKT || bEKR)
 735     {
 736         bCom = TRUE;
 737     }
 738
 739     bDim[XX]  = bX;
 740     bDim[YY]  = bY;
 741     bDim[ZZ]  = bZ;
 742     bDim[DIM] = bNorm;
 743
 744     if (bFP)
 745     {
 746         sprintf(sffmt, "\t%s", gmx_real_fullprecision_pfmt);
 747     }
 748     else
 749     {
 750         sprintf(sffmt, "\t%%g");
 751     }
 752     std::string sffmt6 = gmx::formatString("%s%s%s%s%s%s", sffmt, sffmt, sffmt, sffmt, sffmt, sffmt);
 753
 754     bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xtop, nullptr, topbox,
 755                          bCom && (bOX || bOXT || bOV || bOT || bEKT || bEKR));
 756     sfree(xtop);
 757     if ((bMol || bCV || bCF) && !bTop)
 758     {
 759         gmx_fatal(FARGS, "Need a run input file for option -mol, -cv or -cf");
 760     }
 761
 762     if (bMol)
 763     {
 764         indexfn = ftp2fn(efNDX, NFILE, fnm);
 765     }
 766     else
 767     {
 768         indexfn = ftp2fn_null(efNDX, NFILE, fnm);
 769     }
 770
 771     if (!(bCom && !bMol))
 772     {
 773         ngroups = 1;
 774     }
 775     snew(grpname, ngroups);
 776     snew(isize0, ngroups);
 777     snew(index0, ngroups);
 778     get_index(&(top.atoms), indexfn, ngroups, isize0, index0, grpname);
 779
 780     if (bMol)
 781     {
 782         mols    = &(top.mols);
 783         atndx   = mols->index;
 784         ngroups = isize0[0];
 785         snew(isize, ngroups);
 786         snew(index, ngroups);
 787         for (i = 0; i < ngroups; i++)
 788         {
 789             if (index0[0][i] < 0 || index0[0][i] >= mols->nr)
 790             {
 791                 gmx_fatal(FARGS, "Molecule index (%d) is out of range (%d-%d)", index0[0][i] + 1, 1,
 792                           mols->nr);
 793             }
 794             isize[i] = atndx[index0[0][i] + 1] - atndx[index0[0][i]];
 795             snew(index[i], isize[i]);
 796             for (j = 0; j < isize[i]; j++)
 797             {
 798                 index[i][j] = atndx[index0[0][i]] + j;
 799             }
 800         }
 801     }
 802     else
 803     {
 804         isize = isize0;
 805         index = index0;
 806     }
 807     if (bCom)
 808     {
 809         snew(mass, top.atoms.nr);
 810         for (i = 0; i < top.atoms.nr; i++)
 811         {
 812             mass[i] = top.atoms.atom[i].m;
 813         }
 814     }
 815     else
 816     {
 817         mass = nullptr;
 818     }
 819
 820     flags = 0;
 821     std::string label(output_env_get_xvgr_tlabel(oenv));
 822     if (bOX)
 823     {
 824         flags = flags | TRX_READ_X;
 825         outx  = xvgropen(opt2fn("-ox", NFILE, fnm), bCom ? "Center of mass" : "Coordinate", label,
 826                         "Coordinate (nm)", oenv);
 827         make_legend(outx, ngroups, isize0[0], index0[0], grpname, bCom, bMol, bDim, oenv);
 828     }
 829     if (bOXT)
 830     {
 831         flags      = flags | TRX_READ_X;
 832         status_out = open_trx(opt2fn("-oxt", NFILE, fnm), "w");
 833     }
 834     if (bOV)
 835     {
 836         flags = flags | TRX_READ_V;
 837         outv  = xvgropen(opt2fn("-ov", NFILE, fnm), bCom ? "Center of mass velocity" : "Velocity",
 838                         label, "Velocity (nm/ps)", oenv);
 839         make_legend(outv, ngroups, isize0[0], index0[0], grpname, bCom, bMol, bDim, oenv);
 840     }
 841     if (bOF)
 842     {
 843         flags = flags | TRX_READ_F;
 844         outf  = xvgropen(opt2fn("-of", NFILE, fnm), "Force", label,
 845                         "Force (kJ mol\\S-1\\N nm\\S-1\\N)", oenv);
 846         make_legend(outf, ngroups, isize0[0], index0[0], grpname, bCom, bMol, bDim, oenv);
 847     }
 848     if (bOB)
 849     {
 850         outb = xvgropen(opt2fn("-ob", NFILE, fnm), "Box vector elements", label, "(nm)", oenv);
 851
 852         xvgr_legend(outb, 6, box_leg, oenv);
 853     }
 854     if (bOT)
 855     {
 856         bDum[XX]  = FALSE;
 857         bDum[YY]  = FALSE;
 858         bDum[ZZ]  = FALSE;
 859         bDum[DIM] = TRUE;
 860         flags     = flags | TRX_READ_V;
 861         outt      = xvgropen(opt2fn("-ot", NFILE, fnm), "Temperature", label, "(K)", oenv);
 862         make_legend(outt, ngroups, isize[0], index[0], grpname, bCom, bMol, bDum, oenv);
 863     }
 864     if (bEKT)
 865     {
 866         bDum[XX]  = FALSE;
 867         bDum[YY]  = FALSE;
 868         bDum[ZZ]  = FALSE;
 869         bDum[DIM] = TRUE;
 870         flags     = flags | TRX_READ_V;
 871         outekt    = xvgropen(opt2fn("-ekt", NFILE, fnm), "Center of mass translation", label,
 872                           "Energy (kJ mol\\S-1\\N)", oenv);
 873         make_legend(outekt, ngroups, isize[0], index[0], grpname, bCom, bMol, bDum, oenv);
 874     }
 875     if (bEKR)
 876     {
 877         bDum[XX]  = FALSE;
 878         bDum[YY]  = FALSE;
 879         bDum[ZZ]  = FALSE;
 880         bDum[DIM] = TRUE;
 881         flags     = flags | TRX_READ_X | TRX_READ_V;
 882         outekr    = xvgropen(opt2fn("-ekr", NFILE, fnm), "Center of mass rotation", label,
 883                           "Energy (kJ mol\\S-1\\N)", oenv);
 884         make_legend(outekr, ngroups, isize[0], index[0], grpname, bCom, bMol, bDum, oenv);
 885     }
 886     if (bVD)
 887     {
 888         flags = flags | TRX_READ_V;
 889     }
 890     if (bCV)
 891     {
 892         flags = flags | TRX_READ_X | TRX_READ_V;
 893     }
 894     if (bCF)
 895     {
 896         flags = flags | TRX_READ_X | TRX_READ_F;
 897     }
 898     if ((flags == 0) && !bOB)
 899     {
 900         fprintf(stderr, "Please select one or more output file options\n");
 901         exit(0);
 902     }
 903
 904     read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
 905
 906
 907     if ((bOV || bOF) && fn2ftp(ftp2fn(efTRX, NFILE, fnm)) == efXTC)
 908     {
 909         gmx_fatal(FARGS,
 910                   "Cannot extract velocities or forces since your input XTC file does not contain "
 911                   "them.");
 912     }
 913
 914     if (bCV || bCF)
 915     {
 916         snew(sumx, fr.natoms);
 917     }
 918     if (bCV)
 919     {
 920         snew(sumv, fr.natoms);
 921     }
 922     if (bCF)
 923     {
 924         snew(sumf, fr.natoms);
 925     }
 926     nr_xfr = 0;
 927     nr_vfr = 0;
 928     nr_ffr = 0;
 929
 930     if (bCom && bPBC)
 931     {
 932         gpbc = gmx_rmpbc_init(&top.idef, pbcType, fr.natoms);
 933     }
 934
 935     do
 936     {
 937         time = output_env_conv_time(oenv, fr.time);
 938
 939         if (fr.bX && bNoJump && fr.bBox)
 940         {
 941             if (xp)
 942             {
 943                 remove_jump(fr.box, fr.natoms, xp, fr.x);
 944             }
 945             else
 946             {
 947                 snew(xp, fr.natoms);
 948             }
 949             for (i = 0; i < fr.natoms; i++)
 950             {
 951                 copy_rvec(fr.x[i], xp[i]);
 952             }
 953         }
 954
 955         if (fr.bX && bCom && bPBC)
 956         {
 957             gmx_rmpbc_trxfr(gpbc, &fr);
 958         }
 959
 960         if (bVD && fr.bV)
 961         {
 962             update_histo(isize[0], index[0], fr.v, &nvhisto, &vhisto, binwidth);
 963         }
 964
 965         if (bOX && fr.bX)
 966         {
 967             print_data(outx, time, fr.x, mass, bCom, ngroups, isize, index, bDim, sffmt);
 968         }
 969         if (bOXT && fr.bX)
 970         {
 971             t_trxframe frout = fr;
 972             if (!frout.bAtoms)
 973             {
 974                 frout.atoms  = &top.atoms;
 975                 frout.bAtoms = TRUE;
 976             }
 977             frout.bV = FALSE;
 978             frout.bF = FALSE;
 979             write_trx_x(status_out, &frout, mass, bCom, ngroups, isize, index);
 980         }
 981         if (bOV && fr.bV)
 982         {
 983             print_data(outv, time, fr.v, mass, bCom, ngroups, isize, index, bDim, sffmt);
 984         }
 985         if (bOF && fr.bF)
 986         {
 987             print_data(outf, time, fr.f, nullptr, bCom, ngroups, isize, index, bDim, sffmt);
 988         }
 989         if (bOB && fr.bBox)
 990         {
 991             fprintf(outb, "\t%g", fr.time);
 992             fprintf(outb, sffmt6.c_str(), fr.box[XX][XX], fr.box[YY][YY], fr.box[ZZ][ZZ],
 993                     fr.box[YY][XX], fr.box[ZZ][XX], fr.box[ZZ][YY]);
 994             fprintf(outb, "\n");
 995         }
 996         if (bOT && fr.bV)
 997         {
 998             fprintf(outt, " %g", time);
 999             for (i = 0; i < ngroups; i++)
1000             {
1001                 fprintf(outt, sffmt, temp(fr.v, mass, isize[i], index[i]));
1002             }
1003             fprintf(outt, "\n");
1004         }
1005         if (bEKT && fr.bV)
1006         {
1007             fprintf(outekt, " %g", time);
1008             for (i = 0; i < ngroups; i++)
1009             {
1010                 fprintf(outekt, sffmt, ektrans(fr.v, mass, isize[i], index[i]));
1011             }
1012             fprintf(outekt, "\n");
1013         }
1014         if (bEKR && fr.bX && fr.bV)
1015         {
1016             fprintf(outekr, " %g", time);
1017             for (i = 0; i < ngroups; i++)
1018             {
1019                 fprintf(outekr, sffmt, ekrot(fr.x, fr.v, mass, isize[i], index[i]));
1020             }
1021             fprintf(outekr, "\n");
1022         }
1023         if ((bCV || bCF) && fr.bX
1024             && (ctime < 0 || (fr.time >= ctime * 0.999999 && fr.time <= ctime * 1.000001)))
1025         {
1026             for (i = 0; i < fr.natoms; i++)
1027             {
1028                 rvec_inc(sumx[i], fr.x[i]);
1029             }
1030             nr_xfr++;
1031         }
1032         if (bCV && fr.bV)
1033         {
1034             for (i = 0; i < fr.natoms; i++)
1035             {
1036                 rvec_inc(sumv[i], fr.v[i]);
1037             }
1038             nr_vfr++;
1039         }
1040         if (bCF && fr.bF)
1041         {
1042             for (i = 0; i < fr.natoms; i++)
1043             {
1044                 rvec_inc(sumf[i], fr.f[i]);
1045             }
1046             nr_ffr++;
1047         }
1048
1049     } while (read_next_frame(oenv, status, &fr));
1050
1051     if (gpbc != nullptr)
1052     {
1053         gmx_rmpbc_done(gpbc);
1054     }
1055
1056     /* clean up a bit */
1057     close_trx(status);
1058
1059     if (bOX)
1060     {
1061         xvgrclose(outx);
1062     }
1063     if (bOXT)
1064     {
1065         close_trx(status_out);
1066     }
1067     if (bOV)
1068     {
1069         xvgrclose(outv);
1070     }
1071     if (bOF)
1072     {
1073         xvgrclose(outf);
1074     }
1075     if (bOB)
1076     {
1077         xvgrclose(outb);
1078     }
1079     if (bOT)
1080     {
1081         xvgrclose(outt);
1082     }
1083     if (bEKT)
1084     {
1085         xvgrclose(outekt);
1086     }
1087     if (bEKR)
1088     {
1089         xvgrclose(outekr);
1090     }
1091
1092     if (bVD)
1093     {
1094         print_histo(opt2fn("-vd", NFILE, fnm), nvhisto, vhisto, binwidth, oenv);
1095     }
1096
1097     if (bCV || bCF)
1098     {
1099         if (nr_xfr > 1)
1100         {
1101             if (pbcType != PbcType::No && !bNoJump)
1102             {
1103                 fprintf(stderr,
1104                         "\nWARNING: More than one frame was used for option -cv or -cf\n"
1105                         "If atoms jump across the box you should use the -nojump or -ctime "
1106                         "option\n\n");
1107             }
1108             for (i = 0; i < isize[0]; i++)
1109             {
1110                 svmul(1.0 / nr_xfr, sumx[index[0][i]], sumx[index[0][i]]);
1111             }
1112         }
1113         else if (nr_xfr == 0)
1114         {
1115             fprintf(stderr, "\nWARNING: No coordinate frames found for option -cv or -cf\n\n");
1116         }
1117     }
1118     if (bCV)
1119     {
1120         write_pdb_bfac(opt2fn("-cv", NFILE, fnm), opt2fn("-av", NFILE, fnm), "average velocity",
1121                        &(top.atoms), pbcType, topbox, isize[0], index[0], nr_xfr, sumx, nr_vfr,
1122                        sumv, bDim, scale, oenv);
1123     }
1124     if (bCF)
1125     {
1126         write_pdb_bfac(opt2fn("-cf", NFILE, fnm), opt2fn("-af", NFILE, fnm), "average force",
1127                        &(top.atoms), pbcType, topbox, isize[0], index[0], nr_xfr, sumx, nr_ffr,
1128                        sumf, bDim, scale, oenv);
1129     }
1130
1131     /* view it */
1132     view_all(oenv, NFILE, fnm);
1133
1134     done_top(&top);
1135     // Free index and isize only if they are distinct from index0 and isize0
1136     if (bMol)
1137     {
1138         for (int i = 0; i < ngroups; i++)
1139         {
1140             sfree(index[i]);
1141         }
1142         sfree(index);
1143         sfree(isize);
1144     }
1145     for (int i = 0; i < ngroups; i++)
1146     {
1147         sfree(index0[i]);
1148         sfree(grpname[i]);
1149     }
1150     sfree(index0);
1151     sfree(isize0);
1152     sfree(grpname);
1153     done_frame(&fr);
1154     output_env_done(oenv);
1155
1156     return 0;
1157 }