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