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