2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2019,2020, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
38 * Implements methods from coordinatefile.h
40 * \author Paul Bauer <paul.bauer.q@gmail.com>
41 * \ingroup module_coordinateio
46 #include "coordinatefile.h"
50 #include "gromacs/options.h"
51 #include "gromacs/coordinateio/outputadapters.h"
52 #include "gromacs/coordinateio/requirements.h"
53 #include "gromacs/fileio/trxio.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/topology/mtop_util.h"
56 #include "gromacs/trajectory/trajectoryframe.h"
57 #include "gromacs/utility/exceptions.h"
63 * Get the internal file type from the \p filename.
65 * \param[in] filename Filename of output file.
66 * \throws InvalidInputError When unable to work on an emoty file name.
67 * \returns integer value of file type.
69 static int getFileType(const std::string& filename)
72 if (!filename.empty())
74 filetype = fn2ftp(filename.c_str());
78 GMX_THROW(InvalidInputError("Can not open file with an empty name"));
84 * Get the flag representing the requirements for a given file output.
86 * Also checks if the supplied topology is sufficient through the pointer
89 * \param[in] filetype Internal file type used to check requirements.
90 * \throws InvalidInputError When encountering an invalid file type.
91 * \returns Requirements represent by the bitmask in the return type.
93 static unsigned long getSupportedOutputAdapters(int filetype)
95 unsigned long supportedOutputAdapters = 0;
96 supportedOutputAdapters |= convertFlag(CoordinateFileFlags::Base);
100 supportedOutputAdapters |=
101 (convertFlag(CoordinateFileFlags::RequireForceOutput)
102 | convertFlag(CoordinateFileFlags::RequireVelocityOutput)
103 | convertFlag(CoordinateFileFlags::RequireAtomConnections)
104 | convertFlag(CoordinateFileFlags::RequireAtomInformation)
105 | convertFlag(CoordinateFileFlags::RequireChangedOutputPrecision));
108 supportedOutputAdapters |= (convertFlag(CoordinateFileFlags::RequireAtomConnections)
109 | convertFlag(CoordinateFileFlags::RequireAtomInformation));
112 supportedOutputAdapters |= (convertFlag(CoordinateFileFlags::RequireAtomInformation)
113 | convertFlag(CoordinateFileFlags::RequireVelocityOutput));
116 supportedOutputAdapters |= (convertFlag(CoordinateFileFlags::RequireForceOutput)
117 | convertFlag(CoordinateFileFlags::RequireVelocityOutput));
120 supportedOutputAdapters |= (convertFlag(CoordinateFileFlags::RequireChangedOutputPrecision));
123 default: GMX_THROW(InvalidInputError("Invalid file type"));
125 return supportedOutputAdapters;
129 * Creates a new container object with the user requested IOutputAdapter derived
130 * methods attached to it.
132 * \param[in] requirements Specifications for modules to add.
133 * \param[in] atoms Local copy of atom information to use.
134 * \param[in] sel Selection to use for choosing atoms to write out.
135 * \param[in] abilities Specifications for what the output method can do.
136 * \returns New container for IoutputAdapter derived methods.
138 static OutputAdapterContainer addOutputAdapters(const OutputRequirements& requirements,
140 const Selection& sel,
141 unsigned long abilities)
143 OutputAdapterContainer output(abilities);
145 /* An adapter only gets added if the user has specified a non-default
146 * behaviour. In most cases, this behaviour is passive, meaning that
147 * output gets written if it exists in the input and if the output
150 if (requirements.velocity != ChangeSettingType::PreservedIfPresent)
152 output.addAdapter(std::make_unique<SetVelocities>(requirements.velocity),
153 CoordinateFileFlags::RequireVelocityOutput);
155 if (requirements.force != ChangeSettingType::PreservedIfPresent)
157 output.addAdapter(std::make_unique<SetForces>(requirements.force),
158 CoordinateFileFlags::RequireForceOutput);
160 if (requirements.precision != ChangeFrameInfoType::PreservedIfPresent)
162 output.addAdapter(std::make_unique<SetPrecision>(requirements.prec),
163 CoordinateFileFlags::RequireChangedOutputPrecision);
165 if (requirements.atoms != ChangeAtomsType::PreservedIfPresent)
167 output.addAdapter(std::make_unique<SetAtoms>(requirements.atoms, std::move(atoms)),
168 CoordinateFileFlags::RequireAtomInformation);
170 if (requirements.frameTime != ChangeFrameTimeType::PreservedIfPresent)
172 switch (requirements.frameTime)
174 case (ChangeFrameTimeType::StartTime):
175 output.addAdapter(std::make_unique<SetStartTime>(requirements.startTimeValue),
176 CoordinateFileFlags::RequireNewFrameStartTime);
178 case (ChangeFrameTimeType::TimeStep):
179 output.addAdapter(std::make_unique<SetTimeStep>(requirements.timeStepValue),
180 CoordinateFileFlags::RequireNewFrameTimeStep);
182 case (ChangeFrameTimeType::Both):
183 output.addAdapter(std::make_unique<SetStartTime>(requirements.startTimeValue),
184 CoordinateFileFlags::RequireNewFrameStartTime);
185 output.addAdapter(std::make_unique<SetTimeStep>(requirements.timeStepValue),
186 CoordinateFileFlags::RequireNewFrameTimeStep);
191 if (requirements.box != ChangeFrameInfoType::PreservedIfPresent)
193 output.addAdapter(std::make_unique<SetBox>(requirements.newBox), CoordinateFileFlags::RequireNewBox);
197 output.addAdapter(std::make_unique<OutputSelector>(sel),
198 CoordinateFileFlags::RequireCoordinateSelection);
203 std::unique_ptr<TrajectoryFrameWriter> createTrajectoryFrameWriter(const gmx_mtop_t* top,
204 const Selection& sel,
205 const std::string& filename,
207 OutputRequirements requirements)
210 * Currently the requirements object is expected to be processed and valid,
211 * meaning that e.g. a new box is specified if requested by the option,
212 * or that time values have been set if the corresponding values are set.
213 * This will need to get revisited when the code that builds this object from
214 * the user options gets merged.
216 int filetype = getFileType(filename);
217 unsigned long abilities = getSupportedOutputAdapters(filetype);
219 // first, check if we have a special output format that needs atoms
220 if ((filetype == efPDB) || (filetype == efGRO))
222 if (requirements.atoms == ChangeAtomsType::Never)
225 InconsistentInputError("Can not write to PDB or GRO when"
226 "explicitly turning atom information off"));
228 if (requirements.atoms != ChangeAtomsType::AlwaysFromStructure)
230 requirements.atoms = ChangeAtomsType::Always;
233 OutputAdapterContainer outputAdapters =
234 addOutputAdapters(requirements, std::move(atoms), sel, abilities);
236 TrajectoryFrameWriterPointer trajectoryFrameWriter(
237 new TrajectoryFrameWriter(filename, filetype, sel, top, std::move(outputAdapters)));
238 return trajectoryFrameWriter;
243 * Create a deep copy of a t_trxframe \p input into \p copy
245 * When running the analysis tools and changing values with the
246 * outputadapters, a deep copy of the \p input coordinate frame has to be
247 * created first to ensure that the data is not changed if it is needed for other
248 * tools following with analysis later. Therefore, the data is passed
249 * to \p copy by performing a deep copy first.
251 * The method allocates new storage for coordinates of the x, v, and f arrays
252 * in the new coordinate frame. This means that those arrays need to be free'd
253 * after the frame has been processed and been written to disk.
255 * \param[in] input Reference input coordinate frame.
256 * \param[in,out] copy Pointer to new output frame that will receive the deep copy.
257 * \param[in] xvec Pointer to local coordinate storage vector.
258 * \param[in] vvec Pointer to local velocity storage vector.
259 * \param[in] fvec Pointer to local force storage vector.
260 * \param[in] indexvec Pointer to local index storage vector.
263 deepCopy_t_trxframe(const t_trxframe& input, t_trxframe* copy, RVec* xvec, RVec* vvec, RVec* fvec, int* indexvec)
265 copy->not_ok = input.not_ok;
266 copy->bStep = input.bStep;
267 copy->bTime = input.bTime;
268 copy->bLambda = input.bLambda;
269 copy->bFepState = input.bFepState;
270 copy->bAtoms = input.bAtoms;
271 copy->bPrec = input.bPrec;
275 copy->bBox = input.bBox;
276 copy->bDouble = input.bDouble;
277 copy->natoms = input.natoms;
278 copy->step = input.step;
279 copy->time = input.time;
280 copy->lambda = input.lambda;
281 copy->fep_state = input.fep_state;
284 copy->atoms = input.atoms;
286 copy->prec = input.prec;
289 copy->x = as_rvec_array(xvec);
293 copy->v = as_rvec_array(vvec);
297 copy->f = as_rvec_array(fvec);
301 copy->index = indexvec;
305 copy->index = nullptr;
307 for (int i = 0; i < copy->natoms; i++)
311 copy_rvec(input.x[i], copy->x[i]);
315 copy_rvec(input.v[i], copy->v[i]);
319 copy_rvec(input.f[i], copy->f[i]);
323 copy->index[i] = input.index[i];
326 copy_mat(input.box, copy->box);
327 copy->bPBC = input.bPBC;
328 copy->pbcType = input.pbcType;
332 * Method to open TNG file.
334 * Only need extra method to open this kind of file as it may need access to
335 * a Selection \p sel, if it is valid. Otherwise atom indices will be taken
336 * from the topology \p mtop.
338 * \param[in] name Name of the output file.
339 * \param[in] sel Reference to selection.
340 * \param[in] mtop Pointer to topology, tested before that it is valid.
341 * \todo Those should be methods in a replacement for t_trxstatus instead.
343 static t_trxstatus* openTNG(const std::string& name, const Selection& sel, const gmx_mtop_t* mtop)
345 const char* filemode = "w";
348 GMX_ASSERT(sel.hasOnlyAtoms(), "Can only work with selections consisting out of atoms");
349 return trjtools_gmx_prepare_tng_writing(name.c_str(), filemode[0],
350 nullptr, // infile_, //how to get the input file here?
351 nullptr, sel.atomCount(), mtop, sel.atomIndices(),
356 return trjtools_gmx_prepare_tng_writing(name.c_str(), filemode[0],
357 nullptr, // infile_, //how to get the input file here?
358 nullptr, mtop->natoms, mtop, get_atom_index(mtop),
363 TrajectoryFileOpener::~TrajectoryFileOpener()
365 close_trx(outputFile_);
368 t_trxstatus* TrajectoryFileOpener::outputFile()
370 if (outputFile_ == nullptr)
372 const char* filemode = "w";
375 case (efTNG): outputFile_ = openTNG(outputFileName_, sel_, mtop_); break;
380 case (efG96): outputFile_ = open_trx(outputFileName_.c_str(), filemode); break;
381 default: GMX_THROW(InvalidInputError("Invalid file type"));
387 void TrajectoryFrameWriter::prepareAndWriteFrame(const int framenumber, const t_trxframe& input)
389 if (!outputAdapters_.isEmpty())
392 clear_trxframe(&local, true);
393 localX_.resize(input.natoms);
394 localIndex_.resize(input.natoms);
397 localV_.resize(input.natoms);
401 localF_.resize(input.natoms);
403 deepCopy_t_trxframe(input, &local, localX_.data(), localV_.data(), localF_.data(),
405 for (const auto& outputAdapter : outputAdapters_.getAdapters())
409 outputAdapter->processFrame(framenumber, &local);
412 write_trxframe(file_.outputFile(), &local, nullptr);
416 write_trxframe(file_.outputFile(), const_cast<t_trxframe*>(&input), nullptr);