3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
65 #define TIME_EXPLICIT 0
66 #define TIME_CONTINUE 1
71 #define FLAGS (TRX_READ_X | TRX_READ_V | TRX_READ_F)
73 static void scan_trj_files(char **fnms, int nfiles, real *readtime,
74 real *timestep, atom_id imax,
75 const output_env_t oenv)
77 /* Check start time of all files */
84 for (i = 0; i < nfiles; i++)
86 ok = read_first_frame(oenv, &status, fnms[i], &fr, FLAGS);
89 gmx_fatal(FARGS,"\nCouldn't read frame from file." );
95 fprintf(stderr,"\nWARNING: Couldn't find a time in the frame.\n");
106 if(natoms!=fr.natoms)
107 gmx_fatal(FARGS,"\nDifferent numbers of atoms (%d/%d) in files",
112 if(fr.natoms <= imax)
114 gmx_fatal(FARGS,"\nNot enough atoms (%d) for index group (%d)",
119 ok=read_next_frame(oenv,status,&fr);
122 timestep[i] = fr.time - readtime[i];
137 fprintf(stderr,"\n");
141 static void sort_files(char **fnms, real *settime, int nfile)
147 for (i = 0; i < nfile; i++)
150 for (j = i + 1; j < nfile; j++)
152 if (settime[j] < settime[minidx])
159 timeswap = settime[i];
160 settime[i] = settime[minidx];
161 settime[minidx] = timeswap;
163 fnms[i] = fnms[minidx];
164 fnms[minidx] = chptr;
169 static void edit_files(char **fnms, int nfiles, real *readtime, real *timestep,
170 real *settime, int *cont_type, gmx_bool bSetTime,
171 gmx_bool bSort, const output_env_t oenv)
175 char inputstring[STRLEN], *chptr;
179 fprintf(stderr, "\n\nEnter the new start time (%s) for each file.\n"
180 "There are two special options, both disable sorting:\n\n"
181 "c (continue) - The start time is taken from the end\n"
182 "of the previous file. Use it when your continuation run\n"
183 "restarts with t=0.\n\n"
184 "l (last) - The time in this file will be changed the\n"
185 "same amount as in the previous. Use it when the time in the\n"
186 "new run continues from the end of the previous one,\n"
187 "since this takes possible overlap into account.\n\n",
188 output_env_get_time_unit(oenv));
192 " File Current start (%s) New start (%s)\n"
193 "---------------------------------------------------------\n",
194 output_env_get_time_unit(oenv), output_env_get_time_unit(oenv));
196 for (i = 0; i < nfiles; i++)
198 fprintf(stderr, "%25s %10.3f %s ", fnms[i],
199 output_env_conv_time(oenv, readtime[i]), output_env_get_time_unit(oenv));
203 if (NULL == fgets(inputstring, STRLEN - 1, stdin))
205 gmx_fatal(FARGS,"Error reading user input" );
208 inputstring[strlen(inputstring)-1]=0;
210 if(inputstring[0]=='c' || inputstring[0]=='C')
212 cont_type[i]=TIME_CONTINUE;
217 else if(inputstring[0]=='l' ||
220 cont_type[i]=TIME_LAST;
227 settime[i]=strtod(inputstring,&chptr)*
228 output_env_get_time_invfactor(oenv);
229 if(chptr==inputstring)
231 fprintf(stderr,"'%s' not recognized as a floating point number, 'c' or 'l'. "
232 "Try again: ",inputstring);
235 cont_type[i]=TIME_EXPLICIT;
241 if(cont_type[0]!=TIME_EXPLICIT)
243 cont_type[0]=TIME_EXPLICIT;
249 for(i=0;i<nfiles;i++)
251 settime[i]=readtime[i];
256 fprintf(stderr,"Sorting disabled.\n");
260 sort_files(fnms,settime,nfiles);
262 /* Write out the new order and start times */
263 fprintf(stderr,"\nSummary of files and start times used:\n\n"
264 " File Start time Time step\n"
265 "---------------------------------------------------------\n");
266 for(i=0;i<nfiles;i++)
270 fprintf(stderr,"%25s %10.3f %s %10.3f %s",
272 output_env_conv_time(oenv,settime[i]),output_env_get_time_unit(oenv),
273 output_env_conv_time(oenv,timestep[i]),output_env_get_time_unit(oenv));
275 cont_type[i-1]==TIME_EXPLICIT && settime[i]==settime[i-1] )
276 fprintf(stderr," WARNING: same Start time as previous");
277 fprintf(stderr,"\n");
280 fprintf(stderr,"%25s Continue from last file\n",fnms[i]);
283 fprintf(stderr,"%25s Change by same amount as last file\n",
287 fprintf(stderr,"\n");
289 settime[nfiles]=FLT_MAX;
290 cont_type[nfiles]=TIME_EXPLICIT;
291 readtime[nfiles]=FLT_MAX;
294 static void do_demux(int nset, char *fnms[], char *fnms_out[], int nval,
295 real **value, real *time, real dt_remd, int isize,
296 atom_id index[], real dt, const output_env_t oenv)
298 int i, j, k, natoms, nnn;
299 t_trxstatus **fp_in, **fp_out;
300 gmx_bool bCont, *bSet;
301 real t, first_time = 0;
309 for (i = 0; (i < nset); i++)
311 nnn = read_first_frame(oenv, &(fp_in[i]), fnms[i], &(trx[i]),
316 first_time = trx[i].time;
318 else if (natoms != nnn)
320 gmx_fatal(FARGS,"Trajectory file %s has %d atoms while previous trajs had %d atoms" ,fnms[i],nnn,natoms);
326 else if (t != trx[i].time)
328 gmx_fatal(FARGS,"Trajectory file %s has time %f while previous trajs had time %f",fnms[i],trx[i].time,t);
333 for(i=0; (i<nset); i++)
335 fp_out[i] = open_trx(fnms_out[i],"w");
338 if (gmx_nint(time[k] - t) != 0)
340 gmx_fatal(FARGS,"First time in demuxing table does not match trajectories");
344 while ((k+1 < nval) && ((trx[0].time - time[k+1]) > dt_remd*0.1))
350 fprintf(debug,"trx[0].time = %g, time[k] = %g\n",trx[0].time,time[k]);
352 for(i=0; (i<nset); i++)
356 for(i=0; (i<nset); i++)
358 j = gmx_nint(value[i][k]);
359 range_check(j,0,nset);
362 gmx_fatal(FARGS,"Demuxing the same replica %d twice at time %f",
367 if (dt==0 || bRmod(trx[i].time,first_time,dt))
371 write_trxframe_indexed(fp_out[j],&trx[i],isize,index,NULL);
375 write_trxframe(fp_out[j],&trx[i],NULL);
381 for(i=0; (i<nset); i++)
383 bCont = bCont && read_next_frame(oenv,fp_in[i],&trx[i]);
387 for(i=0; (i<nset); i++)
390 close_trx(fp_out[i]);
394 int gmx_trjcat(int argc, char *argv[])
399 "[TT]trjcat[tt] concatenates several input trajectory files in sorted order. ",
400 "In case of double time frames the one in the later file is used. ",
401 "By specifying [TT]-settime[tt] you will be asked for the start time ",
402 "of each file. The input files are taken from the command line, ",
403 "such that a command like [TT]trjcat -f *.trr -o fixed.trr[tt] should do ",
404 "the trick. Using [TT]-cat[tt], you can simply paste several files ",
405 "together without removal of frames with identical time stamps.[PAR]",
406 "One important option is inferred when the output file is amongst the",
407 "input files. In that case that particular file will be appended to",
408 "which implies you do not need to store double the amount of data.",
409 "Obviously the file to append to has to be the one with lowest starting",
410 "time since one can only append at the end of a file.[PAR]",
411 "If the [TT]-demux[tt] option is given, the N trajectories that are",
412 "read, are written in another order as specified in the [TT].xvg[tt] file.",
413 "The [TT].xvg[tt] file should contain something like:[PAR]",
414 "[TT]0 0 1 2 3 4 5[BR]",
415 "2 1 0 2 3 5 4[tt][BR]",
416 "Where the first number is the time, and subsequent numbers point to",
417 "trajectory indices.",
418 "The frames corresponding to the numbers present at the first line",
419 "are collected into the output trajectory. If the number of frames in",
420 "the trajectory does not match that in the [TT].xvg[tt] file then the program",
421 "tries to be smart. Beware." };
422 static gmx_bool bVels = TRUE;
424 static gmx_bool bCat = FALSE;
425 static gmx_bool bSort = TRUE;
426 static gmx_bool bKeepLast = FALSE;
427 static gmx_bool bKeepLastAppend = FALSE;
428 static gmx_bool bOverwrite = FALSE;
429 static gmx_bool bSetTime = FALSE;
430 static gmx_bool bDeMux;
431 static real begin = -1;
432 static real end = -1;
438 { "-b", FALSE, etTIME,
439 { &begin }, "First time to use (%t)" },
440 { "-e", FALSE, etTIME,
441 { &end }, "Last time to use (%t)" },
442 { "-dt", FALSE, etTIME,
443 { &dt }, "Only write frame when t MOD dt = first time (%t)" },
444 { "-prec", FALSE, etINT,
445 { &prec }, "Precision for [TT].xtc[tt] and [TT].gro[tt] writing in number of decimal places" },
446 { "-vel", FALSE, etBOOL,
447 { &bVels }, "Read and write velocities if possible" },
448 { "-settime", FALSE, etBOOL,
449 { &bSetTime }, "Change starting time interactively" },
450 { "-sort", FALSE, etBOOL,
451 { &bSort }, "Sort trajectory files (not frames)" },
452 { "-keeplast", FALSE, etBOOL,
453 { &bKeepLast }, "Keep overlapping frames at end of trajectory" },
454 { "-overwrite", FALSE, etBOOL,
455 { &bOverwrite }, "Overwrite overlapping frames during appending" },
456 { "-cat", FALSE, etBOOL,
457 { &bCat }, "Do not discard double time frames" } };
458 #define npargs asize(pa)
459 int ftpin, i, frame, frame_out, step = 0, trjout = 0;
463 t_trxframe fr, frout;
464 char **fnms, **fnms_out, *in_file, *out_file;
466 t_trxstatus *trxout = NULL;
467 gmx_bool bNewFile, bIndex, bWrite;
468 int earliersteps, nfile_in, nfile_out, *cont_type, last_ok_step;
469 real *readtime, *timest, *settime;
470 real first_time = 0, lasttime = NOTSET, last_ok_t = -1, timestep;
471 real last_frame_time, searchtime;
473 atom_id *index = NULL, imax;
475 real **val = NULL, *t = NULL, dt_remd;
482 { efTRX, "-f", NULL, ffRDMULT },
483 { efTRO, "-o", NULL, ffWRMULT },
484 { efNDX, "-n", "index", ffOPTRD },
485 { efXVG, "-demux", "remd", ffOPTRD } };
487 #define NFILE asize(fnm)
489 CopyRight(stderr, argv[0]);
490 parse_common_args(&argc, argv, PCA_BE_NICE | PCA_TIME_UNIT, NFILE, fnm,
491 asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
493 bIndex = ftp2bSet(efNDX, NFILE, fnm);
494 bDeMux = ftp2bSet(efXVG, NFILE, fnm);
495 bSort = bSort && !bDeMux;
500 printf("Select group for output\n");
501 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
504 for (i = 1; i < isize; i++)
506 imax = max(imax, index[i]);
513 val = read_xvg_time(opt2fn("-demux", NFILE, fnm), TRUE,
514 opt2parg_bSet("-b", npargs, pa), begin,
515 opt2parg_bSet("-e", npargs, pa), end, 1, &nset, &n,
517 printf("Read %d sets of %d points, dt = %g\n\n", nset, n, dt_remd);
520 fprintf(debug, "Dump of replica_index.xvg\n");
521 for (i = 0; (i < n); i++)
523 fprintf(debug, "%10g", t[i]);
524 for (j = 0; (j < nset); j++)
526 fprintf(debug, " %3d", gmx_nint(val[j][i]));
528 fprintf(debug, "\n");
532 /* prec is in nr of decimal places, xtcprec is a multiplication factor: */
534 for (i = 0; i < prec; i++)
539 nfile_in = opt2fns(&fnms, "-f", NFILE, fnm);
542 gmx_fatal(FARGS,"No input files!" );
545 if (bDeMux && (nfile_in != nset))
547 gmx_fatal(FARGS,"You have specified %d files and %d entries in the demux table",nfile_in,nset);
550 nfile_out = opt2fns(&fnms_out,"-o",NFILE,fnm);
553 gmx_fatal(FARGS,"No output files!");
555 if ((nfile_out > 1) && !bDeMux)
557 gmx_fatal(FARGS,"Don't know what to do with more than 1 output file if not demultiplexing");
559 else if (bDeMux && (nfile_out != nset) && (nfile_out != 1))
561 gmx_fatal(FARGS,"Number of output files should be 1 or %d (#input files), not %d",nset,nfile_out);
565 if (nfile_out != nset)
567 char *buf = strdup(fnms_out[0]);
569 for(i=0; (i<nset); i++)
571 snew(fnms_out[i],strlen(buf)+32);
572 sprintf(fnms_out[i],"%d_%s",i,buf);
576 do_demux(nfile_in,fnms,fnms_out,n,val,t,dt_remd,isize,index,dt,oenv);
580 snew(readtime,nfile_in+1);
581 snew(timest,nfile_in+1);
582 scan_trj_files(fnms,nfile_in,readtime,timest,imax,oenv);
584 snew(settime,nfile_in+1);
585 snew(cont_type,nfile_in+1);
586 edit_files(fnms,nfile_in,readtime,timest,settime,cont_type,bSetTime,bSort,
589 /* Check whether the output file is amongst the input files
590 * This has to be done after sorting etc.
592 out_file = fnms_out[0];
594 for(i=0; ((i<nfile_in) && (n_append==-1)); i++)
596 if (strcmp(fnms[i],out_file) == 0)
602 fprintf(stderr,"Will append to %s rather than creating a new file\n",
605 else if (n_append != -1)
607 gmx_fatal(FARGS,"Can only append to the first file which is %s (not %s)",
612 /* Not checking input format, could be dangerous :-) */
613 /* Not checking output format, equally dangerous :-) */
617 /* the default is not to change the time at all,
618 * but this is overridden by the edit_files routine
624 trxout = open_trx(out_file,"w");
625 memset(&frout,0,sizeof(frout));
631 if (!read_first_frame(oenv,&status,out_file,&fr,FLAGS))
632 gmx_fatal(FARGS,"Reading first frame from %s",out_file);
634 stfio=trx_get_fileio(status);
635 if (!bKeepLast && !bOverwrite)
637 fprintf(stderr, "\n\nWARNING: Appending without -overwrite implies -keeplast "
638 "between the first two files. \n"
639 "If the trajectories have an overlap and have not been written binary \n"
640 "reproducible this will produce an incorrect trajectory!\n\n");
642 /* Fails if last frame is incomplete
643 * We can't do anything about it without overwriting
645 if (gmx_fio_getftp(stfio) == efXTC)
648 xdr_xtc_get_last_frame_time(gmx_fio_getfp(stfio),
649 gmx_fio_getxdr(stfio),
654 gmx_fatal(FARGS,"Error reading last frame. Maybe seek not supported." );
659 while (read_next_frame(oenv,status,&fr))
663 bKeepLastAppend = TRUE;
665 trxout = open_trx(out_file,"a");
669 if (gmx_fio_getftp(stfio) != efXTC) {
670 gmx_fatal(FARGS,"Overwrite only supported for XTC." );
673 xdr_xtc_get_last_frame_time(gmx_fio_getfp(stfio),
674 gmx_fio_getxdr(stfio),
678 gmx_fatal(FARGS,"Error reading last frame. Maybe seek not supported." );
680 /* xtc_seek_time broken for trajectories containing only 1 or 2 frames
681 * or when seek time = 0 */
682 if (nfile_in>1 && settime[1]<last_frame_time+timest[0]*0.5)
684 /* Jump to one time-frame before the start of next
686 searchtime = settime[1]-timest[0]*1.25;
690 searchtime = last_frame_time;
692 if (xtc_seek_time(stfio,searchtime,fr.natoms,TRUE))
694 gmx_fatal(FARGS,"Error seeking to append position.");
696 read_next_frame(oenv,status,&fr);
697 if (fabs(searchtime - fr.time) > timest[0]*0.5)
699 gmx_fatal(FARGS,"Error seeking: attempted to seek to %f but got %f.",
703 fpos = gmx_fio_ftell(stfio);
705 trxout = open_trx(out_file,"r+");
706 if (gmx_fio_seek(trx_get_fileio(trxout),fpos)) {
707 gmx_fatal(FARGS,"Error seeking to append position.");
710 printf("\n Will append after %f \n",lasttime);
713 /* Lets stitch up some files */
714 timestep = timest[0];
715 for(i=n_append+1; (i<nfile_in); i++)
719 /* set the next time from the last frame in previous file */
724 if(cont_type[i]==TIME_CONTINUE)
727 begin += 0.5*timestep;
728 settime[i]=frout.time;
729 cont_type[i]=TIME_EXPLICIT;
731 else if(cont_type[i]==TIME_LAST)
734 begin += 0.5*timestep;
736 /* Or, if the time in the next part should be changed by the
737 * same amount, start at half a timestep from the last time
738 * so we dont repeat frames.
740 /* I don't understand the comment above, but for all the cases
741 * I tried the code seems to work properly. B. Hess 2008-4-2.
744 /* Or, if time is set explicitly, we check for overlap/gap */
745 if(cont_type[i]==TIME_EXPLICIT)
747 if( ( i < nfile_in ) &&
748 ( frout.time < settime[i]-1.5*timestep ) )
750 fprintf(stderr, "WARNING: Frames around t=%f %s have a different "
751 "spacing than the rest,\n"
752 "might be a gap or overlap that couldn't be corrected "
753 "automatically.\n",output_env_conv_time(oenv,frout.time),
754 output_env_get_time_unit(oenv));
759 /* if we don't have a timestep in the current file, use the old one */
760 if ( timest[i] != 0 )
762 timestep = timest[i];
764 read_first_frame(oenv,&status,fnms[i],&fr,FLAGS);
768 fprintf(stderr,"\nWARNING: Couldn't find a time in the frame.\n");
771 if(cont_type[i]==TIME_EXPLICIT)
773 t_corr=settime[i]-fr.time;
775 /* t_corr is the amount we want to change the time.
776 * If the user has chosen not to change the time for
777 * this part of the trajectory t_corr remains at
778 * the value it had in the last part, changing this
779 * by the same amount.
780 * If no value was given for the first trajectory part
781 * we let the time start at zero, see the edit_files routine.
787 if (lasttime != NOTSET)
789 printf("lasttime %g\n", lasttime);
794 /* copy the input frame to the output frame */
796 /* set the new time by adding the correct calculated above */
797 frout.time += t_corr;
798 /* quit if we have reached the end of what should be written */
799 if((end > 0) && (frout.time > end+GMX_REAL_EPS))
805 /* determine if we should write this frame (dt is handled elsewhere) */
806 if (bCat) /* write all frames of all files */
810 else if ( bKeepLast || (bKeepLastAppend && i==1))
811 /* write till last frame of this traj
812 and skip first frame(s) of next traj */
814 bWrite = ( frout.time > lasttime+0.5*timestep );
816 else /* write till first frame of next traj */
818 bWrite = ( frout.time < settime[i+1]-0.5*timestep );
821 if( bWrite && (frout.time >= begin) )
825 first_time = frout.time;
826 lasttime = frout.time;
827 if (dt==0 || bRmod(frout.time,first_time,dt))
830 last_ok_t=frout.time;
833 fprintf(stderr,"\nContinue writing frames from %s t=%g %s, "
835 fnms[i],output_env_conv_time(oenv,frout.time),output_env_get_time_unit(oenv),
842 write_trxframe_indexed(trxout,&frout,isize,index,
847 write_trxframe(trxout,&frout,NULL);
849 if ( ((frame % 10) == 0) || (frame < 10) )
851 fprintf(stderr," -> frame %6d time %8.3f %s \r",
852 frame_out,output_env_conv_time(oenv,frout.time),output_env_get_time_unit(oenv));
856 } while( read_next_frame(oenv,status,&fr));
866 fprintf(stderr,"\nLast frame written was %d, time %f %s\n",
867 frame,output_env_conv_time(oenv,last_ok_t),output_env_get_time_unit(oenv));