src/tools/gmx_trjcat.c

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