2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2018,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.
36 #ifndef GROMACS_RESTRAINT_RESTRAINTPOTENTIAL_H
37 #define GROMACS_RESTRAINT_RESTRAINTPOTENTIAL_H
40 * \defgroup module_restraint MD restraints
42 * \brief Apply restraints during MD integration.
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.
49 * \brief Declare generic interface for restraint implementations.
51 * \author M. Eric Irrgang <ericirrgang@gmail.com>
54 * \ingroup module_restraint
62 #include "gromacs/math/vectypes.h"
64 // TODO: Get from a header once the public API settles down a bit.
67 class SessionResources;
74 * \brief Provide a vector type name with a more stable interface than RVec and a more stable
75 * implementation than vec3<>.
77 * \ingroup module_restraint
80 * Type alias is used at namespace level
82 using Vector = ::gmx::RVec;
85 * \brief Structure to hold the results of IRestraintPotential::evaluate().
87 * \ingroup module_restraint
89 class PotentialPointData
93 * \brief Initialize a new data structure.
95 PotentialPointData() : PotentialPointData{ Vector(0., 0., 0.), real(0.0) } {}
98 * \brief Initialize from an argument list
100 * \param f Force vector.
101 * \param e Energy value.
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.
106 PotentialPointData(const Vector& f, const real e) : force(f), energy(e) {}
109 * \brief Force vector calculated for first position.
114 * \brief Potential energy calculated for this interaction.
120 * \brief Interface for Restraint potentials.
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.
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$
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.
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.
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.
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
152 * \ingroup module_restraint
154 class IRestraintPotential
157 virtual ~IRestraintPotential() = default;
160 * \brief Calculate a force vector according to two input positions at a given time.
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.
168 * \todo The virtual function call should be replaced by a (table of) function objects retrieved before the run.
170 virtual PotentialPointData evaluate(Vector r1, Vector r2, double t) = 0;
174 * \brief Call-back hook for restraint implementations.
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.
182 * \param v position of the first site
183 * \param v0 position of the second site
184 * \param t simulation time
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
192 * \todo: Provide gmxapi facility for plugin restraints to wrap themselves
193 * with a default implementation to let this class be pure virtual.
195 virtual void update(gmx::Vector v, gmx::Vector v0, double t)
204 * \brief Find out what sites this restraint is configured to act on.
207 virtual std::vector<int> sites() const = 0;
210 * \brief Allow Session-mediated interaction with other resources or workflow elements.
212 * \param resources temporary access to the resources provided by the session for additional configuration.
214 * A module implements this method to receive a handle to resources configured for this particular workflow
218 * \todo This should be more general than the RestraintPotential interface.
220 virtual void bindSession(gmxapi::SessionResources* resources) { (void)resources; }
223 } // end namespace gmx
225 #endif // GMX_PULLING_PULLPOTENTIAL_H