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