src/gromacs/restraint/restraintpotential.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2018,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
  36 #ifndef GROMACS_RESTRAINT_RESTRAINTPOTENTIAL_H
  37 #define GROMACS_RESTRAINT_RESTRAINTPOTENTIAL_H
  38
  39 /*!
  40  * \defgroup module_restraint MD restraints
  41  * \inpublicapi
  42  * \brief Apply restraints during MD integration.
  43  *
  44  * The classes here are available through the public API. To write a restraint
  45  * plugin, implement gmx::IRestraintPotential with a calculation method that
  46  * produces a PotentialPointData for each site to which MD forces will be applied.
  47  */
  48 /*! \file
  49  * \brief Declare generic interface for restraint implementations.
  50  *
  51  * \author M. Eric Irrgang <ericirrgang@gmail.com>
  52  *
  53  * \inpublicapi
  54  * \ingroup module_restraint
  55  */
  56
  57 #include <functional>
  58 #include <memory>
  59 #include <ostream>
  60 #include <vector>
  61
  62 #include "gromacs/math/vectypes.h"
  63
  64 // TODO: Get from a header once the public API settles down a bit.
  65 namespace gmxapi
  66 {
  67 class SessionResources;
  68 }
  69
  70 namespace gmx
  71 {
  72
  73 /*!
  74  * \brief Provide a vector type name with a more stable interface than RVec and a more stable
  75  * implementation than vec3<>.
  76  *
  77  * \ingroup module_restraint
  78  *
  79  * \internal
  80  * Type alias is used at namespace level
  81  */
  82 using Vector = ::gmx::RVec;
  83
  84 /*!
  85  * \brief Structure to hold the results of IRestraintPotential::evaluate().
  86  *
  87  * \ingroup module_restraint
  88  */
  89 class PotentialPointData
  90 {
  91 public:
  92     /*!
  93      * \brief Initialize a new data structure.
  94      */
  95     PotentialPointData() : PotentialPointData{ Vector(0., 0., 0.), real(0.0) } {}
  96
  97     /*!
  98      * \brief Initialize from an argument list
  99      *
 100      * \param f Force vector.
 101      * \param e Energy value.
 102      *
 103      * Note that if force was calculated as a scalar, it needs to be multiplied by a unit
 104      * vector in the direction to which it should be applied.
 105      */
 106     PotentialPointData(const Vector& f, const real e) : force(f), energy(e) {}
 107
 108     /*!
 109      * \brief Force vector calculated for first position.
 110      */
 111     Vector force;
 112
 113     /*!
 114      * \brief Potential energy calculated for this interaction.
 115      */
 116     real energy;
 117 };
 118
 119 /*!
 120  * \brief Interface for Restraint potentials.
 121  *
 122  * Derive from this interface class to implement a restraint potential. The derived class
 123  * must implement the evaluate() member function that produces a PotentialPointData instance.
 124  * For various convenience functions and to decouple the internal library
 125  * interface from implementations of potentials, it is expected that
 126  * restraints will be programmed by subclassing gmx::RestraintPotential<>
 127  * rather than gmx::IRestraintPotential.
 128  *
 129  * For a set of \f$n\f$ coordinates, generate a force field according to a
 130  * scalar potential that is a fun. \f$F_i = - \nabla_{q_i} \Phi (q_0, q_1, ... q_n; t)\f$
 131  *
 132  * Potentials implemented with these classes may be long ranged and are appropriate
 133  * for only a small number of particles to avoid substantial performance impact.
 134  *
 135  * The indices in the above equation refer to the input and output sites specified
 136  * before the simulation is started. In the current interface, the evaluate() virtual
 137  * function allows an implementer to calculate the energy and force acting at the
 138  * first of a pair of sites. The equal and opposite force is applied at the second site.
 139  *
 140  * The potential function is evaluated with a time argument
 141  * and can be updated during the simulation.
 142  * For non-time-varying potentials, the time argument may still be useful for
 143  * internal optimizations, such as managing cached values.
 144  *
 145  * In the simplest and most common case, pairs of sites (atoms or groups)
 146  * are specified by the user and then, during integration, GROMACS provides
 147  * the current positions of each pair for the restraint potential to be evaluated.
 148  * In such a case, the potential can be implemented by overriding evaluate().
 149  * \todo Template headers can help to build compatible calculation methods with different input
 150  * requirements. For reference, see https://github.com/kassonlab/sample_restraint
 151  *
 152  * \ingroup module_restraint
 153  */
 154 class IRestraintPotential
 155 {
 156 public:
 157     virtual ~IRestraintPotential() = default;
 158
 159     /*!
 160      * \brief Calculate a force vector according to two input positions at a given time.
 161      *
 162      * If not overridden by derived class, returns a zero vector.
 163      * \param r1 position of first site
 164      * \param r2 position of second site
 165      * \param t simulation time in picoseconds
 166      * \return force vector and potential energy to be applied by calling code.
 167      *
 168      * \todo The virtual function call should be replaced by a (table of) function objects retrieved before the run.
 169      */
 170     virtual PotentialPointData evaluate(Vector r1, Vector r2, double t) = 0;
 171
 172
 173     /*!
 174      * \brief Call-back hook for restraint implementations.
 175      *
 176      * An update function to be called on the simulation master rank/thread periodically
 177      * by the Restraint framework.
 178      * Receives the same input as the evaluate() method, but is only called on the master
 179      * rank of a simulation to allow implementation code to be thread-safe without knowing
 180      * anything about the domain decomposition.
 181      *
 182      * \param v position of the first site
 183      * \param v0 position of the second site
 184      * \param t simulation time
 185      *
 186      * \internal
 187      * We give the definition here because we don't want plugins to have to link against
 188      * libgromacs right now (complicated header maintenance and no API stability guarantees).
 189      * But once we've had plugin restraints wrap themselves in a Restraint template,
 190      * we can set update = 0
 191      *
 192      * \todo: Provide gmxapi facility for plugin restraints to wrap themselves
 193      * with a default implementation to let this class be pure virtual.
 194      */
 195     virtual void update(gmx::Vector v, gmx::Vector v0, double t)
 196     {
 197         (void)v;
 198         (void)v0;
 199         (void)t;
 200     }
 201
 202
 203     /*!
 204      * \brief Find out what sites this restraint is configured to act on.
 205      * \return
 206      */
 207     virtual std::vector<int> sites() const = 0;
 208
 209     /*!
 210      * \brief Allow Session-mediated interaction with other resources or workflow elements.
 211      *
 212      * \param resources temporary access to the resources provided by the session for additional configuration.
 213      *
 214      * A module implements this method to receive a handle to resources configured for this particular workflow
 215      * element.
 216      *
 217      * \internal
 218      * \todo This should be more general than the RestraintPotential interface.
 219      */
 220     virtual void bindSession(gmxapi::SessionResources* resources) { (void)resources; }
 221 };
 222
 223 } // end namespace gmx
 224
 225 #endif // GMX_PULLING_PULLPOTENTIAL_H