src/tools/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                            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,mtot;
  99     rvec tmp;
 100     double sum[DIM];
 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,bool bCom,
 147                        int ngrps,int isize[],atom_id **index,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,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,bool bCom,bool bMol,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                            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     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_x;
 458         for(i=0; i<isize; i++)
 459         {
 460             svmul(scale,x[index[i]],x[index[i]]);
 461         }
 462         scale = 1.0/nfr_v;
 463         for(i=0; i<isize; i++)
 464         {
 465             svmul(scale,sum[index[i]],sum[index[i]]);
 466         }
 467
 468         fp=xvgropen(xname,title,"Atom","",oenv);
 469         for(i=0; i<isize; i++)
 470         {
 471             fprintf(fp,"%-5d  %10.3f  %10.3f  %10.3f\n",i,
 472                     sum[index[i]][XX],sum[index[i]][YY],sum[index[i]][ZZ]);
 473         }
 474         ffclose(fp);
 475         max  = 0;
 476         maxi = 0;
 477         for(i=0; i<isize; i++)
 478         {
 479             len2 = 0;
 480             for(m=0; m<DIM; m++)
 481             {
 482                 if (bDim[m] || bDim[DIM])
 483                 {
 484                     len2 += sqr(sum[index[i]][m]);
 485                 }
 486             }
 487             if (len2 > max)
 488             {
 489                 max  = len2;
 490                 maxi = index[i];
 491             }
 492         }
 493         if (scale_factor != 0)
 494         {
 495             scale = scale_factor;
 496         }
 497         else
 498         {
 499             if (max == 0)
 500             {
 501                 scale = 1;
 502             }
 503             else
 504             {
 505                 scale = 10.0/sqrt(max);
 506             }
 507         }
 508
 509         printf("Maximum %s is %g on atom %d %s, res. %s %d\n",
 510                title,sqrt(max)/nfr_v,maxi+1,*(atoms->atomname[maxi]),
 511                *(atoms->resinfo[atoms->atom[maxi].resind].name),
 512                atoms->resinfo[atoms->atom[maxi].resind].nr);
 513
 514         if (atoms->pdbinfo == NULL)
 515         {
 516             snew(atoms->pdbinfo,atoms->nr);
 517         }
 518         if (onedim == -1)
 519         {
 520             for(i=0; i<isize; i++)
 521             {
 522                 len2 = 0;
 523                 for(m=0; m<DIM; m++)
 524                 {
 525                     if (bDim[m] || bDim[DIM])
 526                     {
 527                         len2 += sqr(sum[index[i]][m]);
 528                     }
 529                 }
 530                 atoms->pdbinfo[index[i]].bfac = sqrt(len2)*scale;
 531             }
 532         }
 533         else
 534         {
 535             for(i=0; i<isize; i++)
 536             {
 537                 atoms->pdbinfo[index[i]].bfac = sum[index[i]][onedim]*scale;
 538             }
 539         }
 540         write_sto_conf_indexed(fname,title,atoms,x,NULL,ePBC,box,isize,index);
 541     }
 542 }
 543
 544 static void update_histo(int gnx,atom_id index[],rvec v[],
 545                          int *nhisto,int **histo,real binwidth)
 546 {
 547     int  i,m,in,nnn;
 548     real vn,vnmax;
 549
 550     if (histo == NULL)
 551     {
 552         vnmax = 0;
 553         for(i=0; (i<gnx); i++)
 554         {
 555             vn = norm(v[index[i]]);
 556             vnmax = max(vn,vnmax);
 557         }
 558         vnmax *= 2;
 559         *nhisto = 1+(vnmax/binwidth);
 560         snew(*histo,*nhisto);
 561     }
 562     for(i=0; (i<gnx); i++)
 563     {
 564         vn = norm(v[index[i]]);
 565         in = vn/binwidth;
 566         if (in >= *nhisto)
 567         {
 568             nnn = in+100;
 569             fprintf(stderr,"Extending histogram from %d to %d\n",*nhisto,nnn);
 570
 571             srenew(*histo,nnn);
 572             for(m=*nhisto; (m<nnn); m++)
 573             {
 574                 (*histo)[m] = 0;
 575             }
 576             *nhisto = nnn;
 577         }
 578         (*histo)[in]++;
 579     }
 580 }
 581
 582 static void print_histo(const char *fn,int nhisto,int histo[],real binwidth,
 583                         const output_env_t oenv)
 584 {
 585     FILE *fp;
 586     int i;
 587
 588     fp = xvgropen(fn,"Velocity distribution","V (nm/ps)","arbitrary units",
 589                   oenv);
 590     for(i=0; (i<nhisto); i++)
 591     {
 592         fprintf(fp,"%10.3e  %10d\n",i*binwidth,histo[i]);
 593     }
 594     ffclose(fp);
 595 }
 596
 597 int gmx_traj(int argc,char *argv[])
 598 {
 599     const char *desc[] = {
 600         "g_traj plots coordinates, velocities, forces and/or the box.",
 601         "With [TT]-com[tt] the coordinates, velocities and forces are",
 602         "calculated for the center of mass of each group.",
 603         "When [TT]-mol[tt] is set, the numbers in the index file are",
 604         "interpreted as molecule numbers and the same procedure as with",
 605         "[TT]-com[tt] is used for each molecule.[PAR]",
 606         "Option [TT]-ot[tt] plots the temperature of each group,",
 607         "provided velocities are present in the trajectory file.",
 608         "No corrections are made for constrained degrees of freedom!",
 609         "This implies [TT]-com[tt].[PAR]",
 610         "Options [TT]-ekt[tt] and [TT]-ekr[tt] plot the translational and",
 611         "rotational kinetic energy of each group,",
 612         "provided velocities are present in the trajectory file.",
 613         "This implies [TT]-com[tt].[PAR]",
 614         "Options [TT]-cv[tt] and [TT]-cf[tt] write the average velocities",
 615         "and average forces as temperature factors to a pdb file with",
 616         "the average coordinates. The temperature factors are scaled such",
 617         "that the maximum is 10. The scaling can be changed with the option",
 618         "[TT]-scale[tt]. To get the velocities or forces of one",
 619         "frame set both [TT]-b[tt] and [TT]-e[tt] to the time of",
 620         "desired frame. When averaging over frames you might need to use",
 621         "the [TT]-nojump[tt] option to obtain the correct average coordinates.",
 622         "If you select either of these option the average force and velocity",
 623         "for each atom are written to an xvg file as well",
 624         "(specified with [TT]-av[tt] or [TT]-af[tt]).[PAR]",
 625         "Option [TT]-vd[tt] computes a velocity distribution, i.e. the",
 626         "norm of the vector is plotted. In addition in the same graph",
 627         "the kinetic energy distribution is given."
 628     };
 629     static bool bMol=FALSE,bCom=FALSE,bPBC=TRUE,bNoJump=FALSE;
 630     static bool bX=TRUE,bY=TRUE,bZ=TRUE,bNorm=FALSE,bFP=FALSE;
 631     static int  ngroups=1;
 632     static real scale=0,binwidth=1;
 633     t_pargs pa[] = {
 634         { "-com", FALSE, etBOOL, {&bCom},
 635           "Plot data for the com of each group" },
 636         { "-pbc", FALSE, etBOOL, {&bPBC},
 637           "Make molecules whole for COM" },
 638         { "-mol", FALSE, etBOOL, {&bMol},
 639           "Index contains molecule numbers iso atom numbers" },
 640         { "-nojump", FALSE, etBOOL, {&bNoJump},
 641           "Remove jumps of atoms across the box" },
 642         { "-x", FALSE, etBOOL, {&bX},
 643           "Plot X-component" },
 644         { "-y", FALSE, etBOOL, {&bY},
 645           "Plot Y-component" },
 646         { "-z", FALSE, etBOOL, {&bZ},
 647           "Plot Z-component" },
 648         { "-ng",       FALSE, etINT, {&ngroups},
 649           "Number of groups to consider" },
 650         { "-len", FALSE, etBOOL, {&bNorm},
 651           "Plot vector length" },
 652         { "-fp", FALSE, etBOOL, {&bFP},
 653           "Full precision output" },
 654         { "-bin", FALSE, etREAL, {&binwidth},
 655           "Binwidth for velocity histogram (nm/ps)" },
 656         { "-scale", FALSE, etREAL, {&scale},
 657           "Scale factor for pdb output, 0 is autoscale" }
 658     };
 659     FILE       *outx=NULL,*outv=NULL,*outf=NULL,*outb=NULL,*outt=NULL;
 660     FILE       *outekt=NULL,*outekr=NULL;
 661     t_topology top;
 662     int        ePBC;
 663     real       *mass,time;
 664     char       title[STRLEN];
 665     const char *indexfn;
 666     t_trxframe fr,frout;
 667     int        flags,nvhisto=0,*vhisto=NULL;
 668     rvec       *xtop,*xp=NULL;
 669     rvec       *sumxv=NULL,*sumv=NULL,*sumxf=NULL,*sumf=NULL;
 670     matrix     topbox;
 671     t_trxstatus *status;
 672     t_trxstatus *status_out=NULL;
 673     gmx_rmpbc_t  gpbc=NULL;
 674     int        i,j,n;
 675     int        nr_xfr,nr_vfr,nr_ffr;
 676     char       **grpname;
 677     int        *isize0,*isize;
 678     atom_id    **index0,**index;
 679     atom_id    *atndx;
 680     t_block    *mols;
 681     bool       bTop,bOX,bOXT,bOV,bOF,bOB,bOT,bEKT,bEKR,bCV,bCF;
 682     bool       bDim[4],bDum[4],bVD;
 683     char       *sffmt,sffmt6[1024];
 684     const char *box_leg[6] = { "XX", "YY", "ZZ", "YX", "ZX", "ZY" };
 685     output_env_t oenv;
 686
 687     t_filenm fnm[] = {
 688         { efTRX, "-f", NULL, ffREAD },
 689         { efTPS, NULL, NULL, ffREAD },
 690         { efNDX, NULL, NULL, ffOPTRD },
 691         { efXVG, "-ox", "coord.xvg", ffOPTWR },
 692         { efTRX, "-oxt","coord.xtc", ffOPTWR },
 693         { efXVG, "-ov", "veloc.xvg", ffOPTWR },
 694         { efXVG, "-of", "force.xvg", ffOPTWR },
 695         { efXVG, "-ob", "box.xvg",   ffOPTWR },
 696         { efXVG, "-ot", "temp.xvg",  ffOPTWR },
 697         { efXVG, "-ekt","ektrans.xvg", ffOPTWR },
 698         { efXVG, "-ekr","ekrot.xvg", ffOPTWR },
 699         { efXVG, "-vd", "veldist.xvg", ffOPTWR },
 700         { efPDB, "-cv", "veloc.pdb", ffOPTWR },
 701         { efPDB, "-cf", "force.pdb", ffOPTWR },
 702         { efXVG, "-av", "all_veloc.xvg", ffOPTWR },
 703         { efXVG, "-af", "all_force.xvg", ffOPTWR }
 704     };
 705 #define NFILE asize(fnm)
 706
 707     CopyRight(stderr,argv[0]);
 708     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     if (bMol)
 713     {
 714         fprintf(stderr,"Interpreting indexfile entries as molecules.\n"
 715                 "Using center of mass.\n");
 716     }
 717
 718     bOX  = opt2bSet("-ox",NFILE,fnm);
 719     bOXT = opt2bSet("-oxt",NFILE,fnm);
 720     bOV  = opt2bSet("-ov",NFILE,fnm);
 721     bOF  = opt2bSet("-of",NFILE,fnm);
 722     bOB  = opt2bSet("-ob",NFILE,fnm);
 723     bOT  = opt2bSet("-ot",NFILE,fnm);
 724     bEKT = opt2bSet("-ekt",NFILE,fnm);
 725     bEKR = opt2bSet("-ekr",NFILE,fnm);
 726     bCV  = opt2bSet("-cv",NFILE,fnm) || opt2bSet("-av",NFILE,fnm);
 727     bCF  = opt2bSet("-cf",NFILE,fnm) || opt2bSet("-af",NFILE,fnm);
 728     bVD  = opt2bSet("-vd",NFILE,fnm) || opt2parg_bSet("-bin",asize(pa),pa);
 729     if (bMol || bOT || bEKT || bEKR)
 730     {
 731         bCom = TRUE;
 732     }
 733
 734     bDim[XX] = bX;
 735     bDim[YY] = bY;
 736     bDim[ZZ] = bZ;
 737     bDim[DIM] = bNorm;
 738
 739     if (bFP)
 740     {
 741         sffmt = "\t" gmx_real_fullprecision_pfmt;
 742     }
 743     else
 744     {
 745         sffmt = "\t%g";
 746     }
 747     sprintf(sffmt6,"%s%s%s%s%s%s",sffmt,sffmt,sffmt,sffmt,sffmt,sffmt);
 748
 749     bTop = read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,
 750                          &xtop,NULL,topbox,
 751                          bCom && (bOX || bOXT || bOV || bOT || bEKT || bEKR));
 752     sfree(xtop);
 753     if ((bMol || bCV || bCF) && !bTop)
 754     {
 755         gmx_fatal(FARGS,"Need a run input file for option -mol, -cv or -cf");
 756     }
 757
 758     if (bMol)
 759     {
 760         indexfn = ftp2fn(efNDX,NFILE,fnm);
 761     }
 762     else
 763     {
 764         indexfn = ftp2fn_null(efNDX,NFILE,fnm);
 765     }
 766
 767     if (!(bCom && !bMol))
 768     {
 769         ngroups = 1;
 770     }
 771     snew(grpname,ngroups);
 772     snew(isize0,ngroups);
 773     snew(index0,ngroups);
 774     get_index(&(top.atoms),indexfn,ngroups,isize0,index0,grpname);
 775
 776     if (bMol)
 777     {
 778         mols=&(top.mols);
 779         atndx = mols->index;
 780         ngroups = isize0[0];
 781         snew(isize,ngroups);
 782         snew(index,ngroups);
 783         for (i=0; i<ngroups; i++)
 784         {
 785             if (index0[0][i] < 0 || index0[0][i] >= mols->nr)
 786             {
 787                 gmx_fatal(FARGS,"Molecule index (%d) is out of range (%d-%d)",
 788                           index0[0][i]+1,1,mols->nr);
 789             }
 790             isize[i] = atndx[index0[0][i]+1] - atndx[index0[0][i]];
 791             snew(index[i],isize[i]);
 792             for(j=0; j<isize[i]; j++)
 793             {
 794                 index[i][j] = atndx[index0[0][i]] + j;
 795             }
 796         }
 797     }
 798     else
 799     {
 800         isize = isize0;
 801         index = index0;
 802     }
 803     if (bCom)
 804     {
 805         snew(mass,top.atoms.nr);
 806         for(i=0; i<top.atoms.nr; i++)
 807         {
 808             mass[i] = top.atoms.atom[i].m;
 809         }
 810     }
 811     else
 812     {
 813         mass = NULL;
 814     }
 815
 816     flags = 0;
 817     if (bOX)
 818     {
 819         flags = flags | TRX_READ_X;
 820         outx = xvgropen(opt2fn("-ox",NFILE,fnm),
 821                         bCom ? "Center of mass" : "Coordinate",
 822                         output_env_get_xvgr_tlabel(oenv),"Coordinate (nm)",oenv);
 823         make_legend(outx,ngroups,isize0[0],index0[0],grpname,bCom,bMol,bDim,oenv);
 824     }
 825     if (bOXT)
 826     {
 827         flags = flags | TRX_READ_X;
 828         status_out = open_trx(opt2fn("-oxt",NFILE,fnm),"w");
 829     }
 830     if (bOV)
 831     {
 832         flags = flags | TRX_READ_V;
 833         outv = xvgropen(opt2fn("-ov",NFILE,fnm),
 834                         bCom ? "Center of mass velocity" : "Velocity",
 835                         output_env_get_xvgr_tlabel(oenv),"Velocity (nm/ps)",oenv);
 836         make_legend(outv,ngroups,isize0[0],index0[0],grpname,bCom,bMol,bDim,oenv);
 837     }
 838     if (bOF) {
 839         flags = flags | TRX_READ_F;
 840         outf = xvgropen(opt2fn("-of",NFILE,fnm),"Force",
 841                         output_env_get_xvgr_tlabel(oenv),"Force (kJ mol\\S-1\\N nm\\S-1\\N)",
 842                         oenv);
 843         make_legend(outf,ngroups,isize0[0],index0[0],grpname,bCom,bMol,bDim,oenv);
 844     }
 845     if (bOB)
 846     {
 847         outb = xvgropen(opt2fn("-ob",NFILE,fnm),"Box vector elements",
 848                         output_env_get_xvgr_tlabel(oenv),"(nm)",oenv);
 849
 850         xvgr_legend(outb,6,box_leg,oenv);
 851     }
 852     if (bOT)
 853     {
 854         bDum[XX] = FALSE;
 855         bDum[YY] = FALSE;
 856         bDum[ZZ] = FALSE;
 857         bDum[DIM] = TRUE;
 858         flags = flags | TRX_READ_V;
 859         outt = xvgropen(opt2fn("-ot",NFILE,fnm),"Temperature",
 860                         output_env_get_xvgr_tlabel(oenv),"(K)", oenv);
 861         make_legend(outt,ngroups,isize[0],index[0],grpname,bCom,bMol,bDum,oenv);
 862     }
 863     if (bEKT)
 864     {
 865         bDum[XX] = FALSE;
 866         bDum[YY] = FALSE;
 867         bDum[ZZ] = FALSE;
 868         bDum[DIM] = TRUE;
 869         flags = flags | TRX_READ_V;
 870         outekt = xvgropen(opt2fn("-ekt",NFILE,fnm),"Center of mass translation",
 871                           output_env_get_xvgr_tlabel(oenv),"Energy (kJ mol\\S-1\\N)",oenv);
 872         make_legend(outekt,ngroups,isize[0],index[0],grpname,bCom,bMol,bDum,oenv);
 873     }
 874     if (bEKR)
 875     {
 876         bDum[XX] = FALSE;
 877         bDum[YY] = FALSE;
 878         bDum[ZZ] = FALSE;
 879         bDum[DIM] = TRUE;
 880         flags = flags | TRX_READ_X | TRX_READ_V;
 881         outekr = xvgropen(opt2fn("-ekr",NFILE,fnm),"Center of mass rotation",
 882                           output_env_get_xvgr_tlabel(oenv),"Energy (kJ mol\\S-1\\N)",oenv);
 883         make_legend(outekr,ngroups,isize[0],index[0],grpname,bCom,bMol,bDum,oenv);
 884     }
 885     if (bVD)
 886     {
 887         flags = flags | TRX_READ_V;
 888     }
 889     if (bCV)
 890     {
 891         flags = flags | TRX_READ_X | TRX_READ_V;
 892     }
 893     if (bCF)
 894     {
 895         flags = flags | TRX_READ_X | TRX_READ_F;
 896     }
 897     if ((flags == 0) && !bOB)
 898     {
 899         fprintf(stderr,"Please select one or more output file options\n");
 900         exit(0);
 901     }
 902
 903     read_first_frame(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&fr,flags);
 904
 905     if (bCV)
 906     {
 907         snew(sumxv,fr.natoms);
 908         snew(sumv,fr.natoms);
 909     }
 910     if (bCF)
 911     {
 912         snew(sumxf,fr.natoms);
 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)
1012         {
1013             if (fr.bX)
1014             {
1015                 for(i=0; i<fr.natoms; i++)
1016                 {
1017                     rvec_inc(sumxv[i],fr.x[i]);
1018                 }
1019                 nr_xfr++;
1020             }
1021             if (fr.bV)
1022             {
1023                 for(i=0; i<fr.natoms; i++)
1024                     rvec_inc(sumv[i],fr.v[i]);
1025                 nr_vfr++;
1026             }
1027         }
1028         if (bCF)
1029         {
1030             if (fr.bX)
1031             {
1032                 for(i=0; i<fr.natoms; i++)
1033                 {
1034                     rvec_inc(sumxf[i],fr.x[i]);
1035                 }
1036                 nr_xfr++;
1037             }
1038             if (fr.bF)
1039             {
1040                 for(i=0; i<fr.natoms; i++)
1041                 {
1042                     rvec_inc(sumf[i],fr.f[i]);
1043                 }
1044                 nr_ffr++;
1045             }
1046         }
1047
1048     } while(read_next_frame(oenv,status,&fr));
1049
1050     gmx_rmpbc_done(gpbc);
1051
1052     /* clean up a bit */
1053     close_trj(status);
1054
1055     if (bOX) ffclose(outx);
1056     if (bOXT) close_trx(status_out);
1057     if (bOV) ffclose(outv);
1058     if (bOF) ffclose(outf);
1059     if (bOB) ffclose(outb);
1060     if (bOT) ffclose(outt);
1061     if (bEKT) ffclose(outekt);
1062     if (bEKR) ffclose(outekr);
1063
1064     if (bVD)
1065     {
1066         print_histo(opt2fn("-vd",NFILE,fnm),nvhisto,vhisto,binwidth,oenv);
1067     }
1068
1069     if ((bCV || bCF) && (nr_vfr>1 || nr_ffr>1) && !bNoJump)
1070     {
1071         fprintf(stderr,"WARNING: More than one frame was used for option -cv or -cf\n"
1072                 "If atoms jump across the box you should use the -nojump option\n");
1073     }
1074
1075     if (bCV)
1076     {
1077         write_pdb_bfac(opt2fn("-cv",NFILE,fnm),
1078                        opt2fn("-av",NFILE,fnm),"average velocity",&(top.atoms),
1079                        ePBC,topbox,isize[0],index[0],nr_xfr,sumxv,
1080                        nr_vfr,sumv,bDim,scale,oenv);
1081     }
1082     if (bCF)
1083     {
1084         write_pdb_bfac(opt2fn("-cf",NFILE,fnm),
1085                        opt2fn("-af",NFILE,fnm),"average force",&(top.atoms),
1086                        ePBC,topbox,isize[0],index[0],nr_xfr,sumxf,
1087                        nr_ffr,sumf,bDim,scale,oenv);
1088     }
1089
1090     /* view it */
1091     view_all(oenv,NFILE, fnm);
1092
1093     thanx(stderr);
1094
1095     return 0;
1096 }