?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixed make_edi.c
[gromacs/rigid-bodies.git] / src / tools / gmx_trjcat.c
blob164e7b156361c58a04e1fca5b8e35f681d29efa6
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.2.0
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.
14
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.
19 *
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.
26 *
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.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * Green Red Orange Magenta Azure Cyan Skyblue
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
39 #include <string.h>
40 #include <math.h>
41 #include "macros.h"
42 #include "sysstuff.h"
43 #include "smalloc.h"
44 #include "typedefs.h"
45 #include "copyrite.h"
46 #include "gmxfio.h"
47 #include "tpxio.h"
48 #include "trnio.h"
49 #include "statutil.h"
50 #include "futil.h"
51 #include "pdbio.h"
52 #include "confio.h"
53 #include "names.h"
54 #include "index.h"
55 #include "vec.h"
56 #include "xtcio.h"
57 #include "do_fit.h"
58 #include "rmpbc.h"
59 #include "wgms.h"
60 #include "magic.h"
61 #include "pbc.h"
62 #include "xvgr.h"
63 #include "xdrf.h"
64 #include "gmx_ana.h"
65
66 #define TIME_EXPLICIT 0
67 #define TIME_CONTINUE 1
68 #define TIME_LAST 2
69 #ifndef FLT_MAX
70 #define FLT_MAX 1e36
71 #endif
72 #define FLAGS (TRX_READ_X | TRX_READ_V | TRX_READ_F)
73
74 static void scan_trj_files(char **fnms, int nfiles, real *readtime,
75 real *timestep, atom_id imax,
76 const output_env_t oenv)
77 {
78 /* Check start time of all files */
79 int i, status, natoms = 0;
80 real t;
81 t_trxframe fr;
82 bool ok;
83
84 for (i = 0; i < nfiles; i++)
85 {
86 ok = read_first_frame(oenv, &status, fnms[i], &fr, FLAGS);
87
88 if (!ok)
89 gmx_fatal(FARGS,"\nCouldn't read frame from file." );
90 if(fr.bTime)
91 readtime[i]=fr.time;
92 else
93 {
94 readtime[i]=0;
95 fprintf(stderr,"\nWARNING: Couldn't find a time in the frame.\n");
96 }
97
98 if(i==0)
99 {
100 natoms=fr.natoms;
101 }
102 else
103 {
104 if (imax==NO_ATID)
105 {
106 if(natoms!=fr.natoms)
107 gmx_fatal(FARGS,"\nDifferent numbers of atoms (%d/%d) in files",
108 natoms,fr.natoms);
109 }
110 else
111 {
112 if(fr.natoms <= imax)
113 {
114 gmx_fatal(FARGS,"\nNot enough atoms (%d) for index group (%d)",
115 fr.natoms,imax);
116 }
117 }
118 }
119 ok=read_next_frame(oenv,status,&fr);
120 if(ok && fr.bTime)
121 {
122 timestep[i] = fr.time - readtime[i];
123 }
124 else
125 {
126 timestep[i] = 0;
127 }
128
129 close_trj(status);
130 }
131 fprintf(stderr,"\n");
132
133 sfree(fr.x);
134 }
135
136 static void sort_files(char **fnms, real *settime, int nfile)
137 {
138 int i, j, minidx;
139 real timeswap;
140 char *chptr;
141
142 for (i = 0; i < nfile; i++)
143 {
144 minidx = i;
145 for (j = i + 1; j < nfile; j++)
146 {
147 if (settime[j] < settime[minidx])
148 {
149 minidx = j;
150 }
151 }
152 if (minidx != i)
153 {
154 timeswap = settime[i];
155 settime[i] = settime[minidx];
156 settime[minidx] = timeswap;
157 chptr = fnms[i];
158 fnms[i] = fnms[minidx];
159 fnms[minidx] = chptr;
160 }
161 }
162 }
163
164 static void edit_files(char **fnms, int nfiles, real *readtime, real *timestep,
165 real *settime, int *cont_type, bool bSetTime,
166 bool bSort, const output_env_t oenv)
167 {
168 int i;
169 bool ok;
170 char inputstring[STRLEN], *chptr;
171
172 if (bSetTime)
173 {
174 fprintf(stderr, "\n\nEnter the new start time (%s) for each file.\n"
175 "There are two special options, both disable sorting:\n\n"
176 "c (continue) - The start time is taken from the end\n"
177 "of the previous file. Use it when your continuation run\n"
178 "restarts with t=0.\n\n"
179 "l (last) - The time in this file will be changed the\n"
180 "same amount as in the previous. Use it when the time in the\n"
181 "new run continues from the end of the previous one,\n"
182 "since this takes possible overlap into account.\n\n",
183 output_env_get_time_unit(oenv));
184
185 fprintf(
186 stderr,
187 " File Current start (%s) New start (%s)\n"
188 "---------------------------------------------------------\n",
189 output_env_get_time_unit(oenv), output_env_get_time_unit(oenv));
190
191 for (i = 0; i < nfiles; i++)
192 {
193 fprintf(stderr, "%25s %10.3f %s ", fnms[i],
194 output_env_conv_time(oenv, readtime[i]), output_env_get_time_unit(oenv));
195 ok = FALSE;
196 do
197 {
198 if (NULL == fgets(inputstring, STRLEN - 1, stdin))
199 {
200 gmx_fatal(FARGS,"Error reading user input" );
201 }
202
203 inputstring[strlen(inputstring)-1]=0;
204
205 if(inputstring[0]=='c' || inputstring[0]=='C')
206 {
207 cont_type[i]=TIME_CONTINUE;
208 bSort=FALSE;
209 ok=TRUE;
210 settime[i]=FLT_MAX;
211 }
212 else if(inputstring[0]=='l' ||
213 inputstring[0]=='L')
214 {
215 cont_type[i]=TIME_LAST;
216 bSort=FALSE;
217 ok=TRUE;
218 settime[i]=FLT_MAX;
219 }
220 else
221 {
222 settime[i]=strtod(inputstring,&chptr)*
223 output_env_get_time_invfactor(oenv);
224 if(chptr==inputstring)
225 {
226 fprintf(stderr,"'%s' not recognized as a floating point number, 'c' or 'l'. "
227 "Try again: ",inputstring);
228 }
229 else {
230 cont_type[i]=TIME_EXPLICIT;
231 ok=TRUE;
232 }
233 }
234 } while (!ok);
235 }
236 if(cont_type[0]!=TIME_EXPLICIT)
237 {
238 cont_type[0]=TIME_EXPLICIT;
239 settime[0]=0;
240 }
241 }
242 else
243 {
244 for(i=0;i<nfiles;i++)
245 {
246 settime[i]=readtime[i];
247 }
248 }
249 if(!bSort)
250 {
251 fprintf(stderr,"Sorting disabled.\n");
252 }
253 else
254 {
255 sort_files(fnms,settime,nfiles);
256 }
257 /* Write out the new order and start times */
258 fprintf(stderr,"\nSummary of files and start times used:\n\n"
259 " File Start time Time step\n"
260 "---------------------------------------------------------\n");
261 for(i=0;i<nfiles;i++)
262 switch(cont_type[i])
263 {
264 case TIME_EXPLICIT:
265 fprintf(stderr,"%25s %10.3f %s %10.3f %s",
266 fnms[i],
267 output_env_conv_time(oenv,settime[i]),output_env_get_time_unit(oenv),
268 output_env_conv_time(oenv,timestep[i]),output_env_get_time_unit(oenv));
269 if ( i>0 &&
270 cont_type[i-1]==TIME_EXPLICIT && settime[i]==settime[i-1] )
271 fprintf(stderr," WARNING: same Start time as previous");
272 fprintf(stderr,"\n");
273 break;
274 case TIME_CONTINUE:
275 fprintf(stderr,"%25s Continue from last file\n",fnms[i]);
276 break;
277 case TIME_LAST:
278 fprintf(stderr,"%25s Change by same amount as last file\n",
279 fnms[i]);
280 break;
281 }
282 fprintf(stderr,"\n");
283
284 settime[nfiles]=FLT_MAX;
285 cont_type[nfiles]=TIME_EXPLICIT;
286 readtime[nfiles]=FLT_MAX;
287 }
288
289 static void do_demux(int nset, char *fnms[], char *fnms_out[], int nval,
290 real **value, real *time, real dt_remd, int isize,
291 atom_id index[], real dt, const output_env_t oenv)
292 {
293 int i, j, k, natoms, nnn;
294 int *fp_in, *fp_out;
295 bool bCont, *bSet;
296 real t, first_time = 0;
297 t_trxframe *trx;
298
299 snew(fp_in,nset);
300 snew(trx,nset);
301 snew(bSet,nset);
302 natoms = -1;
303 t = -1;
304 for (i = 0; (i < nset); i++)
305 {
306 nnn = read_first_frame(oenv, &(fp_in[i]), fnms[i], &(trx[i]),
307 TRX_NEED_X);
308 if (natoms == -1)
309 {
310 natoms = nnn;
311 first_time = trx[i].time;
312 }
313 else if (natoms != nnn)
314 {
315 gmx_fatal(FARGS,"Trajectory file %s has %d atoms while previous trajs had %d atoms" ,fnms[i],nnn,natoms);
316 }
317 if (t == -1)
318 {
319 t = trx[i].time;
320 }
321 else if (t != trx[i].time)
322 {
323 gmx_fatal(FARGS,"Trajectory file %s has time %f while previous trajs had time %f",fnms[i],trx[i].time,t);
324 }
325 }
326
327 snew(fp_out,nset);
328 for(i=0; (i<nset); i++)
329 {
330 fp_out[i] = open_trx(fnms_out[i],"w");
331 }
332 k = 0;
333 if (gmx_nint(time[k] - t) != 0)
334 {
335 gmx_fatal(FARGS,"First time in demuxing table does not match trajectories");
336 }
337 do
338 {
339 while ((k+1 < nval) && ((trx[0].time - time[k+1]) > dt_remd*0.1))
340 {
341 k++;
342 }
343 if (debug)
344 {
345 fprintf(debug,"trx[0].time = %g, time[k] = %g\n",trx[0].time,time[k]);
346 }
347 for(i=0; (i<nset); i++)
348 {
349 bSet[i] = FALSE;
350 }
351 for(i=0; (i<nset); i++)
352 {
353 j = gmx_nint(value[i][k]);
354 range_check(j,0,nset);
355 if (bSet[j])
356 {
357 gmx_fatal(FARGS,"Demuxing the same replica %d twice at time %f",
358 j,trx[0].time);
359 }
360 bSet[j] = TRUE;
361
362 if (dt==0 || bRmod(trx[i].time,first_time,dt))
363 {
364 if (index)
365 {
366 write_trxframe_indexed(fp_out[j],&trx[i],isize,index,NULL);
367 }
368 else
369 {
370 write_trxframe(fp_out[j],&trx[i],NULL);
371 }
372 }
373 }
374
375 bCont = (k < nval);
376 for(i=0; (i<nset); i++)
377 {
378 bCont = bCont && read_next_frame(oenv,fp_in[i],&trx[i]);
379 }
380 } while (bCont);
381
382 for(i=0; (i<nset); i++)
383 {
384 close_trx(fp_in[i]);
385 close_trx(fp_out[i]);
386 }
387 }
388
389 int gmx_trjcat(int argc, char *argv[])
390 {
391 const char
392 *desc[] =
393 {
394 "trjcat concatenates several input trajectory files in sorted order. ",
395 "In case of double time frames the one in the later file is used. ",
396 "By specifying [TT]-settime[tt] you will be asked for the start time ",
397 "of each file. The input files are taken from the command line, ",
398 "such that a command like [TT]trjcat -o fixed.trr *.trr[tt] should do ",
399 "the trick. Using [TT]-cat[tt] you can simply paste several files ",
400 "together without removal of frames with identical time stamps.[PAR]",
401 "One important option is inferred when the output file is amongst the",
402 "input files. In that case that particular file will be appended to",
403 "which implies you do not need to store double the amount of data.",
404 "Obviously the file to append to has to be the one with lowest starting",
405 "time since one can only append at the end of a file.[PAR]",
406 "If the [TT]-demux[tt] option is given, the N trajectories that are",
407 "read, are written in another order as specified in the xvg file.",
408 "The xvg file should contain something like:[PAR]",
409 "0 0 1 2 3 4 5[BR]",
410 "2 1 0 2 3 5 4[BR]",
411 "Where the first number is the time, and subsequent numbers point to",
412 "trajectory indices.",
413 "The frames corresponding to the numbers present at the first line",
414 "are collected into the output trajectory. If the number of frames in",
415 "the trajectory does not match that in the xvg file then the program",
416 "tries to be smart. Beware." };
417 static bool bVels = TRUE;
418 static int prec = 3;
419 static bool bCat = FALSE;
420 static bool bSort = TRUE;
421 static bool bKeepLast = FALSE;
422 static bool bKeepLastAppend = FALSE;
423 static bool bOverwrite = FALSE;
424 static bool bSetTime = FALSE;
425 static bool bDeMux;
426 static real begin = -1;
427 static real end = -1;
428 static real dt = 0;
429
430 t_pargs
431 pa[] =
432 {
433 { "-b", FALSE, etTIME,
434 { &begin }, "First time to use (%t)" },
435 { "-e", FALSE, etTIME,
436 { &end }, "Last time to use (%t)" },
437 { "-dt", FALSE, etTIME,
438 { &dt }, "Only write frame when t MOD dt = first time (%t)" },
439 { "-prec", FALSE, etINT,
440 { &prec }, "Precision for .xtc and .gro writing in number of decimal places" },
441 { "-vel", FALSE, etBOOL,
442 { &bVels }, "Read and write velocities if possible" },
443 { "-settime", FALSE, etBOOL,
444 { &bSetTime }, "Change starting time interactively" },
445 { "-sort", FALSE, etBOOL,
446 { &bSort }, "Sort trajectory files (not frames)" },
447 { "-keeplast", FALSE, etBOOL,
448 { &bKeepLast }, "keep overlapping frames at end of trajectory" },
449 { "-overwrite", FALSE, etBOOL,
450 { &bOverwrite }, "overwrite overlapping frames during appending" },
451 { "-cat", FALSE, etBOOL,
452 { &bCat }, "do not discard double time frames" } };
453 #define npargs asize(pa)
454 int status, ftpin, i, frame, frame_out, step = 0, trjout = 0;
455 rvec *x, *v;
456 real xtcpr, t_corr;
457 t_trxframe fr, frout;
458 char **fnms, **fnms_out, *in_file, *out_file;
459 int n_append;
460 int trxout = -1;
461 bool bNewFile, bIndex, bWrite;
462 int earliersteps, nfile_in, nfile_out, *cont_type, last_ok_step;
463 real *readtime, *timest, *settime;
464 real first_time = 0, lasttime = NOTSET, last_ok_t = -1, timestep;
465 real last_frame_time, searchtime;
466 int isize, j;
467 atom_id *index = NULL, imax;
468 char *grpname;
469 real **val = NULL, *t = NULL, dt_remd;
470 int n, nset;
471 bool bOK;
472 off_t fpos;
473 output_env_t oenv;
474 t_filenm fnm[] =
475 {
476 { efTRX, "-f", NULL, ffRDMULT },
477 { efTRO, "-o", NULL, ffWRMULT },
478 { efNDX, "-n", "index", ffOPTRD },
479 { efXVG, "-demux", "remd", ffOPTRD } };
480
481 #define NFILE asize(fnm)
482
483 CopyRight(stderr, argv[0]);
484 parse_common_args(&argc, argv, PCA_BE_NICE | PCA_TIME_UNIT, NFILE, fnm,
485 asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
486
487 bIndex = ftp2bSet(efNDX, NFILE, fnm);
488 bDeMux = ftp2bSet(efXVG, NFILE, fnm);
489 bSort = bSort && !bDeMux;
490
491 imax = NO_ATID;
492 if (bIndex)
493 {
494 printf("Select group for output\n");
495 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
496 /* scan index */
497 imax = index[0];
498 for (i = 1; i < isize; i++)
499 {
500 imax = max(imax, index[i]);
501 }
502 }
503 if (bDeMux)
504 {
505 nset = 0;
506 dt_remd = 0;
507 val = read_xvg_time(opt2fn("-demux", NFILE, fnm), TRUE,
508 opt2parg_bSet("-b", npargs, pa), begin,
509 opt2parg_bSet("-e", npargs, pa), end, 1, &nset, &n,
510 &dt_remd, &t);
511 printf("Read %d sets of %d points, dt = %g\n\n", nset, n, dt_remd);
512 if (debug)
513 {
514 fprintf(debug, "Dump of replica_index.xvg\n");
515 for (i = 0; (i < n); i++)
516 {
517 fprintf(debug, "%10g", t[i]);
518 for (j = 0; (j < nset); j++)
519 {
520 fprintf(debug, " %3d", gmx_nint(val[j][i]));
521 }
522 fprintf(debug, "\n");
523 }
524 }
525 }
526 /* prec is in nr of decimal places, xtcprec is a multiplication factor: */
527 xtcpr = 1;
528 for (i = 0; i < prec; i++)
529 {
530 xtcpr *= 10;
531 }
532
533 nfile_in = opt2fns(&fnms, "-f", NFILE, fnm);
534 if (!nfile_in)
535 {
536 gmx_fatal(FARGS,"No input files!" );
537 }
538
539 if (bDeMux && (nfile_in != nset))
540 {
541 gmx_fatal(FARGS,"You have specified %d files and %d entries in the demux table",nfile_in,nset);
542 }
543
544 nfile_out = opt2fns(&fnms_out,"-o",NFILE,fnm);
545 if (!nfile_out)
546 {
547 gmx_fatal(FARGS,"No output files!");
548 }
549 if ((nfile_out > 1) && !bDeMux)
550 {
551 gmx_fatal(FARGS,"Don't know what to do with more than 1 output file if not demultiplexing");
552 }
553 else if (bDeMux && (nfile_out != nset) && (nfile_out != 1))
554 {
555 gmx_fatal(FARGS,"Number of output files should be 1 or %d (#input files), not %d",nset,nfile_out);
556 }
557 if (bDeMux)
558 {
559 if (nfile_out != nset)
560 {
561 char *buf = strdup(fnms_out[0]);
562 snew(fnms_out,nset);
563 for(i=0; (i<nset); i++)
564 {
565 snew(fnms_out[i],strlen(buf)+32);
566 sprintf(fnms_out[i],"%d_%s",i,buf);
567 }
568 }
569 do_demux(nfile_in,fnms,fnms_out,n,val,t,dt_remd,isize,index,dt,oenv);
570 }
571 else
572 {
573 snew(readtime,nfile_in+1);
574 snew(timest,nfile_in+1);
575 scan_trj_files(fnms,nfile_in,readtime,timest,imax,oenv);
576
577 snew(settime,nfile_in+1);
578 snew(cont_type,nfile_in+1);
579 edit_files(fnms,nfile_in,readtime,timest,settime,cont_type,bSetTime,bSort,
580 oenv);
581
582 /* Check whether the output file is amongst the input files
583 * This has to be done after sorting etc.
584 */
585 out_file = fnms_out[0];
586 n_append = -1;
587 for(i=0; ((i<nfile_in) && (n_append==-1)); i++)
588 {
589 if (strcmp(fnms[i],out_file) == 0)
590 {
591 n_append = i;
592 }
593 }
594 if (n_append == 0) {
595 fprintf(stderr,"Will append to %s rather than creating a new file\n",
596 out_file);
597 }
598 else if (n_append != -1)
599 {
600 gmx_fatal(FARGS,"Can only append to the first file which is %s (not %s)",
601 fnms[0],out_file);
602 }
603 earliersteps=0;
604
605 /* Not checking input format, could be dangerous :-) */
606 /* Not checking output format, equally dangerous :-) */
607
608 frame=-1;
609 frame_out=-1;
610 /* the default is not to change the time at all,
611 * but this is overridden by the edit_files routine
612 */
613 t_corr=0;
614
615 if (n_append == -1)
616 {
617 trxout = open_trx(out_file,"w");
618 memset(&frout,0,sizeof(frout));
619 }
620 else
621 {
622 if (!read_first_frame(oenv,&status,out_file,&fr,FLAGS))
623 gmx_fatal(FARGS,"Reading first frame from %s",out_file);
624 if (!bKeepLast && !bOverwrite)
625 {
626 fprintf(stderr, "\n\nWARNING: Appending without -overwrite implies -keeplast "
627 "between the first two files. \n"
628 "If the trajectories have an overlap and have not been written binary \n"
629 "reproducible this will produce an incorrect trajectory!\n\n");
630
631 /* Fails if last frame is incomplete
632 * We can't do anything about it without overwriting
633 * */
634 if (gmx_fio_getftp(status) == efXTC)
635 {
636 lasttime = xdr_xtc_get_last_frame_time(gmx_fio_getfp(status),gmx_fio_getxdr(status),fr.natoms,&bOK);
637 fr.time = lasttime;
638 if (!bOK)
639 {
640 gmx_fatal(FARGS,"Error reading last frame. Maybe seek not supported." );
641 }
642 }
643 else
644 {
645 while (read_next_frame(oenv,status,&fr))
646 ;
647 lasttime = fr.time;
648 }
649 bKeepLastAppend = TRUE;
650 close_trj(status);
651 trxout = open_trx(out_file,"a");
652 }
653 else if (bOverwrite)
654 {
655 if (gmx_fio_getftp(status) != efXTC) {
656 gmx_fatal(FARGS,"Overwrite only supported for XTC." );
657 }
658 last_frame_time = xdr_xtc_get_last_frame_time(gmx_fio_getfp(status),gmx_fio_getxdr(status),fr.natoms,&bOK);
659 if (!bOK)
660 {
661 gmx_fatal(FARGS,"Error reading last frame. Maybe seek not supported." );
662 }
663 /* xtc_seek_time broken for trajectories containing only 1 or 2 frames
664 * or when seek time = 0 */
665 if (nfile_in>1 && settime[1]<last_frame_time+timest[0]*0.5)
666 {
667 /* Jump to one time-frame before the start of next
668 * trajectory file */
669 searchtime = settime[1]-timest[0]*1.25;
670 }
671 else
672 {
673 searchtime = last_frame_time;
674 }
675 if (xtc_seek_time(searchtime,status,fr.natoms))
676 {
677 gmx_fatal(FARGS,"Error seeking to append position.");
678 }
679 read_next_frame(oenv,status,&fr);
680 if (fabs(searchtime - fr.time) > timest[0]*0.5)
681 {
682 gmx_fatal(FARGS,"Error seeking: attempted to seek to %f but got %f.",
683 searchtime,fr.time);
684 }
685 lasttime = fr.time;
686 fpos = gmx_fio_ftell(status);
687 close_trj(status);
688 trxout = open_trx(out_file,"r+");
689 if (gmx_fio_seek(trxout,fpos)) {
690 gmx_fatal(FARGS,"Error seeking to append position.");
691 }
692 }
693 printf("\n Will append after %f \n",lasttime);
694 frout = fr;
695 }
696 /* Lets stitch up some files */
697 timestep = timest[0];
698 for(i=n_append+1; (i<nfile_in); i++)
699 {
700 /* Open next file */
701
702 /* set the next time from the last frame in previous file */
703 if (i > 0)
704 {
705 if (frame_out >= 0)
706 {
707 if(cont_type[i]==TIME_CONTINUE)
708 {
709 begin =frout.time;
710 begin += 0.5*timestep;
711 settime[i]=frout.time;
712 cont_type[i]=TIME_EXPLICIT;
713 }
714 else if(cont_type[i]==TIME_LAST)
715 {
716 begin=frout.time;
717 begin += 0.5*timestep;
718 }
719 /* Or, if the time in the next part should be changed by the
720 * same amount, start at half a timestep from the last time
721 * so we dont repeat frames.
722 */
723 /* I don't understand the comment above, but for all the cases
724 * I tried the code seems to work properly. B. Hess 2008-4-2.
725 */
726 }
727 /* Or, if time is set explicitly, we check for overlap/gap */
728 if(cont_type[i]==TIME_EXPLICIT)
729 {
730 if( ( i < nfile_in ) &&
731 ( frout.time < settime[i]-1.5*timestep ) )
732 {
733 fprintf(stderr, "WARNING: Frames around t=%f %s have a different "
734 "spacing than the rest,\n"
735 "might be a gap or overlap that couldn't be corrected "
736 "automatically.\n",output_env_conv_time(oenv,frout.time),
737 output_env_get_time_unit(oenv));
738 }
739 }
740 }
741
742 /* if we don't have a timestep in the current file, use the old one */
743 if ( timest[i] != 0 )
744 {
745 timestep = timest[i];
746 }
747 read_first_frame(oenv,&status,fnms[i],&fr,FLAGS);
748 if(!fr.bTime)
749 {
750 fr.time=0;
751 fprintf(stderr,"\nWARNING: Couldn't find a time in the frame.\n");
752 }
753
754 if(cont_type[i]==TIME_EXPLICIT)
755 {
756 t_corr=settime[i]-fr.time;
757 }
758 /* t_corr is the amount we want to change the time.
759 * If the user has chosen not to change the time for
760 * this part of the trajectory t_corr remains at
761 * the value it had in the last part, changing this
762 * by the same amount.
763 * If no value was given for the first trajectory part
764 * we let the time start at zero, see the edit_files routine.
765 */
766
767 bNewFile=TRUE;
768
769 printf("\n");
770 if (lasttime != NOTSET)
771 {
772 printf("lasttime %g\n", lasttime);
773 }
774
775 do
776 {
777 /* copy the input frame to the output frame */
778 frout=fr;
779 /* set the new time by adding the correct calculated above */
780 frout.time += t_corr;
781 /* quit if we have reached the end of what should be written */
782 if((end > 0) && (frout.time > end+GMX_REAL_EPS))
783 {
784 i=nfile_in;
785 break;
786 }
787
788 /* determine if we should write this frame (dt is handled elsewhere) */
789 if (bCat) /* write all frames of all files */
790 {
791 bWrite = TRUE;
792 }
793 else if ( bKeepLast || (bKeepLastAppend && i==1))
794 /* write till last frame of this traj
795 and skip first frame(s) of next traj */
796 {
797 bWrite = ( frout.time > lasttime+0.5*timestep );
798 }
799 else /* write till first frame of next traj */
800 {
801 bWrite = ( frout.time < settime[i+1]-0.5*timestep );
802 }
803
804 if( bWrite && (frout.time >= begin) )
805 {
806 frame++;
807 if (frame_out == -1)
808 first_time = frout.time;
809 lasttime = frout.time;
810 if (dt==0 || bRmod(frout.time,first_time,dt))
811 {
812 frame_out++;
813 last_ok_t=frout.time;
814 if(bNewFile)
815 {
816 fprintf(stderr,"\nContinue writing frames from %s t=%g %s, "
817 "frame=%d \n",
818 fnms[i],output_env_conv_time(oenv,frout.time),output_env_get_time_unit(oenv),
819 frame);
820 bNewFile=FALSE;
821 }
822
823 if (bIndex)
824 {
825 write_trxframe_indexed(trxout,&frout,isize,index,NULL);
826 }
827 else
828 {
829 write_trxframe(trxout,&frout,NULL);
830 }
831 if ( ((frame % 10) == 0) || (frame < 10) )
832 {
833 fprintf(stderr," -> frame %6d time %8.3f %s \r",
834 frame_out,output_env_conv_time(oenv,frout.time),output_env_get_time_unit(oenv));
835 }
836 }
837 }
838 } while( read_next_frame(oenv,status,&fr));
839
840 close_trj(status);
841
842 earliersteps+=step;
843 }
844 if (trxout >= 0)
845 {
846 close_trx(trxout);
847 }
848 fprintf(stderr,"\nLast frame written was %d, time %f %s\n",
849 frame,output_env_conv_time(oenv,last_ok_t),output_env_get_time_unit(oenv));
850 }
851 thanx(stderr);
852
853 return 0;
854 }