src/gromacs/tools/trjcat.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "trjcat.h"
  41
  42 #include <cmath>
  43 #include <cstdlib>
  44 #include <cstring>
  45
  46 #include <algorithm>
  47 #include <string>
  48
  49 #include "gromacs/commandline/pargs.h"
  50 #include "gromacs/fileio/confio.h"
  51 #include "gromacs/fileio/gmxfio.h"
  52 #include "gromacs/fileio/pdbio.h"
  53 #include "gromacs/fileio/tngio.h"
  54 #include "gromacs/fileio/trxio.h"
  55 #include "gromacs/fileio/xtcio.h"
  56 #include "gromacs/fileio/xvgr.h"
  57 #include "gromacs/math/vec.h"
  58 #include "gromacs/mdtypes/md_enums.h"
  59 #include "gromacs/topology/index.h"
  60 #include "gromacs/trajectory/trajectoryframe.h"
  61 #include "gromacs/utility/arrayref.h"
  62 #include "gromacs/utility/arraysize.h"
  63 #include "gromacs/utility/fatalerror.h"
  64 #include "gromacs/utility/futil.h"
  65 #include "gromacs/utility/smalloc.h"
  66 #include "gromacs/utility/stringutil.h"
  67
  68 #define TIME_EXPLICIT 0
  69 #define TIME_CONTINUE 1
  70 #define TIME_LAST 2
  71 #ifndef FLT_MAX
  72 #    define FLT_MAX 1e36
  73 #endif
  74 #define FLAGS (TRX_READ_X | TRX_READ_V | TRX_READ_F)
  75
  76 static void scan_trj_files(gmx::ArrayRef<const std::string> files,
  77                            real*                            readtime,
  78                            real*                            timestep,
  79                            int                              imax,
  80                            const gmx_output_env_t*          oenv)
  81 {
  82     /* Check start time of all files */
  83     int          natoms = 0;
  84     t_trxstatus* status;
  85     t_trxframe   fr;
  86     bool         ok;
  87
  88     for (gmx::index i = 0; i < files.ssize(); i++)
  89     {
  90         ok = read_first_frame(oenv, &status, files[i].c_str(), &fr, FLAGS);
  91
  92         if (!ok)
  93         {
  94             gmx_fatal(FARGS, "\nCouldn't read frame from file.");
  95         }
  96         if (fr.bTime)
  97         {
  98             readtime[i] = fr.time;
  99         }
 100         else
 101         {
 102             readtime[i] = 0;
 103             fprintf(stderr, "\nWARNING: Couldn't find a time in the frame.\n");
 104         }
 105
 106         if (i == 0)
 107         {
 108             natoms = fr.natoms;
 109         }
 110         else
 111         {
 112             if (imax == -1)
 113             {
 114                 if (natoms != fr.natoms)
 115                 {
 116                     gmx_fatal(FARGS, "\nDifferent numbers of atoms (%d/%d) in files", natoms, fr.natoms);
 117                 }
 118             }
 119             else
 120             {
 121                 if (fr.natoms <= imax)
 122                 {
 123                     gmx_fatal(FARGS, "\nNot enough atoms (%d) for index group (%d)", fr.natoms, imax);
 124                 }
 125             }
 126         }
 127         ok = read_next_frame(oenv, status, &fr);
 128         if (ok && fr.bTime)
 129         {
 130             timestep[i] = fr.time - readtime[i];
 131         }
 132         else
 133         {
 134             timestep[i] = 0;
 135         }
 136
 137         close_trx(status);
 138         if (fr.bX)
 139         {
 140             sfree(fr.x);
 141         }
 142         if (fr.bV)
 143         {
 144             sfree(fr.v);
 145         }
 146         if (fr.bF)
 147         {
 148             sfree(fr.f);
 149         }
 150     }
 151     fprintf(stderr, "\n");
 152 }
 153
 154 static void sort_files(gmx::ArrayRef<std::string> files, real* settime)
 155 {
 156     for (gmx::index i = 0; i < files.ssize(); i++)
 157     {
 158         gmx::index minidx = i;
 159         for (gmx::index j = i + 1; j < files.ssize(); j++)
 160         {
 161             if (settime[j] < settime[minidx])
 162             {
 163                 minidx = j;
 164             }
 165         }
 166         if (minidx != i)
 167         {
 168             real timeswap   = settime[i];
 169             settime[i]      = settime[minidx];
 170             settime[minidx] = timeswap;
 171             std::swap(files[i], files[minidx]);
 172         }
 173     }
 174 }
 175
 176 static void edit_files(gmx::ArrayRef<std::string> files,
 177                        real*                      readtime,
 178                        real*                      timestep,
 179                        real*                      settime,
 180                        int*                       cont_type,
 181                        gmx_bool                   bSetTime,
 182                        gmx_bool                   bSort,
 183                        const gmx_output_env_t*    oenv)
 184 {
 185     gmx_bool ok;
 186     char     inputstring[STRLEN], *chptr;
 187
 188     auto timeUnit = output_env_get_time_unit(oenv);
 189     if (bSetTime)
 190     {
 191         fprintf(stderr,
 192                 "\n\nEnter the new start time (%s) for each file.\n"
 193                 "There are two special options, both disable sorting:\n\n"
 194                 "c (continue) - The start time is taken from the end\n"
 195                 "of the previous file. Use it when your continuation run\n"
 196                 "restarts with t=0.\n\n"
 197                 "l (last) - The time in this file will be changed the\n"
 198                 "same amount as in the previous. Use it when the time in the\n"
 199                 "new run continues from the end of the previous one,\n"
 200                 "since this takes possible overlap into account.\n\n",
 201                 timeUnit.c_str());
 202
 203         fprintf(stderr,
 204                 "          File             Current start (%s)  New start (%s)\n"
 205                 "---------------------------------------------------------\n",
 206                 timeUnit.c_str(), timeUnit.c_str());
 207
 208         for (gmx::index i = 0; i < files.ssize(); i++)
 209         {
 210             fprintf(stderr, "%25s   %10.3f %s          ", files[i].c_str(),
 211                     output_env_conv_time(oenv, readtime[i]), timeUnit.c_str());
 212             ok = FALSE;
 213             do
 214             {
 215                 if (nullptr == fgets(inputstring, STRLEN - 1, stdin))
 216                 {
 217                     gmx_fatal(FARGS, "Error reading user input");
 218                 }
 219
 220                 inputstring[std::strlen(inputstring) - 1] = 0;
 221
 222                 if (inputstring[0] == 'c' || inputstring[0] == 'C')
 223                 {
 224                     cont_type[i] = TIME_CONTINUE;
 225                     bSort        = FALSE;
 226                     ok           = TRUE;
 227                     settime[i]   = FLT_MAX;
 228                 }
 229                 else if (inputstring[0] == 'l' || inputstring[0] == 'L')
 230                 {
 231                     cont_type[i] = TIME_LAST;
 232                     bSort        = FALSE;
 233                     ok           = TRUE;
 234                     settime[i]   = FLT_MAX;
 235                 }
 236                 else
 237                 {
 238                     settime[i] = strtod(inputstring, &chptr) * output_env_get_time_invfactor(oenv);
 239                     if (chptr == inputstring)
 240                     {
 241                         fprintf(stderr,
 242                                 "'%s' not recognized as a floating point number, 'c' or 'l'. "
 243                                 "Try again: ",
 244                                 inputstring);
 245                     }
 246                     else
 247                     {
 248                         cont_type[i] = TIME_EXPLICIT;
 249                         ok           = TRUE;
 250                     }
 251                 }
 252             } while (!ok);
 253         }
 254         if (cont_type[0] != TIME_EXPLICIT)
 255         {
 256             cont_type[0] = TIME_EXPLICIT;
 257             settime[0]   = 0;
 258         }
 259     }
 260     else
 261     {
 262         for (gmx::index i = 0; i < files.ssize(); i++)
 263         {
 264             settime[i] = readtime[i];
 265         }
 266     }
 267     if (!bSort)
 268     {
 269         fprintf(stderr, "Sorting disabled.\n");
 270     }
 271     else
 272     {
 273         sort_files(files, settime);
 274     }
 275     /* Write out the new order and start times */
 276     fprintf(stderr,
 277             "\nSummary of files and start times used:\n\n"
 278             "          File                Start time       Time step\n"
 279             "---------------------------------------------------------\n");
 280     for (gmx::index i = 0; i < files.ssize(); i++)
 281     {
 282         switch (cont_type[i])
 283         {
 284             case TIME_EXPLICIT:
 285                 fprintf(stderr, "%25s   %10.3f %s   %10.3f %s", files[i].c_str(),
 286                         output_env_conv_time(oenv, settime[i]), timeUnit.c_str(),
 287                         output_env_conv_time(oenv, timestep[i]), timeUnit.c_str());
 288                 if (i > 0 && cont_type[i - 1] == TIME_EXPLICIT && settime[i] == settime[i - 1])
 289                 {
 290                     fprintf(stderr, " WARNING: same Start time as previous");
 291                 }
 292                 fprintf(stderr, "\n");
 293                 break;
 294             case TIME_CONTINUE:
 295                 fprintf(stderr, "%25s        Continue from last file\n", files[i].c_str());
 296                 break;
 297             case TIME_LAST:
 298                 fprintf(stderr, "%25s        Change by same amount as last file\n", files[i].c_str());
 299                 break;
 300         }
 301     }
 302     fprintf(stderr, "\n");
 303
 304     settime[files.size()]   = FLT_MAX;
 305     cont_type[files.size()] = TIME_EXPLICIT;
 306     readtime[files.size()]  = FLT_MAX;
 307 }
 308
 309 static void do_demux(gmx::ArrayRef<const std::string> inFiles,
 310                      gmx::ArrayRef<const std::string> outFiles,
 311                      int                              nval,
 312                      real**                           value,
 313                      real*                            time,
 314                      real                             dt_remd,
 315                      int                              isize,
 316                      int                              index[],
 317                      real                             dt,
 318                      const gmx_output_env_t*          oenv)
 319 {
 320     int           k, natoms;
 321     t_trxstatus **fp_in, **fp_out;
 322     gmx_bool      bCont, *bSet;
 323     real          t, first_time = 0;
 324     t_trxframe*   trx;
 325
 326     snew(fp_in, inFiles.size());
 327     snew(trx, inFiles.size());
 328     snew(bSet, inFiles.size());
 329     natoms = -1;
 330     t      = -1;
 331     for (gmx::index i = 0; i < inFiles.ssize(); i++)
 332     {
 333         read_first_frame(oenv, &(fp_in[i]), inFiles[i].c_str(), &(trx[i]), TRX_NEED_X);
 334         if (natoms == -1)
 335         {
 336             natoms     = trx[i].natoms;
 337             first_time = trx[i].time;
 338         }
 339         else if (natoms != trx[i].natoms)
 340         {
 341             gmx_fatal(FARGS, "Trajectory file %s has %d atoms while previous trajs had %d atoms",
 342                       inFiles[i].c_str(), trx[i].natoms, natoms);
 343         }
 344         if (t == -1)
 345         {
 346             t = trx[i].time;
 347         }
 348         else if (t != trx[i].time)
 349         {
 350             gmx_fatal(FARGS, "Trajectory file %s has time %f while previous trajs had time %f",
 351                       inFiles[i].c_str(), trx[i].time, t);
 352         }
 353     }
 354
 355     snew(fp_out, inFiles.size());
 356     for (gmx::index i = 0; i < inFiles.ssize(); i++)
 357     {
 358         fp_out[i] = open_trx(outFiles[i].c_str(), "w");
 359     }
 360     k = 0;
 361     if (std::round(time[k] - t) != 0)
 362     {
 363         gmx_fatal(FARGS, "First time in demuxing table does not match trajectories");
 364     }
 365     do
 366     {
 367         while ((k + 1 < nval) && ((trx[0].time - time[k + 1]) > dt_remd * 0.1))
 368         {
 369             k++;
 370         }
 371         if (debug)
 372         {
 373             fprintf(debug, "trx[0].time = %g, time[k] = %g\n", trx[0].time, time[k]);
 374         }
 375         for (gmx::index i = 0; i < inFiles.ssize(); i++)
 376         {
 377             bSet[i] = FALSE;
 378         }
 379         for (gmx::index i = 0; i < inFiles.ssize(); i++)
 380         {
 381             int j = gmx::roundToInt(value[i][k]);
 382             range_check(j, 0, inFiles.size());
 383             if (bSet[j])
 384             {
 385                 gmx_fatal(FARGS, "Demuxing the same replica %d twice at time %f", j, trx[0].time);
 386             }
 387             bSet[j] = TRUE;
 388
 389             if (dt == 0 || bRmod(trx[i].time, first_time, dt))
 390             {
 391                 if (index)
 392                 {
 393                     write_trxframe_indexed(fp_out[j], &trx[i], isize, index, nullptr);
 394                 }
 395                 else
 396                 {
 397                     write_trxframe(fp_out[j], &trx[i], nullptr);
 398                 }
 399             }
 400         }
 401
 402         bCont = (k < nval);
 403         for (gmx::index i = 0; i < inFiles.ssize(); i++)
 404         {
 405             bCont = bCont && read_next_frame(oenv, fp_in[i], &trx[i]);
 406         }
 407     } while (bCont);
 408
 409     for (gmx::index i = 0; i < inFiles.ssize(); i++)
 410     {
 411         close_trx(fp_in[i]);
 412         close_trx(fp_out[i]);
 413     }
 414 }
 415
 416 int gmx_trjcat(int argc, char* argv[])
 417 {
 418     const char* desc[] = {
 419         "[THISMODULE] concatenates several input trajectory files in sorted order. ",
 420         "In case of double time frames the one in the later file is used. ",
 421         "By specifying [TT]-settime[tt] you will be asked for the start time ",
 422         "of each file. The input files are taken from the command line, ",
 423         "such that a command like [TT]gmx trjcat -f *.trr -o fixed.trr[tt] should do ",
 424         "the trick. Using [TT]-cat[tt], you can simply paste several files ",
 425         "together without removal of frames with identical time stamps.[PAR]",
 426         "One important option is inferred when the output file is amongst the",
 427         "input files. In that case that particular file will be appended to",
 428         "which implies you do not need to store double the amount of data.",
 429         "Obviously the file to append to has to be the one with lowest starting",
 430         "time since one can only append at the end of a file.[PAR]",
 431         "If the [TT]-demux[tt] option is given, the N trajectories that are",
 432         "read, are written in another order as specified in the [REF].xvg[ref] file.",
 433         "The [REF].xvg[ref] file should contain something like::",
 434         "",
 435         "    0  0  1  2  3  4  5",
 436         "    2  1  0  2  3  5  4",
 437         "",
 438         "The first number is the time, and subsequent numbers point to",
 439         "trajectory indices.",
 440         "The frames corresponding to the numbers present at the first line",
 441         "are collected into the output trajectory. If the number of frames in",
 442         "the trajectory does not match that in the [REF].xvg[ref] file then the program",
 443         "tries to be smart. Beware."
 444     };
 445     static gmx_bool bCat            = FALSE;
 446     static gmx_bool bSort           = TRUE;
 447     static gmx_bool bKeepLast       = FALSE;
 448     static gmx_bool bKeepLastAppend = FALSE;
 449     static gmx_bool bOverwrite      = FALSE;
 450     static gmx_bool bSetTime        = FALSE;
 451     static gmx_bool bDeMux;
 452     static real     begin = -1;
 453     static real     end   = -1;
 454     static real     dt    = 0;
 455
 456     t_pargs pa[] = {
 457         { "-b", FALSE, etTIME, { &begin }, "First time to use (%t)" },
 458         { "-e", FALSE, etTIME, { &end }, "Last time to use (%t)" },
 459         { "-dt", FALSE, etTIME, { &dt }, "Only write frame when t MOD dt = first time (%t)" },
 460         { "-settime", FALSE, etBOOL, { &bSetTime }, "Change starting time interactively" },
 461         { "-sort", FALSE, etBOOL, { &bSort }, "Sort trajectory files (not frames)" },
 462         { "-keeplast",
 463           FALSE,
 464           etBOOL,
 465           { &bKeepLast },
 466           "Keep overlapping frames at end of trajectory" },
 467         { "-overwrite",
 468           FALSE,
 469           etBOOL,
 470           { &bOverwrite },
 471           "Overwrite overlapping frames during appending" },
 472         { "-cat", FALSE, etBOOL, { &bCat }, "Do not discard double time frames" }
 473     };
 474 #define npargs asize(pa)
 475     int               ftpin, i, frame, frame_out;
 476     t_trxstatus *     status, *trxout = nullptr;
 477     real              t_corr;
 478     t_trxframe        fr, frout;
 479     int               n_append;
 480     gmx_bool          bNewFile, bIndex, bWrite;
 481     int*              cont_type;
 482     real *            readtime, *timest, *settime;
 483     real              first_time = 0, lasttime = 0, last_ok_t = -1, timestep;
 484     gmx_bool          lastTimeSet = FALSE;
 485     real              last_frame_time, searchtime;
 486     int               isize = 0, j;
 487     int *             index = nullptr, imax;
 488     char*             grpname;
 489     real **           val = nullptr, *t = nullptr, dt_remd;
 490     int               n, nset, ftpout = -1, prevEndStep = 0, filetype;
 491     gmx_off_t         fpos;
 492     gmx_output_env_t* oenv;
 493     t_filenm          fnm[] = { { efTRX, "-f", nullptr, ffRDMULT },
 494                        { efTRO, "-o", nullptr, ffWRMULT },
 495                        { efNDX, "-n", "index", ffOPTRD },
 496                        { efXVG, "-demux", "remd", ffOPTRD } };
 497
 498 #define NFILE asize(fnm)
 499
 500     if (!parse_common_args(&argc, argv, PCA_TIME_UNIT, NFILE, fnm, asize(pa), pa, asize(desc), desc,
 501                            0, nullptr, &oenv))
 502     {
 503         return 0;
 504     }
 505     fprintf(stdout,
 506             "Note that major changes are planned in future for "
 507             "trjcat, to improve usability and utility.");
 508
 509     auto timeUnit = output_env_get_time_unit(oenv);
 510
 511     bIndex = ftp2bSet(efNDX, NFILE, fnm);
 512     bDeMux = ftp2bSet(efXVG, NFILE, fnm);
 513     bSort  = bSort && !bDeMux;
 514
 515     imax = -1;
 516     if (bIndex)
 517     {
 518         printf("Select group for output\n");
 519         rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 520         /* scan index */
 521         imax = index[0];
 522         for (i = 1; i < isize; i++)
 523         {
 524             imax = std::max(imax, index[i]);
 525         }
 526     }
 527     if (bDeMux)
 528     {
 529         nset    = 0;
 530         dt_remd = 0;
 531         val     = read_xvg_time(opt2fn("-demux", NFILE, fnm), TRUE, opt2parg_bSet("-b", npargs, pa),
 532                             begin, opt2parg_bSet("-e", npargs, pa), end, 1, &nset, &n, &dt_remd, &t);
 533         printf("Read %d sets of %d points, dt = %g\n\n", nset, n, dt_remd);
 534         if (debug)
 535         {
 536             fprintf(debug, "Dump of replica_index.xvg\n");
 537             for (i = 0; (i < n); i++)
 538             {
 539                 fprintf(debug, "%10g", t[i]);
 540                 for (j = 0; (j < nset); j++)
 541                 {
 542                     fprintf(debug, "  %3d", static_cast<int>(std::round(val[j][i])));
 543                 }
 544                 fprintf(debug, "\n");
 545             }
 546         }
 547     }
 548
 549     gmx::ArrayRef<const std::string> inFiles = opt2fns("-f", NFILE, fnm);
 550     if (inFiles.empty())
 551     {
 552         gmx_fatal(FARGS, "No input files!");
 553     }
 554
 555     if (bDeMux && ssize(inFiles) != nset)
 556     {
 557         gmx_fatal(FARGS, "You have specified %td files and %d entries in the demux table",
 558                   inFiles.ssize(), nset);
 559     }
 560
 561     ftpin = fn2ftp(inFiles[0].c_str());
 562
 563     if (ftpin != efTRR && ftpin != efXTC && ftpin != efTNG)
 564     {
 565         gmx_fatal(FARGS, "gmx trjcat can only handle binary trajectory formats (trr, xtc, tng)");
 566     }
 567
 568     for (const std::string& inFile : inFiles)
 569     {
 570         if (ftpin != fn2ftp(inFile.c_str()))
 571         {
 572             gmx_fatal(FARGS, "All input files must be of the same (trr, xtc or tng) format");
 573         }
 574     }
 575
 576     gmx::ArrayRef<const std::string> outFiles = opt2fns("-o", NFILE, fnm);
 577     if (outFiles.empty())
 578     {
 579         gmx_fatal(FARGS, "No output files!");
 580     }
 581     if ((outFiles.size() > 1) && !bDeMux)
 582     {
 583         gmx_fatal(FARGS,
 584                   "Don't know what to do with more than 1 output file if  not demultiplexing");
 585     }
 586     else if (bDeMux && ssize(outFiles) != nset && outFiles.size() != 1)
 587     {
 588         gmx_fatal(FARGS, "Number of output files should be 1 or %d (#input files), not %td", nset,
 589                   outFiles.ssize());
 590     }
 591     if (bDeMux)
 592     {
 593         auto outFilesDemux = gmx::copyOf(outFiles);
 594         if (gmx::ssize(outFilesDemux) != nset)
 595         {
 596             std::string name = outFilesDemux[0];
 597             outFilesDemux.resize(nset);
 598             for (i = 0; (i < nset); i++)
 599             {
 600                 outFilesDemux[0] = gmx::formatString("%d_%s", i, name.c_str());
 601             }
 602         }
 603         do_demux(inFiles, outFilesDemux, n, val, t, dt_remd, isize, index, dt, oenv);
 604     }
 605     else
 606     {
 607         snew(readtime, inFiles.size() + 1);
 608         snew(timest, inFiles.size() + 1);
 609         scan_trj_files(inFiles, readtime, timest, imax, oenv);
 610
 611         snew(settime, inFiles.size() + 1);
 612         snew(cont_type, inFiles.size() + 1);
 613         auto inFilesEdited = gmx::copyOf(inFiles);
 614         edit_files(inFilesEdited, readtime, timest, settime, cont_type, bSetTime, bSort, oenv);
 615
 616         /* Check whether the output file is amongst the input files
 617          * This has to be done after sorting etc.
 618          */
 619         const char* out_file = outFiles[0].c_str();
 620         ftpout               = fn2ftp(out_file);
 621         n_append             = -1;
 622         for (size_t i = 0; i < inFilesEdited.size() && n_append == -1; i++)
 623         {
 624             if (std::strcmp(inFilesEdited[i].c_str(), out_file) == 0)
 625             {
 626                 n_append = i;
 627             }
 628         }
 629         if (n_append == 0)
 630         {
 631             fprintf(stderr, "Will append to %s rather than creating a new file\n", out_file);
 632         }
 633         else if (n_append != -1)
 634         {
 635             gmx_fatal(FARGS, "Can only append to the first file which is %s (not %s)",
 636                       inFilesEdited[0].c_str(), out_file);
 637         }
 638
 639         /* Not checking input format, could be dangerous :-) */
 640         /* Not checking output format, equally dangerous :-) */
 641
 642         frame     = -1;
 643         frame_out = -1;
 644         /* the default is not to change the time at all,
 645          * but this is overridden by the edit_files routine
 646          */
 647         t_corr = 0;
 648
 649         if (n_append == -1)
 650         {
 651             if (ftpout == efTNG)
 652             {
 653                 if (ftpout != ftpin)
 654                 {
 655                     gmx_fatal(FARGS, "When writing TNG the input file format must also be TNG");
 656                 }
 657                 if (bIndex)
 658                 {
 659                     trxout = trjtools_gmx_prepare_tng_writing(
 660                             out_file, 'w', nullptr, inFilesEdited[0].c_str(), isize, nullptr,
 661                             gmx::arrayRefFromArray(index, isize), grpname);
 662                 }
 663                 else
 664                 {
 665                     trxout = trjtools_gmx_prepare_tng_writing(
 666                             out_file, 'w', nullptr, inFilesEdited[0].c_str(), -1, nullptr, {}, nullptr);
 667                 }
 668             }
 669             else
 670             {
 671                 trxout = open_trx(out_file, "w");
 672             }
 673             std::memset(&frout, 0, sizeof(frout));
 674         }
 675         else
 676         {
 677             t_fileio* stfio;
 678
 679             if (!read_first_frame(oenv, &status, out_file, &fr, FLAGS))
 680             {
 681                 gmx_fatal(FARGS, "Reading first frame from %s", out_file);
 682             }
 683
 684             stfio = trx_get_fileio(status);
 685             if (!bKeepLast && !bOverwrite)
 686             {
 687                 fprintf(stderr,
 688                         "\n\nWARNING: Appending without -overwrite implies -keeplast "
 689                         "between the first two files. \n"
 690                         "If the trajectories have an overlap and have not been written binary \n"
 691                         "reproducible this will produce an incorrect trajectory!\n\n");
 692
 693                 filetype = gmx_fio_getftp(stfio);
 694                 /* Fails if last frame is incomplete
 695                  * We can't do anything about it without overwriting
 696                  * */
 697                 if (filetype == efXTC || filetype == efTNG)
 698                 {
 699                     lasttime = trx_get_time_of_final_frame(status);
 700                     fr.time  = lasttime;
 701                 }
 702                 else
 703                 {
 704                     while (read_next_frame(oenv, status, &fr)) {}
 705                     lasttime = fr.time;
 706                 }
 707                 lastTimeSet     = TRUE;
 708                 bKeepLastAppend = TRUE;
 709                 close_trx(status);
 710                 trxout = open_trx(out_file, "a");
 711             }
 712             else if (bOverwrite)
 713             {
 714                 if (gmx_fio_getftp(stfio) != efXTC)
 715                 {
 716                     gmx_fatal(FARGS, "Overwrite only supported for XTC.");
 717                 }
 718                 last_frame_time = trx_get_time_of_final_frame(status);
 719
 720                 /* xtc_seek_time broken for trajectories containing only 1 or 2 frames
 721                  *     or when seek time = 0 */
 722                 if (inFilesEdited.size() > 1 && settime[1] < last_frame_time + timest[0] * 0.5)
 723                 {
 724                     /* Jump to one time-frame before the start of next
 725                      *  trajectory file */
 726                     searchtime = settime[1] - timest[0] * 1.25;
 727                 }
 728                 else
 729                 {
 730                     searchtime = last_frame_time;
 731                 }
 732                 if (xtc_seek_time(stfio, searchtime, fr.natoms, TRUE))
 733                 {
 734                     gmx_fatal(FARGS, "Error seeking to append position.");
 735                 }
 736                 read_next_frame(oenv, status, &fr);
 737                 if (std::abs(searchtime - fr.time) > timest[0] * 0.5)
 738                 {
 739                     gmx_fatal(FARGS, "Error seeking: attempted to seek to %f but got %f.",
 740                               searchtime, fr.time);
 741                 }
 742                 lasttime    = fr.time;
 743                 lastTimeSet = TRUE;
 744                 fpos        = gmx_fio_ftell(stfio);
 745                 close_trx(status);
 746                 trxout = open_trx(out_file, "r+");
 747                 if (gmx_fio_seek(trx_get_fileio(trxout), fpos))
 748                 {
 749                     gmx_fatal(FARGS, "Error seeking to append position.");
 750                 }
 751             }
 752             if (lastTimeSet)
 753             {
 754                 printf("\n Will append after %f \n", lasttime);
 755             }
 756             frout = fr;
 757         }
 758         /* Lets stitch up some files */
 759         timestep = timest[0];
 760         for (size_t i = n_append + 1; i < inFilesEdited.size(); i++)
 761         {
 762             /* Open next file */
 763
 764             /* set the next time from the last frame in previous file */
 765             if (i > 0)
 766             {
 767                 /* When writing TNG the step determine which frame to write. Use an
 768                  * offset to be able to increase steps properly when changing files. */
 769                 if (ftpout == efTNG)
 770                 {
 771                     prevEndStep = frout.step;
 772                 }
 773
 774                 if (frame_out >= 0)
 775                 {
 776                     if (cont_type[i] == TIME_CONTINUE)
 777                     {
 778                         begin = frout.time;
 779                         begin += 0.5 * timestep;
 780                         settime[i]   = frout.time;
 781                         cont_type[i] = TIME_EXPLICIT;
 782                     }
 783                     else if (cont_type[i] == TIME_LAST)
 784                     {
 785                         begin = frout.time;
 786                         begin += 0.5 * timestep;
 787                     }
 788                     /* Or, if the time in the next part should be changed by the
 789                      * same amount, start at half a timestep from the last time
 790                      * so we dont repeat frames.
 791                      */
 792                     /* I don't understand the comment above, but for all the cases
 793                      * I tried the code seems to work properly. B. Hess 2008-4-2.
 794                      */
 795                 }
 796                 /* Or, if time is set explicitly, we check for overlap/gap */
 797                 if (cont_type[i] == TIME_EXPLICIT)
 798                 {
 799                     if (i < inFilesEdited.size() && frout.time < settime[i] - 1.5 * timestep)
 800                     {
 801                         fprintf(stderr,
 802                                 "WARNING: Frames around t=%f %s have a different "
 803                                 "spacing than the rest,\n"
 804                                 "might be a gap or overlap that couldn't be corrected "
 805                                 "automatically.\n",
 806                                 output_env_conv_time(oenv, frout.time), timeUnit.c_str());
 807                     }
 808                 }
 809             }
 810
 811             /* if we don't have a timestep in the current file, use the old one */
 812             if (timest[i] != 0)
 813             {
 814                 timestep = timest[i];
 815             }
 816             read_first_frame(oenv, &status, inFilesEdited[i].c_str(), &fr, FLAGS);
 817             if (!fr.bTime)
 818             {
 819                 fr.time = 0;
 820                 fprintf(stderr, "\nWARNING: Couldn't find a time in the frame.\n");
 821             }
 822
 823             if (cont_type[i] == TIME_EXPLICIT)
 824             {
 825                 t_corr = settime[i] - fr.time;
 826             }
 827             /* t_corr is the amount we want to change the time.
 828              * If the user has chosen not to change the time for
 829              * this part of the trajectory t_corr remains at
 830              * the value it had in the last part, changing this
 831              * by the same amount.
 832              * If no value was given for the first trajectory part
 833              * we let the time start at zero, see the edit_files routine.
 834              */
 835
 836             bNewFile = TRUE;
 837
 838             if (!lastTimeSet)
 839             {
 840                 lasttime    = 0;
 841                 lastTimeSet = true;
 842             }
 843             printf("\n");
 844             printf("lasttime %g\n", lasttime);
 845
 846             do
 847             {
 848                 /* copy the input frame to the output frame */
 849                 frout = fr;
 850                 /* set the new time by adding the correct calculated above */
 851                 frout.time += t_corr;
 852                 if (ftpout == efTNG)
 853                 {
 854                     frout.step += prevEndStep;
 855                 }
 856                 /* quit if we have reached the end of what should be written */
 857                 if ((end > 0) && (frout.time > end + GMX_REAL_EPS))
 858                 {
 859                     i = inFilesEdited.size();
 860                     break;
 861                 }
 862
 863                 /* determine if we should write this frame (dt is handled elsewhere) */
 864                 if (bCat) /* write all frames of all files */
 865                 {
 866                     bWrite = TRUE;
 867                 }
 868                 else if (bKeepLast || (bKeepLastAppend && i == 1))
 869                 /* write till last frame of this traj
 870                    and skip first frame(s) of next traj */
 871                 {
 872                     bWrite = (frout.time > lasttime + 0.5 * timestep);
 873                 }
 874                 else /* write till first frame of next traj */
 875                 {
 876                     bWrite = (frout.time < settime[i + 1] - 0.5 * timestep);
 877                 }
 878
 879                 if (bWrite && (frout.time >= begin))
 880                 {
 881                     frame++;
 882                     if (frame_out == -1)
 883                     {
 884                         first_time = frout.time;
 885                     }
 886                     lasttime    = frout.time;
 887                     lastTimeSet = TRUE;
 888                     if (dt == 0 || bRmod(frout.time, first_time, dt))
 889                     {
 890                         frame_out++;
 891                         last_ok_t = frout.time;
 892                         if (bNewFile)
 893                         {
 894                             fprintf(stderr,
 895                                     "\nContinue writing frames from %s t=%g %s, "
 896                                     "frame=%d      \n",
 897                                     inFilesEdited[i].c_str(),
 898                                     output_env_conv_time(oenv, frout.time), timeUnit.c_str(), frame);
 899                             bNewFile = FALSE;
 900                         }
 901
 902                         if (bIndex)
 903                         {
 904                             write_trxframe_indexed(trxout, &frout, isize, index, nullptr);
 905                         }
 906                         else
 907                         {
 908                             write_trxframe(trxout, &frout, nullptr);
 909                         }
 910                         if (((frame % 10) == 0) || (frame < 10))
 911                         {
 912                             fprintf(stderr, " ->  frame %6d time %8.3f %s     \r", frame_out,
 913                                     output_env_conv_time(oenv, frout.time), timeUnit.c_str());
 914                             fflush(stderr);
 915                         }
 916                     }
 917                 }
 918             } while (read_next_frame(oenv, status, &fr));
 919
 920             close_trx(status);
 921         }
 922         if (trxout)
 923         {
 924             close_trx(trxout);
 925         }
 926         fprintf(stderr, "\nLast frame written was %d, time %f %s\n", frame,
 927                 output_env_conv_time(oenv, last_ok_t), timeUnit.c_str());
 928     }
 929
 930     return 0;
 931 }