python_packaging/sample_restraint/src/cpp/ensemblepotential.h

   1 #ifndef HARMONICRESTRAINT_ENSEMBLEPOTENTIAL_H
   2 #define HARMONICRESTRAINT_ENSEMBLEPOTENTIAL_H
   3
   4 /*! \file
   5  * \brief Provide restrained ensemble MD potential for GROMACS plugin.
   6  *
   7  * The restraint implemented here uses a facility provided by gmxapi to perform averaging of some
   8  * array data across an ensemble of simulations. Simpler pair restraints can use less of this
   9  * example code.
  10  *
  11  * Contains a lot of boiler plate that is being generalized and migrate out of this file, but other
  12  * pair restraints can be implemented by following the example in this and ``ensemblepotential.cpp``.
  13  * The ``CMakeLists.txt`` file will need to be updated if you add additional source files, and
  14  * ``src/pythonmodule/export_plugin.cpp`` will need to be updated if you add or change the name of
  15  * potentials.
  16  *
  17  * \author M. Eric Irrgang <ericirrgang@gmail.com>
  18  */
  19
  20 #include <array>
  21 #include <memory>
  22 #include <mutex>
  23 #include <vector>
  24
  25 #include "gmxapi/gromacsfwd.h"
  26 #include "gmxapi/session.h"
  27 #include "gmxapi/md/mdmodule.h"
  28
  29 #include "gromacs/restraint/restraintpotential.h"
  30 #include "gromacs/utility/basedefinitions.h"
  31 #include "gromacs/utility/real.h"
  32
  33 #include "sessionresources.h"
  34
  35 namespace plugin
  36 {
  37
  38 // Histogram for a single restrained pair.
  39 using PairHist = std::vector<double>;
  40
  41 struct ensemble_input_param_type
  42 {
  43     /// distance histogram parameters
  44     size_t nBins{ 0 };
  45     double binWidth{ 0. };
  46
  47     /// Flat-bottom potential boundaries.
  48     double minDist{ 0 };
  49     double maxDist{ 0 };
  50
  51     /// Experimental reference distribution.
  52     PairHist experimental{};
  53
  54     /// Number of samples to store during each window.
  55     unsigned int nSamples{ 0 };
  56     double       samplePeriod{ 0 };
  57
  58     /// Number of windows to use for smoothing histogram updates.
  59     unsigned int nWindows{ 0 };
  60
  61     /// Harmonic force coefficient
  62     double k{ 0 };
  63     /// Smoothing factor: width of Gaussian interpolation for histogram
  64     double sigma{ 0 };
  65 };
  66
  67 // \todo We should be able to automate a lot of the parameter setting stuff
  68 // by having the developer specify a map of parameter names and the corresponding type, but that could get tricky.
  69 // The statically compiled fast parameter structure would be generated with a recursive variadic template
  70 // the way a tuple is. ref https://eli.thegreenplace.net/2014/variadic-templates-in-c/
  71
  72 std::unique_ptr<ensemble_input_param_type> makeEnsembleParams(size_t                     nbins,
  73                                                               double                     binWidth,
  74                                                               double                     minDist,
  75                                                               double                     maxDist,
  76                                                               const std::vector<double>& experimental,
  77                                                               unsigned int               nSamples,
  78                                                               double       samplePeriod,
  79                                                               unsigned int nWindows,
  80                                                               double       k,
  81                                                               double       sigma);
  82
  83 /*!
  84  * \brief a residue-pair bias calculator for use in restrained-ensemble simulations.
  85  *
  86  * Applies a force between two sites according to the difference between an experimentally observed
  87  * site pair distance distribution and the distance distribution observed earlier in the simulation
  88  * trajectory. The sampled distribution is averaged from the previous `nwindows` histograms from all
  89  * ensemble members. Each window contains a histogram populated with `nsamples` distances recorded
  90  * at `sample_period` step intervals.
  91  *
  92  * \internal
  93  * During a the window_update_period steps of a window, the potential applied is a harmonic function
  94  * of the difference between the sampled and experimental histograms. At the beginning of the
  95  * window, this difference is found and a Gaussian blur is applied.
  96  */
  97 class EnsemblePotential
  98 {
  99 public:
 100     using input_param_type = ensemble_input_param_type;
 101
 102     /* No default constructor. Parameters must be provided. */
 103     EnsemblePotential() = delete;
 104
 105     /*!
 106      * \brief Constructor called by the wrapper code to produce a new instance.
 107      *
 108      * This constructor is called once per simulation per GROMACS process. Note that until
 109      * gmxapi 0.0.8 there is only one instance per simulation in a thread-MPI simulation.
 110      *
 111      * \param params
 112      */
 113     explicit EnsemblePotential(const input_param_type& params);
 114
 115     /*!
 116      * \brief Deprecated constructor taking a parameter list.
 117      *
 118      * \param nbins
 119      * \param binWidth
 120      * \param minDist
 121      * \param maxDist
 122      * \param experimental
 123      * \param nSamples
 124      * \param samplePeriod
 125      * \param nWindows
 126      * \param k
 127      * \param sigma
 128      */
 129     EnsemblePotential(size_t       nbins,
 130                       double       binWidth,
 131                       double       minDist,
 132                       double       maxDist,
 133                       PairHist     experimental,
 134                       unsigned int nSamples,
 135                       double       samplePeriod,
 136                       unsigned int nWindows,
 137                       double       k,
 138                       double       sigma);
 139
 140     /*!
 141      * \brief Evaluates the pair restraint potential.
 142      *
 143      * In parallel simulations, the gmxapi framework does not make guarantees about where or
 144      * how many times this function is called. It should be simple and stateless; it should not
 145      * update class member data (see ``ensemblepotential.cpp``. For a more controlled API hook
 146      * and to manage state in the object, use ``callback()``.
 147      *
 148      * \param v position of the site for which force is being calculated.
 149      * \param v0 reference site (other member of the pair).
 150      * \param t current simulation time (ps).
 151      * \return container for force and potential energy data.
 152      */
 153     // Implementation note for the future: If dispatching this virtual function is not fast
 154     // enough, the compiler may be able to better optimize a free
 155     // function that receives the current restraint as an argument.
 156     gmx::PotentialPointData calculate(gmx::Vector v, gmx::Vector v0, gmx_unused double t);
 157
 158     /*!
 159      * \brief An update function to be called on the simulation master rank/thread periodically by the Restraint framework.
 160      *
 161      * Defining this function in a plugin potential is optional. If the function is defined,
 162      * the restraint framework calls this function (on the first rank only in a parallel simulation) before calling calculate().
 163      *
 164      * The callback may use resources provided by the Session in the callback to perform updates
 165      * to the local or global state of an ensemble of simulations. Future gmxapi releases will
 166      * include additional optimizations, allowing call-back frequency to be expressed, and more
 167      * general Session resources, as well as more flexible call signatures.
 168      */
 169     void callback(gmx::Vector v, gmx::Vector v0, double t, const Resources& resources);
 170
 171 private:
 172     /// Width of bins (distance) in histogram
 173     size_t nBins_;
 174     double binWidth_;
 175
 176     /// Flat-bottom potential boundaries.
 177     double minDist_;
 178     double maxDist_;
 179     /// Smoothed historic distribution for this restraint. An element of the array of restraints in this simulation.
 180     // Was `hij` in earlier code.
 181     PairHist histogram_;
 182     PairHist experimental_;
 183
 184     /// Number of samples to store during each window.
 185     unsigned int nSamples_;
 186     unsigned int currentSample_;
 187     double       samplePeriod_;
 188     double       nextSampleTime_;
 189     /// Accumulated list of samples during a new window.
 190     std::vector<double> distanceSamples_;
 191
 192     /// Number of windows to use for smoothing histogram updates.
 193     size_t nWindows_;
 194     size_t currentWindow_;
 195     double windowStartTime_;
 196     double nextWindowUpdateTime_;
 197     /// The history of nwindows histograms for this restraint.
 198     std::vector<std::unique_ptr<plugin::Matrix<double>>> windows_;
 199
 200     /// Harmonic force coefficient
 201     double k_;
 202     /// Smoothing factor: width of Gaussian interpolation for histogram
 203     double sigma_;
 204 };
 205
 206 /*!
 207  * \brief Use EnsemblePotential to implement a RestraintPotential
 208  *
 209  * This is boiler plate that will be templated and moved.
 210  */
 211 class EnsembleRestraint : public ::gmx::IRestraintPotential, private EnsemblePotential
 212 {
 213 public:
 214     using EnsemblePotential::input_param_type;
 215
 216     EnsembleRestraint(std::vector<int>           sites,
 217                       const input_param_type&    params,
 218                       std::shared_ptr<Resources> resources) :
 219         EnsemblePotential(params), sites_{ std::move(sites) }, resources_{ std::move(resources) }
 220     {
 221     }
 222
 223     ~EnsembleRestraint() override = default;
 224
 225     /*!
 226      * \brief Implement required interface of gmx::IRestraintPotential
 227      *
 228      * \return list of configured site indices.
 229      *
 230      * \todo remove to template header
 231      * \todo abstraction of site references
 232      */
 233     std::vector<int> sites() const override { return sites_; }
 234
 235     /*!
 236      * \brief Implement the interface gmx::IRestraintPotential
 237      *
 238      * Dispatch to calculate() method.
 239      *
 240      * \param r1 coordinate of first site
 241      * \param r2 reference coordinate (second site)
 242      * \param t simulation time
 243      * \return calculated force and energy
 244      *
 245      * \todo remove to template header.
 246      */
 247     gmx::PotentialPointData evaluate(gmx::Vector r1, gmx::Vector r2, double t) override
 248     {
 249         return calculate(r1, r2, t);
 250     };
 251
 252     /*!
 253      * \brief An update function to be called on the simulation master rank/thread periodically by the Restraint framework.
 254      *
 255      * Implements optional override of gmx::IRestraintPotential::update
 256      *
 257      * This boilerplate will disappear into the Restraint template in an upcoming gmxapi release.
 258      */
 259     void update(gmx::Vector v, gmx::Vector v0, double t) override
 260     {
 261         // Todo: use a callback period to mostly bypass this and avoid excessive mutex locking.
 262         callback(v, v0, t, *resources_);
 263     };
 264
 265     /*!
 266      * \brief Implement the binding protocol that allows access to Session resources.
 267      *
 268      * The client receives a non-owning pointer to the session and cannot extent the life of the
 269      * session. In the future we can use a more formal handle mechanism.
 270      *
 271      * \param session pointer to the current session
 272      */
 273     void bindSession(gmxapi::SessionResources* session) override
 274     {
 275         resources_->setSession(session);
 276     }
 277
 278     void setResources(std::unique_ptr<Resources>&& resources) { resources_ = std::move(resources); }
 279
 280 private:
 281     std::vector<int> sites_;
 282     //        double callbackPeriod_;
 283     //        double nextCallback_;
 284     std::shared_ptr<Resources> resources_;
 285 };
 286
 287
 288 // Important: Just declare the template instantiation here for client code.
 289 // We will explicitly instantiate a definition in the .cpp file where the input_param_type is defined.
 290 extern template class RestraintModule<EnsembleRestraint>;
 291
 292 } // end namespace plugin
 293
 294 #endif // HARMONICRESTRAINT_ENSEMBLEPOTENTIAL_H