src/gromacs/gmxana/gmx_traj.c

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