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, 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 }