2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
68 #define TIME_EXPLICIT 0
69 #define TIME_CONTINUE 1
74 #define FLAGS (TRX_READ_X | TRX_READ_V | TRX_READ_F)
76 static void scan_trj_files(char **fnms, int nfiles, real *readtime,
77 real *timestep, atom_id imax,
78 const output_env_t oenv)
80 /* Check start time of all files */
87 for (i = 0; i < nfiles; i++)
89 ok = read_first_frame(oenv, &status, fnms[i], &fr, FLAGS);
92 gmx_fatal(FARGS,"\nCouldn't read frame from file." );
98 fprintf(stderr,"\nWARNING: Couldn't find a time in the frame.\n");
109 if(natoms!=fr.natoms)
110 gmx_fatal(FARGS,"\nDifferent numbers of atoms (%d/%d) in files",
115 if(fr.natoms <= imax)
117 gmx_fatal(FARGS,"\nNot enough atoms (%d) for index group (%d)",
122 ok=read_next_frame(oenv,status,&fr);
125 timestep[i] = fr.time - readtime[i];
140 fprintf(stderr,"\n");
144 static void sort_files(char **fnms, real *settime, int nfile)
150 for (i = 0; i < nfile; i++)
153 for (j = i + 1; j < nfile; j++)
155 if (settime[j] < settime[minidx])
162 timeswap = settime[i];
163 settime[i] = settime[minidx];
164 settime[minidx] = timeswap;
166 fnms[i] = fnms[minidx];
167 fnms[minidx] = chptr;
172 static void edit_files(char **fnms, int nfiles, real *readtime, real *timestep,
173 real *settime, int *cont_type, gmx_bool bSetTime,
174 gmx_bool bSort, const output_env_t oenv)
178 char inputstring[STRLEN], *chptr;
182 fprintf(stderr, "\n\nEnter the new start time (%s) for each file.\n"
183 "There are two special options, both disable sorting:\n\n"
184 "c (continue) - The start time is taken from the end\n"
185 "of the previous file. Use it when your continuation run\n"
186 "restarts with t=0.\n\n"
187 "l (last) - The time in this file will be changed the\n"
188 "same amount as in the previous. Use it when the time in the\n"
189 "new run continues from the end of the previous one,\n"
190 "since this takes possible overlap into account.\n\n",
191 output_env_get_time_unit(oenv));
195 " File Current start (%s) New start (%s)\n"
196 "---------------------------------------------------------\n",
197 output_env_get_time_unit(oenv), output_env_get_time_unit(oenv));
199 for (i = 0; i < nfiles; i++)
201 fprintf(stderr, "%25s %10.3f %s ", fnms[i],
202 output_env_conv_time(oenv, readtime[i]), output_env_get_time_unit(oenv));
206 if (NULL == fgets(inputstring, STRLEN - 1, stdin))
208 gmx_fatal(FARGS,"Error reading user input" );
211 inputstring[strlen(inputstring)-1]=0;
213 if(inputstring[0]=='c' || inputstring[0]=='C')
215 cont_type[i]=TIME_CONTINUE;
220 else if(inputstring[0]=='l' ||
223 cont_type[i]=TIME_LAST;
230 settime[i]=strtod(inputstring,&chptr)*
231 output_env_get_time_invfactor(oenv);
232 if(chptr==inputstring)
234 fprintf(stderr,"'%s' not recognized as a floating point number, 'c' or 'l'. "
235 "Try again: ",inputstring);
238 cont_type[i]=TIME_EXPLICIT;
244 if(cont_type[0]!=TIME_EXPLICIT)
246 cont_type[0]=TIME_EXPLICIT;
252 for(i=0;i<nfiles;i++)
254 settime[i]=readtime[i];
259 fprintf(stderr,"Sorting disabled.\n");
263 sort_files(fnms,settime,nfiles);
265 /* Write out the new order and start times */
266 fprintf(stderr,"\nSummary of files and start times used:\n\n"
267 " File Start time Time step\n"
268 "---------------------------------------------------------\n");
269 for(i=0;i<nfiles;i++)
273 fprintf(stderr,"%25s %10.3f %s %10.3f %s",
275 output_env_conv_time(oenv,settime[i]),output_env_get_time_unit(oenv),
276 output_env_conv_time(oenv,timestep[i]),output_env_get_time_unit(oenv));
278 cont_type[i-1]==TIME_EXPLICIT && settime[i]==settime[i-1] )
279 fprintf(stderr," WARNING: same Start time as previous");
280 fprintf(stderr,"\n");
283 fprintf(stderr,"%25s Continue from last file\n",fnms[i]);
286 fprintf(stderr,"%25s Change by same amount as last file\n",
290 fprintf(stderr,"\n");
292 settime[nfiles]=FLT_MAX;
293 cont_type[nfiles]=TIME_EXPLICIT;
294 readtime[nfiles]=FLT_MAX;
297 static void do_demux(int nset, char *fnms[], char *fnms_out[], int nval,
298 real **value, real *time, real dt_remd, int isize,
299 atom_id index[], real dt, const output_env_t oenv)
301 int i, j, k, natoms, nnn;
302 t_trxstatus **fp_in, **fp_out;
303 gmx_bool bCont, *bSet;
304 real t, first_time = 0;
312 for (i = 0; (i < nset); i++)
314 nnn = read_first_frame(oenv, &(fp_in[i]), fnms[i], &(trx[i]),
319 first_time = trx[i].time;
321 else if (natoms != nnn)
323 gmx_fatal(FARGS,"Trajectory file %s has %d atoms while previous trajs had %d atoms" ,fnms[i],nnn,natoms);
329 else if (t != trx[i].time)
331 gmx_fatal(FARGS,"Trajectory file %s has time %f while previous trajs had time %f",fnms[i],trx[i].time,t);
336 for(i=0; (i<nset); i++)
338 fp_out[i] = open_trx(fnms_out[i],"w");
341 if (gmx_nint(time[k] - t) != 0)
343 gmx_fatal(FARGS,"First time in demuxing table does not match trajectories");
347 while ((k+1 < nval) && ((trx[0].time - time[k+1]) > dt_remd*0.1))
353 fprintf(debug,"trx[0].time = %g, time[k] = %g\n",trx[0].time,time[k]);
355 for(i=0; (i<nset); i++)
359 for(i=0; (i<nset); i++)
361 j = gmx_nint(value[i][k]);
362 range_check(j,0,nset);
365 gmx_fatal(FARGS,"Demuxing the same replica %d twice at time %f",
370 if (dt==0 || bRmod(trx[i].time,first_time,dt))
374 write_trxframe_indexed(fp_out[j],&trx[i],isize,index,NULL);
378 write_trxframe(fp_out[j],&trx[i],NULL);
384 for(i=0; (i<nset); i++)
386 bCont = bCont && read_next_frame(oenv,fp_in[i],&trx[i]);
390 for(i=0; (i<nset); i++)
393 close_trx(fp_out[i]);
397 int gmx_trjcat(int argc, char *argv[])
402 "[TT]trjcat[tt] concatenates several input trajectory files in sorted order. ",
403 "In case of double time frames the one in the later file is used. ",
404 "By specifying [TT]-settime[tt] you will be asked for the start time ",
405 "of each file. The input files are taken from the command line, ",
406 "such that a command like [TT]trjcat -f *.trr -o fixed.trr[tt] should do ",
407 "the trick. Using [TT]-cat[tt], you can simply paste several files ",
408 "together without removal of frames with identical time stamps.[PAR]",
409 "One important option is inferred when the output file is amongst the",
410 "input files. In that case that particular file will be appended to",
411 "which implies you do not need to store double the amount of data.",
412 "Obviously the file to append to has to be the one with lowest starting",
413 "time since one can only append at the end of a file.[PAR]",
414 "If the [TT]-demux[tt] option is given, the N trajectories that are",
415 "read, are written in another order as specified in the [TT].xvg[tt] file.",
416 "The [TT].xvg[tt] file should contain something like:[PAR]",
417 "[TT]0 0 1 2 3 4 5[BR]",
418 "2 1 0 2 3 5 4[tt][BR]",
419 "Where the first number is the time, and subsequent numbers point to",
420 "trajectory indices.",
421 "The frames corresponding to the numbers present at the first line",
422 "are collected into the output trajectory. If the number of frames in",
423 "the trajectory does not match that in the [TT].xvg[tt] file then the program",
424 "tries to be smart. Beware." };
425 static gmx_bool bVels = TRUE;
427 static gmx_bool bCat = FALSE;
428 static gmx_bool bSort = TRUE;
429 static gmx_bool bKeepLast = FALSE;
430 static gmx_bool bKeepLastAppend = FALSE;
431 static gmx_bool bOverwrite = FALSE;
432 static gmx_bool bSetTime = FALSE;
433 static gmx_bool bDeMux;
434 static real begin = -1;
435 static real end = -1;
441 { "-b", FALSE, etTIME,
442 { &begin }, "First time to use (%t)" },
443 { "-e", FALSE, etTIME,
444 { &end }, "Last time to use (%t)" },
445 { "-dt", FALSE, etTIME,
446 { &dt }, "Only write frame when t MOD dt = first time (%t)" },
447 { "-prec", FALSE, etINT,
448 { &prec }, "Precision for [TT].xtc[tt] and [TT].gro[tt] writing in number of decimal places" },
449 { "-vel", FALSE, etBOOL,
450 { &bVels }, "Read and write velocities if possible" },
451 { "-settime", FALSE, etBOOL,
452 { &bSetTime }, "Change starting time interactively" },
453 { "-sort", FALSE, etBOOL,
454 { &bSort }, "Sort trajectory files (not frames)" },
455 { "-keeplast", FALSE, etBOOL,
456 { &bKeepLast }, "Keep overlapping frames at end of trajectory" },
457 { "-overwrite", FALSE, etBOOL,
458 { &bOverwrite }, "Overwrite overlapping frames during appending" },
459 { "-cat", FALSE, etBOOL,
460 { &bCat }, "Do not discard double time frames" } };
461 #define npargs asize(pa)
462 int ftpin, i, frame, frame_out, step = 0, trjout = 0;
466 t_trxframe fr, frout;
467 char **fnms, **fnms_out, *in_file, *out_file;
469 t_trxstatus *trxout = NULL;
470 gmx_bool bNewFile, bIndex, bWrite;
471 int earliersteps, nfile_in, nfile_out, *cont_type, last_ok_step;
472 real *readtime, *timest, *settime;
473 real first_time = 0, lasttime = NOTSET, last_ok_t = -1, timestep;
474 real last_frame_time, searchtime;
476 atom_id *index = NULL, imax;
478 real **val = NULL, *t = NULL, dt_remd;
485 { efTRX, "-f", NULL, ffRDMULT },
486 { efTRO, "-o", NULL, ffWRMULT },
487 { efNDX, "-n", "index", ffOPTRD },
488 { efXVG, "-demux", "remd", ffOPTRD } };
490 #define NFILE asize(fnm)
492 CopyRight(stderr, argv[0]);
493 parse_common_args(&argc, argv, PCA_BE_NICE | PCA_TIME_UNIT, NFILE, fnm,
494 asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
496 bIndex = ftp2bSet(efNDX, NFILE, fnm);
497 bDeMux = ftp2bSet(efXVG, NFILE, fnm);
498 bSort = bSort && !bDeMux;
503 printf("Select group for output\n");
504 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
507 for (i = 1; i < isize; i++)
509 imax = max(imax, index[i]);
516 val = read_xvg_time(opt2fn("-demux", NFILE, fnm), TRUE,
517 opt2parg_bSet("-b", npargs, pa), begin,
518 opt2parg_bSet("-e", npargs, pa), end, 1, &nset, &n,
520 printf("Read %d sets of %d points, dt = %g\n\n", nset, n, dt_remd);
523 fprintf(debug, "Dump of replica_index.xvg\n");
524 for (i = 0; (i < n); i++)
526 fprintf(debug, "%10g", t[i]);
527 for (j = 0; (j < nset); j++)
529 fprintf(debug, " %3d", gmx_nint(val[j][i]));
531 fprintf(debug, "\n");
535 /* prec is in nr of decimal places, xtcprec is a multiplication factor: */
537 for (i = 0; i < prec; i++)
542 nfile_in = opt2fns(&fnms, "-f", NFILE, fnm);
545 gmx_fatal(FARGS,"No input files!" );
548 if (bDeMux && (nfile_in != nset))
550 gmx_fatal(FARGS,"You have specified %d files and %d entries in the demux table",nfile_in,nset);
553 nfile_out = opt2fns(&fnms_out,"-o",NFILE,fnm);
556 gmx_fatal(FARGS,"No output files!");
558 if ((nfile_out > 1) && !bDeMux)
560 gmx_fatal(FARGS,"Don't know what to do with more than 1 output file if not demultiplexing");
562 else if (bDeMux && (nfile_out != nset) && (nfile_out != 1))
564 gmx_fatal(FARGS,"Number of output files should be 1 or %d (#input files), not %d",nset,nfile_out);
568 if (nfile_out != nset)
570 char *buf = strdup(fnms_out[0]);
572 for(i=0; (i<nset); i++)
574 snew(fnms_out[i],strlen(buf)+32);
575 sprintf(fnms_out[i],"%d_%s",i,buf);
579 do_demux(nfile_in,fnms,fnms_out,n,val,t,dt_remd,isize,index,dt,oenv);
583 snew(readtime,nfile_in+1);
584 snew(timest,nfile_in+1);
585 scan_trj_files(fnms,nfile_in,readtime,timest,imax,oenv);
587 snew(settime,nfile_in+1);
588 snew(cont_type,nfile_in+1);
589 edit_files(fnms,nfile_in,readtime,timest,settime,cont_type,bSetTime,bSort,
592 /* Check whether the output file is amongst the input files
593 * This has to be done after sorting etc.
595 out_file = fnms_out[0];
597 for(i=0; ((i<nfile_in) && (n_append==-1)); i++)
599 if (strcmp(fnms[i],out_file) == 0)
605 fprintf(stderr,"Will append to %s rather than creating a new file\n",
608 else if (n_append != -1)
610 gmx_fatal(FARGS,"Can only append to the first file which is %s (not %s)",
615 /* Not checking input format, could be dangerous :-) */
616 /* Not checking output format, equally dangerous :-) */
620 /* the default is not to change the time at all,
621 * but this is overridden by the edit_files routine
627 trxout = open_trx(out_file,"w");
628 memset(&frout,0,sizeof(frout));
634 if (!read_first_frame(oenv,&status,out_file,&fr,FLAGS))
635 gmx_fatal(FARGS,"Reading first frame from %s",out_file);
637 stfio=trx_get_fileio(status);
638 if (!bKeepLast && !bOverwrite)
640 fprintf(stderr, "\n\nWARNING: Appending without -overwrite implies -keeplast "
641 "between the first two files. \n"
642 "If the trajectories have an overlap and have not been written binary \n"
643 "reproducible this will produce an incorrect trajectory!\n\n");
645 /* Fails if last frame is incomplete
646 * We can't do anything about it without overwriting
648 if (gmx_fio_getftp(stfio) == efXTC)
651 xdr_xtc_get_last_frame_time(gmx_fio_getfp(stfio),
652 gmx_fio_getxdr(stfio),
657 gmx_fatal(FARGS,"Error reading last frame. Maybe seek not supported." );
662 while (read_next_frame(oenv,status,&fr))
666 bKeepLastAppend = TRUE;
668 trxout = open_trx(out_file,"a");
672 if (gmx_fio_getftp(stfio) != efXTC) {
673 gmx_fatal(FARGS,"Overwrite only supported for XTC." );
676 xdr_xtc_get_last_frame_time(gmx_fio_getfp(stfio),
677 gmx_fio_getxdr(stfio),
681 gmx_fatal(FARGS,"Error reading last frame. Maybe seek not supported." );
683 /* xtc_seek_time broken for trajectories containing only 1 or 2 frames
684 * or when seek time = 0 */
685 if (nfile_in>1 && settime[1]<last_frame_time+timest[0]*0.5)
687 /* Jump to one time-frame before the start of next
689 searchtime = settime[1]-timest[0]*1.25;
693 searchtime = last_frame_time;
695 if (xtc_seek_time(stfio,searchtime,fr.natoms,TRUE))
697 gmx_fatal(FARGS,"Error seeking to append position.");
699 read_next_frame(oenv,status,&fr);
700 if (fabs(searchtime - fr.time) > timest[0]*0.5)
702 gmx_fatal(FARGS,"Error seeking: attempted to seek to %f but got %f.",
706 fpos = gmx_fio_ftell(stfio);
708 trxout = open_trx(out_file,"r+");
709 if (gmx_fio_seek(trx_get_fileio(trxout),fpos)) {
710 gmx_fatal(FARGS,"Error seeking to append position.");
713 printf("\n Will append after %f \n",lasttime);
716 /* Lets stitch up some files */
717 timestep = timest[0];
718 for(i=n_append+1; (i<nfile_in); i++)
722 /* set the next time from the last frame in previous file */
727 if(cont_type[i]==TIME_CONTINUE)
730 begin += 0.5*timestep;
731 settime[i]=frout.time;
732 cont_type[i]=TIME_EXPLICIT;
734 else if(cont_type[i]==TIME_LAST)
737 begin += 0.5*timestep;
739 /* Or, if the time in the next part should be changed by the
740 * same amount, start at half a timestep from the last time
741 * so we dont repeat frames.
743 /* I don't understand the comment above, but for all the cases
744 * I tried the code seems to work properly. B. Hess 2008-4-2.
747 /* Or, if time is set explicitly, we check for overlap/gap */
748 if(cont_type[i]==TIME_EXPLICIT)
750 if( ( i < nfile_in ) &&
751 ( frout.time < settime[i]-1.5*timestep ) )
753 fprintf(stderr, "WARNING: Frames around t=%f %s have a different "
754 "spacing than the rest,\n"
755 "might be a gap or overlap that couldn't be corrected "
756 "automatically.\n",output_env_conv_time(oenv,frout.time),
757 output_env_get_time_unit(oenv));
762 /* if we don't have a timestep in the current file, use the old one */
763 if ( timest[i] != 0 )
765 timestep = timest[i];
767 read_first_frame(oenv,&status,fnms[i],&fr,FLAGS);
771 fprintf(stderr,"\nWARNING: Couldn't find a time in the frame.\n");
774 if(cont_type[i]==TIME_EXPLICIT)
776 t_corr=settime[i]-fr.time;
778 /* t_corr is the amount we want to change the time.
779 * If the user has chosen not to change the time for
780 * this part of the trajectory t_corr remains at
781 * the value it had in the last part, changing this
782 * by the same amount.
783 * If no value was given for the first trajectory part
784 * we let the time start at zero, see the edit_files routine.
790 if (lasttime != NOTSET)
792 printf("lasttime %g\n", lasttime);
797 /* copy the input frame to the output frame */
799 /* set the new time by adding the correct calculated above */
800 frout.time += t_corr;
801 /* quit if we have reached the end of what should be written */
802 if((end > 0) && (frout.time > end+GMX_REAL_EPS))
808 /* determine if we should write this frame (dt is handled elsewhere) */
809 if (bCat) /* write all frames of all files */
813 else if ( bKeepLast || (bKeepLastAppend && i==1))
814 /* write till last frame of this traj
815 and skip first frame(s) of next traj */
817 bWrite = ( frout.time > lasttime+0.5*timestep );
819 else /* write till first frame of next traj */
821 bWrite = ( frout.time < settime[i+1]-0.5*timestep );
824 if( bWrite && (frout.time >= begin) )
828 first_time = frout.time;
829 lasttime = frout.time;
830 if (dt==0 || bRmod(frout.time,first_time,dt))
833 last_ok_t=frout.time;
836 fprintf(stderr,"\nContinue writing frames from %s t=%g %s, "
838 fnms[i],output_env_conv_time(oenv,frout.time),output_env_get_time_unit(oenv),
845 write_trxframe_indexed(trxout,&frout,isize,index,
850 write_trxframe(trxout,&frout,NULL);
852 if ( ((frame % 10) == 0) || (frame < 10) )
854 fprintf(stderr," -> frame %6d time %8.3f %s \r",
855 frame_out,output_env_conv_time(oenv,frout.time),output_env_get_time_unit(oenv));
859 } while( read_next_frame(oenv,status,&fr));
869 fprintf(stderr,"\nLast frame written was %d, time %f %s\n",
870 frame,output_env_conv_time(oenv,last_ok_t),output_env_get_time_unit(oenv));