src/tools/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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43 #include <string.h>
  44 #include "statutil.h"
  45 #include "sysstuff.h"
  46 #include "typedefs.h"
  47 #include "smalloc.h"
  48 #include "macros.h"
  49 #include "vec.h"
  50 #include "pbc.h"
  51 #include "copyrite.h"
  52 #include "futil.h"
  53 #include "statutil.h"
  54 #include "index.h"
  55 #include "mshift.h"
  56 #include "xvgr.h"
  57 #include "tpxio.h"
  58 #include "rmpbc.h"
  59 #include "physics.h"
  60 #include "nrjac.h"
  61 #include "confio.h"
  62 #include "gmx_ana.h"
  63
  64
  65 static void low_print_data(FILE *fp,real time,rvec x[],int n,atom_id *index,
  66                            gmx_bool bDim[],const char *sffmt)
  67 {
  68     int i,ii,d;
  69
  70     fprintf(fp," %g",time);
  71     for(i=0; i<n; i++)
  72     {
  73         if (index != NULL)
  74         {
  75             ii = index[i];
  76         }
  77         else
  78         {
  79             ii = i;
  80         }
  81         for(d=0; d<DIM; d++)
  82         {
  83             if (bDim[d])
  84             {
  85                 fprintf(fp,sffmt,x[ii][d]);
  86             }
  87         }
  88         if (bDim[DIM])
  89         {
  90             fprintf(fp,sffmt,norm(x[ii]));
  91         }
  92     }
  93     fprintf(fp,"\n");
  94 }
  95
  96 static void average_data(rvec x[],rvec xav[],real *mass,
  97                          int ngrps,int isize[],atom_id **index)
  98 {
  99     int  g,i,ind,d;
 100     real m;
 101     rvec tmp;
 102     double sum[DIM],mtot;
 103
 104     for(g=0; g<ngrps; g++)
 105     {
 106         clear_dvec(sum);
 107         clear_rvec(xav[g]);
 108         mtot = 0;
 109         for(i=0; i<isize[g]; i++)
 110         {
 111             ind = index[g][i];
 112             if (mass != NULL)
 113             {
 114                 m = mass[ind];
 115                 svmul(m,x[ind],tmp);
 116                 for(d=0; d<DIM; d++)
 117                 {
 118                     sum[d] += tmp[d];
 119                 }
 120                 mtot += m;
 121             }
 122             else
 123             {
 124                 for(d=0; d<DIM; d++)
 125                 {
 126                     sum[d] += x[ind][d];
 127                 }
 128             }
 129         }
 130         if (mass != NULL)
 131         {
 132             for(d=0; d<DIM; d++)
 133             {
 134                 xav[g][d] = sum[d]/mtot;
 135             }
 136         }
 137         else
 138         {
 139             /* mass=NULL, so these are forces: sum only (do not average) */
 140             for(d=0; d<DIM; d++)
 141             {
 142                 xav[g][d] = sum[d];
 143             }
 144         }
 145     }
 146 }
 147
 148 static void print_data(FILE *fp,real time,rvec x[],real *mass,gmx_bool bCom,
 149                        int ngrps,int isize[],atom_id **index,gmx_bool bDim[],
 150                        const char *sffmt)
 151 {
 152     static rvec *xav=NULL;
 153
 154     if (bCom)
 155     {
 156         if (xav==NULL)
 157         {
 158             snew(xav,ngrps);
 159         }
 160         average_data(x,xav,mass,ngrps,isize,index);
 161         low_print_data(fp,time,xav,ngrps,NULL,bDim,sffmt);
 162     }
 163     else
 164     {
 165         low_print_data(fp,time,x,isize[0],index[0],bDim,sffmt);
 166     }
 167 }
 168
 169 static void write_trx_x(t_trxstatus *status,t_trxframe *fr,real *mass,gmx_bool bCom,
 170                         int ngrps,int isize[],atom_id **index)
 171 {
 172     static rvec *xav=NULL;
 173     static t_atoms *atoms=NULL;
 174     t_trxframe fr_av;
 175     int i;
 176
 177     fr->bV = FALSE;
 178     fr->bF = FALSE;
 179     if (bCom)
 180     {
 181         if (xav==NULL)
 182         {
 183             snew(xav,ngrps);
 184             snew(atoms,1);
 185             *atoms = *fr->atoms;
 186             snew(atoms->atom,ngrps);
 187             atoms->nr = ngrps;
 188             /* Note that atom and residue names will be the ones
 189              * of the first atom in each group.
 190              */
 191             for(i=0; i<ngrps; i++)
 192             {
 193                 atoms->atom[i]     = fr->atoms->atom[index[i][0]];
 194                 atoms->atomname[i] = fr->atoms->atomname[index[i][0]];
 195             }
 196         }
 197         average_data(fr->x,xav,mass,ngrps,isize,index);
 198         fr_av = *fr;
 199         fr_av.natoms = ngrps;
 200         fr_av.atoms  = atoms;
 201         fr_av.x      = xav;
 202         write_trxframe(status,&fr_av,NULL);
 203     }
 204     else
 205     {
 206         write_trxframe_indexed(status,fr,isize[0],index[0],NULL);
 207     }
 208 }
 209
 210 static void make_legend(FILE *fp,int ngrps,int isize,atom_id index[],
 211                         char **name,gmx_bool bCom,gmx_bool bMol,gmx_bool bDim[],
 212                         const output_env_t oenv)
 213 {
 214     char **leg;
 215     const char *dimtxt[] = { " X", " Y", " Z", "" };
 216     int n,i,j,d;
 217
 218     if (bCom)
 219     {
 220         n = ngrps;
 221     }
 222     else
 223     {
 224         n = isize;
 225     }
 226
 227     snew(leg,4*n);
 228     j=0;
 229     for(i=0; i<n; i++)
 230     {
 231         for(d=0; d<=DIM; d++)
 232         {
 233             if (bDim[d])
 234             {
 235                 snew(leg[j],STRLEN);
 236                 if (bMol)
 237                 {
 238                     sprintf(leg[j],"mol %d%s",index[i]+1,dimtxt[d]);
 239                 }
 240                 else if (bCom)
 241                 {
 242                     sprintf(leg[j],"%s%s",name[i],dimtxt[d]);
 243                 }
 244                 else
 245                 {
 246                     sprintf(leg[j],"atom %d%s",index[i]+1,dimtxt[d]);
 247                 }
 248                 j++;
 249             }
 250         }
 251     }
 252     xvgr_legend(fp,j,(const char**)leg,oenv);
 253
 254     for(i=0; i<j; i++)
 255     {
 256         sfree(leg[i]);
 257     }
 258     sfree(leg);
 259 }
 260
 261 static real ekrot(rvec x[],rvec v[],real mass[],int isize,atom_id index[])
 262 {
 263     static real **TCM=NULL,**L;
 264     double tm,m0,lxx,lxy,lxz,lyy,lyz,lzz,ekrot;
 265     rvec   a0,ocm;
 266     dvec   dx,b0;
 267     dvec   xcm,vcm,acm;
 268     int    i,j,m,n;
 269
 270     if (TCM == NULL)
 271     {
 272         snew(TCM,DIM);
 273         for(i=0; i<DIM; i++)
 274         {
 275             snew(TCM[i],DIM);
 276         }
 277         snew(L,DIM);
 278         for(i=0; i<DIM; i++)
 279         {
 280             snew(L[i],DIM);
 281         }
 282     }
 283
 284     clear_dvec(xcm);
 285     clear_dvec(vcm);
 286     clear_dvec(acm);
 287     tm=0.0;
 288     for(i=0; i<isize; i++)
 289     {
 290         j = index[i];
 291         m0 = mass[j];
 292         tm += m0;
 293         cprod(x[j],v[j],a0);
 294         for(m=0; (m<DIM); m++)
 295         {
 296             xcm[m] += m0*x[j][m]; /* c.o.m. position */
 297             vcm[m] += m0*v[j][m]; /* c.o.m. velocity */
 298             acm[m] += m0*a0[m];   /* rotational velocity around c.o.m. */
 299         }
 300     }
 301     dcprod(xcm,vcm,b0);
 302     for(m=0; (m<DIM); m++)
 303     {
 304         xcm[m] /= tm;
 305         vcm[m] /= tm;
 306         acm[m] -= b0[m]/tm;
 307   }
 308
 309     lxx=lxy=lxz=lyy=lyz=lzz=0.0;
 310     for(i=0; i<isize; i++)
 311     {
 312         j = index[i];
 313         m0 = mass[j];
 314         for(m=0; m<DIM; m++)
 315         {
 316             dx[m] = x[j][m] - xcm[m];
 317         }
 318         lxx += dx[XX]*dx[XX]*m0;
 319         lxy += dx[XX]*dx[YY]*m0;
 320         lxz += dx[XX]*dx[ZZ]*m0;
 321         lyy += dx[YY]*dx[YY]*m0;
 322         lyz += dx[YY]*dx[ZZ]*m0;
 323         lzz += dx[ZZ]*dx[ZZ]*m0;
 324     }
 325
 326     L[XX][XX] =  lyy + lzz;
 327     L[YY][XX] = -lxy;
 328     L[ZZ][XX] = -lxz;
 329     L[XX][YY] = -lxy;
 330     L[YY][YY] =  lxx + lzz;
 331     L[ZZ][YY] = -lyz;
 332     L[XX][ZZ] = -lxz;
 333     L[YY][ZZ] = -lyz;
 334     L[ZZ][ZZ] =  lxx + lyy;
 335
 336     m_inv_gen(L,DIM,TCM);
 337
 338     /* Compute omega (hoeksnelheid) */
 339     clear_rvec(ocm);
 340     ekrot = 0;
 341     for(m=0; m<DIM; m++)
 342     {
 343         for(n=0; n<DIM; n++)
 344         {
 345             ocm[m] += TCM[m][n]*acm[n];
 346         }
 347         ekrot += 0.5*ocm[m]*acm[m];
 348     }
 349
 350     return ekrot;
 351 }
 352
 353 static real ektrans(rvec v[],real mass[],int isize,atom_id index[])
 354 {
 355     dvec   mvcom;
 356     double mtot;
 357     int    i,j,d;
 358
 359     clear_dvec(mvcom);
 360     mtot = 0;
 361     for(i=0; i<isize; i++)
 362     {
 363         j = index[i];
 364         for(d=0; d<DIM; d++)
 365         {
 366             mvcom[d] += mass[j]*v[j][d];
 367         }
 368         mtot += mass[j];
 369     }
 370
 371     return dnorm2(mvcom)/(mtot*2);
 372 }
 373
 374 static real temp(rvec v[],real mass[],int isize,atom_id index[])
 375 {
 376     double ekin2;
 377     int  i,j;
 378
 379     ekin2 = 0;
 380     for(i=0; i<isize; i++)
 381     {
 382         j = index[i];
 383         ekin2 += mass[j]*norm2(v[j]);
 384     }
 385
 386     return ekin2/(3*isize*BOLTZ);
 387 }
 388
 389 static void remove_jump(matrix box,int natoms,rvec xp[],rvec x[])
 390 {
 391     rvec hbox;
 392     int d,i,m;
 393
 394     for(d=0; d<DIM; d++)
 395     {
 396         hbox[d] = 0.5*box[d][d];
 397     }
 398     for(i=0; i<natoms; i++)
 399     {
 400         for(m=DIM-1; m>=0; m--)
 401         {
 402             while (x[i][m]-xp[i][m] <= -hbox[m])
 403             {
 404                 for(d=0; d<=m; d++)
 405                 {
 406                     x[i][d] += box[m][d];
 407                 }
 408             }
 409             while (x[i][m]-xp[i][m] > hbox[m])
 410             {
 411                 for(d=0; d<=m; d++)
 412                 {
 413                     x[i][d] -= box[m][d];
 414                 }
 415             }
 416         }
 417     }
 418 }
 419
 420 static void write_pdb_bfac(const char *fname,const char *xname,
 421                            const char *title,t_atoms *atoms,int ePBC,matrix box,
 422                            int isize,atom_id *index,int nfr_x,rvec *x,
 423                            int nfr_v,rvec *sum,
 424                            gmx_bool bDim[],real scale_factor,
 425                            const output_env_t oenv)
 426 {
 427     FILE    *fp;
 428     real    max,len2,scale;
 429     atom_id maxi;
 430     int     i,m,onedim;
 431     gmx_bool    bOne;
 432
 433     if ((nfr_x == 0) || (nfr_v == 0))
 434     {
 435         fprintf(stderr,"No frames found for %s, will not write %s\n",
 436                 title,fname);
 437     }
 438     else
 439     {
 440         fprintf(stderr,"Used %d frames for %s\n",nfr_x,"coordinates");
 441         fprintf(stderr,"Used %d frames for %s\n",nfr_v,title);
 442         onedim = -1;
 443         if (!bDim[DIM])
 444         {
 445             m = 0;
 446             for(i=0; i<DIM; i++)
 447             {
 448                 if (bDim[i])
 449                 {
 450                     onedim = i;
 451                     m++;
 452                 }
 453             }
 454             if (m != 1)
 455             {
 456                 onedim = -1;
 457             }
 458         }
 459         scale = 1.0/nfr_v;
 460         for(i=0; i<isize; i++)
 461         {
 462             svmul(scale,sum[index[i]],sum[index[i]]);
 463         }
 464
 465         fp=xvgropen(xname,title,"Atom","",oenv);
 466         for(i=0; i<isize; i++)
 467         {
 468             fprintf(fp,"%-5d  %10.3f  %10.3f  %10.3f\n",1+i,
 469                     sum[index[i]][XX],sum[index[i]][YY],sum[index[i]][ZZ]);
 470         }
 471         ffclose(fp);
 472         max  = 0;
 473         maxi = 0;
 474         for(i=0; i<isize; i++)
 475         {
 476             len2 = 0;
 477             for(m=0; m<DIM; m++)
 478             {
 479                 if (bDim[m] || bDim[DIM])
 480                 {
 481                     len2 += sqr(sum[index[i]][m]);
 482                 }
 483             }
 484             if (len2 > max)
 485             {
 486                 max  = len2;
 487                 maxi = index[i];
 488             }
 489         }
 490         if (scale_factor != 0)
 491         {
 492             scale = scale_factor;
 493         }
 494         else
 495         {
 496             if (max == 0)
 497             {
 498                 scale = 1;
 499             }
 500             else
 501             {
 502                 scale = 10.0/sqrt(max);
 503             }
 504         }
 505
 506         printf("Maximum %s is %g on atom %d %s, res. %s %d\n",
 507                title,sqrt(max),maxi+1,*(atoms->atomname[maxi]),
 508                *(atoms->resinfo[atoms->atom[maxi].resind].name),
 509                atoms->resinfo[atoms->atom[maxi].resind].nr);
 510
 511         if (atoms->pdbinfo == NULL)
 512         {
 513             snew(atoms->pdbinfo,atoms->nr);
 514         }
 515         if (onedim == -1)
 516         {
 517             for(i=0; i<isize; i++)
 518             {
 519                 len2 = 0;
 520                 for(m=0; m<DIM; m++)
 521                 {
 522                     if (bDim[m] || bDim[DIM])
 523                     {
 524                         len2 += sqr(sum[index[i]][m]);
 525                     }
 526                 }
 527                 atoms->pdbinfo[index[i]].bfac = sqrt(len2)*scale;
 528             }
 529         }
 530         else
 531         {
 532             for(i=0; i<isize; i++)
 533             {
 534                 atoms->pdbinfo[index[i]].bfac = sum[index[i]][onedim]*scale;
 535             }
 536         }
 537         write_sto_conf_indexed(fname,title,atoms,x,NULL,ePBC,box,isize,index);
 538     }
 539 }
 540
 541 static void update_histo(int gnx,atom_id index[],rvec v[],
 542                          int *nhisto,int **histo,real binwidth)
 543 {
 544     int  i,m,in,nnn;
 545     real vn,vnmax;
 546
 547     if (*histo == NULL)
 548     {
 549         vnmax = 0;
 550         for(i=0; (i<gnx); i++)
 551         {
 552             vn = norm(v[index[i]]);
 553             vnmax = max(vn,vnmax);
 554         }
 555         vnmax *= 2;
 556         *nhisto = 1+(vnmax/binwidth);
 557         snew(*histo,*nhisto);
 558     }
 559     for(i=0; (i<gnx); i++)
 560     {
 561         vn = norm(v[index[i]]);
 562         in = vn/binwidth;
 563         if (in >= *nhisto)
 564         {
 565             nnn = in+100;
 566             fprintf(stderr,"Extending histogram from %d to %d\n",*nhisto,nnn);
 567
 568             srenew(*histo,nnn);
 569             for(m=*nhisto; (m<nnn); m++)
 570             {
 571                 (*histo)[m] = 0;
 572             }
 573             *nhisto = nnn;
 574         }
 575         (*histo)[in]++;
 576     }
 577 }
 578
 579 static void print_histo(const char *fn,int nhisto,int histo[],real binwidth,
 580                         const output_env_t oenv)
 581 {
 582     FILE *fp;
 583     int i;
 584
 585     fp = xvgropen(fn,"Velocity distribution","V (nm/ps)","arbitrary units",
 586                   oenv);
 587     for(i=0; (i<nhisto); i++)
 588     {
 589         fprintf(fp,"%10.3e  %10d\n",i*binwidth,histo[i]);
 590     }
 591     ffclose(fp);
 592 }
 593
 594 int gmx_traj(int argc,char *argv[])
 595 {
 596     const char *desc[] = {
 597         "[TT]g_traj[tt] plots coordinates, velocities, forces and/or the box.",
 598         "With [TT]-com[tt] the coordinates, velocities and forces are",
 599         "calculated for the center of mass of each group.",
 600         "When [TT]-mol[tt] is set, the numbers in the index file are",
 601         "interpreted as molecule numbers and the same procedure as with",
 602         "[TT]-com[tt] is used for each molecule.[PAR]",
 603         "Option [TT]-ot[tt] plots the temperature of each group,",
 604         "provided velocities are present in the trajectory file.",
 605         "No corrections are made for constrained degrees of freedom!",
 606         "This implies [TT]-com[tt].[PAR]",
 607         "Options [TT]-ekt[tt] and [TT]-ekr[tt] plot the translational and",
 608         "rotational kinetic energy of each group,",
 609         "provided velocities are present in the trajectory file.",
 610         "This implies [TT]-com[tt].[PAR]",
 611         "Options [TT]-cv[tt] and [TT]-cf[tt] write the average velocities",
 612         "and average forces as temperature factors to a [TT].pdb[tt] file with",
 613         "the average coordinates or the coordinates at [TT]-ctime[tt].",
 614         "The temperature factors are scaled such that the maximum is 10.",
 615         "The scaling can be changed with the option [TT]-scale[tt].",
 616         "To get the velocities or forces of one",
 617         "frame set both [TT]-b[tt] and [TT]-e[tt] to the time of",
 618         "desired frame. When averaging over frames you might need to use",
 619         "the [TT]-nojump[tt] option to obtain the correct average coordinates.",
 620         "If you select either of these option the average force and velocity",
 621         "for each atom are written to an [TT].xvg[tt] file as well",
 622         "(specified with [TT]-av[tt] or [TT]-af[tt]).[PAR]",
 623         "Option [TT]-vd[tt] computes a velocity distribution, i.e. the",
 624         "norm of the vector is plotted. In addition in the same graph",
 625         "the kinetic energy distribution is given."
 626     };
 627     static gmx_bool bMol=FALSE,bCom=FALSE,bPBC=TRUE,bNoJump=FALSE;
 628     static gmx_bool bX=TRUE,bY=TRUE,bZ=TRUE,bNorm=FALSE,bFP=FALSE;
 629     static int  ngroups=1;
 630     static real ctime=-1,scale=0,binwidth=1;
 631     t_pargs pa[] = {
 632         { "-com", FALSE, etBOOL, {&bCom},
 633           "Plot data for the com of each group" },
 634         { "-pbc", FALSE, etBOOL, {&bPBC},
 635           "Make molecules whole for COM" },
 636         { "-mol", FALSE, etBOOL, {&bMol},
 637           "Index contains molecule numbers iso atom numbers" },
 638         { "-nojump", FALSE, etBOOL, {&bNoJump},
 639           "Remove jumps of atoms across the box" },
 640         { "-x", FALSE, etBOOL, {&bX},
 641           "Plot X-component" },
 642         { "-y", FALSE, etBOOL, {&bY},
 643           "Plot Y-component" },
 644         { "-z", FALSE, etBOOL, {&bZ},
 645           "Plot Z-component" },
 646         { "-ng",       FALSE, etINT, {&ngroups},
 647           "Number of groups to consider" },
 648         { "-len", FALSE, etBOOL, {&bNorm},
 649           "Plot vector length" },
 650         { "-fp", FALSE, etBOOL, {&bFP},
 651           "Full precision output" },
 652         { "-bin", FALSE, etREAL, {&binwidth},
 653           "Binwidth for velocity histogram (nm/ps)" },
 654         { "-ctime", FALSE, etREAL, {&ctime},
 655           "Use frame at this time for x in [TT]-cv[tt] and [TT]-cf[tt] instead of the average x" },
 656         { "-scale", FALSE, etREAL, {&scale},
 657           "Scale factor for [TT].pdb[tt] 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       *sumx=NULL,*sumv=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     gmx_bool       bTop,bOX,bOXT,bOV,bOF,bOB,bOT,bEKT,bEKR,bCV,bCF;
 682     gmx_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 || bCF)
 906     {
 907         snew(sumx,fr.natoms);
 908     }
 909     if (bCV)
 910     {
 911         snew(sumv,fr.natoms);
 912     }
 913     if (bCF)
 914     {
 915         snew(sumf,fr.natoms);
 916     }
 917     nr_xfr = 0;
 918     nr_vfr = 0;
 919     nr_ffr = 0;
 920
 921     if (bCom && bPBC)
 922     {
 923         gpbc = gmx_rmpbc_init(&top.idef,ePBC,fr.natoms,fr.box);
 924     }
 925
 926     do
 927     {
 928         time = output_env_conv_time(oenv,fr.time);
 929
 930         if (fr.bX && bNoJump && fr.bBox)
 931         {
 932             if (xp)
 933             {
 934                 remove_jump(fr.box,fr.natoms,xp,fr.x);
 935             }
 936             else
 937             {
 938                 snew(xp,fr.natoms);
 939             }
 940             for(i=0; i<fr.natoms; i++)
 941             {
 942                 copy_rvec(fr.x[i],xp[i]);
 943             }
 944         }
 945
 946         if (fr.bX && bCom && bPBC)
 947         {
 948             gmx_rmpbc_trxfr(gpbc,&fr);
 949         }
 950
 951         if (bVD && fr.bV)
 952         {
 953             update_histo(isize[0],index[0],fr.v,&nvhisto,&vhisto,binwidth);
 954         }
 955
 956         if (bOX && fr.bX)
 957         {
 958             print_data(outx,time,fr.x,mass,bCom,ngroups,isize,index,bDim,sffmt);
 959         }
 960         if (bOXT && fr.bX)
 961         {
 962             frout = fr;
 963             if (!frout.bAtoms)
 964             {
 965                 frout.atoms  = &top.atoms;
 966                 frout.bAtoms = TRUE;
 967             }
 968             write_trx_x(status_out,&frout,mass,bCom,ngroups,isize,index);
 969         }
 970         if (bOV && fr.bV)
 971         {
 972             print_data(outv,time,fr.v,mass,bCom,ngroups,isize,index,bDim,sffmt);
 973         }
 974         if (bOF && fr.bF)
 975         {
 976             print_data(outf,time,fr.f,NULL,bCom,ngroups,isize,index,bDim,sffmt);
 977         }
 978         if (bOB && fr.bBox)
 979         {
 980             fprintf(outb,"\t%g",fr.time);
 981             fprintf(outb,sffmt6,
 982                     fr.box[XX][XX],fr.box[YY][YY],fr.box[ZZ][ZZ],
 983                     fr.box[YY][XX],fr.box[ZZ][XX],fr.box[ZZ][YY]);
 984             fprintf(outb,"\n");
 985         }
 986         if (bOT && fr.bV)
 987         {
 988             fprintf(outt," %g",time);
 989             for(i=0; i<ngroups; i++)
 990             {
 991                 fprintf(outt,sffmt,temp(fr.v,mass,isize[i],index[i]));
 992             }
 993             fprintf(outt,"\n");
 994         }
 995         if (bEKT && fr.bV)
 996         {
 997             fprintf(outekt," %g",time);
 998             for(i=0; i<ngroups; i++)
 999             {
1000                 fprintf(outekt,sffmt,ektrans(fr.v,mass,isize[i],index[i]));
1001             }
1002             fprintf(outekt,"\n");
1003         }
1004         if (bEKR && fr.bX && fr.bV)
1005         {
1006             fprintf(outekr," %g",time);
1007             for(i=0; i<ngroups; i++)
1008             {
1009                 fprintf(outekr,sffmt,ekrot(fr.x,fr.v,mass,isize[i],index[i]));
1010             }
1011             fprintf(outekr,"\n");
1012         }
1013         if ((bCV || bCF) && fr.bX &&
1014             (ctime < 0 || (fr.time >= ctime*0.999999 &&
1015                            fr.time <= ctime*1.000001)))
1016         {
1017             for(i=0; i<fr.natoms; i++)
1018             {
1019                 rvec_inc(sumx[i],fr.x[i]);
1020             }
1021             nr_xfr++;
1022         }
1023         if (bCV && fr.bV)
1024         {
1025             for(i=0; i<fr.natoms; i++)
1026             {
1027                 rvec_inc(sumv[i],fr.v[i]);
1028             }
1029             nr_vfr++;
1030         }
1031         if (bCF && fr.bF)
1032         {
1033             for(i=0; i<fr.natoms; i++)
1034             {
1035                 rvec_inc(sumf[i],fr.f[i]);
1036             }
1037             nr_ffr++;
1038         }
1039
1040     } while(read_next_frame(oenv,status,&fr));
1041
1042     if (gpbc != NULL)
1043     {
1044         gmx_rmpbc_done(gpbc);
1045     }
1046
1047     /* clean up a bit */
1048     close_trj(status);
1049
1050     if (bOX) ffclose(outx);
1051     if (bOXT) close_trx(status_out);
1052     if (bOV) ffclose(outv);
1053     if (bOF) ffclose(outf);
1054     if (bOB) ffclose(outb);
1055     if (bOT) ffclose(outt);
1056     if (bEKT) ffclose(outekt);
1057     if (bEKR) ffclose(outekr);
1058
1059     if (bVD)
1060     {
1061         print_histo(opt2fn("-vd",NFILE,fnm),nvhisto,vhisto,binwidth,oenv);
1062     }
1063
1064     if (bCV || bCF)
1065     {
1066         if (nr_xfr > 1)
1067         {
1068             if (ePBC != epbcNONE && !bNoJump)
1069             {
1070                 fprintf(stderr,"\nWARNING: More than one frame was used for option -cv or -cf\n"
1071                         "If atoms jump across the box you should use the -nojump or -ctime option\n\n");
1072             }
1073             for(i=0; i<isize[0]; i++)
1074             {
1075                 svmul(1.0/nr_xfr,sumx[index[0][i]],sumx[index[0][i]]);
1076             }
1077         }
1078         else if (nr_xfr == 0)
1079         {
1080             fprintf(stderr,"\nWARNING: No coordinate frames found for option -cv or -cf\n\n");
1081         }
1082     }
1083     if (bCV)
1084     {
1085         write_pdb_bfac(opt2fn("-cv",NFILE,fnm),
1086                        opt2fn("-av",NFILE,fnm),"average velocity",&(top.atoms),
1087                        ePBC,topbox,isize[0],index[0],nr_xfr,sumx,
1088                        nr_vfr,sumv,bDim,scale,oenv);
1089     }
1090     if (bCF)
1091     {
1092         write_pdb_bfac(opt2fn("-cf",NFILE,fnm),
1093                        opt2fn("-af",NFILE,fnm),"average force",&(top.atoms),
1094                        ePBC,topbox,isize[0],index[0],nr_xfr,sumx,
1095                        nr_ffr,sumf,bDim,scale,oenv);
1096     }
1097
1098     /* view it */
1099     view_all(oenv,NFILE, fnm);
1100
1101     thanx(stderr);
1102
1103     return 0;
1104 }