src/tools/g_traj.c

   1 /*
   2  * $Id$
   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 #include <math.h>
  37 #include <string.h>
  38 #include "statutil.h"
  39 #include "sysstuff.h"
  40 #include "typedefs.h"
  41 #include "smalloc.h"
  42 #include "macros.h"
  43 #include "vec.h"
  44 #include "pbc.h"
  45 #include "copyrite.h"
  46 #include "futil.h"
  47 #include "statutil.h"
  48 #include "index.h"
  49 #include "mshift.h"
  50 #include "xvgr.h"
  51 #include "tpxio.h"
  52 #include "rmpbc.h"
  53 #include "physics.h"
  54 #include "nrjac.h"
  55 #include "confio.h"
  56
  57 static void low_print_data(FILE *fp,real time,rvec x[],int n,atom_id *index,
  58                            bool bDim[])
  59 {
  60   int i,d;
  61
  62   fprintf(fp," %g",time);
  63   for(i=0; i<n; i++) {
  64     for(d=0; d<DIM; d++)
  65       if (bDim[d]) {
  66         if (index)
  67           fprintf(fp,"\t%g",x[index[i]][d]);
  68         else
  69           fprintf(fp,"\t%g",x[i][d]);
  70       }
  71     if (bDim[DIM]) {
  72       if (index)
  73           fprintf(fp,"\t%g",norm(x[index[i]]));
  74         else
  75           fprintf(fp,"\t%g",norm(x[i]));
  76     }
  77   }
  78   fprintf(fp,"\n");
  79 }
  80
  81 static void average_data(rvec x[],rvec xav[],real *mass,
  82                          int ngrps,int isize[],atom_id **index)
  83 {
  84   int  g,i,ind,d;
  85   real m,mtot;
  86   rvec tmp;
  87   double sum[DIM];
  88
  89   for(g=0; g<ngrps; g++) {
  90     for(d=0; d<DIM; d++)
  91       sum[d] = 0;
  92     clear_rvec(xav[g]);
  93     mtot = 0;
  94     for(i=0; i<isize[g]; i++) {
  95       ind = index[g][i];
  96       if (mass) {
  97         m = mass[ind];
  98         svmul(m,x[ind],tmp);
  99         for(d=0; d<DIM; d++)
 100           sum[d] += tmp[d];
 101         mtot += m;
 102       } else
 103         for(d=0; d<DIM; d++)
 104           sum[d] += x[ind][d];
 105     }
 106     if (mass) {
 107       for(d=0; d<DIM; d++)
 108         xav[g][d] = sum[d]/mtot;
 109     } else {
 110       for(d=0; d<DIM; d++)
 111         xav[g][d] = sum[d]/isize[g];
 112     }
 113   }
 114 }
 115
 116 static void print_data(FILE *fp,real time,rvec x[],real *mass,bool bCom,
 117                        int ngrps,int isize[],atom_id **index,bool bDim[])
 118 {
 119   static rvec *xav=NULL;
 120
 121   if (bCom) {
 122     if (xav==NULL)
 123       snew(xav,ngrps);
 124     average_data(x,xav,mass,ngrps,isize,index);
 125     low_print_data(fp,time,xav,ngrps,NULL,bDim);
 126   } else
 127     low_print_data(fp,time,x,isize[0],index[0],bDim);
 128 }
 129
 130 static void make_legend(FILE *fp,int ngrps,int isize,atom_id index[],
 131                         char **name,bool bCom,bool bMol,bool bDim[])
 132 {
 133   char **leg;
 134   char *dimtxt[] = { " X", " Y", " Z", "" };
 135   int n,i,j,d;
 136
 137   if (bCom)
 138     n = ngrps;
 139   else
 140     n = isize;
 141
 142   snew(leg,4*n);
 143   j=0;
 144   for(i=0; i<n; i++)
 145     for(d=0; d<=DIM; d++)
 146       if (bDim[d]) {
 147         snew(leg[j],STRLEN);
 148         if (bMol)
 149           sprintf(leg[j],"mol %d%s",index[i]+1,dimtxt[d]);
 150         else if (bCom)
 151           sprintf(leg[j],"%s%s",name[i],dimtxt[d]);
 152         else
 153           sprintf(leg[j],"atom %d%s",index[i]+1,dimtxt[d]);
 154         j++;
 155       }
 156   xvgr_legend(fp,j,leg);
 157   for(i=0; i<j; i++)
 158     sfree(leg[i]);
 159   sfree(leg);
 160 }
 161
 162 static real ekrot(rvec x[],rvec v[],real mass[],int isize,atom_id index[])
 163 {
 164   static real **TCM=NULL,**L;
 165   real   tm,m0,lxx,lxy,lxz,lyy,lyz,lzz,ekrot;
 166   rvec   dx,a0,ocm;
 167   rvec   xcm,vcm,acm;
 168   int    i,j,m,n;
 169
 170   if (TCM == NULL) {
 171     snew(TCM,DIM);
 172     for(i=0; i<DIM; i++)
 173       snew(TCM[i],DIM);
 174     snew(L,DIM);
 175     for(i=0; i<DIM; i++)
 176       snew(L[i],DIM);
 177   }
 178
 179   clear_rvec(xcm);
 180   clear_rvec(vcm);
 181   clear_rvec(acm);
 182   tm=0.0;
 183   for(i=0; i<isize; i++) {
 184     j = index[i];
 185     m0=mass[j];
 186     tm+=m0;
 187     oprod(x[j],v[j],a0);
 188     for(m=0; (m<DIM); m++) {
 189       xcm[m]+=m0*x[j][m]; /* c.o.m. position */
 190       vcm[m]+=m0*v[j][m]; /* c.o.m. velocity */
 191       acm[m]+=m0*a0[m];   /* rotational velocity around c.o.m. */
 192     }
 193   }
 194   oprod(xcm,vcm,a0);
 195   for(m=0; (m<DIM); m++) {
 196     xcm[m]/=tm;
 197     vcm[m]/=tm;
 198     acm[m]-=a0[m]/tm;
 199   }
 200
 201   lxx=lxy=lxz=lyy=lyz=lzz=0.0;
 202   for(i=0; i<isize; i++) {
 203     j = index[i];
 204     m0=mass[j];
 205     for(m=0; (m<DIM); m++)
 206       dx[m]=x[j][m]-xcm[m];
 207     lxx+=dx[XX]*dx[XX]*m0;
 208     lxy+=dx[XX]*dx[YY]*m0;
 209     lxz+=dx[XX]*dx[ZZ]*m0;
 210     lyy+=dx[YY]*dx[YY]*m0;
 211     lyz+=dx[YY]*dx[ZZ]*m0;
 212     lzz+=dx[ZZ]*dx[ZZ]*m0;
 213   }
 214
 215   L[XX][XX]=lyy+lzz;
 216   L[YY][XX]=-lxy;
 217   L[ZZ][XX]=-lxz;
 218   L[XX][YY]=-lxy;
 219   L[YY][YY]=lxx+lzz;
 220   L[ZZ][YY]=-lyz;
 221   L[XX][ZZ]=-lxz;
 222   L[YY][ZZ]=-lyz;
 223   L[ZZ][ZZ]=lxx+lyy;
 224
 225   m_inv_gen(L,DIM,TCM);
 226
 227   /* Compute omega (hoeksnelheid) */
 228   clear_rvec(ocm);
 229   ekrot=0;
 230   for(m=0; (m<DIM); m++) {
 231     for(n=0; (n<DIM); n++)
 232       ocm[m]+=TCM[m][n]*acm[n];
 233     ekrot+=0.5*ocm[m]*acm[m];
 234   }
 235   return ekrot;
 236 }
 237
 238 static real ektrans(rvec v[],real mass[],int isize,atom_id index[])
 239 {
 240   rvec   mvcom;
 241   real   mtot=0;
 242   int    i,j,d;
 243
 244   clear_rvec(mvcom);
 245   for(i=0; i<isize; i++) {
 246     j = index[i];
 247     for(d=0; d<DIM; d++)
 248       mvcom[d] += mass[j]*v[j][d];
 249     mtot += mass[j];
 250   }
 251
 252   return norm2(mvcom)/(mtot*2);
 253 }
 254
 255 static real temp(rvec v[],real mass[],int isize,atom_id index[])
 256 {
 257   real ekin2=0;
 258   int  i,j;
 259
 260   for(i=0; i<isize; i++) {
 261     j = index[i];
 262     ekin2 += mass[j]*norm2(v[j]);
 263   }
 264
 265   return ekin2/(3*isize*BOLTZ);
 266 }
 267
 268 static void remove_jump(matrix box,int natoms,rvec xp[],rvec x[])
 269 {
 270   rvec hbox;
 271   int d,i,m;
 272
 273   for(d=0; d<DIM; d++)
 274     hbox[d] = 0.5*box[d][d];
 275   for(i=0; i<natoms; i++)
 276     for(m=DIM-1; m>=0; m--) {
 277       while (x[i][m]-xp[i][m] <= -hbox[m])
 278         for(d=0; d<=m; d++)
 279           x[i][d] += box[m][d];
 280       while (x[i][m]-xp[i][m] > hbox[m])
 281         for(d=0; d<=m; d++)
 282           x[i][d] -= box[m][d];
 283     }
 284 }
 285
 286 static void write_pdb_bfac(char *fname,char *title,t_atoms *atoms,matrix box,
 287                            int isize,atom_id *index,int nfr,rvec *x,rvec *sum,
 288                            bool bDim[],real scale_factor)
 289 {
 290   FILE    *fp;
 291   real    max,len2,scale;
 292   atom_id maxi;
 293   int     i,m,onedim;
 294   bool    bOne;
 295
 296   if (nfr == 0) {
 297     fprintf(stderr,"No frames found for %s, will not write %s\n",title,fname);
 298   } else {
 299     fprintf(stderr,"Used %d frames for %s\n",nfr,title);
 300     onedim = -1;
 301     if (!bDim[DIM]) {
 302       m = 0;
 303       for(i=0; i<DIM; i++)
 304         if (bDim[i]) {
 305           onedim = i;
 306           m++;
 307         }
 308       if (m != 1)
 309         onedim = -1;
 310     }
 311     scale = 1.0/nfr;
 312     for(i=0; i<isize; i++)
 313       svmul(scale,x[index[i]],x[index[i]]);
 314     max  = 0;
 315     maxi = 0;
 316     for(i=0; i<isize; i++) {
 317       len2 = 0;
 318       for(m=0; m<DIM; m++)
 319         if (bDim[m] || bDim[DIM])
 320           len2 += sqr(sum[index[i]][m]);
 321       if (len2 > max) {
 322         max  = len2;
 323         maxi = index[i];
 324       }
 325     }
 326     if (scale_factor != 0) {
 327       scale = scale_factor;
 328     } else {
 329       if (max == 0)
 330         scale = 1;
 331       else
 332         scale = 10.0/sqrt(max);
 333     }
 334
 335     fprintf(stdout,"Maximum %s is %g on atom %d %s, res. %s %d\n",
 336             title,sqrt(max)/nfr,maxi+1,*(atoms->atomname[maxi]),
 337             *(atoms->resname[atoms->atom[maxi].resnr]),
 338             atoms->atom[maxi].resnr+1);
 339
 340     if (atoms->pdbinfo == NULL)
 341       snew(atoms->pdbinfo,atoms->nr);
 342     if (onedim == -1) {
 343       for(i=0; i<isize; i++) {
 344         len2 = 0;
 345         for(m=0; m<DIM; m++)
 346           if (bDim[m] || bDim[DIM])
 347             len2 += sqr(sum[index[i]][m]);
 348         atoms->pdbinfo[index[i]].bfac = sqrt(len2)*scale;
 349       }
 350     } else {
 351       for(i=0; i<isize; i++)
 352         atoms->pdbinfo[index[i]].bfac = sum[index[i]][onedim]*scale;
 353     }
 354     write_sto_conf_indexed(fname,title,atoms,x,NULL,box,isize,index);
 355   }
 356 }
 357
 358 int gmx_traj(int argc,char *argv[])
 359 {
 360   static char *desc[] = {
 361     "g_traj plots coordinates, velocities, forces and/or the box.",
 362     "With [TT]-com[tt] the coordinates, velocities and forces are",
 363     "calculated for the center of mass of each group.",
 364     "When [TT]-mol[tt] is set, the numbers in the index file are",
 365     "interpreted as molecule numbers and the same procedure as with",
 366     "[TT]-com[tt] is used for each molecule.[PAR]",
 367     "Option [TT]-ot[tt] plots the temperature of each group,",
 368     "provided velocities are present in the trajectory file.",
 369     "No corrections are made for constrained degrees of freedom!",
 370     "This implies [TT]-com[tt].[PAR]",
 371     "Options [TT]-ekt[tt] and [TT]-ekr[tt] plot the translational and",
 372     "rotational kinetic energy of each group,",
 373     "provided velocities are present in the trajectory file.",
 374     "This implies [TT]-com[tt].[PAR]",
 375     "Options [TT]-cv[tt] and [TT]-cf[tt] write the average velocities",
 376     "and average forces as temperature factors to a pdb file with",
 377     "the average coordinates. The temperature factors are scaled such",
 378     "that the maximum is 10. The scaling can be changed with the option",
 379     "[TT]-scale[tt]. To get the velocities or forces of one",
 380     "frame set both [TT]-b[tt] and [TT]-e[tt] to the time of",
 381     "desired frame. When averaging over frames you might need to use",
 382     "the [TT]-nojump[tt] option to obtain the correct average coordinates."
 383   };
 384   static bool bMol=FALSE,bCom=FALSE,bNoJump=FALSE;
 385   static bool bX=TRUE,bY=TRUE,bZ=TRUE,bNorm=FALSE;
 386   static real scale=0;
 387   t_pargs pa[] = {
 388     { "-com", FALSE, etBOOL, {&bCom},
 389       "Plot data for the com of each group" },
 390     { "-mol", FALSE, etBOOL, {&bMol},
 391       "Index contains molecule numbers iso atom numbers" },
 392     { "-nojump", FALSE, etBOOL, {&bNoJump},
 393       "Remove jumps of atoms across the box" },
 394     { "-x", FALSE, etBOOL, {&bX},
 395       "Plot X-component" },
 396     { "-y", FALSE, etBOOL, {&bY},
 397       "Plot Y-component" },
 398     { "-z", FALSE, etBOOL, {&bZ},
 399       "Plot Z-component" },
 400     { "-len", FALSE, etBOOL, {&bNorm},
 401       "Plot vector length" },
 402     { "-scale", FALSE, etREAL, {&scale},
 403       "Scale factor for pdb output, 0 is autoscale" }
 404
 405   };
 406   FILE       *outx=NULL,*outv=NULL,*outf=NULL,*outb=NULL,*outt=NULL;
 407   FILE       *outekt=NULL,*outekr=NULL;
 408   t_topology top;
 409   real       *mass,time;
 410   char       title[STRLEN],*indexfn;
 411   t_trxframe fr;
 412   int        flags;
 413   rvec       *xtop,*xp=NULL;
 414   rvec       *sumxv=NULL,*sumv=NULL,*sumxf=NULL,*sumf=NULL;
 415   matrix     topbox;
 416   int        status;
 417   int        i,j,n;
 418   int        ngrps,nr_vfr,nr_ffr;
 419   char       **grpname;
 420   int        *isize0,*isize;
 421   atom_id    **index0,**index;
 422   atom_id    *a,*atndx;
 423   t_block    *mols;
 424   bool       bTop,bOX,bOV,bOF,bOB,bOT,bEKT,bEKR,bCV,bCF,bDim[4],bDum[4];
 425   char       *box_leg[6] = { "XX", "YY", "ZZ", "YX", "ZX", "ZY" };
 426
 427   t_filenm fnm[] = {
 428     { efTRX, "-f", NULL, ffREAD },
 429     { efTPS, NULL, NULL, ffREAD },
 430     { efNDX, NULL, NULL, ffOPTRD },
 431     { efXVG, "-ox", "coord.xvg", ffOPTWR },
 432     { efXVG, "-ov", "veloc.xvg", ffOPTWR },
 433     { efXVG, "-of", "force.xvg", ffOPTWR },
 434     { efXVG, "-ob", "box.xvg",   ffOPTWR },
 435     { efXVG, "-ot", "temp.xvg",  ffOPTWR },
 436     { efXVG, "-ekt","ektrans.xvg", ffOPTWR },
 437     { efXVG, "-ekr","ekrot.xvg", ffOPTWR },
 438     { efPDB, "-cv", "veloc.pdb", ffOPTWR },
 439     { efPDB, "-cf", "force.pdb", ffOPTWR }
 440   };
 441 #define NFILE asize(fnm)
 442
 443   CopyRight(stderr,argv[0]);
 444   parse_common_args(&argc,argv,
 445                     PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE,
 446                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 447
 448   if (bMol)
 449     fprintf(stderr,"Interpreting indexfile entries as molecules.\n"
 450             "Using center of mass.\n");
 451
 452   bOX  = opt2bSet("-ox",NFILE,fnm);
 453   bOV  = opt2bSet("-ov",NFILE,fnm);
 454   bOF  = opt2bSet("-of",NFILE,fnm);
 455   bOB  = opt2bSet("-ob",NFILE,fnm);
 456   bOT  = opt2bSet("-ot",NFILE,fnm);
 457   bEKT = opt2bSet("-ekt",NFILE,fnm);
 458   bEKR = opt2bSet("-ekr",NFILE,fnm);
 459   bCV  = opt2bSet("-cv",NFILE,fnm);
 460   bCF  = opt2bSet("-cf",NFILE,fnm);
 461   if (bMol || bOT || bEKT || bEKR)
 462     bCom = TRUE;
 463
 464   bDim[XX] = bX;
 465   bDim[YY] = bY;
 466   bDim[ZZ] = bZ;
 467   bDim[DIM] = bNorm;
 468
 469   bTop = read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&xtop,NULL,topbox,
 470                        bCom && (bOX || bOV || bOT || bEKT || bEKR));
 471   sfree(xtop);
 472   if ((bMol || bCV || bCF) && !bTop)
 473     fatal_error(0,"Need a run input file for option -mol, -cv or -cf");
 474
 475
 476   if (bMol)
 477     indexfn = ftp2fn(efNDX,NFILE,fnm);
 478   else
 479     indexfn = ftp2fn_null(efNDX,NFILE,fnm);
 480
 481   if (bCom && !bMol) {
 482     fprintf(stderr,"How many groups do you want to analyze? ");
 483     scanf("%d",&ngrps);
 484   } else
 485     ngrps = 1;
 486   snew(grpname,ngrps);
 487   snew(isize0,ngrps);
 488   snew(index0,ngrps);
 489   get_index(&(top.atoms),indexfn,ngrps,isize0,index0,grpname);
 490
 491   if (bMol) {
 492     mols=&(top.blocks[ebMOLS]);
 493     a = mols->a;
 494     atndx = mols->index;
 495     ngrps = isize0[0];
 496     snew(isize,ngrps);
 497     snew(index,ngrps);
 498     for (i=0; i<ngrps; i++) {
 499       isize[i] = atndx[index0[0][i]+1] - atndx[index0[0][i]];
 500       snew(index[i],isize[i]);
 501       for(j=0; j<isize[i]; j++)
 502         index[i][j] = a[atndx[index0[0][i]]+j];
 503     }
 504   } else {
 505     isize = isize0;
 506     index = index0;
 507   }
 508   if (bCom) {
 509     snew(mass,top.atoms.nr);
 510     for(i=0; i<top.atoms.nr; i++)
 511       mass[i] = top.atoms.atom[i].m;
 512   } else
 513     mass = NULL;
 514
 515   flags = 0;
 516   if (bOX) {
 517     flags = flags | TRX_READ_X;
 518     outx = xvgropen(opt2fn("-ox",NFILE,fnm),
 519                     bCom ? "Center of mass" : "Coordinate",
 520                     xvgr_tlabel(),"Coordinate (nm)");
 521     make_legend(outx,ngrps,isize0[0],index0[0],grpname,bCom,bMol,bDim);
 522   }
 523   if (bOV) {
 524     flags = flags | TRX_READ_V;
 525     outv = xvgropen(opt2fn("-ov",NFILE,fnm),
 526                     bCom ? "Center of mass velocity" : "Velocity",
 527                     xvgr_tlabel(),"Velocity (nm/ps)");
 528    make_legend(outv,ngrps,isize0[0],index0[0],grpname,bCom,bMol,bDim);
 529   }
 530   if (bOF) {
 531     flags = flags | TRX_READ_F;
 532     outf = xvgropen(opt2fn("-of",NFILE,fnm),"Force",
 533                     xvgr_tlabel(),"Force (kJ mol\\S-1\\N nm\\S-1\\N)");
 534     make_legend(outf,ngrps,isize0[0],index0[0],grpname,bCom,bMol,bDim);
 535   }
 536   if (bOB) {
 537     outb = xvgropen(opt2fn("-ob",NFILE,fnm),"Box vector elements",
 538                     xvgr_tlabel(),"(nm)");
 539
 540     xvgr_legend(outb,6,box_leg);
 541   }
 542   if (bOT) {
 543     bDum[XX] = FALSE;
 544     bDum[YY] = FALSE;
 545     bDum[ZZ] = FALSE;
 546     bDum[DIM] = TRUE;
 547     flags = flags | TRX_READ_V;
 548     outt = xvgropen(opt2fn("-ot",NFILE,fnm),"Temperature",xvgr_tlabel(),"(K)");
 549     make_legend(outt,ngrps,isize[0],index[0],grpname,bCom,bMol,bDum);
 550   }
 551   if (bEKT) {
 552     bDum[XX] = FALSE;
 553     bDum[YY] = FALSE;
 554     bDum[ZZ] = FALSE;
 555     bDum[DIM] = TRUE;
 556     flags = flags | TRX_READ_V;
 557     outekt = xvgropen(opt2fn("-ekt",NFILE,fnm),"Center of mass translation",
 558                       xvgr_tlabel(),"Energy (kJ mol\\S-1\\N)");
 559     make_legend(outekt,ngrps,isize[0],index[0],grpname,bCom,bMol,bDum);
 560   }
 561   if (bEKR) {
 562     bDum[XX] = FALSE;
 563     bDum[YY] = FALSE;
 564     bDum[ZZ] = FALSE;
 565     bDum[DIM] = TRUE;
 566     flags = flags | TRX_READ_X | TRX_READ_V;
 567     outekr = xvgropen(opt2fn("-ekr",NFILE,fnm),"Center of mass rotation",
 568                       xvgr_tlabel(),"Energy (kJ mol\\S-1\\N)");
 569     make_legend(outekr,ngrps,isize[0],index[0],grpname,bCom,bMol,bDum);
 570   }
 571   if (bCV)
 572     flags = flags | TRX_READ_X | TRX_READ_V;
 573   if (bCF)
 574     flags = flags | TRX_READ_X | TRX_READ_F;
 575   if (flags == 0 && !bOB) {
 576     fprintf(stderr,"Please select one or more output file options\n");
 577     exit(0);
 578   }
 579
 580   read_first_frame(&status,ftp2fn(efTRX,NFILE,fnm),&fr,flags);
 581
 582   if (bCV) {
 583     snew(sumxv,fr.natoms);
 584     snew(sumv,fr.natoms);
 585   }
 586   if (bCF) {
 587     snew(sumxf,fr.natoms);
 588     snew(sumf,fr.natoms);
 589   }
 590   nr_vfr = 0;
 591   nr_ffr = 0;
 592
 593   do {
 594     time = convert_time(fr.time);
 595
 596     if (fr.bX && bNoJump && fr.bBox) {
 597       if (xp)
 598         remove_jump(fr.box,fr.natoms,xp,fr.x);
 599       else
 600         snew(xp,fr.natoms);
 601       for(i=0; i<fr.natoms; i++)
 602         copy_rvec(fr.x[i],xp[i]);
 603     }
 604
 605     if (fr.bX && bCom)
 606       rm_pbc(&(top.idef),fr.natoms,fr.box,fr.x,fr.x);
 607
 608     if (bOX && fr.bX)
 609       print_data(outx,time,fr.x,mass,bCom,ngrps,isize,index,bDim);
 610     if (bOV && fr.bV)
 611       print_data(outv,time,fr.v,mass,bCom,ngrps,isize,index,bDim);
 612     if (bOF && fr.bF)
 613       print_data(outf,time,fr.f,NULL,bCom,ngrps,isize,index,bDim);
 614     if (bOB && fr.bBox)
 615       fprintf(outb,"\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n",fr.time,
 616               fr.box[XX][XX],fr.box[YY][YY],fr.box[ZZ][ZZ],
 617               fr.box[YY][XX],fr.box[ZZ][XX],fr.box[ZZ][YY]);
 618     if (bOT && fr.bV) {
 619       fprintf(outt," %g",time);
 620       for(i=0; i<ngrps; i++)
 621         fprintf(outt,"\t%g",temp(fr.v,mass,isize[i],index[i]));
 622       fprintf(outt,"\n");
 623     }
 624     if (bEKT && fr.bV) {
 625       fprintf(outekt," %g",time);
 626       for(i=0; i<ngrps; i++)
 627         fprintf(outekt,"\t%g",ektrans(fr.v,mass,isize[i],index[i]));
 628       fprintf(outekt,"\n");
 629     }
 630     if (bEKR && fr.bX && fr.bV) {
 631       fprintf(outekr," %g",time);
 632       for(i=0; i<ngrps; i++)
 633         fprintf(outekr,"\t%g",ekrot(fr.x,fr.v,mass,isize[i],index[i]));
 634       fprintf(outekr,"\n");
 635     }
 636     if (bCV && fr.bX && fr.bV) {
 637       for(i=0; i<fr.natoms; i++)
 638         rvec_inc(sumxv[i],fr.x[i]);
 639       for(i=0; i<fr.natoms; i++)
 640         rvec_inc(sumv[i],fr.v[i]);
 641       nr_vfr++;
 642     }
 643     if (bCF && fr.bX && fr.bF) {
 644       for(i=0; i<fr.natoms; i++)
 645         rvec_inc(sumxf[i],fr.x[i]);
 646       for(i=0; i<fr.natoms; i++)
 647         rvec_inc(sumf[i],fr.f[i]);
 648       nr_ffr++;
 649     }
 650
 651   } while(read_next_frame(status,&fr));
 652
 653
 654   /* clean up a bit */
 655   close_trj(status);
 656
 657   if (bOX) fclose(outx);
 658   if (bOV) fclose(outv);
 659   if (bOF) fclose(outf);
 660   if (bOB) fclose(outb);
 661   if (bOT) fclose(outt);
 662   if (bEKT) fclose(outekt);
 663   if (bEKR) fclose(outekr);
 664
 665   if ((bCV || bCF) && (nr_vfr>1 || nr_ffr>1) && !bNoJump)
 666     fprintf(stderr,"WARNING: More than one frame was used for option -cv or -cf\n"
 667             "If atoms jump across the box you should use the -nojump option\n");
 668
 669   if (bCV)
 670     write_pdb_bfac(opt2fn("-cv",NFILE,fnm),"average velocity",&(top.atoms),
 671                    topbox,isize[0],index[0],nr_vfr,sumxv,sumv,bDim,scale);
 672   if (bCF)
 673     write_pdb_bfac(opt2fn("-cf",NFILE,fnm),"average force",&(top.atoms),
 674                    topbox,isize[0],index[0],nr_ffr,sumxf,sumf,bDim,scale);
 675
 676   /* view it */
 677   view_all(NFILE, fnm);
 678
 679   thanx(stderr);
 680
 681   return 0;
 682 }
 683