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