src/tools/trjconv.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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_trjconv_c = "$Id$";
  30
  31 #include <string.h>
  32 #include <math.h>
  33 #include <unistd.h>
  34 #include "macros.h"
  35 #include "sysstuff.h"
  36 #include "smalloc.h"
  37 #include "typedefs.h"
  38 #include "copyrite.h"
  39 #include "gmxfio.h"
  40 #include "tpxio.h"
  41 #include "trnio.h"
  42 #include "statutil.h"
  43 #include "futil.h"
  44 #include "pdbio.h"
  45 #include "confio.h"
  46 #include "names.h"
  47 #include "rdgroup.h"
  48 #include "vec.h"
  49 #include "xtcio.h"
  50 #include "do_fit.h"
  51 #include "rmpbc.h"
  52 #include "wgms.h"
  53 #include "magic.h"
  54 #include "binio.h"
  55 #include "pbc.h"
  56
  57 static void center_x(rvec x[],matrix box,int n,atom_id index[])
  58 {
  59   int i,m,ai;
  60   rvec cmin,cmax,dx;
  61
  62   if (n==0) return;
  63   copy_rvec(x[index[0]],cmin);
  64   copy_rvec(x[index[0]],cmax);
  65   for(i=0; (i<n); i++) {
  66     ai=index[i];
  67     for(m=0; (m<DIM); m++) {
  68       if (x[ai][m] < cmin[m])
  69         cmin[m]=x[ai][m];
  70       else if (x[ai][m] > cmax[m])
  71         cmax[m]=x[ai][m];
  72     }
  73   }
  74   for(m=0; (m<DIM); m++) {
  75     dx[m]=-(box[m][m]-(cmin[m]+cmax[m]))*0.5;
  76   }
  77   for(i=0; (i<n); i++) {
  78     ai=index[i];
  79     rvec_dec(x[ai],dx);
  80   }
  81 }
  82
  83
  84 void check_trn(char *fn)
  85 {
  86   if ((fn2ftp(fn) != efTRJ)  && (fn2ftp(fn) != efTRR))
  87     fatal_error(0,"%s is not a trj file, exiting\n",fn);
  88 }
  89
  90 #ifndef _win_
  91 void do_trunc(char *fn, real t0)
  92 {
  93   int          in;
  94   FILE         *fp;
  95   bool         bStop,bOK;
  96   t_trnheader  sh;
  97   long         fpos;
  98   char         yesno[256];
  99   int          j;
 100   real         t=0;
 101
 102   if (t0 == -1)
 103     fatal_error(0,"You forgot to set the truncation time");
 104
 105   /* Check whether this is a .trj file */
 106   check_trn(fn);
 107
 108   in   = open_trn(fn,"r");
 109   fp   = fio_getfp(in);
 110   if (fp == NULL) {
 111     fprintf(stderr,"Sorry, can not trunc %s, truncation of this filetype is not supported\n",fn);
 112     close_trn(in);
 113   } else {
 114     j    = 0;
 115     fpos = fio_ftell(in);
 116     bStop= FALSE;
 117     while (!bStop && fread_trnheader(in,&sh,&bOK)) {
 118       fread_htrn(in,&sh,NULL,NULL,NULL,NULL);
 119       fpos=ftell(fp);
 120       t=sh.t;
 121       if (t>=t0) {
 122         fseek(fp,fpos,SEEK_SET);
 123         bStop=TRUE;
 124       }
 125     }
 126     if (bStop) {
 127       fprintf(stderr,"Do you REALLY want to truncate this trajectory (%s) at:\n"
 128               "frame %d, time %g, bytes %ld ??? (type YES if so)\n",
 129               fn,j,t,fpos);
 130       scanf("%s",yesno);
 131       if (strcmp(yesno,"YES") == 0) {
 132         fprintf(stderr,"Once again, I'm gonna DO this...\n");
 133         close_trn(in);
 134         truncate(fn,fpos);
 135       }
 136       else {
 137         fprintf(stderr,"Ok, I'll forget about it\n");
 138       }
 139     }
 140     else {
 141       fprintf(stderr,"Already at end of file (t=%g)...\n",t);
 142       close_trn(in);
 143     }
 144   }
 145 }
 146 #endif
 147
 148 bool bRmod(double a,double b)
 149 {
 150   int iq;
 151   double tol = 1e-6;
 152
 153   iq = ((1.0+tol)*a)/b;
 154
 155   if (fabs(a-b*iq) <= tol*a)
 156     return TRUE;
 157   else
 158     return FALSE;
 159 }
 160
 161 int main(int argc,char *argv[])
 162 {
 163   static char *desc[] = {
 164     "trjconv can convert trajectory files in many ways:[BR]",
 165     "[BB]1.[bb] from one format to another[BR]",
 166     "[BB]2.[bb] select a subset of atoms[BR]",
 167     "[BB]3.[bb] remove periodicity from molecules[BR]",
 168     "[BB]4.[bb] keep multimeric molecules together[BR]",
 169     "[BB]5.[bb] center atoms in the box[BR]",
 170     "[BB]6.[bb] fit atoms to reference structure[BR]",
 171     "[BB]7.[bb] remove duplicate frames[BR]",
 172     "[BB]8.[bb] reduce the number of frames[BR]",
 173     "[BB]9.[bb] change the timestamps of the frames (e.g. t0 and delta-t)",
 174     "[PAR]",
 175     "The program [TT]trjcat[tt] can concatenate multiple trajectory files.",
 176     "[PAR]",
 177     "Currently seven formats are supported for input and output:",
 178     "[TT].xtc[tt], [TT].trr[tt], [TT].trj[tt], [TT].gro[tt], [TT].g96[tt],",
 179     "[TT].pdb[tt] and [TT].g87[tt].",
 180     "The file formats are detected from the file extension.",
 181     "For [TT].gro[tt] and [TT].xtc[tt] files the output precision ",
 182     "can be given as a number of ",
 183     "decimal places. Note that velocities are only supported in ",
 184     "[TT].trr[tt], [TT].trj[tt], [TT].gro[tt] and [TT].g96[tt] files.[PAR]",
 185     "The option [TT]-app[tt] can be used to",
 186     "append output to an existing trajectory file.",
 187     "No checks are performed to ensure integrity",
 188     "of the resulting combined trajectory file.",
 189     "[TT].pdb[tt] files with all frames concatenated can be viewed with",
 190     "[TT]rasmol -nmrpdb[tt].[PAR]",
 191     "It is possible to select part of your trajectory and write it out",
 192     "to a new trajectory file in order to save disk space, e.g. for leaving",
 193     "out the water from a trajectory of a protein in water.",
 194     "[BB]ALWAYS[bb] put the original trajectory on tape!",
 195     "We recommend to use the portable [TT].xtc[tt] format for your analysis",
 196     "to save disk space and to have portable files.[PAR]",
 197     "There are two options for fitting the trajectory to a reference",
 198     "either for essential dynamics analysis or for whatever.",
 199     "The first option is just plain fitting to a reference structure",
 200     "in the structure file, the second option is a progressive fit",
 201     "in which the first timeframe is fitted to the reference structure ",
 202     "in the structure file to obtain and each subsequent timeframe is ",
 203     "fitted to the previously fitted structure. This way a continuous",
 204     "trajectory is generated, which might not be the case when using the",
 205     "regular fit method, e.g. when your protein undergoes large",
 206     "conformational transitions.[PAR]",
 207     "The option [TT]-pbc[tt] sets the type of periodic boundary condition",
 208     "treatment. [TT]whole[tt] makes broken molecules whole (a run input",
 209     "file is required). [TT]-pbc[tt] is changed form [TT]none[tt] to",
 210     "[TT]whole[tt] when [TT]-fit[tt] or [TT]-pfit[tt] is set.",
 211     "[TT]inbox[tt] puts all the atoms in the box.",
 212     "[TT]nojump[tt] checks if atoms jump across the box and then puts",
 213     "them back. This has the effect that all molecules",
 214     "will remain whole (provided they were whole in the initial",
 215     "conformation), note that this ensures a continuous trajectory but",
 216     "molecules may diffuse out of the box. The starting configuration",
 217     "for this procedure is taken from the structure file, if one is",
 218     "supplied, otherwise it is the first frame.",
 219     "Use [TT]-center[tt] to put the system in the center of the box.",
 220     "This is especially useful for multimeric proteins, since this",
 221     "procedure will ensure the subunits stay together in the trajectory",
 222     "(due to PBC, they might be separated), providing they were together",
 223     "in the initial conformation.[PAR]",
 224     "With the option [TT]-dt[tt] it is possible to reduce the number of ",
 225     "frames in the output. This option relies on the accuracy of the times",
 226     "in your input trajectory, so if these are inaccurate use the",
 227     "[TT]-timestep[tt]",
 228     "option to modify the time (this can be done simultaneously).[PAR]",
 229     "Using [TT]-trunc[tt] trjconv can truncate [TT].trj[tt] in place, i.e.",
 230     "without copying the file. This is useful when a run has crashed",
 231     "during disk I/O (one more disk full), or when two contiguous",
 232     "trajectories must be concatenated without have double frames.[PAR]",
 233     "Also the option [TT]-checkdouble[tt] may be used to remove all",
 234     "duplicate frames from such a concatenated trajectory, this is done",
 235     "by ignoring all frames with a time smaller than or equal to the previous",
 236     "frame. [TT]trjcat[tt] is more suitable for concatenating trajectory",
 237     "files.[PAR]",
 238     "The option [TT]-dump[tt] can be used to extract a frame at or near",
 239     "one specific time from your trajectory."
 240   };
 241
 242   static char *pbc_opt[] = { NULL, "none", "whole", "inbox", "nojump", NULL };
 243
 244   static bool  bAppend=FALSE,bSeparate=FALSE,bVels=TRUE;
 245   static bool  bCenter=FALSE,bFit=FALSE,bPFit=FALSE,bBox=TRUE;
 246   static bool  bCheckDouble=FALSE;
 247   static int   skip_nr=1,prec=3;
 248   static real  tzero=0.0,delta_t=0.0,timestep=0.0,ttrunc=-1,tdump=-1;
 249   static rvec  newbox = {0,0,0}, shift = {0,0,0};
 250   static char  *exec_command=NULL;
 251
 252   t_pargs pa[] = {
 253     { "-pbc", FALSE,  etENUM, {pbc_opt},
 254       "PBC treatment" },
 255     { "-center", FALSE,  etBOOL, {&bCenter},
 256       "Center atoms in box" },
 257     { "-box", FALSE, etRVEC, {&newbox},
 258       "Size for new cubic box (default: read from input)" },
 259     { "-shift", FALSE, etRVEC, {&shift},
 260       "All coordinates will be shifted by framenr*shift" },
 261     { "-fit", FALSE,  etBOOL, {&bFit},
 262       "Fit molecule to ref structure in the structure file" },
 263     { "-pfit", FALSE,  etBOOL, {&bPFit},
 264       "Progressive fit, to the previous fitted structure" },
 265     { "-prec", FALSE,  etINT,  {&prec},
 266       "Precision for .xtc and .gro writing in number of decimal places" },
 267     { "-vel", FALSE, etBOOL, {&bVels},
 268       "Read and write velocities if possible" },
 269     { "-skip", FALSE,  etINT, {&skip_nr},
 270       "Only write every nr-th frame" },
 271     { "-dt", FALSE,  etREAL, {&delta_t},
 272       "Only write frame when t MOD dt = first time" },
 273     { "-t0", FALSE,  etREAL, {&tzero},
 274       "Starting time for trajectory"
 275       "(default: don't change)"},
 276 #ifndef _win_
 277     { "-trunc", FALSE, etREAL, {&ttrunc},
 278       "Truncate input trj file after this amount of ps" },
 279 #endif
 280     { "-dump", FALSE, etREAL, {&tdump},
 281       "Dump frame nearest specified time" },
 282     { "-g87box", FALSE,  etBOOL, {&bBox},
 283       "Write a box for .g87" },
 284     { "-exec", FALSE,  etSTR, {&exec_command},
 285       "Execute command for every output frame with the frame number "
 286       "as argument" },
 287     { "-timestep", FALSE,  etREAL, {&timestep},
 288       "Change time step between frames" },
 289     { "-app", FALSE,  etBOOL, {&bAppend},
 290       "Append output"},
 291     { "-sep", FALSE,  etBOOL, {&bSeparate},
 292       "Write each frame to a separate .gro or .pdb file"},
 293     { "-checkdouble", FALSE, etBOOL, {&bCheckDouble},
 294       "Only write frames with time larger than previous frame" }
 295   };
 296
 297   FILE         *out=NULL;
 298   int          trjout=0;
 299   int          status,ftp,ftpin,file_nr;
 300   rvec         *x,*xn,*xout,*v,*vn,*vout;
 301   rvec         *xp,x_shift;
 302   real         xtcpr, lambda,*w_rls=NULL;
 303   matrix       box;
 304   int          m,i,d,frame,outframe,natoms=0,nout,nre,step;
 305 #define SKIP 10
 306   t_topology   top;
 307   t_atoms      *atoms=NULL,useatoms;
 308   atom_id      *index;
 309   char         *grpname;
 310   int          ifit,irms;
 311   atom_id      *ind_fit,*ind_rms;
 312   char         *gn_fit,*gn_rms;
 313   real         t,pt,tshift=0,t0=-1,dt=0.001;
 314   bool         bPBC,bInBox,bNoJump;
 315   bool         bCopy,bDoIt,bIndex,bTDump,bSetTime,bTPS=FALSE,bDTset=FALSE;
 316   bool         bExec,bTimeStep=FALSE,bDumpFrame=FALSE,bToldYouOnce=FALSE;
 317   bool         bHaveNextFrame,bHaveX=FALSE,bHaveV=FALSE,bSetBox;
 318   char         *grpnm;
 319   char         *top_file,*in_file,*out_file,out_file2[256],*charpt;
 320   char         top_title[256],title[256],command[256],filemode[5];
 321   int          xdr=0;
 322
 323   t_filenm fnm[] = {
 324     { efTRX, "-f",  NULL, ffREAD },
 325     { efTRX, "-o", "trajout", ffWRITE },
 326     { efTPS, NULL,  NULL, ffOPTRD },
 327     { efNDX, NULL,  NULL, ffOPTRD }
 328   };
 329 #define NFILE asize(fnm)
 330
 331   CopyRight(stderr,argv[0]);
 332   parse_common_args(&argc,argv,PCA_CAN_TIME,TRUE,
 333                     NFILE,fnm,asize(pa),pa,asize(desc),desc,
 334                     0,NULL);
 335
 336   top_file=ftp2fn(efTPS,NFILE,fnm);
 337
 338   /* Check command line */
 339   in_file=opt2fn("-f",NFILE,fnm);
 340   if (ttrunc != -1) {
 341 #ifndef _win_
 342     do_trunc(in_file,ttrunc);
 343 #endif
 344   }
 345   else {
 346     if (bPFit)
 347       bFit = TRUE;
 348     bSetBox   = opt2parg_bSet("-box", asize(pa), pa);
 349     bSetTime  = opt2parg_bSet("-t0", asize(pa), pa);
 350     bExec     = opt2parg_bSet("-exec", asize(pa), pa);
 351     bTimeStep = opt2parg_bSet("-timestep", asize(pa), pa);
 352     bTDump    = opt2parg_bSet("-dump", asize(pa), pa);
 353     bPBC      = (strcmp(pbc_opt[0],"whole") == 0);
 354     bInBox    = (strcmp(pbc_opt[0],"inbox") == 0);
 355     bNoJump   = (strcmp(pbc_opt[0],"nojump") == 0);
 356     if (bFit && (strcmp(pbc_opt[0],"none") == 0))
 357       bPBC = TRUE;
 358
 359     /* prec is in nr of decimal places, xtcprec is a multiplication factor: */
 360     xtcpr=1;
 361     for (i=0; i<prec; i++)
 362       xtcpr*=10;
 363
 364     bIndex=ftp2bSet(efNDX,NFILE,fnm);
 365
 366     /* Determine output type */
 367     out_file=opt2fn("-o",NFILE,fnm);
 368     ftp=fn2ftp(out_file);
 369     fprintf(stderr,"Will write %s: %s\n",ftp2ext(ftp),ftp2desc(ftp));
 370     if (bVels) {
 371       /* check if velocities are possible in input and output files */
 372       ftpin=fn2ftp(in_file);
 373       bVels= ((ftp==efTRR) ||(ftp==efTRJ) || (ftp==efGRO))
 374         && ((ftpin==efTRR) ||(ftpin==efTRJ) || (ftpin==efGRO));
 375     }
 376     if (!bVels) {
 377       bHaveX=TRUE;
 378       bHaveV=FALSE;
 379     }
 380
 381     /* skipping */
 382     if (skip_nr <= 0) {
 383       fprintf(stderr,"The number of frames to skip is <= 0. "
 384               "Writing out all frames.\n\n");
 385       skip_nr = 1;
 386     }
 387
 388     /* Determine whether to read a topology */
 389     bTPS = (ftp2bSet(efTPS,NFILE,fnm) ||
 390             bPBC || bFit || (ftp == efGRO) || (ftp == efPDB));
 391
 392     /* Determine if when can read index groups */
 393     bIndex = (bIndex || bTPS);
 394
 395     if (bTPS) {
 396       read_tps_conf(top_file,top_title,&top,&xp,NULL,box,bFit);
 397       atoms=&top.atoms;
 398       /* top_title is only used for gro and pdb,
 399        * the header in such a file is top_title t= ...
 400        * to prevent a double t=, remove it from top_title
 401        */
 402       if ((charpt=strstr(top_title," t= ")))
 403           charpt[0]='\0';
 404     }
 405
 406     if (bFit) {
 407       printf("Select group for root least squares fit\n");
 408       get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
 409                 1,&ifit,&ind_fit,&gn_fit);
 410
 411       if (ifit < 3)
 412         fatal_error(0,"Need at least 3 points to fit!\n");
 413     }
 414
 415     if (bIndex) {
 416       printf("Select group for output\n");
 417       get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),
 418                 1,&nout,&index,&grpnm);
 419     }
 420     else {
 421       /* no index file, so read natoms from TRX */
 422       natoms=read_first_x(&status,in_file,&t,&x,box);
 423       close_trj(status);
 424       snew(index,natoms);
 425       for(i=0;i<natoms;i++)
 426         index[i]=i;
 427       nout=natoms;
 428     }
 429
 430     /* if xp was not snew-ed before, do it now */
 431     if (!xp)
 432       snew(xp, natoms);
 433
 434     if (bFit) {
 435       snew(w_rls,atoms->nr);
 436       for(i=0; (i<ifit); i++)
 437         w_rls[ind_fit[i]]=atoms->atom[ind_fit[i]].m;
 438
 439       /* Restore reference structure and set to origin,
 440          store original location (to put structure back) */
 441       rm_pbc(&(top.idef),atoms->nr,box,xp,xp);
 442       copy_rvec(xp[index[0]],x_shift);
 443       reset_x(ifit,ind_fit,nout,index,xp,w_rls);
 444       rvec_dec(x_shift,xp[index[0]]);
 445     }
 446
 447     /* Make atoms struct for output in GRO or PDB files */
 448     if ((ftp == efGRO) || ((ftp == efG96) && bTPS) || (ftp == efPDB)) {
 449       /* get memory for stuff to go in pdb file */
 450       init_t_atoms(&useatoms,atoms->nr,FALSE);
 451       sfree(useatoms.resname);
 452       useatoms.resname=atoms->resname;
 453       for(i=0;(i<nout);i++) {
 454         useatoms.atomname[i]=atoms->atomname[index[i]];
 455         useatoms.atom[i].resnr=atoms->atom[index[i]].resnr;
 456         useatoms.nres=max(useatoms.nres,useatoms.atom[i].resnr+1);
 457       }
 458       useatoms.nr=nout;
 459     }
 460     /* open trj file for reading */
 461     if (bVels)
 462       natoms=read_first_x_or_v(&status,in_file,&t,&x,&v,box);
 463     else
 464       natoms=read_first_x     (&status,in_file,&t,&x,   box);
 465     if (bSetTime)
 466       tshift=tzero-t;
 467     else
 468       tzero=t;
 469
 470     if (natoms == 0)
 471       fatal_error(0,"No atoms found in file %s",in_file);
 472
 473     /* open output for writing */
 474     if ((bAppend) && (fexist(out_file))) {
 475       strcpy(filemode,"a");
 476       fprintf(stderr,"APPENDING to existing file %s\n",out_file);
 477     } else
 478       strcpy(filemode,"w");
 479     switch (ftp) {
 480     case efXTC:
 481       xdr = open_xtc(out_file,filemode);
 482       break;
 483     case efG87:
 484       out=ffopen(out_file,filemode);
 485       break;
 486     case efTRR:
 487     case efTRJ:
 488       out=NULL;
 489       trjout = open_tpx(out_file,filemode);
 490       break;
 491     case efGRO:
 492     case efG96:
 493     case efPDB:
 494       if (!bSeparate)
 495         out=ffopen(out_file,filemode);
 496       break;
 497     }
 498
 499     bCopy=FALSE;
 500     if (bIndex)
 501       /* check if index is meaningful */
 502       for(i=0; i<nout; i++) {
 503         if (index[i] >= natoms)
 504           fatal_error(0,"Element %d of index group %s is %d, this is larger than the number of atoms in the trajectory file (%d)",i,grpnm,index[i]+1,natoms);
 505         bCopy = bCopy || (i != index[i]);
 506       }
 507
 508     if (bCopy) {
 509       snew(xn,nout);
 510       xout=xn;
 511       snew(vn,nout);
 512       vout=vn;
 513     } else {
 514       xout=x;
 515       vout=v;
 516     }
 517
 518     if (ftp == efG87)
 519       fprintf(out,"Generated by %s. #atoms=%d, %s BOX is stored in "
 520               "this file.\n",Program(),nout,bBox ? "a" : "NO");
 521
 522     /* Start the big loop over frames */
 523     file_nr  =  0;
 524     frame    =  0;
 525     outframe =  0;
 526     pt       = -666.0;
 527     bDTset   = FALSE;
 528
 529     do {
 530       /* generate new box */
 531       if (bSetBox) {
 532         clear_mat(box);
 533         for (m=0; m<DIM; m++)
 534           box[m][m] = newbox[m];
 535       }
 536
 537       if (bTDump) {
 538         /* determine timestep */
 539         if (t0 == -1) {
 540           t0=t;
 541         } else {
 542           if (!bDTset) {
 543             dt=t-t0;
 544             bDTset=TRUE;
 545           }
 546         }
 547         bDumpFrame = (t >= tdump-(0.5*dt) ) && (t <= tdump+(0.5*dt) );
 548       }
 549
 550       /* determine if an atom jumped across the box and reset it if so */
 551       if (bNoJump && (bTPS || frame!=0)) {
 552         for(i=0; (i<natoms); i++)
 553           for(d=0; (d<DIM); d++)
 554             if ( x[i][d]-xp[i][d] > 0.5*box[d][d] )
 555               x[i][d] -= box[d][d];
 556             else if ( x[i][d]-xp[i][d] < -0.5*box[d][d] )
 557               x[i][d] += box[d][d];
 558       }
 559
 560       if (bPFit) {
 561         /* Now modify the coords according to the flags,
 562            for normal fit, this is only done for output frames */
 563         rm_pbc(&(top.idef),natoms,box,x,x);
 564
 565         reset_x(ifit,ind_fit,nout,index,x,w_rls);
 566         do_fit(natoms,w_rls,xp,x);
 567       }
 568
 569       /* store this set of coordinates for future use */
 570       if (bPFit || bNoJump) {
 571         for(i=0; (i<natoms); i++) {
 572           copy_rvec(x[i],xp[i]);
 573           rvec_inc(x[i],x_shift);
 574         }
 575       }
 576
 577       if ( bCheckDouble && (t<=pt) && !bToldYouOnce ) {
 578         fprintf(stderr,"\nRemoving duplicate frame(s) after t=%g "
 579                 "(t=%g, frame %d)\n",pt,t,frame);
 580         bToldYouOnce=TRUE;
 581       }
 582
 583       if ( ( ( !bTDump && (frame % skip_nr == 0) ) || bDumpFrame  ) &&
 584            ( !bCheckDouble || ( bCheckDouble && (t > pt) ) ) ) {
 585
 586         /* remember time from this frame */
 587         pt = t;
 588
 589         /* calc new time */
 590         if (bTimeStep)
 591           t=tzero+frame*timestep;
 592         else
 593           if (bSetTime)
 594             t=t+tshift;
 595
 596         if (bTDump)
 597           fprintf(stderr,"\nDumping frame at t= %g\n",t);
 598
 599         /* check for writing at each delta_t */
 600         bDoIt=(delta_t == 0);
 601         if (!bDoIt)
 602           bDoIt=bRmod(t-tzero, delta_t);
 603
 604         if (bDoIt || bTDump) {
 605           /* print sometimes */
 606           if (bToldYouOnce) {
 607             bToldYouOnce=FALSE;
 608             fprintf(stderr,"\nContinue writing frames from t=%g, frame=%d\n",
 609                     t,outframe);
 610           }
 611           if ( ((outframe % SKIP) == 0) || (outframe < SKIP) )
 612             fprintf(stderr," ->  frame %6d time %8.3f      \r",outframe,t);
 613
 614           if (!bPFit) {
 615             /* Now modify the coords according to the flags,
 616                for PFit we did this already! */
 617             if (bPBC)
 618               rm_pbc(&(top.idef),natoms,box,x,x);
 619
 620             if (bFit) {
 621               reset_x(ifit,ind_fit,nout,index,x,w_rls);
 622               do_fit(natoms,w_rls,xp,x);
 623               for(i=0; (i<natoms); i++)
 624                 rvec_inc(x[i],x_shift);
 625             }
 626           }
 627
 628           if (bCenter) {
 629             center_x(x,box,nout,index);
 630           }
 631
 632           if (bCopy) {
 633             for(i=0; (i<nout); i++) {
 634               copy_rvec(x[index[i]],xout[i]);
 635               if (bVels) copy_rvec(v[index[i]],vout[i]);
 636             }
 637           }
 638
 639           /* this should make sure all atoms in output are really inside
 640              a rectangular box. Was needed to make movies.
 641              Peter Tieleman, Mon Jul 15 1996
 642              */
 643           if (bInBox)
 644             put_atoms_in_box(nout,box,xout);
 645
 646           if (opt2parg_bSet("-shift",asize(pa),pa))
 647             for(i=0; (i<nout); i++)
 648               for (d=0; (d<DIM); d++)
 649                 xout[i][d] += outframe*shift[d];
 650
 651           if (bVels) {
 652           /* check if we have velocities and/or coordinates,
 653              don't ask me why you can have a frame w/o coords !? */
 654             bHaveV=FALSE;
 655             bHaveX=FALSE;
 656             for (i=0; (i<natoms); i++)
 657               for (d=0; (d<DIM); d++) {
 658                 bHaveV=bHaveV || v[i][d];
 659                 bHaveX=bHaveX || x[i][d];
 660               }
 661           }
 662
 663           switch(ftp) {
 664           case efTRJ:
 665           case efTRR:
 666             fwrite_trn(trjout,frame,t,0,box,
 667                        nout,
 668                        bHaveX ? xout : NULL,
 669                        bHaveV ? vout : NULL,
 670                        NULL);
 671
 672             break;
 673           case efG87:
 674             write_gms(out,nout,xout,bBox ? box : NULL );
 675             break;
 676           case efXTC:
 677             write_xtc(xdr,nout,frame,t,box,xout,xtcpr);
 678             break;
 679           case efGRO:
 680           case efG96:
 681           case efPDB:
 682             sprintf(title,"Generated by trjconv : %s t= %9.5f",top_title,t);
 683             if (bSeparate) {
 684               sprintf(out_file2,"%d_%s",file_nr,out_file);
 685               out=ffopen(out_file2,"w");
 686             }
 687             switch(ftp) {
 688             case efGRO:
 689               write_hconf_p(out,title,&useatoms,prec,xout,bHaveV?vout:NULL,box);
 690               break;
 691             case efPDB:
 692               fprintf(out,"REMARK    GENERATED BY TRJCONV\n");
 693               sprintf(title,"%s t= %9.5f",top_title,t);
 694               write_pdbfile(out,title,&useatoms,xout,box,0,TRUE);
 695               break;
 696             case efG96:
 697               if (bSeparate || bTDump)
 698                 write_g96_conf(out,title,&useatoms,xout,bHaveV?vout:NULL,box,
 699                                0,NULL);
 700               else {
 701                 if (outframe == 0)
 702                   fprintf(out,"TITLE\n%s\nEND\n",title);
 703                 fprintf(out,"TIMESTEP\n%9d%15.9f\nEND\nPOSITIONRED\n",
 704                         frame,t);
 705                 for(i=0; i<nout; i++)
 706                   fprintf(out,"%15.9f%15.9f%15.9f\n",
 707                           xout[i][XX],xout[i][YY],xout[i][ZZ]);
 708                 fprintf(out,"END\nBOX\n%15.9f%15.9f%15.9f\nEND\n",
 709                         box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
 710               }
 711               break;
 712             }
 713             if (bSeparate) {
 714               ffclose(out);
 715               out = NULL;
 716               file_nr++;
 717             }
 718             break;
 719           default:
 720             fatal_error(0,"DHE, ftp=%d\n",ftp);
 721           }
 722
 723           /* execute command */
 724           if (bExec) {
 725             char c[255];
 726             sprintf(c,"%s  %d",exec_command,file_nr-1);
 727             /*fprintf(stderr,"Executing '%s'\n",c);*/
 728             system(c);
 729           }
 730           outframe++;
 731         }
 732       }
 733       frame++;
 734       if (bVels) {
 735         bHaveNextFrame=read_next_x_or_v(status,&t,natoms,x,v,box);
 736       } else
 737         bHaveNextFrame=read_next_x     (status,&t,natoms,x,  box);
 738     } while (!bDumpFrame && bHaveNextFrame);
 739
 740     if ( bTDump && !bDumpFrame )
 741       fprintf(stderr,"\nWARNING no frame found near specified time (%g), "
 742               "trajectory ended at %g\n",tdump,t);
 743     fprintf(stderr,"\n");
 744
 745     close_trj(status);
 746
 747     if (ftp == efXTC)
 748       close_xtc(xdr);
 749     else if (out != NULL)
 750       fclose(out);
 751   }
 752
 753   thanx(stdout);
 754
 755   return 0;
 756 }