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