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