src/gromacs/gmxana/gmx_traj.c

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