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