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