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