2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2019, 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
47 #include "coordinatefile.h"
51 #include "gromacs/coordinateio/outputadapters.h"
52 #include "gromacs/coordinateio/requirements.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/trajectory/trajectoryframe.h"
55 #include "gromacs/utility/exceptions.h"
61 * Get the internal file type from the \p filename.
63 * \param[in] filename Filename of output file.
64 * \throws InvalidInputError When unable to work on an emoty file name.
65 * \returns integer value of file type.
67 static int getFileType(const std::string& filename)
70 if (!filename.empty())
72 filetype = fn2ftp(filename.c_str());
76 GMX_THROW(InvalidInputError("Can not open file with an empty name"));
82 * Get the flag representing the requirements for a given file output.
84 * Also checks if the supplied topology is sufficient through the pointer
87 * \param[in] filetype Internal file type used to check requirements.
88 * \throws InvalidInputError When encountering an invalid file type.
89 * \returns Requirements represent by the bitmask in the return type.
91 static unsigned long getSupportedOutputAdapters(int filetype)
93 unsigned long supportedOutputAdapters = 0;
94 supportedOutputAdapters |= convertFlag(CoordinateFileFlags::Base);
98 supportedOutputAdapters |=
99 (convertFlag(CoordinateFileFlags::RequireForceOutput)
100 | convertFlag(CoordinateFileFlags::RequireVelocityOutput)
101 | convertFlag(CoordinateFileFlags::RequireAtomConnections)
102 | convertFlag(CoordinateFileFlags::RequireAtomInformation)
103 | convertFlag(CoordinateFileFlags::RequireChangedOutputPrecision));
106 supportedOutputAdapters |= (convertFlag(CoordinateFileFlags::RequireAtomConnections)
107 | convertFlag(CoordinateFileFlags::RequireAtomInformation));
110 supportedOutputAdapters |= (convertFlag(CoordinateFileFlags::RequireAtomInformation)
111 | convertFlag(CoordinateFileFlags::RequireVelocityOutput));
114 supportedOutputAdapters |= (convertFlag(CoordinateFileFlags::RequireForceOutput)
115 | convertFlag(CoordinateFileFlags::RequireVelocityOutput));
118 supportedOutputAdapters |= (convertFlag(CoordinateFileFlags::RequireChangedOutputPrecision));
121 default: GMX_THROW(InvalidInputError("Invalid file type"));
123 return supportedOutputAdapters;
127 * Creates a new container object with the user requested IOutputAdapter derived
128 * methods attached to it.
130 * \param[in] requirements Specifications for modules to add.
131 * \param[in] atoms Local copy of atom information to use.
132 * \param[in] sel Selection to use for choosing atoms to write out.
133 * \param[in] abilities Specifications for what the output method can do.
134 * \returns New container for IoutputAdapter derived methods.
136 static OutputAdapterContainer addOutputAdapters(const OutputRequirements& requirements,
138 const Selection& sel,
139 unsigned long abilities)
141 OutputAdapterContainer output(abilities);
143 /* An adapter only gets added if the user has specified a non-default
144 * behaviour. In most cases, this behaviour is passive, meaning that
145 * output gets written if it exists in the input and if the output
148 if (requirements.velocity != ChangeSettingType::PreservedIfPresent)
150 output.addAdapter(std::make_unique<SetVelocities>(requirements.velocity),
151 CoordinateFileFlags::RequireVelocityOutput);
153 if (requirements.force != ChangeSettingType::PreservedIfPresent)
155 output.addAdapter(std::make_unique<SetForces>(requirements.force),
156 CoordinateFileFlags::RequireForceOutput);
158 if (requirements.precision != ChangeFrameInfoType::PreservedIfPresent)
160 output.addAdapter(std::make_unique<SetPrecision>(requirements.prec),
161 CoordinateFileFlags::RequireChangedOutputPrecision);
163 if (requirements.atoms != ChangeAtomsType::PreservedIfPresent)
165 output.addAdapter(std::make_unique<SetAtoms>(requirements.atoms, std::move(atoms)),
166 CoordinateFileFlags::RequireAtomInformation);
168 if (requirements.frameTime != ChangeFrameTimeType::PreservedIfPresent)
170 switch (requirements.frameTime)
172 case (ChangeFrameTimeType::StartTime):
173 output.addAdapter(std::make_unique<SetStartTime>(requirements.startTimeValue),
174 CoordinateFileFlags::RequireNewFrameStartTime);
176 case (ChangeFrameTimeType::TimeStep):
177 output.addAdapter(std::make_unique<SetTimeStep>(requirements.timeStepValue),
178 CoordinateFileFlags::RequireNewFrameTimeStep);
180 case (ChangeFrameTimeType::Both):
181 output.addAdapter(std::make_unique<SetStartTime>(requirements.startTimeValue),
182 CoordinateFileFlags::RequireNewFrameStartTime);
183 output.addAdapter(std::make_unique<SetTimeStep>(requirements.timeStepValue),
184 CoordinateFileFlags::RequireNewFrameTimeStep);
189 if (requirements.box != ChangeFrameInfoType::PreservedIfPresent)
191 output.addAdapter(std::make_unique<SetBox>(requirements.newBox), CoordinateFileFlags::RequireNewBox);
195 output.addAdapter(std::make_unique<OutputSelector>(sel),
196 CoordinateFileFlags::RequireCoordinateSelection);
201 std::unique_ptr<TrajectoryFrameWriter> createTrajectoryFrameWriter(const gmx_mtop_t* top,
202 const Selection& sel,
203 const std::string& filename,
205 OutputRequirements requirements)
208 * Currently the requirements object is expected to be processed and valid,
209 * meaning that e.g. a new box is specified if requested by the option,
210 * or that time values have been set if the corresponding values are set.
211 * This will need to get revisited when the code that builds this object from
212 * the user options gets merged.
214 int filetype = getFileType(filename);
215 unsigned long abilities = getSupportedOutputAdapters(filetype);
217 // first, check if we have a special output format that needs atoms
218 if ((filetype == efPDB) || (filetype == efGRO))
220 if (requirements.atoms == ChangeAtomsType::Never)
223 InconsistentInputError("Can not write to PDB or GRO when"
224 "explicitly turning atom information off"));
226 if (requirements.atoms != ChangeAtomsType::AlwaysFromStructure)
228 requirements.atoms = ChangeAtomsType::Always;
231 OutputAdapterContainer outputAdapters =
232 addOutputAdapters(requirements, std::move(atoms), sel, abilities);
234 TrajectoryFrameWriterPointer trajectoryFrameWriter(
235 new TrajectoryFrameWriter(filename, filetype, sel, top, std::move(outputAdapters)));
236 return trajectoryFrameWriter;
241 * Create a deep copy of a t_trxframe \p input into \p copy
243 * When running the analysis tools and changing values with the
244 * outputadapters, a deep copy of the \p input coordinate frame has to be
245 * created first to ensure that the data is not changed if it is needed for other
246 * tools following with analysis later. Therefore, the data is passed
247 * to \p copy by performing a deep copy first.
249 * The method allocates new storage for coordinates of the x, v, and f arrays
250 * in the new coordinate frame. This means that those arrays need to be free'd
251 * after the frame has been processed and been written to disk.
253 * \param[in] input Reference input coordinate frame.
254 * \param[in,out] copy Pointer to new output frame that will receive the deep copy.
255 * \param[in] xvec Pointer to local coordinate storage vector.
256 * \param[in] vvec Pointer to local velocity storage vector.
257 * \param[in] fvec Pointer to local force storage vector.
258 * \param[in] indexvec Pointer to local index storage vector.
261 deepCopy_t_trxframe(const t_trxframe& input, t_trxframe* copy, RVec* xvec, RVec* vvec, RVec* fvec, int* indexvec)
263 copy->not_ok = input.not_ok;
264 copy->bStep = input.bStep;
265 copy->bTime = input.bTime;
266 copy->bLambda = input.bLambda;
267 copy->bFepState = input.bFepState;
268 copy->bAtoms = input.bAtoms;
269 copy->bPrec = input.bPrec;
273 copy->bBox = input.bBox;
274 copy->bDouble = input.bDouble;
275 copy->natoms = input.natoms;
276 copy->step = input.step;
277 copy->time = input.time;
278 copy->lambda = input.lambda;
279 copy->fep_state = input.fep_state;
282 copy->atoms = input.atoms;
284 copy->prec = input.prec;
287 copy->x = as_rvec_array(xvec);
291 copy->v = as_rvec_array(vvec);
295 copy->f = as_rvec_array(fvec);
299 copy->index = indexvec;
303 copy->index = nullptr;
305 for (int i = 0; i < copy->natoms; i++)
309 copy_rvec(input.x[i], copy->x[i]);
313 copy_rvec(input.v[i], copy->v[i]);
317 copy_rvec(input.f[i], copy->f[i]);
321 copy->index[i] = input.index[i];
324 copy_mat(input.box, copy->box);
325 copy->bPBC = input.bPBC;
326 copy->ePBC = input.ePBC;
330 * Method to open TNG file.
332 * Only need extra method to open this kind of file as it may need access to
333 * a Selection \p sel, if it is valid. Otherwise atom indices will be taken
334 * from the topology \p mtop.
336 * \param[in] name Name of the output file.
337 * \param[in] sel Reference to selection.
338 * \param[in] mtop Pointer to topology, tested before that it is valid.
339 * \todo Those should be methods in a replacement for t_trxstatus instead.
341 static t_trxstatus* openTNG(const std::string& name, const Selection& sel, const gmx_mtop_t* mtop)
343 const char* filemode = "w";
346 GMX_ASSERT(sel.hasOnlyAtoms(), "Can only work with selections consisting out of atoms");
347 return trjtools_gmx_prepare_tng_writing(name.c_str(), filemode[0],
348 nullptr, // infile_, //how to get the input file here?
349 nullptr, sel.atomCount(), mtop, sel.atomIndices(),
354 return trjtools_gmx_prepare_tng_writing(name.c_str(), filemode[0],
355 nullptr, // infile_, //how to get the input file here?
356 nullptr, mtop->natoms, mtop, get_atom_index(mtop),
361 TrajectoryFileOpener::~TrajectoryFileOpener()
363 close_trx(outputFile_);
366 t_trxstatus* TrajectoryFileOpener::outputFile()
368 if (outputFile_ == nullptr)
370 const char* filemode = "w";
373 case (efTNG): outputFile_ = openTNG(outputFileName_, sel_, mtop_); break;
378 case (efG96): outputFile_ = open_trx(outputFileName_.c_str(), filemode); break;
379 default: GMX_THROW(InvalidInputError("Invalid file type"));
385 void TrajectoryFrameWriter::prepareAndWriteFrame(const int framenumber, const t_trxframe& input)
387 if (!outputAdapters_.isEmpty())
390 clear_trxframe(&local, true);
391 localX_.resize(input.natoms);
392 localIndex_.resize(input.natoms);
395 localV_.resize(input.natoms);
399 localF_.resize(input.natoms);
401 deepCopy_t_trxframe(input, &local, localX_.data(), localV_.data(), localF_.data(),
403 for (const auto& outputAdapter : outputAdapters_.getAdapters())
407 outputAdapter->processFrame(framenumber, &local);
410 write_trxframe(file_.outputFile(), &local, nullptr);
414 write_trxframe(file_.outputFile(), const_cast<t_trxframe*>(&input), nullptr);