src/gromacs/coordinateio/coordinatefile.cpp

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