src/gromacs/coordinateio/coordinatefile.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020,2021, 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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35 /*!\internal
  36  * \file
  37  * \brief
  38  * Implements methods from coordinatefile.h
  39  *
  40  * \author Paul Bauer <paul.bauer.q@gmail.com>
  41  * \ingroup module_coordinateio
  42  */
  43
  44 #include "gmxpre.h"
  45
  46 #include "coordinatefile.h"
  47
  48 #include <algorithm>
  49
  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"
  58
  59 namespace gmx
  60 {
  61
  62 /*!\brief
  63  *  Get the internal file type from the \p filename.
  64  *
  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.
  68  */
  69 static int getFileType(const std::string& filename)
  70 {
  71     int filetype = efNR;
  72     if (!filename.empty())
  73     {
  74         filetype = fn2ftp(filename.c_str());
  75     }
  76     else
  77     {
  78         GMX_THROW(InvalidInputError("Can not open file with an empty name"));
  79     }
  80     return filetype;
  81 }
  82
  83 /*!\brief
  84  * Get the flag representing the requirements for a given file output.
  85  *
  86  * Also checks if the supplied topology is sufficient through the pointer
  87  * to \p mtop.
  88  *
  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.
  92  */
  93 static unsigned long getSupportedOutputAdapters(int filetype)
  94 {
  95     unsigned long supportedOutputAdapters = 0;
  96     supportedOutputAdapters |= convertFlag(CoordinateFileFlags::Base);
  97     switch (filetype)
  98     {
  99         case (efTNG):
 100             supportedOutputAdapters |=
 101                     (convertFlag(CoordinateFileFlags::RequireForceOutput)
 102                      | convertFlag(CoordinateFileFlags::RequireVelocityOutput)
 103                      | convertFlag(CoordinateFileFlags::RequireAtomConnections)
 104                      | convertFlag(CoordinateFileFlags::RequireAtomInformation)
 105                      | convertFlag(CoordinateFileFlags::RequireChangedOutputPrecision));
 106             break;
 107         case (efPDB):
 108             supportedOutputAdapters |= (convertFlag(CoordinateFileFlags::RequireAtomConnections)
 109                                         | convertFlag(CoordinateFileFlags::RequireAtomInformation));
 110             break;
 111         case (efGRO):
 112             supportedOutputAdapters |= (convertFlag(CoordinateFileFlags::RequireAtomInformation)
 113                                         | convertFlag(CoordinateFileFlags::RequireVelocityOutput));
 114             break;
 115         case (efTRR):
 116             supportedOutputAdapters |= (convertFlag(CoordinateFileFlags::RequireForceOutput)
 117                                         | convertFlag(CoordinateFileFlags::RequireVelocityOutput));
 118             break;
 119         case (efXTC):
 120             supportedOutputAdapters |= (convertFlag(CoordinateFileFlags::RequireChangedOutputPrecision));
 121             break;
 122         case (efG96): break;
 123         default: GMX_THROW(InvalidInputError("Invalid file type"));
 124     }
 125     return supportedOutputAdapters;
 126 }
 127
 128 /*! \brief
 129  * Creates a new container object with the user requested IOutputAdapter derived
 130  * methods attached to it.
 131  *
 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.
 137  */
 138 static OutputAdapterContainer addOutputAdapters(const OutputRequirements& requirements,
 139                                                 AtomsDataPtr              atoms,
 140                                                 const Selection&          sel,
 141                                                 unsigned long             abilities)
 142 {
 143     OutputAdapterContainer output(abilities);
 144
 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
 148      * type supports it.
 149      */
 150     if (requirements.velocity != ChangeSettingType::PreservedIfPresent)
 151     {
 152         output.addAdapter(std::make_unique<SetVelocities>(requirements.velocity),
 153                           CoordinateFileFlags::RequireVelocityOutput);
 154     }
 155     if (requirements.force != ChangeSettingType::PreservedIfPresent)
 156     {
 157         output.addAdapter(std::make_unique<SetForces>(requirements.force),
 158                           CoordinateFileFlags::RequireForceOutput);
 159     }
 160     if (requirements.precision != ChangeFrameInfoType::PreservedIfPresent)
 161     {
 162         output.addAdapter(std::make_unique<SetPrecision>(requirements.prec),
 163                           CoordinateFileFlags::RequireChangedOutputPrecision);
 164     }
 165     if (requirements.atoms != ChangeAtomsType::PreservedIfPresent)
 166     {
 167         output.addAdapter(std::make_unique<SetAtoms>(requirements.atoms, std::move(atoms)),
 168                           CoordinateFileFlags::RequireAtomInformation);
 169     }
 170     if (requirements.frameTime != ChangeFrameTimeType::PreservedIfPresent)
 171     {
 172         switch (requirements.frameTime)
 173         {
 174             case (ChangeFrameTimeType::StartTime):
 175                 output.addAdapter(std::make_unique<SetStartTime>(requirements.startTimeValue),
 176                                   CoordinateFileFlags::RequireNewFrameStartTime);
 177                 break;
 178             case (ChangeFrameTimeType::TimeStep):
 179                 output.addAdapter(std::make_unique<SetTimeStep>(requirements.timeStepValue),
 180                                   CoordinateFileFlags::RequireNewFrameTimeStep);
 181                 break;
 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);
 187                 break;
 188             default: break;
 189         }
 190     }
 191     if (requirements.box != ChangeFrameInfoType::PreservedIfPresent)
 192     {
 193         output.addAdapter(std::make_unique<SetBox>(requirements.newBox), CoordinateFileFlags::RequireNewBox);
 194     }
 195     if (sel.isValid())
 196     {
 197         output.addAdapter(std::make_unique<OutputSelector>(sel),
 198                           CoordinateFileFlags::RequireCoordinateSelection);
 199     }
 200     return output;
 201 }
 202
 203 std::unique_ptr<TrajectoryFrameWriter> createTrajectoryFrameWriter(const gmx_mtop_t*  top,
 204                                                                    const Selection&   sel,
 205                                                                    const std::string& filename,
 206                                                                    AtomsDataPtr       atoms,
 207                                                                    OutputRequirements requirements)
 208 {
 209     /* TODO
 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.
 215      */
 216     int           filetype  = getFileType(filename);
 217     unsigned long abilities = getSupportedOutputAdapters(filetype);
 218
 219     // first, check if we have a special output format that needs atoms
 220     if ((filetype == efPDB) || (filetype == efGRO))
 221     {
 222         if (requirements.atoms == ChangeAtomsType::Never)
 223         {
 224             GMX_THROW(
 225                     InconsistentInputError("Can not write to PDB or GRO when"
 226                                            "explicitly turning atom information off"));
 227         }
 228         if (requirements.atoms != ChangeAtomsType::AlwaysFromStructure)
 229         {
 230             requirements.atoms = ChangeAtomsType::Always;
 231         }
 232     }
 233     OutputAdapterContainer outputAdapters =
 234             addOutputAdapters(requirements, std::move(atoms), sel, abilities);
 235
 236     TrajectoryFrameWriterPointer trajectoryFrameWriter(
 237             new TrajectoryFrameWriter(filename, filetype, sel, top, std::move(outputAdapters)));
 238     return trajectoryFrameWriter;
 239 }
 240
 241
 242 /*! \brief
 243  * Create a deep copy of a t_trxframe \p input into \p copy
 244  *
 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.
 250  *
 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.
 254  *
 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.
 261  */
 262 static void
 263 deepCopy_t_trxframe(const t_trxframe& input, t_trxframe* copy, RVec* xvec, RVec* vvec, RVec* fvec, int* indexvec)
 264 {
 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;
 272     copy->bX        = input.bX;
 273     copy->bV        = input.bV;
 274     copy->bF        = input.bF;
 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;
 282     if (input.bAtoms)
 283     {
 284         copy->atoms = input.atoms;
 285     }
 286     copy->prec = input.prec;
 287     if (copy->bX)
 288     {
 289         copy->x = as_rvec_array(xvec);
 290     }
 291     if (copy->bV)
 292     {
 293         copy->v = as_rvec_array(vvec);
 294     }
 295     if (copy->bF)
 296     {
 297         copy->f = as_rvec_array(fvec);
 298     }
 299     if (input.index)
 300     {
 301         copy->index = indexvec;
 302     }
 303     else
 304     {
 305         copy->index = nullptr;
 306     }
 307     for (int i = 0; i < copy->natoms; i++)
 308     {
 309         if (copy->bX)
 310         {
 311             copy_rvec(input.x[i], copy->x[i]);
 312         }
 313         if (copy->bV)
 314         {
 315             copy_rvec(input.v[i], copy->v[i]);
 316         }
 317         if (copy->bF)
 318         {
 319             copy_rvec(input.f[i], copy->f[i]);
 320         }
 321         if (input.index)
 322         {
 323             copy->index[i] = input.index[i];
 324         }
 325     }
 326     copy_mat(input.box, copy->box);
 327     copy->bPBC    = input.bPBC;
 328     copy->pbcType = input.pbcType;
 329 }
 330
 331 /*! \brief
 332  * Method to open TNG file.
 333  *
 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.
 337  *
 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.
 342  */
 343 static t_trxstatus* openTNG(const std::string& name, const Selection& sel, const gmx_mtop_t* mtop)
 344 {
 345     const char* filemode = "w";
 346     if (sel.isValid())
 347     {
 348         GMX_ASSERT(sel.hasOnlyAtoms(), "Can only work with selections consisting out of atoms");
 349         return trjtools_gmx_prepare_tng_writing(name.c_str(),
 350                                                 filemode[0],
 351                                                 nullptr, // infile_, //how to get the input file here?
 352                                                 nullptr,
 353                                                 sel.atomCount(),
 354                                                 mtop,
 355                                                 sel.atomIndices(),
 356                                                 sel.name());
 357     }
 358     else
 359     {
 360         return trjtools_gmx_prepare_tng_writing(name.c_str(),
 361                                                 filemode[0],
 362                                                 nullptr, // infile_, //how to get the input file here?
 363                                                 nullptr,
 364                                                 mtop->natoms,
 365                                                 mtop,
 366                                                 get_atom_index(*mtop),
 367                                                 "System");
 368     }
 369 }
 370
 371 TrajectoryFileOpener::~TrajectoryFileOpener()
 372 {
 373     close_trx(outputFile_);
 374 }
 375
 376 t_trxstatus* TrajectoryFileOpener::outputFile()
 377 {
 378     if (outputFile_ == nullptr)
 379     {
 380         const char* filemode = "w";
 381         switch (filetype_)
 382         {
 383             case (efTNG): outputFile_ = openTNG(outputFileName_, sel_, mtop_); break;
 384             case (efPDB):
 385             case (efGRO):
 386             case (efTRR):
 387             case (efXTC):
 388             case (efG96): outputFile_ = open_trx(outputFileName_.c_str(), filemode); break;
 389             default: GMX_THROW(InvalidInputError("Invalid file type"));
 390         }
 391     }
 392     return outputFile_;
 393 }
 394
 395 void TrajectoryFrameWriter::prepareAndWriteFrame(const int framenumber, const t_trxframe& input)
 396 {
 397     if (!outputAdapters_.isEmpty())
 398     {
 399         t_trxframe local;
 400         clear_trxframe(&local, true);
 401         localX_.resize(input.natoms);
 402         localIndex_.resize(input.natoms);
 403         if (input.bV)
 404         {
 405             localV_.resize(input.natoms);
 406         }
 407         if (input.bF)
 408         {
 409             localF_.resize(input.natoms);
 410         }
 411         deepCopy_t_trxframe(
 412                 input, &local, localX_.data(), localV_.data(), localF_.data(), localIndex_.data());
 413         for (const auto& outputAdapter : outputAdapters_.getAdapters())
 414         {
 415             if (outputAdapter)
 416             {
 417                 outputAdapter->processFrame(framenumber, &local);
 418             }
 419         }
 420         write_trxframe(file_.outputFile(), &local, nullptr);
 421     }
 422     else
 423     {
 424         write_trxframe(file_.outputFile(), const_cast<t_trxframe*>(&input), nullptr);
 425     }
 426 }
 427
 428 } // namespace gmx