src/gromacs/restraint/restraintmdmodule_impl.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_RESTRAINTMDMODULE_IMPL_H
  37 #define GROMACS_RESTRAINTMDMODULE_IMPL_H
  38
  39 /*! \libinternal \file
  40  * \brief Implementation details for RestraintMDModule
  41  *
  42  * \author M. Eric Irrgang <ericirrgang@gmail.com>
  43  *
  44  * \ingroup module_restraint
  45  */
  46
  47 #include <iostream>
  48 #include <mutex>
  49
  50 #include "gromacs/domdec/domdec_struct.h"
  51 #include "gromacs/domdec/ga2la.h"
  52 #include "gromacs/gmxlib/network.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/mdtypes/commrec.h"
  55 #include "gromacs/mdtypes/iforceprovider.h"
  56 #include "gromacs/mdtypes/imdmodule.h"
  57 #include "gromacs/mdtypes/imdoutputprovider.h"
  58 #include "gromacs/mdtypes/imdpoptionprovider.h"
  59 #include "gromacs/mdtypes/mdatom.h"
  60 #include "gromacs/pbcutil/pbc.h"
  61 #include "gromacs/restraint/restraintpotential.h"
  62 #include "gromacs/utility/arrayref.h"
  63
  64 namespace gmx
  65 {
  66
  67 /*! \libinternal
  68  * \brief Abstraction for a restraint interaction site.
  69  *
  70  * A restraint may operate on a single atom or some other entity, such as a selection of atoms.
  71  * The Restraint implementation is very independent from how coordinates are provided or what they mean.
  72  *
  73  * First implementation can only represent single atoms in a global context.
  74  *
  75  * Ultimately, this should be replaced with a more universal facility for acting
  76  * on distributed atom data or simple transformations thereof.
  77  *
  78  * \inlibraryapi
  79  * \ingroup module_restraint
  80  */
  81 class Site
  82 {
  83     public:
  84         /*! \brief Construct from global atom indices
  85          *
  86          * \param globalIndex Atom index in the global state (as input to the simulation)
  87          */
  88         explicit Site(int globalIndex) :
  89             index_(globalIndex),
  90             r_(0, 0, 0)
  91         {}
  92
  93         /*!
  94          * \brief Explicitly define copies.
  95          *
  96          * Implicit definition is not possible because of the mutex member,
  97          * and a copy constructor is necessary to use Site in a std::vector. Really, we should make
  98          * a copy point to the same implementation object to reuse its cache.
  99
 100          */
 101         Site(const Site &site) :
 102             index_(site.index_),
 103             r_(site.r_)
 104         {}
 105
 106         /*! \brief Disallow assignment.
 107          *
 108          * Assignment doesn't make sense because it implies that a site's meaning is fuzzy.
 109          * If the definition of a site is changing, just make a new site.
 110          * There's nothing to be gained by reusing one or by creating it uninitialized.
 111          */
 112         Site &operator=(const Site &) = delete;
 113
 114         /*!
 115          * \brief Get the global atom index of an atomic site.
 116          *
 117          * \return global index provided at construction.
 118          *
 119          */
 120         int index() const { return index_; }
 121
 122         /*!
 123          * \brief Get the position of this site at time t.
 124          *
 125          * \param cr Communications record.
 126          * \param nx Number of locally available atoms (size of local atom data arrays)
 127          * \param x Array of locally available atom coordinates.
 128          * \param t the current time.
 129          * \return position vector.
 130          *
 131          * \internal
 132          * By providing the current time, we can cache results in order to use them once per timestep.
 133          * In the long term, we would prefer to also allow client code to preregister interest in a
 134          * position at a given time, or issue "futures".
 135          */
 136         RVec centerOfMass(const t_commrec                &cr,
 137                           size_t                          nx,
 138                           ArrayRef<const RVec>            x,
 139                           double               gmx_unused t)
 140         {
 141             // Center of mass to return for the site. Currently the only form of site
 142             // implemented is as a global atomic coordinate.
 143             gmx::RVec r = {0, 0, 0};
 144             if (DOMAINDECOMP(&cr)) // Domain decomposition
 145             {
 146                 // Get global-to-local indexing structure
 147                 auto crossRef = cr.dd->ga2la;
 148                 GMX_ASSERT(crossRef, "Domain decomposition must provide global/local cross-reference.");
 149                 if (const auto localIndex = crossRef->findHome(index_))
 150                 {
 151                     GMX_ASSERT(localIndex, "Expect not to reach this point if findHome does not find index_.");
 152                     GMX_ASSERT(*localIndex < static_cast<decltype(*localIndex)>(nx),
 153                                "We assume that the local index cannot be larger than the number of atoms.");
 154                     GMX_ASSERT(*localIndex >= 0, "localIndex is a signed type, but is assumed >0.");
 155                     // If atom is local, get its location
 156                     copy_rvec(x[*localIndex], r);
 157                 }
 158                 else
 159                 {
 160                     // Nothing to contribute on this rank. Leave position == [0,0,0].
 161                 }
 162                 // AllReduce across the ranks of the simulation to get the center-of-mass
 163                 // of the site locally available everywhere. For single-atom sites, this
 164                 // is trivial: exactly one rank should have a non-zero position.
 165                 // For future multi-atom selections,
 166                 // we will receive weighted center-of-mass contributions from
 167                 // each rank and combine to get the global center of mass.
 168                 // \todo use generalized "pull group" facility when available.
 169                 std::array<double, 3> buffer {{r[0], r[1], r[2]}};
 170                 // This should be an all-reduce sum, which gmx_sumd appears to be.
 171                 gmx_sumd(3, buffer.data(), &cr);
 172                 r[0] = static_cast<real>(buffer[0]);
 173                 r[1] = static_cast<real>(buffer[1]);
 174                 r[2] = static_cast<real>(buffer[2]);
 175
 176             }   // end domain decomposition branch
 177             else
 178             {
 179                 // No DD so all atoms are local.
 180                 copy_rvec(x[index_], r);
 181                 (void)nx;
 182             }
 183             // Update cache and cache status.
 184             copy_rvec(r, r_);
 185
 186             return r_;
 187         }
 188
 189     private:
 190         /*!
 191          * \brief Global index of the single-atom site.
 192          *
 193          * \todo This class should be a specialization of a more general Site data source.
 194          * \todo use LocalAtomSet
 195          */
 196         const int               index_;
 197
 198         /*!
 199          * \brief Last known value of the center-of-mass.
 200          *
 201          * Updated with centerOfMass().
 202          */
 203         RVec                    r_;
 204 };
 205
 206 /*! \internal
 207  * \brief Provide IForceProvider for RestraintMDModuleImpl
 208  *
 209  * Adapter class from IForceProvider to IRestraintPotential.
 210  * Objects of this type are uniquely owned by instances of RestraintMDModuleImpl. The object will
 211  * dispatch calls to IForceProvider->calculateForces() to the functor managed by RestraintMDModuleImpl.
 212  * \ingroup module_restraint
 213  */
 214 class RestraintForceProvider final : public gmx::IForceProvider
 215 {
 216     public:
 217         /*!
 218          * \brief Can only be constructed when initialized from a restraint.
 219          */
 220         RestraintForceProvider() = delete;
 221
 222         ~RestraintForceProvider() = default;
 223
 224         /*!
 225          * \brief RAII construction with an IRestraintPotential
 226          *
 227          * Note, this object must outlive the pointer that will be provided to ForceProviders.
 228          * \param restraint handle to an object providing restraint potential calculation
 229          * \param sites List of atomic site indices
 230          */
 231         explicit RestraintForceProvider(std::shared_ptr<gmx::IRestraintPotential> restraint,
 232                                         const std::vector<int>                   &sites);
 233
 234         /*!
 235          * \brief Implement the IForceProvider interface.
 236          *
 237          * Update the force array with restraint contribution(s) for local atoms.
 238          *
 239          * RestraintForceProvider is implemented with the assumption that few
 240          * restraints apply to many atoms.
 241          * That is, the number of restraints affecting a large number of atoms is small,
 242          * though there may be several restraints that apply to few atoms each.
 243          * Under this assumption, it is considered computationally inexpensive to iterate
 244          * over restraints in an outer loop and iterate over atoms within each restraint.
 245          * This would be an invalid assumption if, say, several restraints applied
 246          * to an entire membrane or the entire solvent group.
 247          *
 248          * If the assumption causes performance problems, we can look for a good
 249          * way to reduce from several restraints in a single pass or a very
 250          * lightweight way to determine whether a given restraint applies to a given atom.
 251          * There is also the notion in the pulling code of a limited number of
 252          * "pull groups" used by the "pull coordinates".
 253          * The right optimization will depend on how the code is being used.
 254          *
 255          * Call the evaluator(s) for the restraints for the configured sites.
 256          * Forces are applied to atoms in the first and last site listed.
 257          * Intermediate sites are used as reference coordinates when the relevant
 258          * vector between sites is on the order of half a box length or otherwise
 259          * ambiguous in the case of periodic boundary conditions.
 260          */
 261         void calculateForces(const ForceProviderInput &forceProviderInput,
 262                              ForceProviderOutput      *forceProviderOutput) override;
 263
 264     private:
 265         std::shared_ptr<gmx::IRestraintPotential> restraint_;
 266         std::vector<Site> sites_;
 267 };
 268
 269 /*! \internal
 270  * \brief IMDModule implementation for RestraintMDModule.
 271  *
 272  * Provides IMDModule interface.
 273  *
 274  * \ingroup module_restraint
 275  */
 276 class RestraintMDModuleImpl final
 277 {
 278     public:
 279         RestraintMDModuleImpl() = delete;
 280         /*!
 281          * \brief Wrap an object implementing IRestraintPotential
 282          *
 283          * \param restraint handle to restraint to wrap.
 284          * \param sites list of sites for framework to process for restraint force calculator.
 285          */
 286         RestraintMDModuleImpl(std::shared_ptr<gmx::IRestraintPotential> restraint,
 287                               const std::vector<int>&sites);
 288
 289         /*!
 290          * \brief Allow moves.
 291          *
 292          * \{
 293          */
 294         RestraintMDModuleImpl(RestraintMDModuleImpl &&) noexcept            = default;
 295         RestraintMDModuleImpl &operator=(RestraintMDModuleImpl &&) noexcept = default;
 296         /*! \} */
 297
 298         ~RestraintMDModuleImpl();
 299
 300         /*!
 301          * \brief Unused implementation of IMDModule aspect
 302          */
 303         IMdpOptionProvider *mdpOptionProvider();
 304
 305         /*!
 306          * \brief Unused implementation of IMDModule aspect
 307          */
 308         IMDOutputProvider *outputProvider();
 309
 310         /*!
 311          * \brief Implement IMDModule interface.
 312          *
 313          * \param forceProviders force module manager in the force record that will call this.
 314          *
 315          * The calling code must ensure that this object stays alive as long as forceProviders needs
 316          * the RestraintForceProvider, since forceProviders can't. Typically that is the duration of a do_md() call.
 317          */
 318         void initForceProviders(ForceProviders *forceProviders);
 319
 320         //! handle to RestraintForceProvider implementation
 321         std::unique_ptr<RestraintForceProvider>  forceProvider_;
 322 };
 323
 324 }      // end namespace gmx
 325
 326 #endif //GROMACS_RESTRAINTMDMODULE_IMPL_H