src/gromacs/gmxana/gmx_trjconv.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, 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 <cmath>
  40 #include <cstdlib>
  41 #include <cstring>
  42
  43 #include <algorithm>
  44
  45 #include "gromacs/commandline/pargs.h"
  46 #include "gromacs/commandline/viewit.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "gromacs/fileio/g96io.h"
  49 #include "gromacs/fileio/gmxfio.h"
  50 #include "gromacs/fileio/groio.h"
  51 #include "gromacs/fileio/pdbio.h"
  52 #include "gromacs/fileio/tngio.h"
  53 #include "gromacs/fileio/tpxio.h"
  54 #include "gromacs/fileio/trrio.h"
  55 #include "gromacs/fileio/trxio.h"
  56 #include "gromacs/fileio/xtcio.h"
  57 #include "gromacs/fileio/xvgr.h"
  58 #include "gromacs/gmxana/gmx_ana.h"
  59 #include "gromacs/math/do_fit.h"
  60 #include "gromacs/math/functions.h"
  61 #include "gromacs/math/vec.h"
  62 #include "gromacs/mdtypes/md_enums.h"
  63 #include "gromacs/pbcutil/pbc.h"
  64 #include "gromacs/pbcutil/rmpbc.h"
  65 #include "gromacs/topology/index.h"
  66 #include "gromacs/topology/topology.h"
  67 #include "gromacs/trajectory/trajectoryframe.h"
  68 #include "gromacs/utility/arrayref.h"
  69 #include "gromacs/utility/arraysize.h"
  70 #include "gromacs/utility/fatalerror.h"
  71 #include "gromacs/utility/futil.h"
  72 #include "gromacs/utility/smalloc.h"
  73
  74 enum {
  75     euSel, euRect, euTric, euCompact, euNR
  76 };
  77
  78
  79 static void calc_pbc_cluster(int ecenter, int nrefat, t_topology *top, int ePBC,
  80                              rvec x[], int index[], matrix box)
  81 {
  82     int       m, i, j, j0, j1, jj, ai, aj;
  83     int       imin, jmin;
  84     real      fac, min_dist2;
  85     rvec      dx, xtest, box_center;
  86     int       nmol, imol_center;
  87     int      *molind;
  88     gmx_bool *bMol, *bTmp;
  89     rvec     *m_com, *m_shift;
  90     t_pbc     pbc;
  91     int      *cluster;
  92     int      *added;
  93     int       ncluster, nadded;
  94     real      tmp_r2;
  95
  96     calc_box_center(ecenter, box, box_center);
  97
  98     /* Initiate the pbc structure */
  99     std::memset(&pbc, 0, sizeof(pbc));
 100     set_pbc(&pbc, ePBC, box);
 101
 102     /* Convert atom index to molecular */
 103     nmol   = top->mols.nr;
 104     molind = top->mols.index;
 105     snew(bMol, nmol);
 106     snew(m_com, nmol);
 107     snew(m_shift, nmol);
 108     snew(cluster, nmol);
 109     snew(added, nmol);
 110     snew(bTmp, top->atoms.nr);
 111
 112     for (i = 0; (i < nrefat); i++)
 113     {
 114         /* Mark all molecules in the index */
 115         ai       = index[i];
 116         bTmp[ai] = TRUE;
 117         /* Binary search assuming the molecules are sorted */
 118         j0 = 0;
 119         j1 = nmol-1;
 120         while (j0 < j1)
 121         {
 122             if (ai < molind[j0+1])
 123             {
 124                 j1 = j0;
 125             }
 126             else if (ai >= molind[j1])
 127             {
 128                 j0 = j1;
 129             }
 130             else
 131             {
 132                 jj = (j0+j1)/2;
 133                 if (ai < molind[jj+1])
 134                 {
 135                     j1 = jj;
 136                 }
 137                 else
 138                 {
 139                     j0 = jj;
 140                 }
 141             }
 142         }
 143         bMol[j0] = TRUE;
 144     }
 145     /* Double check whether all atoms in all molecules that are marked are part
 146      * of the cluster. Simultaneously compute the center of geometry.
 147      */
 148     min_dist2   = 10*gmx::square(trace(box));
 149     imol_center = -1;
 150     ncluster    = 0;
 151     for (i = 0; i < nmol; i++)
 152     {
 153         for (j = molind[i]; j < molind[i+1]; j++)
 154         {
 155             if (bMol[i] && !bTmp[j])
 156             {
 157                 gmx_fatal(FARGS, "Molecule %d marked for clustering but not atom %d in it - check your index!", i+1, j+1);
 158             }
 159             else if (!bMol[i] && bTmp[j])
 160             {
 161                 gmx_fatal(FARGS, "Atom %d marked for clustering but not molecule %d - this is an internal error...", j+1, i+1);
 162             }
 163             else if (bMol[i])
 164             {
 165                 /* Make molecule whole, move 2nd and higher atom to same periodicity as 1st atom in molecule */
 166                 if (j > molind[i])
 167                 {
 168                     pbc_dx(&pbc, x[j], x[j-1], dx);
 169                     rvec_add(x[j-1], dx, x[j]);
 170                 }
 171                 /* Compute center of geometry of molecule - m_com[i] was zeroed when we did snew() on it! */
 172                 rvec_inc(m_com[i], x[j]);
 173             }
 174         }
 175         if (bMol[i])
 176         {
 177             /* Normalize center of geometry */
 178             fac = 1.0/(molind[i+1]-molind[i]);
 179             for (m = 0; (m < DIM); m++)
 180             {
 181                 m_com[i][m] *= fac;
 182             }
 183             /* Determine which molecule is closest to the center of the box */
 184             pbc_dx(&pbc, box_center, m_com[i], dx);
 185             tmp_r2 = iprod(dx, dx);
 186
 187             if (tmp_r2 < min_dist2)
 188             {
 189                 min_dist2   = tmp_r2;
 190                 imol_center = i;
 191             }
 192             cluster[ncluster++] = i;
 193         }
 194     }
 195     sfree(bTmp);
 196
 197     if (ncluster <= 0)
 198     {
 199         fprintf(stderr, "No molecules selected in the cluster\n");
 200         return;
 201     }
 202     else if (imol_center == -1)
 203     {
 204         fprintf(stderr, "No central molecules could be found\n");
 205         return;
 206     }
 207
 208     nadded            = 0;
 209     added[nadded++]   = imol_center;
 210     bMol[imol_center] = FALSE;
 211
 212     while (nadded < ncluster)
 213     {
 214         /* Find min distance between cluster molecules and those remaining to be added */
 215         min_dist2   = 10*gmx::square(trace(box));
 216         imin        = -1;
 217         jmin        = -1;
 218         /* Loop over added mols */
 219         for (i = 0; i < nadded; i++)
 220         {
 221             ai = added[i];
 222             /* Loop over all mols */
 223             for (j = 0; j < ncluster; j++)
 224             {
 225                 aj = cluster[j];
 226                 /* check those remaining to be added */
 227                 if (bMol[aj])
 228                 {
 229                     pbc_dx(&pbc, m_com[aj], m_com[ai], dx);
 230                     tmp_r2 = iprod(dx, dx);
 231                     if (tmp_r2 < min_dist2)
 232                     {
 233                         min_dist2   = tmp_r2;
 234                         imin        = ai;
 235                         jmin        = aj;
 236                     }
 237                 }
 238             }
 239         }
 240
 241         /* Add the best molecule */
 242         added[nadded++]   = jmin;
 243         bMol[jmin]        = FALSE;
 244         /* Calculate the shift from the ai molecule */
 245         pbc_dx(&pbc, m_com[jmin], m_com[imin], dx);
 246         rvec_add(m_com[imin], dx, xtest);
 247         rvec_sub(xtest, m_com[jmin], m_shift[jmin]);
 248         rvec_inc(m_com[jmin], m_shift[jmin]);
 249
 250         for (j = molind[jmin]; j < molind[jmin+1]; j++)
 251         {
 252             rvec_inc(x[j], m_shift[jmin]);
 253         }
 254         fprintf(stdout, "\rClustering iteration %d of %d...", nadded, ncluster);
 255         fflush(stdout);
 256     }
 257
 258     sfree(added);
 259     sfree(cluster);
 260     sfree(bMol);
 261     sfree(m_com);
 262     sfree(m_shift);
 263
 264     fprintf(stdout, "\n");
 265 }
 266
 267 static void put_molecule_com_in_box(int unitcell_enum, int ecenter,
 268                                     t_block *mols,
 269                                     int natoms, t_atom atom[],
 270                                     int ePBC, matrix box, rvec x[])
 271 {
 272     int     i, j;
 273     int     d;
 274     rvec    com, shift, box_center;
 275     real    m;
 276     double  mtot;
 277     t_pbc   pbc;
 278
 279     calc_box_center(ecenter, box, box_center);
 280     set_pbc(&pbc, ePBC, box);
 281     if (mols->nr <= 0)
 282     {
 283         gmx_fatal(FARGS, "There are no molecule descriptions. I need a .tpr file for this pbc option.");
 284     }
 285     for (i = 0; (i < mols->nr); i++)
 286     {
 287         /* calc COM */
 288         clear_rvec(com);
 289         mtot = 0;
 290         for (j = mols->index[i]; (j < mols->index[i+1] && j < natoms); j++)
 291         {
 292             m = atom[j].m;
 293             for (d = 0; d < DIM; d++)
 294             {
 295                 com[d] += m*x[j][d];
 296             }
 297             mtot += m;
 298         }
 299         /* calculate final COM */
 300         svmul(1.0/mtot, com, com);
 301
 302         /* check if COM is outside box */
 303         gmx::RVec newCom;
 304         copy_rvec(com, newCom);
 305         auto      newComArrayRef = gmx::arrayRefFromArray(&newCom, 1);
 306         switch (unitcell_enum)
 307         {
 308             case euRect:
 309                 put_atoms_in_box(ePBC, box, newComArrayRef);
 310                 break;
 311             case euTric:
 312                 put_atoms_in_triclinic_unitcell(ecenter, box, newComArrayRef);
 313                 break;
 314             case euCompact:
 315                 put_atoms_in_compact_unitcell(ePBC, ecenter, box, newComArrayRef);
 316                 break;
 317         }
 318         rvec_sub(newCom, com, shift);
 319         if (norm2(shift) > 0)
 320         {
 321             if (debug)
 322             {
 323                 fprintf(debug, "\nShifting position of molecule %d "
 324                         "by %8.3f  %8.3f  %8.3f\n", i+1,
 325                         shift[XX], shift[YY], shift[ZZ]);
 326             }
 327             for (j = mols->index[i]; (j < mols->index[i+1] && j < natoms); j++)
 328             {
 329                 rvec_inc(x[j], shift);
 330             }
 331         }
 332     }
 333 }
 334
 335 static void put_residue_com_in_box(int unitcell_enum, int ecenter,
 336                                    int natoms, t_atom atom[],
 337                                    int ePBC, matrix box, rvec x[])
 338 {
 339     int              i, j, res_start, res_end;
 340     int              d, presnr;
 341     real             m;
 342     double           mtot;
 343     rvec             box_center, com, shift;
 344     static const int NOTSET = -12347;
 345     calc_box_center(ecenter, box, box_center);
 346
 347     presnr    = NOTSET;
 348     res_start = 0;
 349     clear_rvec(com);
 350     mtot = 0;
 351     for (i = 0; i < natoms+1; i++)
 352     {
 353         if (i == natoms || (presnr != atom[i].resind && presnr != NOTSET))
 354         {
 355             /* calculate final COM */
 356             res_end = i;
 357             svmul(1.0/mtot, com, com);
 358
 359             /* check if COM is outside box */
 360             gmx::RVec newCom;
 361             copy_rvec(com, newCom);
 362             auto      newComArrayRef = gmx::arrayRefFromArray(&newCom, 1);
 363             switch (unitcell_enum)
 364             {
 365                 case euRect:
 366                     put_atoms_in_box(ePBC, box, newComArrayRef);
 367                     break;
 368                 case euTric:
 369                     put_atoms_in_triclinic_unitcell(ecenter, box, newComArrayRef);
 370                     break;
 371                 case euCompact:
 372                     put_atoms_in_compact_unitcell(ePBC, ecenter, box, newComArrayRef);
 373                     break;
 374             }
 375             rvec_sub(newCom, com, shift);
 376             if (norm2(shift))
 377             {
 378                 if (debug)
 379                 {
 380                     fprintf(debug, "\nShifting position of residue %d (atoms %d-%d) "
 381                             "by %g,%g,%g\n", atom[res_start].resind+1,
 382                             res_start+1, res_end+1, shift[XX], shift[YY], shift[ZZ]);
 383                 }
 384                 for (j = res_start; j < res_end; j++)
 385                 {
 386                     rvec_inc(x[j], shift);
 387                 }
 388             }
 389             clear_rvec(com);
 390             mtot = 0;
 391
 392             /* remember start of new residue */
 393             res_start = i;
 394         }
 395         if (i < natoms)
 396         {
 397             /* calc COM */
 398             m = atom[i].m;
 399             for (d = 0; d < DIM; d++)
 400             {
 401                 com[d] += m*x[i][d];
 402             }
 403             mtot += m;
 404
 405             presnr = atom[i].resind;
 406         }
 407     }
 408 }
 409
 410 static void center_x(int ecenter, rvec x[], matrix box, int n, int nc, int ci[])
 411 {
 412     int  i, m, ai;
 413     rvec cmin, cmax, box_center, dx;
 414
 415     if (nc > 0)
 416     {
 417         copy_rvec(x[ci[0]], cmin);
 418         copy_rvec(x[ci[0]], cmax);
 419         for (i = 0; i < nc; i++)
 420         {
 421             ai = ci[i];
 422             for (m = 0; m < DIM; m++)
 423             {
 424                 if (x[ai][m] < cmin[m])
 425                 {
 426                     cmin[m] = x[ai][m];
 427                 }
 428                 else if (x[ai][m] > cmax[m])
 429                 {
 430                     cmax[m] = x[ai][m];
 431                 }
 432             }
 433         }
 434         calc_box_center(ecenter, box, box_center);
 435         for (m = 0; m < DIM; m++)
 436         {
 437             dx[m] = box_center[m]-(cmin[m]+cmax[m])*0.5;
 438         }
 439
 440         for (i = 0; i < n; i++)
 441         {
 442             rvec_inc(x[i], dx);
 443         }
 444     }
 445 }
 446
 447 static void mk_filenm(char *base, const char *ext, int ndigit, int file_nr,
 448                       char out_file[])
 449 {
 450     char nbuf[128];
 451     int  nd = 0, fnr;
 452
 453     std::strcpy(out_file, base);
 454     fnr = file_nr;
 455     do
 456     {
 457         fnr /= 10;
 458         nd++;
 459     }
 460     while (fnr > 0);
 461
 462     if (nd < ndigit)
 463     {
 464         std::strncat(out_file, "00000000000", ndigit-nd);
 465     }
 466     sprintf(nbuf, "%d.", file_nr);
 467     std::strcat(out_file, nbuf);
 468     std::strcat(out_file, ext);
 469 }
 470
 471 static void check_trr(const char *fn)
 472 {
 473     if (fn2ftp(fn) != efTRR)
 474     {
 475         gmx_fatal(FARGS, "%s is not a trajectory file, exiting\n", fn);
 476     }
 477 }
 478
 479 static void do_trunc(const char *fn, real t0)
 480 {
 481     t_fileio        *in;
 482     FILE            *fp;
 483     gmx_bool         bStop, bOK;
 484     gmx_trr_header_t sh;
 485     gmx_off_t        fpos;
 486     char             yesno[256];
 487     int              j;
 488     real             t = 0;
 489
 490     if (t0 == -1)
 491     {
 492         gmx_fatal(FARGS, "You forgot to set the truncation time");
 493     }
 494
 495     /* Check whether this is a .trr file */
 496     check_trr(fn);
 497
 498     in   = gmx_trr_open(fn, "r");
 499     fp   = gmx_fio_getfp(in);
 500     if (fp == nullptr)
 501     {
 502         fprintf(stderr, "Sorry, can not trunc %s, truncation of this filetype is not supported\n", fn);
 503         gmx_trr_close(in);
 504     }
 505     else
 506     {
 507         j     = 0;
 508         fpos  = gmx_fio_ftell(in);
 509         bStop = FALSE;
 510         while (!bStop && gmx_trr_read_frame_header(in, &sh, &bOK))
 511         {
 512             gmx_trr_read_frame_data(in, &sh, nullptr, nullptr, nullptr, nullptr);
 513             fpos = gmx_ftell(fp);
 514             t    = sh.t;
 515             if (t >= t0)
 516             {
 517                 gmx_fseek(fp, fpos, SEEK_SET);
 518                 bStop = TRUE;
 519             }
 520         }
 521         if (bStop)
 522         {
 523             fprintf(stderr, "Do you REALLY want to truncate this trajectory (%s) at:\n"
 524                     "frame %d, time %g, bytes %ld ??? (type YES if so)\n",
 525                     fn, j, t, (long int)fpos);
 526             if (1 != scanf("%s", yesno))
 527             {
 528                 gmx_fatal(FARGS, "Error reading user input");
 529             }
 530             if (std::strcmp(yesno, "YES") == 0)
 531             {
 532                 fprintf(stderr, "Once again, I'm gonna DO this...\n");
 533                 gmx_trr_close(in);
 534                 if (0 != gmx_truncate(fn, fpos))
 535                 {
 536                     gmx_fatal(FARGS, "Error truncating file %s", fn);
 537                 }
 538             }
 539             else
 540             {
 541                 fprintf(stderr, "Ok, I'll forget about it\n");
 542             }
 543         }
 544         else
 545         {
 546             fprintf(stderr, "Already at end of file (t=%g)...\n", t);
 547             gmx_trr_close(in);
 548         }
 549     }
 550 }
 551
 552 /*! \brief Read a full molecular topology if useful and available.
 553  *
 554  * If the input trajectory file is not in TNG format, and the output
 555  * file is in TNG format, then we want to try to read a full topology
 556  * (if available), so that we can write molecule information to the
 557  * output file. The full topology provides better molecule information
 558  * than is available from the normal t_topology data used by GROMACS
 559  * tools.
 560  *
 561  * Also, the t_topology is only read under (different) particular
 562  * conditions. If both apply, then a .tpr file might be read
 563  * twice. Trying to fix this redundancy while trjconv is still an
 564  * all-purpose tool does not seem worthwhile.
 565  *
 566  * Because of the way gmx_prepare_tng_writing is implemented, the case
 567  * where the input TNG file has no molecule information will never
 568  * lead to an output TNG file having molecule information. Since
 569  * molecule information will generally be present if the input TNG
 570  * file was written by a GROMACS tool, this seems like reasonable
 571  * behaviour. */
 572 static gmx_mtop_t *read_mtop_for_tng(const char *tps_file,
 573                                      const char *input_file,
 574                                      const char *output_file)
 575 {
 576     gmx_mtop_t *mtop = nullptr;
 577
 578     if (fn2bTPX(tps_file) &&
 579         efTNG != fn2ftp(input_file) &&
 580         efTNG == fn2ftp(output_file))
 581     {
 582         int temp_natoms = -1;
 583         snew(mtop, 1);
 584         read_tpx(tps_file, nullptr, nullptr, &temp_natoms,
 585                  nullptr, nullptr, mtop);
 586     }
 587
 588     return mtop;
 589 }
 590
 591 int gmx_trjconv(int argc, char *argv[])
 592 {
 593     const char *desc[] = {
 594         "[THISMODULE] can convert trajectory files in many ways:",
 595         "",
 596         "* from one format to another",
 597         "* select a subset of atoms",
 598         "* change the periodicity representation",
 599         "* keep multimeric molecules together",
 600         "* center atoms in the box",
 601         "* fit atoms to reference structure",
 602         "* reduce the number of frames",
 603         "* change the timestamps of the frames ([TT]-t0[tt] and [TT]-timestep[tt])",
 604         "* cut the trajectory in small subtrajectories according",
 605         "  to information in an index file. This allows subsequent analysis of",
 606         "  the subtrajectories that could, for example, be the result of a",
 607         "  cluster analysis. Use option [TT]-sub[tt].",
 608         "  This assumes that the entries in the index file are frame numbers and",
 609         "  dumps each group in the index file to a separate trajectory file.",
 610         "* select frames within a certain range of a quantity given",
 611         "  in an [REF].xvg[ref] file.",
 612         "",
 613         "[gmx-trjcat] is better suited for concatenating multiple trajectory files.",
 614         "[PAR]",
 615
 616         "The following formats are supported for input and output:",
 617         "[REF].xtc[ref], [REF].trr[ref], [REF].gro[ref], [TT].g96[tt]",
 618         "and [REF].pdb[ref].",
 619         "The file formats are detected from the file extension.",
 620         "The precision of [REF].xtc[ref] and [REF].gro[ref] output is taken from the",
 621         "input file for [REF].xtc[ref], [REF].gro[ref] and [REF].pdb[ref],",
 622         "and from the [TT]-ndec[tt] option for other input formats. The precision",
 623         "is always taken from [TT]-ndec[tt], when this option is set.",
 624         "All other formats have fixed precision. [REF].trr[ref]",
 625         "output can be single or double precision, depending on the precision",
 626         "of the [THISMODULE] binary.",
 627         "Note that velocities are only supported in",
 628         "[REF].trr[ref], [REF].gro[ref] and [TT].g96[tt] files.[PAR]",
 629
 630         "Option [TT]-sep[tt] can be used to write every frame to a separate",
 631         "[TT].gro, .g96[tt] or [REF].pdb[ref] file. By default, all frames all written to one file.",
 632         "[REF].pdb[ref] files with all frames concatenated can be viewed with",
 633         "[TT]rasmol -nmrpdb[tt].[PAR]",
 634
 635         "It is possible to select part of your trajectory and write it out",
 636         "to a new trajectory file in order to save disk space, e.g. for leaving",
 637         "out the water from a trajectory of a protein in water.",
 638         "[BB]ALWAYS[bb] put the original trajectory on tape!",
 639         "We recommend to use the portable [REF].xtc[ref] format for your analysis",
 640         "to save disk space and to have portable files.[PAR]",
 641
 642         "There are two options for fitting the trajectory to a reference",
 643         "either for essential dynamics analysis, etc.",
 644         "The first option is just plain fitting to a reference structure",
 645         "in the structure file. The second option is a progressive fit",
 646         "in which the first timeframe is fitted to the reference structure ",
 647         "in the structure file to obtain and each subsequent timeframe is ",
 648         "fitted to the previously fitted structure. This way a continuous",
 649         "trajectory is generated, which might not be the case when using the",
 650         "regular fit method, e.g. when your protein undergoes large",
 651         "conformational transitions.[PAR]",
 652
 653         "Option [TT]-pbc[tt] sets the type of periodic boundary condition",
 654         "treatment:",
 655         "",
 656         " * [TT]mol[tt] puts the center of mass of molecules in the box,",
 657         "   and requires a run input file to be supplied with [TT]-s[tt].",
 658         " * [TT]res[tt] puts the center of mass of residues in the box.",
 659         " * [TT]atom[tt] puts all the atoms in the box.",
 660         " * [TT]nojump[tt] checks if atoms jump across the box and then puts",
 661         "   them back. This has the effect that all molecules",
 662         "   will remain whole (provided they were whole in the initial",
 663         "   conformation). [BB]Note[bb] that this ensures a continuous trajectory but",
 664         "   molecules may diffuse out of the box. The starting configuration",
 665         "   for this procedure is taken from the structure file, if one is",
 666         "   supplied, otherwise it is the first frame.",
 667         " * [TT]cluster[tt] clusters all the atoms in the selected index",
 668         "   such that they are all closest to the center of mass of the cluster,",
 669         "   which is iteratively updated. [BB]Note[bb] that this will only give meaningful",
 670         "   results if you in fact have a cluster. Luckily that can be checked",
 671         "   afterwards using a trajectory viewer. Note also that if your molecules",
 672         "   are broken this will not work either.",
 673         " * [TT]whole[tt] only makes broken molecules whole.",
 674         "",
 675
 676         "Option [TT]-ur[tt] sets the unit cell representation for options",
 677         "[TT]mol[tt], [TT]res[tt] and [TT]atom[tt] of [TT]-pbc[tt].",
 678         "All three options give different results for triclinic boxes and",
 679         "identical results for rectangular boxes.",
 680         "[TT]rect[tt] is the ordinary brick shape.",
 681         "[TT]tric[tt] is the triclinic unit cell.",
 682         "[TT]compact[tt] puts all atoms at the closest distance from the center",
 683         "of the box. This can be useful for visualizing e.g. truncated octahedra",
 684         "or rhombic dodecahedra. The center for options [TT]tric[tt] and [TT]compact[tt]",
 685         "is [TT]tric[tt] (see below), unless the option [TT]-boxcenter[tt]",
 686         "is set differently.[PAR]",
 687
 688         "Option [TT]-center[tt] centers the system in the box. The user can",
 689         "select the group which is used to determine the geometrical center.",
 690         "Option [TT]-boxcenter[tt] sets the location of the center of the box",
 691         "for options [TT]-pbc[tt] and [TT]-center[tt]. The center options are:",
 692         "[TT]tric[tt]: half of the sum of the box vectors,",
 693         "[TT]rect[tt]: half of the box diagonal,",
 694         "[TT]zero[tt]: zero.",
 695         "Use option [TT]-pbc mol[tt] in addition to [TT]-center[tt] when you",
 696         "want all molecules in the box after the centering.[PAR]",
 697
 698         "Option [TT]-box[tt] sets the size of the new box. This option only works",
 699         "for leading dimensions and is thus generally only useful for rectangular boxes.",
 700         "If you want to modify only some of the dimensions, e.g. when reading from",
 701         "a trajectory, you can use -1 for those dimensions that should stay the same",
 702
 703         "It is not always possible to use combinations of [TT]-pbc[tt],",
 704         "[TT]-fit[tt], [TT]-ur[tt] and [TT]-center[tt] to do exactly what",
 705         "you want in one call to [THISMODULE]. Consider using multiple",
 706         "calls, and check out the GROMACS website for suggestions.[PAR]",
 707
 708         "With [TT]-dt[tt], it is possible to reduce the number of ",
 709         "frames in the output. This option relies on the accuracy of the times",
 710         "in your input trajectory, so if these are inaccurate use the",
 711         "[TT]-timestep[tt] option to modify the time (this can be done",
 712         "simultaneously). For making smooth movies, the program [gmx-filter]",
 713         "can reduce the number of frames while using low-pass frequency",
 714         "filtering, this reduces aliasing of high frequency motions.[PAR]",
 715
 716         "Using [TT]-trunc[tt] [THISMODULE] can truncate [REF].trr[ref] in place, i.e.",
 717         "without copying the file. This is useful when a run has crashed",
 718         "during disk I/O (i.e. full disk), or when two contiguous",
 719         "trajectories must be concatenated without having double frames.[PAR]",
 720
 721         "Option [TT]-dump[tt] can be used to extract a frame at or near",
 722         "one specific time from your trajectory, but only works reliably",
 723         "if the time interval between frames is uniform.[PAR]",
 724
 725         "Option [TT]-drop[tt] reads an [REF].xvg[ref] file with times and values.",
 726         "When options [TT]-dropunder[tt] and/or [TT]-dropover[tt] are set,",
 727         "frames with a value below and above the value of the respective options",
 728         "will not be written."
 729     };
 730
 731     int         pbc_enum;
 732     enum
 733     {
 734         epSel,
 735         epNone,
 736         epComMol,
 737         epComRes,
 738         epComAtom,
 739         epNojump,
 740         epCluster,
 741         epWhole,
 742         epNR
 743     };
 744     const char *pbc_opt[epNR + 1] =
 745     {
 746         nullptr, "none", "mol", "res", "atom", "nojump", "cluster", "whole",
 747         nullptr
 748     };
 749
 750     int         unitcell_enum;
 751     const char *unitcell_opt[euNR+1] =
 752     { nullptr, "rect", "tric", "compact", nullptr };
 753
 754     enum
 755     {
 756         ecSel, ecTric, ecRect, ecZero, ecNR
 757     };
 758     const char *center_opt[ecNR+1] =
 759     { nullptr, "tric", "rect", "zero", nullptr };
 760     int         ecenter;
 761
 762     int         fit_enum;
 763     enum
 764     {
 765         efSel, efNone, efFit, efFitXY, efReset, efResetXY, efPFit, efNR
 766     };
 767     const char *fit[efNR + 1] =
 768     {
 769         nullptr, "none", "rot+trans", "rotxy+transxy", "translation", "transxy",
 770         "progressive", nullptr
 771     };
 772
 773     static gmx_bool  bSeparate     = FALSE, bVels = TRUE, bForce = FALSE, bCONECT = FALSE;
 774     static gmx_bool  bCenter       = FALSE;
 775     static int       skip_nr       = 1, ndec = 3, nzero = 0;
 776     static real      tzero         = 0, delta_t = 0, timestep = 0, ttrunc = -1, tdump = -1, split_t = 0;
 777     static rvec      newbox        = {0, 0, 0}, shift = {0, 0, 0}, trans = {0, 0, 0};
 778     static char     *exec_command  = nullptr;
 779     static real      dropunder     = 0, dropover = 0;
 780     static gmx_bool  bRound        = FALSE;
 781
 782     t_pargs
 783         pa[] =
 784     {
 785         { "-skip", FALSE, etINT,
 786           { &skip_nr }, "Only write every nr-th frame" },
 787         { "-dt", FALSE, etTIME,
 788           { &delta_t },
 789           "Only write frame when t MOD dt = first time (%t)" },
 790         { "-round", FALSE, etBOOL,
 791           { &bRound }, "Round measurements to nearest picosecond"},
 792         { "-dump", FALSE, etTIME,
 793           { &tdump }, "Dump frame nearest specified time (%t)" },
 794         { "-t0", FALSE, etTIME,
 795           { &tzero },
 796           "Starting time (%t) (default: don't change)" },
 797         { "-timestep", FALSE, etTIME,
 798           { &timestep },
 799           "Change time step between input frames (%t)" },
 800         { "-pbc", FALSE, etENUM,
 801           { pbc_opt },
 802           "PBC treatment (see help text for full description)" },
 803         { "-ur", FALSE, etENUM,
 804           { unitcell_opt }, "Unit-cell representation" },
 805         { "-center", FALSE, etBOOL,
 806           { &bCenter }, "Center atoms in box" },
 807         { "-boxcenter", FALSE, etENUM,
 808           { center_opt }, "Center for -pbc and -center" },
 809         { "-box", FALSE, etRVEC,
 810           { newbox },
 811           "Size for new cubic box (default: read from input)" },
 812         { "-trans", FALSE, etRVEC,
 813           { trans },
 814           "All coordinates will be translated by trans. This "
 815           "can advantageously be combined with -pbc mol -ur "
 816           "compact." },
 817         { "-shift", FALSE, etRVEC,
 818           { shift },
 819           "All coordinates will be shifted by framenr*shift" },
 820         { "-fit", FALSE, etENUM,
 821           { fit },
 822           "Fit molecule to ref structure in the structure file" },
 823         { "-ndec", FALSE, etINT,
 824           { &ndec },
 825           "Precision for .xtc and .gro writing in number of "
 826           "decimal places" },
 827         { "-vel", FALSE, etBOOL,
 828           { &bVels }, "Read and write velocities if possible" },
 829         { "-force", FALSE, etBOOL,
 830           { &bForce }, "Read and write forces if possible" },
 831         { "-trunc", FALSE, etTIME,
 832           { &ttrunc },
 833           "Truncate input trajectory file after this time (%t)" },
 834         { "-exec", FALSE, etSTR,
 835           { &exec_command },
 836           "Execute command for every output frame with the "
 837           "frame number as argument" },
 838         { "-split", FALSE, etTIME,
 839           { &split_t },
 840           "Start writing new file when t MOD split = first "
 841           "time (%t)" },
 842         { "-sep", FALSE, etBOOL,
 843           { &bSeparate },
 844           "Write each frame to a separate .gro, .g96 or .pdb "
 845           "file" },
 846         { "-nzero", FALSE, etINT,
 847           { &nzero },
 848           "If the -sep flag is set, use these many digits "
 849           "for the file numbers and prepend zeros as needed" },
 850         { "-dropunder", FALSE, etREAL,
 851           { &dropunder }, "Drop all frames below this value" },
 852         { "-dropover", FALSE, etREAL,
 853           { &dropover }, "Drop all frames above this value" },
 854         { "-conect", FALSE, etBOOL,
 855           { &bCONECT },
 856           "Add conect records when writing [REF].pdb[ref] files. Useful "
 857           "for visualization of non-standard molecules, e.g. "
 858           "coarse grained ones" }
 859     };
 860 #define NPA asize(pa)
 861
 862     FILE             *out    = nullptr;
 863     t_trxstatus      *trxout = nullptr;
 864     t_trxstatus      *trxin;
 865     int               file_nr;
 866     t_trxframe        fr, frout;
 867     int               flags;
 868     rvec             *xmem  = nullptr, *vmem = nullptr, *fmem = nullptr;
 869     rvec             *xp    = nullptr, x_shift, hbox;
 870     real             *w_rls = nullptr;
 871     int               m, i, d, frame, outframe, natoms, nout, ncent, newstep = 0, model_nr;
 872 #define SKIP 10
 873     t_topology        top;
 874     gmx_mtop_t       *mtop  = nullptr;
 875     gmx_conect        gc    = nullptr;
 876     int               ePBC  = -1;
 877     t_atoms          *atoms = nullptr, useatoms;
 878     matrix            top_box;
 879     int              *index = nullptr, *cindex = nullptr;
 880     char             *grpnm = nullptr;
 881     int              *frindex, nrfri;
 882     char             *frname;
 883     int               ifit, my_clust = -1;
 884     int              *ind_fit;
 885     char             *gn_fit;
 886     t_cluster_ndx    *clust           = nullptr;
 887     t_trxstatus     **clust_status    = nullptr;
 888     int              *clust_status_id = nullptr;
 889     int               ntrxopen        = 0;
 890     int              *nfwritten       = nullptr;
 891     int               ndrop           = 0, ncol, drop0 = 0, drop1 = 0, dropuse = 0;
 892     double          **dropval;
 893     real              tshift = 0, dt = -1, prec;
 894     gmx_bool          bFit, bPFit, bReset;
 895     int               nfitdim;
 896     gmx_rmpbc_t       gpbc = nullptr;
 897     gmx_bool          bRmPBC, bPBCWhole, bPBCcomRes, bPBCcomMol, bPBCcomAtom, bPBC, bNoJump, bCluster;
 898     gmx_bool          bCopy, bDoIt, bIndex, bTDump, bSetTime, bTPS = FALSE, bDTset = FALSE;
 899     gmx_bool          bExec, bTimeStep = FALSE, bDumpFrame = FALSE, bSetPrec, bNeedPrec;
 900     gmx_bool          bHaveFirstFrame, bHaveNextFrame, bSetBox, bSetUR, bSplit = FALSE;
 901     gmx_bool          bSubTraj = FALSE, bDropUnder = FALSE, bDropOver = FALSE, bTrans = FALSE;
 902     gmx_bool          bWriteFrame, bSplitHere;
 903     const char       *top_file, *in_file, *out_file = nullptr;
 904     char              out_file2[256], *charpt;
 905     char             *outf_base = nullptr;
 906     const char       *outf_ext  = nullptr;
 907     char              top_title[256], title[256], timestr[32], stepstr[32], filemode[5];
 908     gmx_output_env_t *oenv;
 909
 910     t_filenm          fnm[] = {
 911         { efTRX, "-f",   nullptr,      ffREAD  },
 912         { efTRO, "-o",   nullptr,      ffWRITE },
 913         { efTPS, nullptr,   nullptr,      ffOPTRD },
 914         { efNDX, nullptr,   nullptr,      ffOPTRD },
 915         { efNDX, "-fr",  "frames",  ffOPTRD },
 916         { efNDX, "-sub", "cluster", ffOPTRD },
 917         { efXVG, "-drop", "drop",    ffOPTRD }
 918     };
 919 #define NFILE asize(fnm)
 920
 921     if (!parse_common_args(&argc, argv,
 922                            PCA_CAN_BEGIN | PCA_CAN_END | PCA_CAN_VIEW |
 923                            PCA_TIME_UNIT,
 924                            NFILE, fnm, NPA, pa, asize(desc), desc,
 925                            0, nullptr, &oenv))
 926     {
 927         return 0;
 928     }
 929
 930     top_file = ftp2fn(efTPS, NFILE, fnm);
 931
 932     /* Check command line */
 933     in_file = opt2fn("-f", NFILE, fnm);
 934
 935     if (ttrunc != -1)
 936     {
 937         do_trunc(in_file, ttrunc);
 938     }
 939     else
 940     {
 941         /* mark active cmdline options */
 942         bSetBox    = opt2parg_bSet("-box", NPA, pa);
 943         bSetTime   = opt2parg_bSet("-t0", NPA, pa);
 944         bSetPrec   = opt2parg_bSet("-ndec", NPA, pa);
 945         bSetUR     = opt2parg_bSet("-ur", NPA, pa);
 946         bExec      = opt2parg_bSet("-exec", NPA, pa);
 947         bTimeStep  = opt2parg_bSet("-timestep", NPA, pa);
 948         bTDump     = opt2parg_bSet("-dump", NPA, pa);
 949         bDropUnder = opt2parg_bSet("-dropunder", NPA, pa);
 950         bDropOver  = opt2parg_bSet("-dropover", NPA, pa);
 951         bTrans     = opt2parg_bSet("-trans", NPA, pa);
 952         bSplit     = (split_t != 0);
 953
 954         /* parse enum options */
 955         fit_enum      = nenum(fit);
 956         bFit          = (fit_enum == efFit || fit_enum == efFitXY);
 957         bReset        = (fit_enum == efReset || fit_enum == efResetXY);
 958         bPFit         = fit_enum == efPFit;
 959         pbc_enum      = nenum(pbc_opt);
 960         bPBCWhole     = pbc_enum == epWhole;
 961         bPBCcomRes    = pbc_enum == epComRes;
 962         bPBCcomMol    = pbc_enum == epComMol;
 963         bPBCcomAtom   = pbc_enum == epComAtom;
 964         bNoJump       = pbc_enum == epNojump;
 965         bCluster      = pbc_enum == epCluster;
 966         bPBC          = pbc_enum != epNone;
 967         unitcell_enum = nenum(unitcell_opt);
 968         ecenter       = nenum(center_opt) - ecTric;
 969
 970         /* set and check option dependencies */
 971         if (bPFit)
 972         {
 973             bFit = TRUE;        /* for pfit, fit *must* be set */
 974         }
 975         if (bFit)
 976         {
 977             bReset = TRUE;       /* for fit, reset *must* be set */
 978         }
 979         nfitdim = 0;
 980         if (bFit || bReset)
 981         {
 982             nfitdim = (fit_enum == efFitXY || fit_enum == efResetXY) ? 2 : 3;
 983         }
 984         bRmPBC = bFit || bPBCWhole || bPBCcomRes || bPBCcomMol;
 985
 986         if (bSetUR)
 987         {
 988             if (!(bPBCcomRes || bPBCcomMol ||  bPBCcomAtom))
 989             {
 990                 fprintf(stderr,
 991                         "WARNING: Option for unitcell representation (-ur %s)\n"
 992                         "         only has effect in combination with -pbc %s, %s or %s.\n"
 993                         "         Ingoring unitcell representation.\n\n",
 994                         unitcell_opt[0], pbc_opt[2], pbc_opt[3], pbc_opt[4]);
 995             }
 996         }
 997         if (bFit && bPBC)
 998         {
 999             gmx_fatal(FARGS, "PBC condition treatment does not work together with rotational fit.\n"
1000                       "Please do the PBC condition treatment first and then run trjconv in a second step\n"
1001                       "for the rotational fit.\n"
1002                       "First doing the rotational fit and then doing the PBC treatment gives incorrect\n"
1003                       "results!");
1004         }
1005
1006         /* ndec is in nr of decimal places, prec is a multiplication factor: */
1007         prec = 1;
1008         for (i = 0; i < ndec; i++)
1009         {
1010             prec *= 10;
1011         }
1012
1013         bIndex = ftp2bSet(efNDX, NFILE, fnm);
1014
1015
1016         /* Determine output type */
1017         out_file = opt2fn("-o", NFILE, fnm);
1018         int ftp  = fn2ftp(out_file);
1019         fprintf(stderr, "Will write %s: %s\n", ftp2ext(ftp), ftp2desc(ftp));
1020         bNeedPrec = (ftp == efXTC || ftp == efGRO);
1021         int ftpin = fn2ftp(in_file);
1022         if (bVels)
1023         {
1024             /* check if velocities are possible in input and output files */
1025             bVels = (ftp == efTRR || ftp == efGRO ||
1026                      ftp == efG96 || ftp == efTNG)
1027                 && (ftpin == efTRR || ftpin == efGRO ||
1028                     ftpin == efG96 || ftpin == efTNG || ftpin == efCPT);
1029         }
1030         if (bSeparate || bSplit)
1031         {
1032             outf_ext = std::strrchr(out_file, '.');
1033             if (outf_ext == nullptr)
1034             {
1035                 gmx_fatal(FARGS, "Output file name '%s' does not contain a '.'", out_file);
1036             }
1037             outf_base = gmx_strdup(out_file);
1038             outf_base[outf_ext - out_file] = '\0';
1039         }
1040
1041         bSubTraj = opt2bSet("-sub", NFILE, fnm);
1042         if (bSubTraj)
1043         {
1044             if ((ftp != efXTC) && (ftp != efTRR))
1045             {
1046                 /* It seems likely that other trajectory file types
1047                  * could work here. */
1048                 gmx_fatal(FARGS, "Can only use the sub option with output file types "
1049                           "xtc and trr");
1050             }
1051             clust = cluster_index(nullptr, opt2fn("-sub", NFILE, fnm));
1052
1053             /* Check for number of files disabled, as FOPEN_MAX is not the correct
1054              * number to check for. In my linux box it is only 16.
1055              */
1056             if (0 && (clust->clust->nr > FOPEN_MAX-4))
1057             {
1058                 gmx_fatal(FARGS, "Can not open enough (%d) files to write all the"
1059                           " trajectories.\ntry splitting the index file in %d parts.\n"
1060                           "FOPEN_MAX = %d",
1061                           clust->clust->nr, 1+clust->clust->nr/FOPEN_MAX, FOPEN_MAX);
1062             }
1063             gmx_warning("The -sub option could require as many open output files as there are\n"
1064                         "index groups in the file (%d). If you get I/O errors opening new files,\n"
1065                         "try reducing the number of index groups in the file, and perhaps\n"
1066                         "using trjconv -sub several times on different chunks of your index file.\n",
1067                         clust->clust->nr);
1068
1069             snew(clust_status, clust->clust->nr);
1070             snew(clust_status_id, clust->clust->nr);
1071             snew(nfwritten, clust->clust->nr);
1072             for (i = 0; (i < clust->clust->nr); i++)
1073             {
1074                 clust_status[i]    = nullptr;
1075                 clust_status_id[i] = -1;
1076             }
1077             bSeparate = bSplit = FALSE;
1078         }
1079         /* skipping */
1080         if (skip_nr <= 0)
1081         {
1082         }
1083
1084         mtop = read_mtop_for_tng(top_file, in_file, out_file);
1085
1086         /* Determine whether to read a topology */
1087         bTPS = (ftp2bSet(efTPS, NFILE, fnm) ||
1088                 bRmPBC || bReset || bPBCcomMol || bCluster ||
1089                 (ftp == efGRO) || (ftp == efPDB) || bCONECT);
1090
1091         /* Determine if when can read index groups */
1092         bIndex = (bIndex || bTPS);
1093
1094         if (bTPS)
1095         {
1096             read_tps_conf(top_file, &top, &ePBC, &xp, nullptr, top_box,
1097                           bReset || bPBCcomRes);
1098             std::strncpy(top_title, *top.name, 255);
1099             top_title[255] = '\0';
1100             atoms          = &top.atoms;
1101
1102             if (0 == top.mols.nr && (bCluster || bPBCcomMol))
1103             {
1104                 gmx_fatal(FARGS, "Option -pbc %s requires a .tpr file for the -s option", pbc_opt[pbc_enum]);
1105             }
1106
1107             /* top_title is only used for gro and pdb,
1108              * the header in such a file is top_title, followed by
1109              * t= ... and/or step= ...
1110              * to prevent double t= or step=, remove it from top_title.
1111              * From GROMACS-2018 we only write t/step when the frame actually
1112              * has a valid time/step, so we need to check for both separately.
1113              */
1114             if ((charpt = std::strstr(top_title, " t= ")))
1115             {
1116                 charpt[0] = '\0';
1117             }
1118             if ((charpt = std::strstr(top_title, " step= ")))
1119             {
1120                 charpt[0] = '\0';
1121             }
1122
1123             if (bCONECT)
1124             {
1125                 gc = gmx_conect_generate(&top);
1126             }
1127             if (bRmPBC)
1128             {
1129                 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1130             }
1131         }
1132
1133         /* get frame number index */
1134         frindex = nullptr;
1135         if (opt2bSet("-fr", NFILE, fnm))
1136         {
1137             printf("Select groups of frame number indices:\n");
1138             rd_index(opt2fn("-fr", NFILE, fnm), 1, &nrfri, (int **)&frindex, &frname);
1139             if (debug)
1140             {
1141                 for (i = 0; i < nrfri; i++)
1142                 {
1143                     fprintf(debug, "frindex[%4d]=%4d\n", i, frindex[i]);
1144                 }
1145             }
1146         }
1147
1148         /* get index groups etc. */
1149         if (bReset)
1150         {
1151             printf("Select group for %s fit\n",
1152                    bFit ? "least squares" : "translational");
1153             get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1154                       1, &ifit, &ind_fit, &gn_fit);
1155
1156             if (bFit)
1157             {
1158                 if (ifit < 2)
1159                 {
1160                     gmx_fatal(FARGS, "Need at least 2 atoms to fit!\n");
1161                 }
1162                 else if (ifit == 3)
1163                 {
1164                     fprintf(stderr, "WARNING: fitting with only 2 atoms is not unique\n");
1165                 }
1166             }
1167         }
1168         else if (bCluster)
1169         {
1170             printf("Select group for clustering\n");
1171             get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1172                       1, &ifit, &ind_fit, &gn_fit);
1173         }
1174
1175         if (bIndex)
1176         {
1177             if (bCenter)
1178             {
1179                 printf("Select group for centering\n");
1180                 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1181                           1, &ncent, &cindex, &grpnm);
1182             }
1183             printf("Select group for output\n");
1184             get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
1185                       1, &nout, &index, &grpnm);
1186         }
1187         else
1188         {
1189             /* no index file, so read natoms from TRX */
1190             if (!read_first_frame(oenv, &trxin, in_file, &fr, TRX_DONT_SKIP))
1191             {
1192                 gmx_fatal(FARGS, "Could not read a frame from %s", in_file);
1193             }
1194             natoms = fr.natoms;
1195             close_trx(trxin);
1196             sfree(fr.x);
1197             snew(index, natoms);
1198             for (i = 0; i < natoms; i++)
1199             {
1200                 index[i] = i;
1201             }
1202             nout = natoms;
1203             if (bCenter)
1204             {
1205                 ncent  = nout;
1206                 cindex = index;
1207             }
1208         }
1209
1210         if (bReset)
1211         {
1212             snew(w_rls, atoms->nr);
1213             for (i = 0; (i < ifit); i++)
1214             {
1215                 w_rls[ind_fit[i]] = atoms->atom[ind_fit[i]].m;
1216             }
1217
1218             /* Restore reference structure and set to origin,
1219                store original location (to put structure back) */
1220             if (bRmPBC)
1221             {
1222                 gmx_rmpbc(gpbc, top.atoms.nr, top_box, xp);
1223             }
1224             copy_rvec(xp[index[0]], x_shift);
1225             reset_x_ndim(nfitdim, ifit, ind_fit, atoms->nr, nullptr, xp, w_rls);
1226             rvec_dec(x_shift, xp[index[0]]);
1227         }
1228         else
1229         {
1230             clear_rvec(x_shift);
1231         }
1232
1233         if (bDropUnder || bDropOver)
1234         {
1235             /* Read the .xvg file with the drop values */
1236             fprintf(stderr, "\nReading drop file ...");
1237             ndrop = read_xvg(opt2fn("-drop", NFILE, fnm), &dropval, &ncol);
1238             fprintf(stderr, " %d time points\n", ndrop);
1239             if (ndrop == 0 || ncol < 2)
1240             {
1241                 gmx_fatal(FARGS, "Found no data points in %s",
1242                           opt2fn("-drop", NFILE, fnm));
1243             }
1244             drop0 = 0;
1245             drop1 = 0;
1246         }
1247
1248         /* Make atoms struct for output in GRO or PDB files */
1249         if ((ftp == efGRO) || ((ftp == efG96) && bTPS) || (ftp == efPDB))
1250         {
1251             /* get memory for stuff to go in .pdb file, and initialize
1252              * the pdbinfo structure part if the input has it.
1253              */
1254             init_t_atoms(&useatoms, atoms->nr, atoms->havePdbInfo);
1255             sfree(useatoms.resinfo);
1256             useatoms.resinfo = atoms->resinfo;
1257             for (i = 0; (i < nout); i++)
1258             {
1259                 useatoms.atomname[i] = atoms->atomname[index[i]];
1260                 useatoms.atom[i]     = atoms->atom[index[i]];
1261                 if (atoms->havePdbInfo)
1262                 {
1263                     useatoms.pdbinfo[i]  = atoms->pdbinfo[index[i]];
1264                 }
1265                 useatoms.nres        = std::max(useatoms.nres, useatoms.atom[i].resind+1);
1266             }
1267             useatoms.nr = nout;
1268         }
1269         /* select what to read */
1270         if (ftp == efTRR)
1271         {
1272             flags = TRX_READ_X;
1273         }
1274         else
1275         {
1276             flags = TRX_NEED_X;
1277         }
1278         if (bVels)
1279         {
1280             flags = flags | TRX_READ_V;
1281         }
1282         if (bForce)
1283         {
1284             flags = flags | TRX_READ_F;
1285         }
1286
1287         /* open trx file for reading */
1288         bHaveFirstFrame = read_first_frame(oenv, &trxin, in_file, &fr, flags);
1289         if (fr.bPrec)
1290         {
1291             fprintf(stderr, "\nPrecision of %s is %g (nm)\n", in_file, 1/fr.prec);
1292         }
1293         if (bNeedPrec)
1294         {
1295             if (bSetPrec || !fr.bPrec)
1296             {
1297                 fprintf(stderr, "\nSetting output precision to %g (nm)\n", 1/prec);
1298             }
1299             else
1300             {
1301                 fprintf(stderr, "Using output precision of %g (nm)\n", 1/prec);
1302             }
1303         }
1304
1305         if (bHaveFirstFrame)
1306         {
1307             if (bTDump)
1308             {
1309                 // Determine timestep (assuming constant spacing for now) if we
1310                 // need to dump frames based on time. This is required so we do not
1311                 // skip the first frame if that was the one that should have been dumped
1312                 double firstFrameTime = fr.time;
1313                 if (read_next_frame(oenv, trxin, &fr))
1314                 {
1315                     dt     = fr.time - firstFrameTime;
1316                     bDTset = TRUE;
1317                     if (dt <= 0)
1318                     {
1319                         fprintf(stderr, "Warning: Frame times are not incrementing - will dump first frame.\n");
1320                     }
1321                 }
1322                 // Now close and reopen so we are at first frame again
1323                 close_trx(trxin);
1324                 done_frame(&fr);
1325                 // Reopen at first frame (We already know it exists if we got here)
1326                 read_first_frame(oenv, &trxin, in_file, &fr, flags);
1327             }
1328
1329             set_trxframe_ePBC(&fr, ePBC);
1330             natoms = fr.natoms;
1331
1332             if (bSetTime)
1333             {
1334                 tshift = tzero-fr.time;
1335             }
1336             else
1337             {
1338                 tzero = fr.time;
1339             }
1340
1341             bCopy = FALSE;
1342             if (bIndex)
1343             {
1344                 /* check if index is meaningful */
1345                 for (i = 0; i < nout; i++)
1346                 {
1347                     if (index[i] >= natoms)
1348                     {
1349                         gmx_fatal(FARGS,
1350                                   "Index[%d] %d is larger than the number of atoms in the\n"
1351                                   "trajectory file (%d). There is a mismatch in the contents\n"
1352                                   "of your -f, -s and/or -n files.", i, index[i]+1, natoms);
1353                     }
1354                     bCopy = bCopy || (i != index[i]);
1355                 }
1356             }
1357
1358             /* open output for writing */
1359             std::strcpy(filemode, "w");
1360             switch (ftp)
1361             {
1362                 case efTNG:
1363                     trxout = trjtools_gmx_prepare_tng_writing(out_file,
1364                                                               filemode[0],
1365                                                               trxin,
1366                                                               nullptr,
1367                                                               nout,
1368                                                               mtop,
1369                                                               index,
1370                                                               grpnm);
1371                     break;
1372                 case efXTC:
1373                 case efTRR:
1374                     out = nullptr;
1375                     if (!bSplit && !bSubTraj)
1376                     {
1377                         trxout = open_trx(out_file, filemode);
1378                     }
1379                     break;
1380                 case efGRO:
1381                 case efG96:
1382                 case efPDB:
1383                     if (( !bSeparate && !bSplit ) && !bSubTraj)
1384                     {
1385                         out = gmx_ffopen(out_file, filemode);
1386                     }
1387                     break;
1388                 default:
1389                     gmx_incons("Illegal output file format");
1390             }
1391
1392             if (bCopy)
1393             {
1394                 snew(xmem, nout);
1395                 if (bVels)
1396                 {
1397                     snew(vmem, nout);
1398                 }
1399                 if (bForce)
1400                 {
1401                     snew(fmem, nout);
1402                 }
1403             }
1404
1405             /* Start the big loop over frames */
1406             file_nr  =  0;
1407             frame    =  0;
1408             outframe =  0;
1409             model_nr =  0;
1410
1411             /* Main loop over frames */
1412             do
1413             {
1414                 if (!fr.bStep)
1415                 {
1416                     /* set the step */
1417                     fr.step = newstep;
1418                     newstep++;
1419                 }
1420                 if (bSubTraj)
1421                 {
1422                     /*if (frame >= clust->clust->nra)
1423                        gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file\n");*/
1424                     if (frame > clust->maxframe)
1425                     {
1426                         my_clust = -1;
1427                     }
1428                     else
1429                     {
1430                         my_clust = clust->inv_clust[frame];
1431                     }
1432                     if ((my_clust < 0) || (my_clust >= clust->clust->nr) ||
1433                         (my_clust == -1))
1434                     {
1435                         my_clust = -1;
1436                     }
1437                 }
1438
1439                 if (bSetBox)
1440                 {
1441                     /* generate new box */
1442                     if (fr.bBox == FALSE)
1443                     {
1444                         clear_mat(fr.box);
1445                     }
1446                     for (m = 0; m < DIM; m++)
1447                     {
1448                         if (newbox[m] >= 0)
1449                         {
1450                             fr.box[m][m] = newbox[m];
1451                         }
1452                         else
1453                         {
1454                             if (fr.bBox == FALSE)
1455                             {
1456                                 gmx_fatal(FARGS, "Cannot preserve a box that does not exist.\n");
1457                             }
1458                         }
1459                     }
1460                 }
1461
1462                 if (bTrans)
1463                 {
1464                     for (i = 0; i < natoms; i++)
1465                     {
1466                         rvec_inc(fr.x[i], trans);
1467                     }
1468                 }
1469
1470                 if (bTDump)
1471                 {
1472                     // If we could not read two frames or times are not incrementing
1473                     // we have almost no idea what to do,
1474                     // but dump the first frame so output is not broken.
1475                     if (dt <= 0 || !bDTset)
1476                     {
1477                         bDumpFrame = true;
1478                     }
1479                     else
1480                     {
1481                         // Dump the frame if we are less than half a frame time
1482                         // below it. This will also ensure we at least dump a
1483                         // somewhat reasonable frame if the spacing is unequal
1484                         // and we have overrun the frame time. Once we dump one
1485                         // frame based on time we quit, so it does not matter
1486                         // that this might be true for all subsequent frames too.
1487                         bDumpFrame = (fr.time > tdump-0.5*dt);
1488                     }
1489                 }
1490                 else
1491                 {
1492                     bDumpFrame = FALSE;
1493                 }
1494
1495                 /* determine if an atom jumped across the box and reset it if so */
1496                 if (bNoJump && (bTPS || frame != 0))
1497                 {
1498                     for (d = 0; d < DIM; d++)
1499                     {
1500                         hbox[d] = 0.5*fr.box[d][d];
1501                     }
1502                     for (i = 0; i < natoms; i++)
1503                     {
1504                         if (bReset)
1505                         {
1506                             rvec_dec(fr.x[i], x_shift);
1507                         }
1508                         for (m = DIM-1; m >= 0; m--)
1509                         {
1510                             if (hbox[m] > 0)
1511                             {
1512                                 while (fr.x[i][m]-xp[i][m] <= -hbox[m])
1513                                 {
1514                                     for (d = 0; d <= m; d++)
1515                                     {
1516                                         fr.x[i][d] += fr.box[m][d];
1517                                     }
1518                                 }
1519                                 while (fr.x[i][m]-xp[i][m] > hbox[m])
1520                                 {
1521                                     for (d = 0; d <= m; d++)
1522                                     {
1523                                         fr.x[i][d] -= fr.box[m][d];
1524                                     }
1525                                 }
1526                             }
1527                         }
1528                     }
1529                 }
1530                 else if (bCluster)
1531                 {
1532                     calc_pbc_cluster(ecenter, ifit, &top, ePBC, fr.x, ind_fit, fr.box);
1533                 }
1534
1535                 if (bPFit)
1536                 {
1537                     /* Now modify the coords according to the flags,
1538                        for normal fit, this is only done for output frames */
1539                     if (bRmPBC)
1540                     {
1541                         gmx_rmpbc_trxfr(gpbc, &fr);
1542                     }
1543
1544                     reset_x_ndim(nfitdim, ifit, ind_fit, natoms, nullptr, fr.x, w_rls);
1545                     do_fit(natoms, w_rls, xp, fr.x);
1546                 }
1547
1548                 /* store this set of coordinates for future use */
1549                 if (bPFit || bNoJump)
1550                 {
1551                     if (xp == nullptr)
1552                     {
1553                         snew(xp, natoms);
1554                     }
1555                     for (i = 0; (i < natoms); i++)
1556                     {
1557                         copy_rvec(fr.x[i], xp[i]);
1558                         rvec_inc(fr.x[i], x_shift);
1559                     }
1560                 }
1561
1562                 if (frindex)
1563                 {
1564                     /* see if we have a frame from the frame index group */
1565                     for (i = 0; i < nrfri && !bDumpFrame; i++)
1566                     {
1567                         bDumpFrame = frame == frindex[i];
1568                     }
1569                 }
1570                 if (debug && bDumpFrame)
1571                 {
1572                     fprintf(debug, "dumping %d\n", frame);
1573                 }
1574
1575                 bWriteFrame =
1576                     ( ( !bTDump && !frindex && frame % skip_nr == 0 ) || bDumpFrame );
1577
1578                 if (bWriteFrame && (bDropUnder || bDropOver))
1579                 {
1580                     while (dropval[0][drop1] < fr.time && drop1+1 < ndrop)
1581                     {
1582                         drop0 = drop1;
1583                         drop1++;
1584                     }
1585                     if (std::abs(dropval[0][drop0] - fr.time)
1586                         < std::abs(dropval[0][drop1] - fr.time))
1587                     {
1588                         dropuse = drop0;
1589                     }
1590                     else
1591                     {
1592                         dropuse = drop1;
1593                     }
1594                     if ((bDropUnder && dropval[1][dropuse] < dropunder) ||
1595                         (bDropOver && dropval[1][dropuse] > dropover))
1596                     {
1597                         bWriteFrame = FALSE;
1598                     }
1599                 }
1600
1601                 if (bWriteFrame)
1602                 {
1603                     /* We should avoid modifying the input frame,
1604                      * but since here we don't have the output frame yet,
1605                      * we introduce a temporary output frame time variable.
1606                      */
1607                     real frout_time;
1608
1609                     frout_time = fr.time;
1610
1611                     /* calc new time */
1612                     if (bTimeStep)
1613                     {
1614                         frout_time = tzero + frame*timestep;
1615                     }
1616                     else
1617                     if (bSetTime)
1618                     {
1619                         frout_time += tshift;
1620                     }
1621
1622                     if (bTDump)
1623                     {
1624                         fprintf(stderr, "\nDumping frame at t= %g %s\n",
1625                                 output_env_conv_time(oenv, frout_time), output_env_get_time_unit(oenv).c_str());
1626                     }
1627
1628                     /* check for writing at each delta_t */
1629                     bDoIt = (delta_t == 0);
1630                     if (!bDoIt)
1631                     {
1632                         if (!bRound)
1633                         {
1634                             bDoIt = bRmod(frout_time, tzero, delta_t);
1635                         }
1636                         else
1637                         {
1638                             /* round() is not C89 compatible, so we do this:  */
1639                             bDoIt = bRmod(std::floor(frout_time+0.5), std::floor(tzero+0.5),
1640                                           std::floor(delta_t+0.5));
1641                         }
1642                     }
1643
1644                     if (bDoIt || bTDump)
1645                     {
1646                         /* print sometimes */
1647                         if ( ((outframe % SKIP) == 0) || (outframe < SKIP) )
1648                         {
1649                             fprintf(stderr, " ->  frame %6d time %8.3f      \r",
1650                                     outframe, output_env_conv_time(oenv, frout_time));
1651                             fflush(stderr);
1652                         }
1653
1654                         if (!bPFit)
1655                         {
1656                             /* Now modify the coords according to the flags,
1657                                for PFit we did this already! */
1658
1659                             if (bRmPBC)
1660                             {
1661                                 gmx_rmpbc_trxfr(gpbc, &fr);
1662                             }
1663
1664                             if (bReset)
1665                             {
1666                                 reset_x_ndim(nfitdim, ifit, ind_fit, natoms, nullptr, fr.x, w_rls);
1667                                 if (bFit)
1668                                 {
1669                                     do_fit_ndim(nfitdim, natoms, w_rls, xp, fr.x);
1670                                 }
1671                                 if (!bCenter)
1672                                 {
1673                                     for (i = 0; i < natoms; i++)
1674                                     {
1675                                         rvec_inc(fr.x[i], x_shift);
1676                                     }
1677                                 }
1678                             }
1679
1680                             if (bCenter)
1681                             {
1682                                 center_x(ecenter, fr.x, fr.box, natoms, ncent, cindex);
1683                             }
1684                         }
1685
1686                         auto positionsArrayRef = gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec *>(fr.x), natoms);
1687                         if (bPBCcomAtom)
1688                         {
1689                             switch (unitcell_enum)
1690                             {
1691                                 case euRect:
1692                                     put_atoms_in_box(ePBC, fr.box, positionsArrayRef);
1693                                     break;
1694                                 case euTric:
1695                                     put_atoms_in_triclinic_unitcell(ecenter, fr.box, positionsArrayRef);
1696                                     break;
1697                                 case euCompact:
1698                                     put_atoms_in_compact_unitcell(ePBC, ecenter, fr.box,
1699                                                                   positionsArrayRef);
1700                                     break;
1701                             }
1702                         }
1703                         if (bPBCcomRes)
1704                         {
1705                             put_residue_com_in_box(unitcell_enum, ecenter,
1706                                                    natoms, atoms->atom, ePBC, fr.box, fr.x);
1707                         }
1708                         if (bPBCcomMol)
1709                         {
1710                             put_molecule_com_in_box(unitcell_enum, ecenter,
1711                                                     &top.mols,
1712                                                     natoms, atoms->atom, ePBC, fr.box, fr.x);
1713                         }
1714                         /* Copy the input trxframe struct to the output trxframe struct */
1715                         frout        = fr;
1716                         frout.time   = frout_time;
1717                         frout.bV     = (frout.bV && bVels);
1718                         frout.bF     = (frout.bF && bForce);
1719                         frout.natoms = nout;
1720                         if (bNeedPrec && (bSetPrec || !fr.bPrec))
1721                         {
1722                             frout.bPrec = TRUE;
1723                             frout.prec  = prec;
1724                         }
1725                         if (bCopy)
1726                         {
1727                             frout.x = xmem;
1728                             if (frout.bV)
1729                             {
1730                                 frout.v = vmem;
1731                             }
1732                             if (frout.bF)
1733                             {
1734                                 frout.f = fmem;
1735                             }
1736                             for (i = 0; i < nout; i++)
1737                             {
1738                                 copy_rvec(fr.x[index[i]], frout.x[i]);
1739                                 if (frout.bV)
1740                                 {
1741                                     copy_rvec(fr.v[index[i]], frout.v[i]);
1742                                 }
1743                                 if (frout.bF)
1744                                 {
1745                                     copy_rvec(fr.f[index[i]], frout.f[i]);
1746                                 }
1747                             }
1748                         }
1749
1750                         if (opt2parg_bSet("-shift", NPA, pa))
1751                         {
1752                             for (i = 0; i < nout; i++)
1753                             {
1754                                 for (d = 0; d < DIM; d++)
1755                                 {
1756                                     frout.x[i][d] += outframe*shift[d];
1757                                 }
1758                             }
1759                         }
1760
1761                         if (!bRound)
1762                         {
1763                             bSplitHere = bSplit && bRmod(frout.time, tzero, split_t);
1764                         }
1765                         else
1766                         {
1767                             /* round() is not C89 compatible, so we do this: */
1768                             bSplitHere = bSplit && bRmod(std::floor(frout.time+0.5),
1769                                                          std::floor(tzero+0.5),
1770                                                          std::floor(split_t+0.5));
1771                         }
1772                         if (bSeparate || bSplitHere)
1773                         {
1774                             mk_filenm(outf_base, ftp2ext(ftp), nzero, file_nr, out_file2);
1775                         }
1776
1777                         switch (ftp)
1778                         {
1779                             case efTNG:
1780                                 write_tng_frame(trxout, &frout);
1781                                 // TODO when trjconv behaves better: work how to read and write lambda
1782                                 break;
1783                             case efTRR:
1784                             case efXTC:
1785                                 if (bSplitHere)
1786                                 {
1787                                     if (trxout)
1788                                     {
1789                                         close_trx(trxout);
1790                                     }
1791                                     trxout = open_trx(out_file2, filemode);
1792                                 }
1793                                 if (bSubTraj)
1794                                 {
1795                                     if (my_clust != -1)
1796                                     {
1797                                         char buf[STRLEN];
1798                                         if (clust_status_id[my_clust] == -1)
1799                                         {
1800                                             sprintf(buf, "%s.%s", clust->grpname[my_clust], ftp2ext(ftp));
1801                                             clust_status[my_clust]    = open_trx(buf, "w");
1802                                             clust_status_id[my_clust] = 1;
1803                                             ntrxopen++;
1804                                         }
1805                                         else if (clust_status_id[my_clust] == -2)
1806                                         {
1807                                             gmx_fatal(FARGS, "File %s.xtc should still be open (%d open .xtc files)\n" "in order to write frame %d. my_clust = %d",
1808                                                       clust->grpname[my_clust], ntrxopen, frame,
1809                                                       my_clust);
1810                                         }
1811                                         write_trxframe(clust_status[my_clust], &frout, gc);
1812                                         nfwritten[my_clust]++;
1813                                         if (nfwritten[my_clust] ==
1814                                             (clust->clust->index[my_clust+1]-
1815                                              clust->clust->index[my_clust]))
1816                                         {
1817                                             close_trx(clust_status[my_clust]);
1818                                             clust_status[my_clust]    = nullptr;
1819                                             clust_status_id[my_clust] = -2;
1820                                             ntrxopen--;
1821                                             if (ntrxopen < 0)
1822                                             {
1823                                                 gmx_fatal(FARGS, "Less than zero open .xtc files!");
1824                                             }
1825                                         }
1826                                     }
1827                                 }
1828                                 else
1829                                 {
1830                                     write_trxframe(trxout, &frout, gc);
1831                                 }
1832                                 break;
1833                             case efGRO:
1834                             case efG96:
1835                             case efPDB:
1836                                 // Only add a generator statement if title is empty,
1837                                 // to avoid multiple generated-by statements from various programs
1838                                 if (std::strlen(top_title) == 0)
1839                                 {
1840                                     sprintf(top_title, "Generated by trjconv");
1841                                 }
1842                                 if (frout.bTime)
1843                                 {
1844                                     sprintf(timestr, " t= %9.5f", frout.time);
1845                                 }
1846                                 else
1847                                 {
1848                                     std::strcpy(timestr, "");
1849                                 }
1850                                 if (frout.bStep)
1851                                 {
1852                                     sprintf(stepstr, " step= %" GMX_PRId64, frout.step);
1853                                 }
1854                                 else
1855                                 {
1856                                     std::strcpy(stepstr, "");
1857                                 }
1858                                 snprintf(title, 256, "%s%s%s", top_title, timestr, stepstr);
1859                                 if (bSeparate || bSplitHere)
1860                                 {
1861                                     out = gmx_ffopen(out_file2, "w");
1862                                 }
1863                                 switch (ftp)
1864                                 {
1865                                     case efGRO:
1866                                         write_hconf_p(out, title, &useatoms,
1867                                                       frout.x, frout.bV ? frout.v : nullptr, frout.box);
1868                                         break;
1869                                     case efPDB:
1870                                         fprintf(out, "REMARK    GENERATED BY TRJCONV\n");
1871                                         /* if reading from pdb, we want to keep the original
1872                                            model numbering else we write the output frame
1873                                            number plus one, because model 0 is not allowed in pdb */
1874                                         if (ftpin == efPDB && fr.bStep && fr.step > model_nr)
1875                                         {
1876                                             model_nr = fr.step;
1877                                         }
1878                                         else
1879                                         {
1880                                             model_nr++;
1881                                         }
1882                                         write_pdbfile(out, title, &useatoms, frout.x,
1883                                                       frout.ePBC, frout.box, ' ', model_nr, gc, TRUE);
1884                                         break;
1885                                     case efG96:
1886                                         const char *outputTitle = "";
1887                                         if (bSeparate || bTDump)
1888                                         {
1889                                             outputTitle = title;
1890                                             if (bTPS)
1891                                             {
1892                                                 frout.bAtoms = TRUE;
1893                                             }
1894                                             frout.atoms  = &useatoms;
1895                                             frout.bStep  = FALSE;
1896                                             frout.bTime  = FALSE;
1897                                         }
1898                                         else
1899                                         {
1900                                             if (outframe == 0)
1901                                             {
1902                                                 outputTitle = title;
1903                                             }
1904                                             frout.bAtoms = FALSE;
1905                                             frout.bStep  = TRUE;
1906                                             frout.bTime  = TRUE;
1907                                         }
1908                                         write_g96_conf(out, outputTitle, &frout, -1, nullptr);
1909                                 }
1910                                 if (bSeparate || bSplitHere)
1911                                 {
1912                                     gmx_ffclose(out);
1913                                     out = nullptr;
1914                                 }
1915                                 break;
1916                             default:
1917                                 gmx_fatal(FARGS, "DHE, ftp=%d\n", ftp);
1918                         }
1919                         if (bSeparate || bSplitHere)
1920                         {
1921                             file_nr++;
1922                         }
1923
1924                         /* execute command */
1925                         if (bExec)
1926                         {
1927                             char c[255];
1928                             sprintf(c, "%s  %d", exec_command, file_nr-1);
1929                             /*fprintf(stderr,"Executing '%s'\n",c);*/
1930                             if (0 != system(c))
1931                             {
1932                                 gmx_fatal(FARGS, "Error executing command: %s", c);
1933                             }
1934                         }
1935                         outframe++;
1936                     }
1937                 }
1938                 frame++;
1939                 bHaveNextFrame = read_next_frame(oenv, trxin, &fr);
1940             }
1941             while (!(bTDump && bDumpFrame) && bHaveNextFrame);
1942         }
1943
1944         if (!bHaveFirstFrame || (bTDump && !bDumpFrame))
1945         {
1946             fprintf(stderr, "\nWARNING no output, "
1947                     "last frame read at t=%g\n", fr.time);
1948         }
1949         fprintf(stderr, "\n");
1950
1951         close_trx(trxin);
1952         sfree(outf_base);
1953
1954         if (bRmPBC)
1955         {
1956             gmx_rmpbc_done(gpbc);
1957         }
1958
1959         if (trxout)
1960         {
1961             close_trx(trxout);
1962         }
1963         else if (out != nullptr)
1964         {
1965             gmx_ffclose(out);
1966         }
1967         if (bSubTraj)
1968         {
1969             for (i = 0; (i < clust->clust->nr); i++)
1970             {
1971                 if (clust_status_id[i] >= 0)
1972                 {
1973                     close_trx(clust_status[i]);
1974                 }
1975             }
1976         }
1977     }
1978
1979     if (mtop)
1980     {
1981         sfree(mtop);
1982     }
1983     if (bTPS)
1984     {
1985         done_top(&top);
1986     }
1987     sfree(xp);
1988     sfree(xmem);
1989     sfree(vmem);
1990     sfree(fmem);
1991     sfree(grpnm);
1992     sfree(index);
1993     sfree(cindex);
1994     done_frame(&fr);
1995
1996     do_view(oenv, out_file, nullptr);
1997
1998     output_env_done(oenv);
1999     return 0;
2000 }