--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "trjconv.h"
+
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+
+#include <algorithm>
+#include <memory>
+
+#include "gromacs/commandline/pargs.h"
+#include "gromacs/commandline/viewit.h"
+#include "gromacs/fileio/confio.h"
+#include "gromacs/fileio/g96io.h"
+#include "gromacs/fileio/gmxfio.h"
+#include "gromacs/fileio/groio.h"
+#include "gromacs/fileio/pdbio.h"
+#include "gromacs/fileio/tngio.h"
+#include "gromacs/fileio/tpxio.h"
+#include "gromacs/fileio/trrio.h"
+#include "gromacs/fileio/trxio.h"
+#include "gromacs/fileio/xtcio.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/math/do_fit.h"
+#include "gromacs/math/functions.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pbcutil/rmpbc.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/trajectory/trajectoryframe.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/arraysize.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
+
+enum {
+ euSel, euRect, euTric, euCompact, euNR
+};
+
+
+static void calc_pbc_cluster(int ecenter, int nrefat, t_topology *top, int ePBC,
+ rvec x[], const int index[], matrix box)
+{
+ int m, i, j, j0, j1, jj, ai, aj;
+ int imin, jmin;
+ real fac, min_dist2;
+ rvec dx, xtest, box_center;
+ int nmol, imol_center;
+ int *molind;
+ gmx_bool *bMol, *bTmp;
+ rvec *m_com, *m_shift;
+ t_pbc pbc;
+ int *cluster;
+ int *added;
+ int ncluster, nadded;
+ real tmp_r2;
+
+ calc_box_center(ecenter, box, box_center);
+
+ /* Initiate the pbc structure */
+ std::memset(&pbc, 0, sizeof(pbc));
+ set_pbc(&pbc, ePBC, box);
+
+ /* Convert atom index to molecular */
+ nmol = top->mols.nr;
+ molind = top->mols.index;
+ snew(bMol, nmol);
+ snew(m_com, nmol);
+ snew(m_shift, nmol);
+ snew(cluster, nmol);
+ snew(added, nmol);
+ snew(bTmp, top->atoms.nr);
+
+ for (i = 0; (i < nrefat); i++)
+ {
+ /* Mark all molecules in the index */
+ ai = index[i];
+ bTmp[ai] = TRUE;
+ /* Binary search assuming the molecules are sorted */
+ j0 = 0;
+ j1 = nmol-1;
+ while (j0 < j1)
+ {
+ if (ai < molind[j0+1])
+ {
+ j1 = j0;
+ }
+ else if (ai >= molind[j1])
+ {
+ j0 = j1;
+ }
+ else
+ {
+ jj = (j0+j1)/2;
+ if (ai < molind[jj+1])
+ {
+ j1 = jj;
+ }
+ else
+ {
+ j0 = jj;
+ }
+ }
+ }
+ bMol[j0] = TRUE;
+ }
+ /* Double check whether all atoms in all molecules that are marked are part
+ * of the cluster. Simultaneously compute the center of geometry.
+ */
+ min_dist2 = 10*gmx::square(trace(box));
+ imol_center = -1;
+ ncluster = 0;
+ for (i = 0; i < nmol; i++)
+ {
+ for (j = molind[i]; j < molind[i+1]; j++)
+ {
+ if (bMol[i] && !bTmp[j])
+ {
+ gmx_fatal(FARGS, "Molecule %d marked for clustering but not atom %d in it - check your index!", i+1, j+1);
+ }
+ else if (!bMol[i] && bTmp[j])
+ {
+ gmx_fatal(FARGS, "Atom %d marked for clustering but not molecule %d - this is an internal error...", j+1, i+1);
+ }
+ else if (bMol[i])
+ {
+ /* Make molecule whole, move 2nd and higher atom to same periodicity as 1st atom in molecule */
+ if (j > molind[i])
+ {
+ pbc_dx(&pbc, x[j], x[j-1], dx);
+ rvec_add(x[j-1], dx, x[j]);
+ }
+ /* Compute center of geometry of molecule - m_com[i] was zeroed when we did snew() on it! */
+ rvec_inc(m_com[i], x[j]);
+ }
+ }
+ if (bMol[i])
+ {
+ /* Normalize center of geometry */
+ fac = 1.0/(molind[i+1]-molind[i]);
+ for (m = 0; (m < DIM); m++)
+ {
+ m_com[i][m] *= fac;
+ }
+ /* Determine which molecule is closest to the center of the box */
+ pbc_dx(&pbc, box_center, m_com[i], dx);
+ tmp_r2 = iprod(dx, dx);
+
+ if (tmp_r2 < min_dist2)
+ {
+ min_dist2 = tmp_r2;
+ imol_center = i;
+ }
+ cluster[ncluster++] = i;
+ }
+ }
+ sfree(bTmp);
+
+ if (ncluster <= 0)
+ {
+ fprintf(stderr, "No molecules selected in the cluster\n");
+ return;
+ }
+ else if (imol_center == -1)
+ {
+ fprintf(stderr, "No central molecules could be found\n");
+ return;
+ }
+
+ nadded = 0;
+ added[nadded++] = imol_center;
+ bMol[imol_center] = FALSE;
+
+ while (nadded < ncluster)
+ {
+ /* Find min distance between cluster molecules and those remaining to be added */
+ min_dist2 = 10*gmx::square(trace(box));
+ imin = -1;
+ jmin = -1;
+ /* Loop over added mols */
+ for (i = 0; i < nadded; i++)
+ {
+ ai = added[i];
+ /* Loop over all mols */
+ for (j = 0; j < ncluster; j++)
+ {
+ aj = cluster[j];
+ /* check those remaining to be added */
+ if (bMol[aj])
+ {
+ pbc_dx(&pbc, m_com[aj], m_com[ai], dx);
+ tmp_r2 = iprod(dx, dx);
+ if (tmp_r2 < min_dist2)
+ {
+ min_dist2 = tmp_r2;
+ imin = ai;
+ jmin = aj;
+ }
+ }
+ }
+ }
+
+ /* Add the best molecule */
+ added[nadded++] = jmin;
+ bMol[jmin] = FALSE;
+ /* Calculate the shift from the ai molecule */
+ pbc_dx(&pbc, m_com[jmin], m_com[imin], dx);
+ rvec_add(m_com[imin], dx, xtest);
+ rvec_sub(xtest, m_com[jmin], m_shift[jmin]);
+ rvec_inc(m_com[jmin], m_shift[jmin]);
+
+ for (j = molind[jmin]; j < molind[jmin+1]; j++)
+ {
+ rvec_inc(x[j], m_shift[jmin]);
+ }
+ fprintf(stdout, "\rClustering iteration %d of %d...", nadded, ncluster);
+ fflush(stdout);
+ }
+
+ sfree(added);
+ sfree(cluster);
+ sfree(bMol);
+ sfree(m_com);
+ sfree(m_shift);
+
+ fprintf(stdout, "\n");
+}
+
+static void put_molecule_com_in_box(int unitcell_enum, int ecenter,
+ t_block *mols,
+ int natoms, t_atom atom[],
+ int ePBC, matrix box, rvec x[])
+{
+ int i, j;
+ int d;
+ rvec com, shift, box_center;
+ real m;
+ double mtot;
+ t_pbc pbc;
+
+ calc_box_center(ecenter, box, box_center);
+ set_pbc(&pbc, ePBC, box);
+ if (mols->nr <= 0)
+ {
+ gmx_fatal(FARGS, "There are no molecule descriptions. I need a .tpr file for this pbc option.");
+ }
+ for (i = 0; (i < mols->nr); i++)
+ {
+ /* calc COM */
+ clear_rvec(com);
+ mtot = 0;
+ for (j = mols->index[i]; (j < mols->index[i+1] && j < natoms); j++)
+ {
+ m = atom[j].m;
+ for (d = 0; d < DIM; d++)
+ {
+ com[d] += m*x[j][d];
+ }
+ mtot += m;
+ }
+ /* calculate final COM */
+ svmul(1.0/mtot, com, com);
+
+ /* check if COM is outside box */
+ gmx::RVec newCom;
+ copy_rvec(com, newCom);
+ auto newComArrayRef = gmx::arrayRefFromArray(&newCom, 1);
+ switch (unitcell_enum)
+ {
+ case euRect:
+ put_atoms_in_box(ePBC, box, newComArrayRef);
+ break;
+ case euTric:
+ put_atoms_in_triclinic_unitcell(ecenter, box, newComArrayRef);
+ break;
+ case euCompact:
+ put_atoms_in_compact_unitcell(ePBC, ecenter, box, newComArrayRef);
+ break;
+ }
+ rvec_sub(newCom, com, shift);
+ if (norm2(shift) > 0)
+ {
+ if (debug)
+ {
+ fprintf(debug, "\nShifting position of molecule %d "
+ "by %8.3f %8.3f %8.3f\n", i+1,
+ shift[XX], shift[YY], shift[ZZ]);
+ }
+ for (j = mols->index[i]; (j < mols->index[i+1] && j < natoms); j++)
+ {
+ rvec_inc(x[j], shift);
+ }
+ }
+ }
+}
+
+static void put_residue_com_in_box(int unitcell_enum, int ecenter,
+ int natoms, t_atom atom[],
+ int ePBC, matrix box, rvec x[])
+{
+ int i, j, res_start, res_end;
+ int d, presnr;
+ real m;
+ double mtot;
+ rvec box_center, com, shift;
+ static const int NOTSET = -12347;
+ calc_box_center(ecenter, box, box_center);
+
+ presnr = NOTSET;
+ res_start = 0;
+ clear_rvec(com);
+ mtot = 0;
+ for (i = 0; i < natoms+1; i++)
+ {
+ if (i == natoms || (presnr != atom[i].resind && presnr != NOTSET))
+ {
+ /* calculate final COM */
+ res_end = i;
+ svmul(1.0/mtot, com, com);
+
+ /* check if COM is outside box */
+ gmx::RVec newCom;
+ copy_rvec(com, newCom);
+ auto newComArrayRef = gmx::arrayRefFromArray(&newCom, 1);
+ switch (unitcell_enum)
+ {
+ case euRect:
+ put_atoms_in_box(ePBC, box, newComArrayRef);
+ break;
+ case euTric:
+ put_atoms_in_triclinic_unitcell(ecenter, box, newComArrayRef);
+ break;
+ case euCompact:
+ put_atoms_in_compact_unitcell(ePBC, ecenter, box, newComArrayRef);
+ break;
+ }
+ rvec_sub(newCom, com, shift);
+ if (norm2(shift) != 0.0f)
+ {
+ if (debug)
+ {
+ fprintf(debug, "\nShifting position of residue %d (atoms %d-%d) "
+ "by %g,%g,%g\n", atom[res_start].resind+1,
+ res_start+1, res_end+1, shift[XX], shift[YY], shift[ZZ]);
+ }
+ for (j = res_start; j < res_end; j++)
+ {
+ rvec_inc(x[j], shift);
+ }
+ }
+ clear_rvec(com);
+ mtot = 0;
+
+ /* remember start of new residue */
+ res_start = i;
+ }
+ if (i < natoms)
+ {
+ /* calc COM */
+ m = atom[i].m;
+ for (d = 0; d < DIM; d++)
+ {
+ com[d] += m*x[i][d];
+ }
+ mtot += m;
+
+ presnr = atom[i].resind;
+ }
+ }
+}
+
+static void center_x(int ecenter, rvec x[], matrix box, int n, int nc, const int ci[])
+{
+ int i, m, ai;
+ rvec cmin, cmax, box_center, dx;
+
+ if (nc > 0)
+ {
+ copy_rvec(x[ci[0]], cmin);
+ copy_rvec(x[ci[0]], cmax);
+ for (i = 0; i < nc; i++)
+ {
+ ai = ci[i];
+ for (m = 0; m < DIM; m++)
+ {
+ if (x[ai][m] < cmin[m])
+ {
+ cmin[m] = x[ai][m];
+ }
+ else if (x[ai][m] > cmax[m])
+ {
+ cmax[m] = x[ai][m];
+ }
+ }
+ }
+ calc_box_center(ecenter, box, box_center);
+ for (m = 0; m < DIM; m++)
+ {
+ dx[m] = box_center[m]-(cmin[m]+cmax[m])*0.5;
+ }
+
+ for (i = 0; i < n; i++)
+ {
+ rvec_inc(x[i], dx);
+ }
+ }
+}
+
+static void mk_filenm(char *base, const char *ext, int ndigit, int file_nr,
+ char out_file[])
+{
+ char nbuf[128];
+ int nd = 0, fnr;
+
+ std::strcpy(out_file, base);
+ fnr = file_nr;
+ do
+ {
+ fnr /= 10;
+ nd++;
+ }
+ while (fnr > 0);
+
+ if (nd < ndigit)
+ {
+ std::strncat(out_file, "00000000000", ndigit-nd);
+ }
+ sprintf(nbuf, "%d.", file_nr);
+ std::strcat(out_file, nbuf);
+ std::strcat(out_file, ext);
+}
+
+static void check_trr(const char *fn)
+{
+ if (fn2ftp(fn) != efTRR)
+ {
+ gmx_fatal(FARGS, "%s is not a trajectory file, exiting\n", fn);
+ }
+}
+
+static void do_trunc(const char *fn, real t0)
+{
+ t_fileio *in;
+ FILE *fp;
+ gmx_bool bStop, bOK;
+ gmx_trr_header_t sh;
+ gmx_off_t fpos;
+ char yesno[256];
+ int j;
+ real t = 0;
+
+ if (t0 == -1)
+ {
+ gmx_fatal(FARGS, "You forgot to set the truncation time");
+ }
+
+ /* Check whether this is a .trr file */
+ check_trr(fn);
+
+ in = gmx_trr_open(fn, "r");
+ fp = gmx_fio_getfp(in);
+ if (fp == nullptr)
+ {
+ fprintf(stderr, "Sorry, can not trunc %s, truncation of this filetype is not supported\n", fn);
+ gmx_trr_close(in);
+ }
+ else
+ {
+ j = 0;
+ fpos = gmx_fio_ftell(in);
+ bStop = FALSE;
+ while (!bStop && gmx_trr_read_frame_header(in, &sh, &bOK))
+ {
+ gmx_trr_read_frame_data(in, &sh, nullptr, nullptr, nullptr, nullptr);
+ fpos = gmx_ftell(fp);
+ t = sh.t;
+ if (t >= t0)
+ {
+ gmx_fseek(fp, fpos, SEEK_SET);
+ bStop = TRUE;
+ }
+ }
+ if (bStop)
+ {
+ fprintf(stderr, "Do you REALLY want to truncate this trajectory (%s) at:\n"
+ "frame %d, time %g, bytes %ld ??? (type YES if so)\n",
+ fn, j, t, static_cast<long int>(fpos));
+ if (1 != scanf("%s", yesno))
+ {
+ gmx_fatal(FARGS, "Error reading user input");
+ }
+ if (std::strcmp(yesno, "YES") == 0)
+ {
+ fprintf(stderr, "Once again, I'm gonna DO this...\n");
+ gmx_trr_close(in);
+ if (0 != gmx_truncate(fn, fpos))
+ {
+ gmx_fatal(FARGS, "Error truncating file %s", fn);
+ }
+ }
+ else
+ {
+ fprintf(stderr, "Ok, I'll forget about it\n");
+ }
+ }
+ else
+ {
+ fprintf(stderr, "Already at end of file (t=%g)...\n", t);
+ gmx_trr_close(in);
+ }
+ }
+}
+
+/*! \brief Read a full molecular topology if useful and available.
+ *
+ * If the input trajectory file is not in TNG format, and the output
+ * file is in TNG format, then we want to try to read a full topology
+ * (if available), so that we can write molecule information to the
+ * output file. The full topology provides better molecule information
+ * than is available from the normal t_topology data used by GROMACS
+ * tools.
+ *
+ * Also, the t_topology is only read under (different) particular
+ * conditions. If both apply, then a .tpr file might be read
+ * twice. Trying to fix this redundancy while trjconv is still an
+ * all-purpose tool does not seem worthwhile.
+ *
+ * Because of the way gmx_prepare_tng_writing is implemented, the case
+ * where the input TNG file has no molecule information will never
+ * lead to an output TNG file having molecule information. Since
+ * molecule information will generally be present if the input TNG
+ * file was written by a GROMACS tool, this seems like reasonable
+ * behaviour. */
+static std::unique_ptr<gmx_mtop_t>
+read_mtop_for_tng(const char *tps_file,
+ const char *input_file,
+ const char *output_file)
+{
+ std::unique_ptr<gmx_mtop_t> mtop;
+
+ if (fn2bTPX(tps_file) &&
+ efTNG != fn2ftp(input_file) &&
+ efTNG == fn2ftp(output_file))
+ {
+ int temp_natoms = -1;
+ mtop = std::make_unique<gmx_mtop_t>();
+ read_tpx(tps_file, nullptr, nullptr, &temp_natoms,
+ nullptr, nullptr, mtop.get());
+ }
+
+ return mtop;
+}
+
+int gmx_trjconv(int argc, char *argv[])
+{
+ const char *desc[] = {
+ "[THISMODULE] can convert trajectory files in many ways:",
+ "",
+ "* from one format to another",
+ "* select a subset of atoms",
+ "* change the periodicity representation",
+ "* keep multimeric molecules together",
+ "* center atoms in the box",
+ "* fit atoms to reference structure",
+ "* reduce the number of frames",
+ "* change the timestamps of the frames ([TT]-t0[tt] and [TT]-timestep[tt])",
+ "* cut the trajectory in small subtrajectories according",
+ " to information in an index file. This allows subsequent analysis of",
+ " the subtrajectories that could, for example, be the result of a",
+ " cluster analysis. Use option [TT]-sub[tt].",
+ " This assumes that the entries in the index file are frame numbers and",
+ " dumps each group in the index file to a separate trajectory file.",
+ "* select frames within a certain range of a quantity given",
+ " in an [REF].xvg[ref] file.",
+ "",
+ "[gmx-trjcat] is better suited for concatenating multiple trajectory files.",
+ "[PAR]",
+
+ "The following formats are supported for input and output:",
+ "[REF].xtc[ref], [REF].trr[ref], [REF].gro[ref], [TT].g96[tt]",
+ "and [REF].pdb[ref].",
+ "The file formats are detected from the file extension.",
+ "The precision of the [REF].xtc[ref] output is taken from the",
+ "input file for [REF].xtc[ref], [REF].gro[ref] and [REF].pdb[ref],",
+ "and from the [TT]-ndec[tt] option for other input formats. The precision",
+ "is always taken from [TT]-ndec[tt], when this option is set.",
+ "All other formats have fixed precision. [REF].trr[ref]",
+ "output can be single or double precision, depending on the precision",
+ "of the [THISMODULE] binary.",
+ "Note that velocities are only supported in",
+ "[REF].trr[ref], [REF].gro[ref] and [TT].g96[tt] files.[PAR]",
+
+ "Option [TT]-sep[tt] can be used to write every frame to a separate",
+ "[TT].gro, .g96[tt] or [REF].pdb[ref] file. By default, all frames all written to one file.",
+ "[REF].pdb[ref] files with all frames concatenated can be viewed with",
+ "[TT]rasmol -nmrpdb[tt].[PAR]",
+
+ "It is possible to select part of your trajectory and write it out",
+ "to a new trajectory file in order to save disk space, e.g. for leaving",
+ "out the water from a trajectory of a protein in water.",
+ "[BB]ALWAYS[bb] put the original trajectory on tape!",
+ "We recommend to use the portable [REF].xtc[ref] format for your analysis",
+ "to save disk space and to have portable files.[PAR]",
+
+ "There are two options for fitting the trajectory to a reference",
+ "either for essential dynamics analysis, etc.",
+ "The first option is just plain fitting to a reference structure",
+ "in the structure file. The second option is a progressive fit",
+ "in which the first timeframe is fitted to the reference structure ",
+ "in the structure file to obtain and each subsequent timeframe is ",
+ "fitted to the previously fitted structure. This way a continuous",
+ "trajectory is generated, which might not be the case when using the",
+ "regular fit method, e.g. when your protein undergoes large",
+ "conformational transitions.[PAR]",
+
+ "Option [TT]-pbc[tt] sets the type of periodic boundary condition",
+ "treatment:",
+ "",
+ " * [TT]mol[tt] puts the center of mass of molecules in the box,",
+ " and requires a run input file to be supplied with [TT]-s[tt].",
+ " * [TT]res[tt] puts the center of mass of residues in the box.",
+ " * [TT]atom[tt] puts all the atoms in the box.",
+ " * [TT]nojump[tt] checks if atoms jump across the box and then puts",
+ " them back. This has the effect that all molecules",
+ " will remain whole (provided they were whole in the initial",
+ " conformation). [BB]Note[bb] that this ensures a continuous trajectory but",
+ " molecules may diffuse out of the box. The starting configuration",
+ " for this procedure is taken from the structure file, if one is",
+ " supplied, otherwise it is the first frame.",
+ " * [TT]cluster[tt] clusters all the atoms in the selected index",
+ " such that they are all closest to the center of mass of the cluster,",
+ " which is iteratively updated. [BB]Note[bb] that this will only give meaningful",
+ " results if you in fact have a cluster. Luckily that can be checked",
+ " afterwards using a trajectory viewer. Note also that if your molecules",
+ " are broken this will not work either.",
+ " * [TT]whole[tt] only makes broken molecules whole.",
+ "",
+
+ "Option [TT]-ur[tt] sets the unit cell representation for options",
+ "[TT]mol[tt], [TT]res[tt] and [TT]atom[tt] of [TT]-pbc[tt].",
+ "All three options give different results for triclinic boxes and",
+ "identical results for rectangular boxes.",
+ "[TT]rect[tt] is the ordinary brick shape.",
+ "[TT]tric[tt] is the triclinic unit cell.",
+ "[TT]compact[tt] puts all atoms at the closest distance from the center",
+ "of the box. This can be useful for visualizing e.g. truncated octahedra",
+ "or rhombic dodecahedra. The center for options [TT]tric[tt] and [TT]compact[tt]",
+ "is [TT]tric[tt] (see below), unless the option [TT]-boxcenter[tt]",
+ "is set differently.[PAR]",
+
+ "Option [TT]-center[tt] centers the system in the box. The user can",
+ "select the group which is used to determine the geometrical center.",
+ "Option [TT]-boxcenter[tt] sets the location of the center of the box",
+ "for options [TT]-pbc[tt] and [TT]-center[tt]. The center options are:",
+ "[TT]tric[tt]: half of the sum of the box vectors,",
+ "[TT]rect[tt]: half of the box diagonal,",
+ "[TT]zero[tt]: zero.",
+ "Use option [TT]-pbc mol[tt] in addition to [TT]-center[tt] when you",
+ "want all molecules in the box after the centering.[PAR]",
+
+ "Option [TT]-box[tt] sets the size of the new box. This option only works",
+ "for leading dimensions and is thus generally only useful for rectangular boxes.",
+ "If you want to modify only some of the dimensions, e.g. when reading from",
+ "a trajectory, you can use -1 for those dimensions that should stay the same",
+
+ "It is not always possible to use combinations of [TT]-pbc[tt],",
+ "[TT]-fit[tt], [TT]-ur[tt] and [TT]-center[tt] to do exactly what",
+ "you want in one call to [THISMODULE]. Consider using multiple",
+ "calls, and check out the GROMACS website for suggestions.[PAR]",
+
+ "With [TT]-dt[tt], it is possible to reduce the number of ",
+ "frames in the output. This option relies on the accuracy of the times",
+ "in your input trajectory, so if these are inaccurate use the",
+ "[TT]-timestep[tt] option to modify the time (this can be done",
+ "simultaneously). For making smooth movies, the program [gmx-filter]",
+ "can reduce the number of frames while using low-pass frequency",
+ "filtering, this reduces aliasing of high frequency motions.[PAR]",
+
+ "Using [TT]-trunc[tt] [THISMODULE] can truncate [REF].trr[ref] in place, i.e.",
+ "without copying the file. This is useful when a run has crashed",
+ "during disk I/O (i.e. full disk), or when two contiguous",
+ "trajectories must be concatenated without having double frames.[PAR]",
+
+ "Option [TT]-dump[tt] can be used to extract a frame at or near",
+ "one specific time from your trajectory, but only works reliably",
+ "if the time interval between frames is uniform.[PAR]",
+
+ "Option [TT]-drop[tt] reads an [REF].xvg[ref] file with times and values.",
+ "When options [TT]-dropunder[tt] and/or [TT]-dropover[tt] are set,",
+ "frames with a value below and above the value of the respective options",
+ "will not be written."
+ };
+
+ int pbc_enum;
+ enum
+ {
+ epSel,
+ epNone,
+ epComMol,
+ epComRes,
+ epComAtom,
+ epNojump,
+ epCluster,
+ epWhole,
+ epNR
+ };
+ const char *pbc_opt[epNR + 1] =
+ {
+ nullptr, "none", "mol", "res", "atom", "nojump", "cluster", "whole",
+ nullptr
+ };
+
+ int unitcell_enum;
+ const char *unitcell_opt[euNR+1] =
+ { nullptr, "rect", "tric", "compact", nullptr };
+
+ enum
+ {
+ ecSel, ecTric, ecRect, ecZero, ecNR
+ };
+ const char *center_opt[ecNR+1] =
+ { nullptr, "tric", "rect", "zero", nullptr };
+ int ecenter;
+
+ int fit_enum;
+ enum
+ {
+ efSel, efNone, efFit, efFitXY, efReset, efResetXY, efPFit, efNR
+ };
+ const char *fit[efNR + 1] =
+ {
+ nullptr, "none", "rot+trans", "rotxy+transxy", "translation", "transxy",
+ "progressive", nullptr
+ };
+
+ static gmx_bool bSeparate = FALSE, bVels = TRUE, bForce = FALSE, bCONECT = FALSE;
+ static gmx_bool bCenter = FALSE;
+ static int skip_nr = 1, ndec = 3, nzero = 0;
+ static real tzero = 0, delta_t = 0, timestep = 0, ttrunc = -1, tdump = -1, split_t = 0;
+ static rvec newbox = {0, 0, 0}, shift = {0, 0, 0}, trans = {0, 0, 0};
+ static char *exec_command = nullptr;
+ static real dropunder = 0, dropover = 0;
+ static gmx_bool bRound = FALSE;
+
+ t_pargs
+ pa[] =
+ {
+ { "-skip", FALSE, etINT,
+ { &skip_nr }, "Only write every nr-th frame" },
+ { "-dt", FALSE, etTIME,
+ { &delta_t },
+ "Only write frame when t MOD dt = first time (%t)" },
+ { "-round", FALSE, etBOOL,
+ { &bRound }, "Round measurements to nearest picosecond"},
+ { "-dump", FALSE, etTIME,
+ { &tdump }, "Dump frame nearest specified time (%t)" },
+ { "-t0", FALSE, etTIME,
+ { &tzero },
+ "Starting time (%t) (default: don't change)" },
+ { "-timestep", FALSE, etTIME,
+ { ×tep },
+ "Change time step between input frames (%t)" },
+ { "-pbc", FALSE, etENUM,
+ { pbc_opt },
+ "PBC treatment (see help text for full description)" },
+ { "-ur", FALSE, etENUM,
+ { unitcell_opt }, "Unit-cell representation" },
+ { "-center", FALSE, etBOOL,
+ { &bCenter }, "Center atoms in box" },
+ { "-boxcenter", FALSE, etENUM,
+ { center_opt }, "Center for -pbc and -center" },
+ { "-box", FALSE, etRVEC,
+ { newbox },
+ "Size for new cubic box (default: read from input)" },
+ { "-trans", FALSE, etRVEC,
+ { trans },
+ "All coordinates will be translated by trans. This "
+ "can advantageously be combined with -pbc mol -ur "
+ "compact." },
+ { "-shift", FALSE, etRVEC,
+ { shift },
+ "All coordinates will be shifted by framenr*shift" },
+ { "-fit", FALSE, etENUM,
+ { fit },
+ "Fit molecule to ref structure in the structure file" },
+ { "-ndec", FALSE, etINT,
+ { &ndec },
+ "Number of decimal places to write to .xtc output" },
+ { "-vel", FALSE, etBOOL,
+ { &bVels }, "Read and write velocities if possible" },
+ { "-force", FALSE, etBOOL,
+ { &bForce }, "Read and write forces if possible" },
+ { "-trunc", FALSE, etTIME,
+ { &ttrunc },
+ "Truncate input trajectory file after this time (%t)" },
+ { "-exec", FALSE, etSTR,
+ { &exec_command },
+ "Execute command for every output frame with the "
+ "frame number as argument" },
+ { "-split", FALSE, etTIME,
+ { &split_t },
+ "Start writing new file when t MOD split = first "
+ "time (%t)" },
+ { "-sep", FALSE, etBOOL,
+ { &bSeparate },
+ "Write each frame to a separate .gro, .g96 or .pdb "
+ "file" },
+ { "-nzero", FALSE, etINT,
+ { &nzero },
+ "If the -sep flag is set, use these many digits "
+ "for the file numbers and prepend zeros as needed" },
+ { "-dropunder", FALSE, etREAL,
+ { &dropunder }, "Drop all frames below this value" },
+ { "-dropover", FALSE, etREAL,
+ { &dropover }, "Drop all frames above this value" },
+ { "-conect", FALSE, etBOOL,
+ { &bCONECT },
+ "Add conect records when writing [REF].pdb[ref] files. Useful "
+ "for visualization of non-standard molecules, e.g. "
+ "coarse grained ones" }
+ };
+#define NPA asize(pa)
+
+ FILE *out = nullptr;
+ t_trxstatus *trxout = nullptr;
+ t_trxstatus *trxin;
+ int file_nr;
+ t_trxframe fr, frout;
+ int flags;
+ rvec *xmem = nullptr, *vmem = nullptr, *fmem = nullptr;
+ rvec *xp = nullptr, x_shift, hbox;
+ real *w_rls = nullptr;
+ int m, i, d, frame, outframe, natoms, nout, ncent, newstep = 0, model_nr;
+#define SKIP 10
+ t_topology *top = nullptr;
+ gmx_conect gc = nullptr;
+ int ePBC = -1;
+ t_atoms *atoms = nullptr, useatoms;
+ matrix top_box;
+ int *index = nullptr, *cindex = nullptr;
+ char *grpnm = nullptr;
+ int *frindex, nrfri;
+ char *frname;
+ int ifit, my_clust = -1;
+ int *ind_fit;
+ char *gn_fit;
+ t_cluster_ndx *clust = nullptr;
+ t_trxstatus **clust_status = nullptr;
+ int *clust_status_id = nullptr;
+ int ntrxopen = 0;
+ int *nfwritten = nullptr;
+ int ndrop = 0, ncol, drop0 = 0, drop1 = 0, dropuse = 0;
+ double **dropval;
+ real tshift = 0, dt = -1, prec;
+ gmx_bool bFit, bPFit, bReset;
+ int nfitdim;
+ gmx_rmpbc_t gpbc = nullptr;
+ gmx_bool bRmPBC, bPBCWhole, bPBCcomRes, bPBCcomMol, bPBCcomAtom, bPBC, bNoJump, bCluster;
+ gmx_bool bCopy, bDoIt, bIndex, bTDump, bSetTime, bTPS = FALSE, bDTset = FALSE;
+ gmx_bool bExec, bTimeStep = FALSE, bDumpFrame = FALSE, bSetXtcPrec, bNeedPrec;
+ gmx_bool bHaveFirstFrame, bHaveNextFrame, bSetBox, bSetUR, bSplit = FALSE;
+ gmx_bool bSubTraj = FALSE, bDropUnder = FALSE, bDropOver = FALSE, bTrans = FALSE;
+ gmx_bool bWriteFrame, bSplitHere;
+ const char *top_file, *in_file, *out_file = nullptr;
+ char out_file2[256], *charpt;
+ char *outf_base = nullptr;
+ const char *outf_ext = nullptr;
+ char top_title[256], timestr[32], stepstr[32], filemode[5];
+ gmx_output_env_t *oenv;
+
+ t_filenm fnm[] = {
+ { efTRX, "-f", nullptr, ffREAD },
+ { efTRO, "-o", nullptr, ffWRITE },
+ { efTPS, nullptr, nullptr, ffOPTRD },
+ { efNDX, nullptr, nullptr, ffOPTRD },
+ { efNDX, "-fr", "frames", ffOPTRD },
+ { efNDX, "-sub", "cluster", ffOPTRD },
+ { efXVG, "-drop", "drop", ffOPTRD }
+ };
+#define NFILE asize(fnm)
+
+ if (!parse_common_args(&argc, argv,
+ PCA_CAN_BEGIN | PCA_CAN_END | PCA_CAN_VIEW |
+ PCA_TIME_UNIT,
+ NFILE, fnm, NPA, pa, asize(desc), desc,
+ 0, nullptr, &oenv))
+ {
+ return 0;
+ }
+ fprintf(stdout, "Note that major changes are planned in future for "
+ "trjconv, to improve usability and utility.\n");
+
+ top_file = ftp2fn(efTPS, NFILE, fnm);
+
+ /* Check command line */
+ in_file = opt2fn("-f", NFILE, fnm);
+
+ if (ttrunc != -1)
+ {
+ do_trunc(in_file, ttrunc);
+ }
+ else
+ {
+ /* mark active cmdline options */
+ bSetBox = opt2parg_bSet("-box", NPA, pa);
+ bSetTime = opt2parg_bSet("-t0", NPA, pa);
+ bSetXtcPrec = opt2parg_bSet("-ndec", NPA, pa);
+ bSetUR = opt2parg_bSet("-ur", NPA, pa);
+ bExec = opt2parg_bSet("-exec", NPA, pa);
+ bTimeStep = opt2parg_bSet("-timestep", NPA, pa);
+ bTDump = opt2parg_bSet("-dump", NPA, pa);
+ bDropUnder = opt2parg_bSet("-dropunder", NPA, pa);
+ bDropOver = opt2parg_bSet("-dropover", NPA, pa);
+ bTrans = opt2parg_bSet("-trans", NPA, pa);
+ bSplit = (split_t != 0);
+
+ /* parse enum options */
+ fit_enum = nenum(fit);
+ bFit = (fit_enum == efFit || fit_enum == efFitXY);
+ bReset = (fit_enum == efReset || fit_enum == efResetXY);
+ bPFit = fit_enum == efPFit;
+ pbc_enum = nenum(pbc_opt);
+ bPBCWhole = pbc_enum == epWhole;
+ bPBCcomRes = pbc_enum == epComRes;
+ bPBCcomMol = pbc_enum == epComMol;
+ bPBCcomAtom = pbc_enum == epComAtom;
+ bNoJump = pbc_enum == epNojump;
+ bCluster = pbc_enum == epCluster;
+ bPBC = pbc_enum != epNone;
+ unitcell_enum = nenum(unitcell_opt);
+ ecenter = nenum(center_opt) - ecTric;
+
+ /* set and check option dependencies */
+ if (bPFit)
+ {
+ bFit = TRUE; /* for pfit, fit *must* be set */
+ }
+ if (bFit)
+ {
+ bReset = TRUE; /* for fit, reset *must* be set */
+ }
+ nfitdim = 0;
+ if (bFit || bReset)
+ {
+ nfitdim = (fit_enum == efFitXY || fit_enum == efResetXY) ? 2 : 3;
+ }
+ bRmPBC = bFit || bPBCWhole || bPBCcomRes || bPBCcomMol;
+
+ if (bSetUR)
+ {
+ if (!(bPBCcomRes || bPBCcomMol || bPBCcomAtom))
+ {
+ fprintf(stderr,
+ "WARNING: Option for unitcell representation (-ur %s)\n"
+ " only has effect in combination with -pbc %s, %s or %s.\n"
+ " Ingoring unitcell representation.\n\n",
+ unitcell_opt[0], pbc_opt[2], pbc_opt[3], pbc_opt[4]);
+ }
+ }
+ if (bFit && bPBC)
+ {
+ gmx_fatal(FARGS, "PBC condition treatment does not work together with rotational fit.\n"
+ "Please do the PBC condition treatment first and then run trjconv in a second step\n"
+ "for the rotational fit.\n"
+ "First doing the rotational fit and then doing the PBC treatment gives incorrect\n"
+ "results!");
+ }
+
+ /* ndec for XTC writing is in nr of decimal places, prec is a multiplication factor: */
+ prec = 1;
+ for (i = 0; i < ndec; i++)
+ {
+ prec *= 10;
+ }
+
+ bIndex = ftp2bSet(efNDX, NFILE, fnm);
+
+
+ /* Determine output type */
+ out_file = opt2fn("-o", NFILE, fnm);
+ int ftp = fn2ftp(out_file);
+ fprintf(stderr, "Will write %s: %s\n", ftp2ext(ftp), ftp2desc(ftp));
+ bNeedPrec = (ftp == efXTC);
+ int ftpin = fn2ftp(in_file);
+ if (bVels)
+ {
+ /* check if velocities are possible in input and output files */
+ bVels = (ftp == efTRR || ftp == efGRO ||
+ ftp == efG96 || ftp == efTNG)
+ && (ftpin == efTRR || ftpin == efGRO ||
+ ftpin == efG96 || ftpin == efTNG || ftpin == efCPT);
+ }
+ if (bSeparate || bSplit)
+ {
+ outf_ext = std::strrchr(out_file, '.');
+ if (outf_ext == nullptr)
+ {
+ gmx_fatal(FARGS, "Output file name '%s' does not contain a '.'", out_file);
+ }
+ outf_base = gmx_strdup(out_file);
+ outf_base[outf_ext - out_file] = '\0';
+ }
+
+ bSubTraj = opt2bSet("-sub", NFILE, fnm);
+ if (bSubTraj)
+ {
+ if ((ftp != efXTC) && (ftp != efTRR))
+ {
+ /* It seems likely that other trajectory file types
+ * could work here. */
+ gmx_fatal(FARGS, "Can only use the sub option with output file types "
+ "xtc and trr");
+ }
+ clust = cluster_index(nullptr, opt2fn("-sub", NFILE, fnm));
+
+ /* Check for number of files disabled, as FOPEN_MAX is not the correct
+ * number to check for. In my linux box it is only 16.
+ */
+ if (/* DISABLES CODE */ (false))
+ //if (clust->clust->nr > FOPEN_MAX-4)
+ {
+ gmx_fatal(FARGS, "Can not open enough (%d) files to write all the"
+ " trajectories.\ntry splitting the index file in %d parts.\n"
+ "FOPEN_MAX = %d",
+ clust->clust->nr, 1+clust->clust->nr/FOPEN_MAX, FOPEN_MAX);
+ }
+ gmx_warning("The -sub option could require as many open output files as there are\n"
+ "index groups in the file (%d). If you get I/O errors opening new files,\n"
+ "try reducing the number of index groups in the file, and perhaps\n"
+ "using trjconv -sub several times on different chunks of your index file.\n",
+ clust->clust->nr);
+
+ snew(clust_status, clust->clust->nr);
+ snew(clust_status_id, clust->clust->nr);
+ snew(nfwritten, clust->clust->nr);
+ for (i = 0; (i < clust->clust->nr); i++)
+ {
+ clust_status[i] = nullptr;
+ clust_status_id[i] = -1;
+ }
+ bSeparate = bSplit = FALSE;
+ }
+ /* skipping */
+ if (skip_nr <= 0)
+ {
+ gmx_fatal(FARGS, "Argument for -skip (%d) needs to be greater or equal to 1.", skip_nr);
+ }
+
+ std::unique_ptr<gmx_mtop_t> mtop = read_mtop_for_tng(top_file, in_file, out_file);
+
+ /* Determine whether to read a topology */
+ bTPS = (ftp2bSet(efTPS, NFILE, fnm) ||
+ bRmPBC || bReset || bPBCcomMol || bCluster ||
+ (ftp == efGRO) || (ftp == efPDB) || bCONECT);
+
+ /* Determine if when can read index groups */
+ bIndex = (bIndex || bTPS);
+
+ if (bTPS)
+ {
+ snew(top, 1);
+ read_tps_conf(top_file, top, &ePBC, &xp, nullptr, top_box,
+ bReset || bPBCcomRes);
+ std::strncpy(top_title, *top->name, 255);
+ top_title[255] = '\0';
+ atoms = &top->atoms;
+
+ if (0 == top->mols.nr && (bCluster || bPBCcomMol))
+ {
+ gmx_fatal(FARGS, "Option -pbc %s requires a .tpr file for the -s option", pbc_opt[pbc_enum]);
+ }
+
+ /* top_title is only used for gro and pdb,
+ * the header in such a file is top_title, followed by
+ * t= ... and/or step= ...
+ * to prevent double t= or step=, remove it from top_title.
+ * From GROMACS-2018 we only write t/step when the frame actually
+ * has a valid time/step, so we need to check for both separately.
+ */
+ if ((charpt = std::strstr(top_title, " t= ")))
+ {
+ charpt[0] = '\0';
+ }
+ if ((charpt = std::strstr(top_title, " step= ")))
+ {
+ charpt[0] = '\0';
+ }
+
+ if (bCONECT)
+ {
+ gc = gmx_conect_generate(top);
+ }
+ if (bRmPBC)
+ {
+ gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
+ }
+ }
+
+ /* get frame number index */
+ frindex = nullptr;
+ if (opt2bSet("-fr", NFILE, fnm))
+ {
+ printf("Select groups of frame number indices:\n");
+ rd_index(opt2fn("-fr", NFILE, fnm), 1, &nrfri, &frindex, &frname);
+ if (debug)
+ {
+ for (i = 0; i < nrfri; i++)
+ {
+ fprintf(debug, "frindex[%4d]=%4d\n", i, frindex[i]);
+ }
+ }
+ }
+
+ /* get index groups etc. */
+ if (bReset)
+ {
+ printf("Select group for %s fit\n",
+ bFit ? "least squares" : "translational");
+ get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
+ 1, &ifit, &ind_fit, &gn_fit);
+
+ if (bFit)
+ {
+ if (ifit < 2)
+ {
+ gmx_fatal(FARGS, "Need at least 2 atoms to fit!\n");
+ }
+ else if (ifit == 3)
+ {
+ fprintf(stderr, "WARNING: fitting with only 2 atoms is not unique\n");
+ }
+ }
+ }
+ else if (bCluster)
+ {
+ printf("Select group for clustering\n");
+ get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
+ 1, &ifit, &ind_fit, &gn_fit);
+ }
+
+ if (bIndex)
+ {
+ if (bCenter)
+ {
+ printf("Select group for centering\n");
+ get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
+ 1, &ncent, &cindex, &grpnm);
+ }
+ printf("Select group for output\n");
+ get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm),
+ 1, &nout, &index, &grpnm);
+ }
+ else
+ {
+ /* no index file, so read natoms from TRX */
+ if (!read_first_frame(oenv, &trxin, in_file, &fr, TRX_DONT_SKIP))
+ {
+ gmx_fatal(FARGS, "Could not read a frame from %s", in_file);
+ }
+ natoms = fr.natoms;
+ close_trx(trxin);
+ sfree(fr.x);
+ snew(index, natoms);
+ for (i = 0; i < natoms; i++)
+ {
+ index[i] = i;
+ }
+ nout = natoms;
+ if (bCenter)
+ {
+ ncent = nout;
+ cindex = index;
+ }
+ }
+
+ if (bReset)
+ {
+ snew(w_rls, atoms->nr);
+ for (i = 0; (i < ifit); i++)
+ {
+ w_rls[ind_fit[i]] = atoms->atom[ind_fit[i]].m;
+ }
+
+ /* Restore reference structure and set to origin,
+ store original location (to put structure back) */
+ if (bRmPBC)
+ {
+ gmx_rmpbc(gpbc, top->atoms.nr, top_box, xp);
+ }
+ copy_rvec(xp[index[0]], x_shift);
+ reset_x_ndim(nfitdim, ifit, ind_fit, atoms->nr, nullptr, xp, w_rls);
+ rvec_dec(x_shift, xp[index[0]]);
+ }
+ else
+ {
+ clear_rvec(x_shift);
+ }
+
+ if (bDropUnder || bDropOver)
+ {
+ /* Read the .xvg file with the drop values */
+ fprintf(stderr, "\nReading drop file ...");
+ ndrop = read_xvg(opt2fn("-drop", NFILE, fnm), &dropval, &ncol);
+ fprintf(stderr, " %d time points\n", ndrop);
+ if (ndrop == 0 || ncol < 2)
+ {
+ gmx_fatal(FARGS, "Found no data points in %s",
+ opt2fn("-drop", NFILE, fnm));
+ }
+ drop0 = 0;
+ drop1 = 0;
+ }
+
+ /* Make atoms struct for output in GRO or PDB files */
+ if ((ftp == efGRO) || ((ftp == efG96) && bTPS) || (ftp == efPDB))
+ {
+ /* get memory for stuff to go in .pdb file, and initialize
+ * the pdbinfo structure part if the input has it.
+ */
+ init_t_atoms(&useatoms, atoms->nr, atoms->havePdbInfo);
+ sfree(useatoms.resinfo);
+ useatoms.resinfo = atoms->resinfo;
+ for (i = 0; (i < nout); i++)
+ {
+ useatoms.atomname[i] = atoms->atomname[index[i]];
+ useatoms.atom[i] = atoms->atom[index[i]];
+ if (atoms->havePdbInfo)
+ {
+ useatoms.pdbinfo[i] = atoms->pdbinfo[index[i]];
+ }
+ useatoms.nres = std::max(useatoms.nres, useatoms.atom[i].resind+1);
+ }
+ useatoms.nr = nout;
+ }
+ /* select what to read */
+ if (ftp == efTRR)
+ {
+ flags = TRX_READ_X;
+ }
+ else
+ {
+ flags = TRX_NEED_X;
+ }
+ if (bVels)
+ {
+ flags = flags | TRX_READ_V;
+ }
+ if (bForce)
+ {
+ flags = flags | TRX_READ_F;
+ }
+
+ /* open trx file for reading */
+ bHaveFirstFrame = read_first_frame(oenv, &trxin, in_file, &fr, flags);
+ if (fr.bPrec)
+ {
+ fprintf(stderr, "\nPrecision of %s is %g (nm)\n", in_file, 1/fr.prec);
+ }
+ if (bNeedPrec)
+ {
+ if (bSetXtcPrec || !fr.bPrec)
+ {
+ fprintf(stderr, "\nSetting output precision to %g (nm)\n", 1/prec);
+ }
+ else
+ {
+ fprintf(stderr, "Using output precision of %g (nm)\n", 1/prec);
+ }
+ }
+
+ if (bHaveFirstFrame)
+ {
+ if (bTDump)
+ {
+ // Determine timestep (assuming constant spacing for now) if we
+ // need to dump frames based on time. This is required so we do not
+ // skip the first frame if that was the one that should have been dumped
+ double firstFrameTime = fr.time;
+ if (read_next_frame(oenv, trxin, &fr))
+ {
+ dt = fr.time - firstFrameTime;
+ bDTset = TRUE;
+ if (dt <= 0)
+ {
+ fprintf(stderr, "Warning: Frame times are not incrementing - will dump first frame.\n");
+ }
+ }
+ // Now close and reopen so we are at first frame again
+ close_trx(trxin);
+ done_frame(&fr);
+ // Reopen at first frame (We already know it exists if we got here)
+ read_first_frame(oenv, &trxin, in_file, &fr, flags);
+ }
+
+ set_trxframe_ePBC(&fr, ePBC);
+ natoms = fr.natoms;
+
+ if (bSetTime)
+ {
+ tshift = tzero-fr.time;
+ }
+ else
+ {
+ tzero = fr.time;
+ }
+
+ bCopy = FALSE;
+ if (bIndex)
+ {
+ /* check if index is meaningful */
+ for (i = 0; i < nout; i++)
+ {
+ if (index[i] >= natoms)
+ {
+ gmx_fatal(FARGS,
+ "Index[%d] %d is larger than the number of atoms in the\n"
+ "trajectory file (%d). There is a mismatch in the contents\n"
+ "of your -f, -s and/or -n files.", i, index[i]+1, natoms);
+ }
+ bCopy = bCopy || (i != index[i]);
+ }
+ }
+
+ /* open output for writing */
+ std::strcpy(filemode, "w");
+ switch (ftp)
+ {
+ case efTNG:
+ trxout = trjtools_gmx_prepare_tng_writing(out_file,
+ filemode[0],
+ trxin,
+ nullptr,
+ nout,
+ mtop.get(),
+ gmx::arrayRefFromArray(index, nout),
+ grpnm);
+ break;
+ case efXTC:
+ case efTRR:
+ out = nullptr;
+ if (!bSplit && !bSubTraj)
+ {
+ trxout = open_trx(out_file, filemode);
+ }
+ break;
+ case efGRO:
+ case efG96:
+ case efPDB:
+ if (( !bSeparate && !bSplit ) && !bSubTraj)
+ {
+ out = gmx_ffopen(out_file, filemode);
+ }
+ break;
+ default:
+ gmx_incons("Illegal output file format");
+ }
+
+ if (bCopy)
+ {
+ snew(xmem, nout);
+ if (bVels)
+ {
+ snew(vmem, nout);
+ }
+ if (bForce)
+ {
+ snew(fmem, nout);
+ }
+ }
+
+ /* Start the big loop over frames */
+ file_nr = 0;
+ frame = 0;
+ outframe = 0;
+ model_nr = 0;
+
+ /* Main loop over frames */
+ do
+ {
+ if (!fr.bStep)
+ {
+ /* set the step */
+ fr.step = newstep;
+ newstep++;
+ }
+ if (bSubTraj)
+ {
+ /*if (frame >= clust->clust->nra)
+ gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file\n");*/
+ if (frame > clust->maxframe)
+ {
+ my_clust = -1;
+ }
+ else
+ {
+ my_clust = clust->inv_clust[frame];
+ }
+ if ((my_clust < 0) || (my_clust >= clust->clust->nr) ||
+ (my_clust == -1))
+ {
+ my_clust = -1;
+ }
+ }
+
+ if (bSetBox)
+ {
+ /* generate new box */
+ if (!fr.bBox)
+ {
+ clear_mat(fr.box);
+ }
+ for (m = 0; m < DIM; m++)
+ {
+ if (newbox[m] >= 0)
+ {
+ fr.box[m][m] = newbox[m];
+ }
+ else
+ {
+ if (!fr.bBox)
+ {
+ gmx_fatal(FARGS, "Cannot preserve a box that does not exist.\n");
+ }
+ }
+ }
+ }
+
+ if (bTrans)
+ {
+ for (i = 0; i < natoms; i++)
+ {
+ rvec_inc(fr.x[i], trans);
+ }
+ }
+
+ if (bTDump)
+ {
+ // If we could not read two frames or times are not incrementing
+ // we have almost no idea what to do,
+ // but dump the first frame so output is not broken.
+ if (dt <= 0 || !bDTset)
+ {
+ bDumpFrame = true;
+ }
+ else
+ {
+ // Dump the frame if we are less than half a frame time
+ // below it. This will also ensure we at least dump a
+ // somewhat reasonable frame if the spacing is unequal
+ // and we have overrun the frame time. Once we dump one
+ // frame based on time we quit, so it does not matter
+ // that this might be true for all subsequent frames too.
+ bDumpFrame = (fr.time > tdump-0.5*dt);
+ }
+ }
+ else
+ {
+ bDumpFrame = FALSE;
+ }
+
+ /* determine if an atom jumped across the box and reset it if so */
+ if (bNoJump && (bTPS || frame != 0))
+ {
+ for (d = 0; d < DIM; d++)
+ {
+ hbox[d] = 0.5*fr.box[d][d];
+ }
+ for (i = 0; i < natoms; i++)
+ {
+ if (bReset)
+ {
+ rvec_dec(fr.x[i], x_shift);
+ }
+ for (m = DIM-1; m >= 0; m--)
+ {
+ if (hbox[m] > 0)
+ {
+ while (fr.x[i][m]-xp[i][m] <= -hbox[m])
+ {
+ for (d = 0; d <= m; d++)
+ {
+ fr.x[i][d] += fr.box[m][d];
+ }
+ }
+ while (fr.x[i][m]-xp[i][m] > hbox[m])
+ {
+ for (d = 0; d <= m; d++)
+ {
+ fr.x[i][d] -= fr.box[m][d];
+ }
+ }
+ }
+ }
+ }
+ }
+ else if (bCluster)
+ {
+ calc_pbc_cluster(ecenter, ifit, top, ePBC, fr.x, ind_fit, fr.box);
+ }
+
+ if (bPFit)
+ {
+ /* Now modify the coords according to the flags,
+ for normal fit, this is only done for output frames */
+ if (bRmPBC)
+ {
+ gmx_rmpbc_trxfr(gpbc, &fr);
+ }
+
+ reset_x_ndim(nfitdim, ifit, ind_fit, natoms, nullptr, fr.x, w_rls);
+ do_fit(natoms, w_rls, xp, fr.x);
+ }
+
+ /* store this set of coordinates for future use */
+ if (bPFit || bNoJump)
+ {
+ if (xp == nullptr)
+ {
+ snew(xp, natoms);
+ }
+ for (i = 0; (i < natoms); i++)
+ {
+ copy_rvec(fr.x[i], xp[i]);
+ rvec_inc(fr.x[i], x_shift);
+ }
+ }
+
+ if (frindex)
+ {
+ /* see if we have a frame from the frame index group */
+ for (i = 0; i < nrfri && !bDumpFrame; i++)
+ {
+ bDumpFrame = frame == frindex[i];
+ }
+ }
+ if (debug && bDumpFrame)
+ {
+ fprintf(debug, "dumping %d\n", frame);
+ }
+
+ bWriteFrame =
+ ( ( !bTDump && (frindex == nullptr) && frame % skip_nr == 0 ) || bDumpFrame );
+
+ if (bWriteFrame && (bDropUnder || bDropOver))
+ {
+ while (dropval[0][drop1] < fr.time && drop1+1 < ndrop)
+ {
+ drop0 = drop1;
+ drop1++;
+ }
+ if (std::abs(dropval[0][drop0] - fr.time)
+ < std::abs(dropval[0][drop1] - fr.time))
+ {
+ dropuse = drop0;
+ }
+ else
+ {
+ dropuse = drop1;
+ }
+ if ((bDropUnder && dropval[1][dropuse] < dropunder) ||
+ (bDropOver && dropval[1][dropuse] > dropover))
+ {
+ bWriteFrame = FALSE;
+ }
+ }
+
+ if (bWriteFrame)
+ {
+ /* We should avoid modifying the input frame,
+ * but since here we don't have the output frame yet,
+ * we introduce a temporary output frame time variable.
+ */
+ real frout_time;
+
+ frout_time = fr.time;
+
+ /* calc new time */
+ if (bTimeStep)
+ {
+ frout_time = tzero + frame*timestep;
+ }
+ else if (bSetTime)
+ {
+ frout_time += tshift;
+ }
+
+ if (bTDump)
+ {
+ fprintf(stderr, "\nDumping frame at t= %g %s\n",
+ output_env_conv_time(oenv, frout_time), output_env_get_time_unit(oenv).c_str());
+ }
+
+ /* check for writing at each delta_t */
+ bDoIt = (delta_t == 0);
+ if (!bDoIt)
+ {
+ if (!bRound)
+ {
+ bDoIt = bRmod(frout_time, tzero, delta_t);
+ }
+ else
+ {
+ /* round() is not C89 compatible, so we do this: */
+ bDoIt = bRmod(std::floor(frout_time+0.5), std::floor(tzero+0.5),
+ std::floor(delta_t+0.5));
+ }
+ }
+
+ if (bDoIt || bTDump)
+ {
+ /* print sometimes */
+ if ( ((outframe % SKIP) == 0) || (outframe < SKIP) )
+ {
+ fprintf(stderr, " -> frame %6d time %8.3f \r",
+ outframe, output_env_conv_time(oenv, frout_time));
+ fflush(stderr);
+ }
+
+ if (!bPFit)
+ {
+ /* Now modify the coords according to the flags,
+ for PFit we did this already! */
+
+ if (bRmPBC)
+ {
+ gmx_rmpbc_trxfr(gpbc, &fr);
+ }
+
+ if (bReset)
+ {
+ reset_x_ndim(nfitdim, ifit, ind_fit, natoms, nullptr, fr.x, w_rls);
+ if (bFit)
+ {
+ do_fit_ndim(nfitdim, natoms, w_rls, xp, fr.x);
+ }
+ if (!bCenter)
+ {
+ for (i = 0; i < natoms; i++)
+ {
+ rvec_inc(fr.x[i], x_shift);
+ }
+ }
+ }
+
+ if (bCenter)
+ {
+ center_x(ecenter, fr.x, fr.box, natoms, ncent, cindex);
+ }
+ }
+
+ auto positionsArrayRef = gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec *>(fr.x), natoms);
+ if (bPBCcomAtom)
+ {
+ switch (unitcell_enum)
+ {
+ case euRect:
+ put_atoms_in_box(ePBC, fr.box, positionsArrayRef);
+ break;
+ case euTric:
+ put_atoms_in_triclinic_unitcell(ecenter, fr.box, positionsArrayRef);
+ break;
+ case euCompact:
+ put_atoms_in_compact_unitcell(ePBC, ecenter, fr.box,
+ positionsArrayRef);
+ break;
+ }
+ }
+ if (bPBCcomRes)
+ {
+ put_residue_com_in_box(unitcell_enum, ecenter,
+ natoms, atoms->atom, ePBC, fr.box, fr.x);
+ }
+ if (bPBCcomMol)
+ {
+ put_molecule_com_in_box(unitcell_enum, ecenter,
+ &top->mols,
+ natoms, atoms->atom, ePBC, fr.box, fr.x);
+ }
+ /* Copy the input trxframe struct to the output trxframe struct */
+ frout = fr;
+ frout.time = frout_time;
+ frout.bV = (frout.bV && bVels);
+ frout.bF = (frout.bF && bForce);
+ frout.natoms = nout;
+ if (bNeedPrec && (bSetXtcPrec || !fr.bPrec))
+ {
+ frout.bPrec = TRUE;
+ frout.prec = prec;
+ }
+ if (bCopy)
+ {
+ frout.x = xmem;
+ if (frout.bV)
+ {
+ frout.v = vmem;
+ }
+ if (frout.bF)
+ {
+ frout.f = fmem;
+ }
+ for (i = 0; i < nout; i++)
+ {
+ copy_rvec(fr.x[index[i]], frout.x[i]);
+ if (frout.bV)
+ {
+ copy_rvec(fr.v[index[i]], frout.v[i]);
+ }
+ if (frout.bF)
+ {
+ copy_rvec(fr.f[index[i]], frout.f[i]);
+ }
+ }
+ }
+
+ if (opt2parg_bSet("-shift", NPA, pa))
+ {
+ for (i = 0; i < nout; i++)
+ {
+ for (d = 0; d < DIM; d++)
+ {
+ frout.x[i][d] += outframe*shift[d];
+ }
+ }
+ }
+
+ if (!bRound)
+ {
+ bSplitHere = bSplit && bRmod(frout.time, tzero, split_t);
+ }
+ else
+ {
+ /* round() is not C89 compatible, so we do this: */
+ bSplitHere = bSplit && bRmod(std::floor(frout.time+0.5),
+ std::floor(tzero+0.5),
+ std::floor(split_t+0.5));
+ }
+ if (bSeparate || bSplitHere)
+ {
+ mk_filenm(outf_base, ftp2ext(ftp), nzero, file_nr, out_file2);
+ }
+
+ std::string title;
+ switch (ftp)
+ {
+ case efTNG:
+ write_tng_frame(trxout, &frout);
+ // TODO when trjconv behaves better: work how to read and write lambda
+ break;
+ case efTRR:
+ case efXTC:
+ if (bSplitHere)
+ {
+ if (trxout)
+ {
+ close_trx(trxout);
+ }
+ trxout = open_trx(out_file2, filemode);
+ }
+ if (bSubTraj)
+ {
+ if (my_clust != -1)
+ {
+ char buf[STRLEN];
+ if (clust_status_id[my_clust] == -1)
+ {
+ sprintf(buf, "%s.%s", clust->grpname[my_clust], ftp2ext(ftp));
+ clust_status[my_clust] = open_trx(buf, "w");
+ clust_status_id[my_clust] = 1;
+ ntrxopen++;
+ }
+ else if (clust_status_id[my_clust] == -2)
+ {
+ gmx_fatal(FARGS, "File %s.xtc should still be open (%d open .xtc files)\n" "in order to write frame %d. my_clust = %d",
+ clust->grpname[my_clust], ntrxopen, frame,
+ my_clust);
+ }
+ write_trxframe(clust_status[my_clust], &frout, gc);
+ nfwritten[my_clust]++;
+ if (nfwritten[my_clust] ==
+ (clust->clust->index[my_clust+1]-
+ clust->clust->index[my_clust]))
+ {
+ close_trx(clust_status[my_clust]);
+ clust_status[my_clust] = nullptr;
+ clust_status_id[my_clust] = -2;
+ ntrxopen--;
+ if (ntrxopen < 0)
+ {
+ gmx_fatal(FARGS, "Less than zero open .xtc files!");
+ }
+ }
+ }
+ }
+ else
+ {
+ write_trxframe(trxout, &frout, gc);
+ }
+ break;
+ case efGRO:
+ case efG96:
+ case efPDB:
+ // Only add a generator statement if title is empty,
+ // to avoid multiple generated-by statements from various programs
+ if (std::strlen(top_title) == 0)
+ {
+ sprintf(top_title, "Generated by trjconv");
+ }
+ if (frout.bTime)
+ {
+ sprintf(timestr, " t= %9.5f", frout.time);
+ }
+ else
+ {
+ std::strcpy(timestr, "");
+ }
+ if (frout.bStep)
+ {
+ sprintf(stepstr, " step= %" PRId64, frout.step);
+ }
+ else
+ {
+ std::strcpy(stepstr, "");
+ }
+ title = gmx::formatString("%s%s%s", top_title, timestr, stepstr);
+ if (bSeparate || bSplitHere)
+ {
+ out = gmx_ffopen(out_file2, "w");
+ }
+ switch (ftp)
+ {
+ case efGRO:
+ write_hconf_p(out, title.c_str(), &useatoms,
+ frout.x, frout.bV ? frout.v : nullptr, frout.box);
+ break;
+ case efPDB:
+ fprintf(out, "REMARK GENERATED BY TRJCONV\n");
+ /* if reading from pdb, we want to keep the original
+ model numbering else we write the output frame
+ number plus one, because model 0 is not allowed in pdb */
+ if (ftpin == efPDB && fr.bStep && fr.step > model_nr)
+ {
+ model_nr = fr.step;
+ }
+ else
+ {
+ model_nr++;
+ }
+ write_pdbfile(out, title.c_str(), &useatoms, frout.x,
+ frout.ePBC, frout.box, ' ', model_nr, gc);
+ break;
+ case efG96:
+ const char *outputTitle = "";
+ if (bSeparate || bTDump)
+ {
+ outputTitle = title.c_str();
+ if (bTPS)
+ {
+ frout.bAtoms = TRUE;
+ }
+ frout.atoms = &useatoms;
+ frout.bStep = FALSE;
+ frout.bTime = FALSE;
+ }
+ else
+ {
+ if (outframe == 0)
+ {
+ outputTitle = title.c_str();
+ }
+ frout.bAtoms = FALSE;
+ frout.bStep = TRUE;
+ frout.bTime = TRUE;
+ }
+ write_g96_conf(out, outputTitle, &frout, -1, nullptr);
+ }
+ if (bSeparate || bSplitHere)
+ {
+ gmx_ffclose(out);
+ out = nullptr;
+ }
+ break;
+ default:
+ gmx_fatal(FARGS, "DHE, ftp=%d\n", ftp);
+ }
+ if (bSeparate || bSplitHere)
+ {
+ file_nr++;
+ }
+
+ /* execute command */
+ if (bExec)
+ {
+ char c[255];
+ sprintf(c, "%s %d", exec_command, file_nr-1);
+ /*fprintf(stderr,"Executing '%s'\n",c);*/
+ if (0 != system(c))
+ {
+ gmx_fatal(FARGS, "Error executing command: %s", c);
+ }
+ }
+ outframe++;
+ }
+ }
+ frame++;
+ bHaveNextFrame = read_next_frame(oenv, trxin, &fr);
+ }
+ while (!(bTDump && bDumpFrame) && bHaveNextFrame);
+ }
+
+ if (!bHaveFirstFrame || (bTDump && !bDumpFrame))
+ {
+ fprintf(stderr, "\nWARNING no output, "
+ "last frame read at t=%g\n", fr.time);
+ }
+ fprintf(stderr, "\n");
+
+ close_trx(trxin);
+ sfree(outf_base);
+
+ if (bRmPBC)
+ {
+ gmx_rmpbc_done(gpbc);
+ }
+
+ if (trxout)
+ {
+ close_trx(trxout);
+ }
+ else if (out != nullptr)
+ {
+ gmx_ffclose(out);
+ }
+ if (bSubTraj)
+ {
+ for (i = 0; (i < clust->clust->nr); i++)
+ {
+ if (clust_status_id[i] >= 0)
+ {
+ close_trx(clust_status[i]);
+ }
+ }
+ }
+ }
+
+ if (bTPS)
+ {
+ done_top(top);
+ sfree(top);
+ }
+ sfree(xp);
+ sfree(xmem);
+ sfree(vmem);
+ sfree(fmem);
+ sfree(grpnm);
+ sfree(index);
+ sfree(cindex);
+ done_frame(&fr);
+
+ do_view(oenv, out_file, nullptr);
+
+ output_env_done(oenv);
+ return 0;
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff