src/gromacs/restraint/restraintmdmodule.cpp

   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 #include "gmxpre.h"
  37
  38 #include "restraintmdmodule.h"
  39
  40 #include <memory>
  41
  42 #include "gromacs/mdtypes/forceoutput.h"
  43 #include "gromacs/mdtypes/iforceprovider.h"
  44
  45 #include "restraintmdmodule-impl.h"
  46
  47 namespace gmx
  48 {
  49
  50 RestraintForceProvider::RestraintForceProvider(std::shared_ptr<IRestraintPotential> restraint,
  51                                                const std::vector<int>              &sites) :
  52     restraint_ {std::move(restraint)}
  53 {
  54     GMX_ASSERT(restraint_, "Valid RestraintForceProviders wrap non-null restraints.");
  55     GMX_ASSERT(sites_.empty(), "");
  56     for (auto && site : sites)
  57     {
  58         sites_.emplace_back(site);
  59     }
  60     if (sites_.size() < 2)
  61     {
  62         GMX_THROW(InvalidInputError("Restraints require at least two sites to calculate forces."));
  63     }
  64 }
  65
  66 void RestraintForceProvider::calculateForces(const ForceProviderInput &forceProviderInput,
  67                                              ForceProviderOutput     * forceProviderOutput)
  68 {
  69     GMX_ASSERT(restraint_, "Restraint must be initialized.");
  70
  71     const auto &mdatoms = forceProviderInput.mdatoms_;
  72     GMX_ASSERT(mdatoms.homenr >= 0, "number of home atoms must be non-negative.");
  73
  74     const auto &box = forceProviderInput.box_;
  75     GMX_ASSERT(check_box(-1, box) == nullptr, "Invalid box.");
  76     t_pbc       pbc {};
  77     set_pbc(&pbc,
  78             -1,
  79             box);
  80
  81     const auto &x  = forceProviderInput.x_;
  82     const auto &cr = forceProviderInput.cr_;
  83     const auto &t  = forceProviderInput.t_;
  84     // Cooperatively get Cartesian coordinates for center of mass of each site
  85     RVec        r1 = sites_[0].centerOfMass(cr,
  86                                             static_cast<size_t>(mdatoms.homenr),
  87                                             x,
  88                                             t);
  89     // r2 is to be constructed as
  90     // r2 = (site[N] - site[N-1]) + (site_{N-1} - site_{N-2}) + ... + (site_2 - site_1) + site_1
  91     // where the minimum image convention is applied to each path but not to the overall sum.
  92     // It is redundant to pass both r1 and r2 to called code, and potentially confusing, since
  93     // r1 may refer to an actual coordinate in the simulation while r2 may be in an expanded
  94     // Cartesian coordinate system. Called code should not use r1 and r2 to attempt to identify
  95     // sites in the simulation. If we need that functionality, we should do it separately by
  96     // allowing called code to look up atoms by tag or global index.
  97     RVec r2 = {r1[0], r1[1], r1[2]};
  98     rvec dr = {0, 0, 0};
  99     // Build r2 by following a path of difference vectors that are each presumed to be less than
 100     // a half-box apart, in case we are battling periodic boundary conditions along the lines of
 101     // a big molecule in a small box.
 102     for (size_t i = 0; i < sites_.size() - 1; ++i)
 103     {
 104         RVec a = sites_[i].centerOfMass(cr,
 105                                         static_cast<size_t>(mdatoms.homenr),
 106                                         x,
 107                                         t);
 108         RVec b = sites_[i + 1].centerOfMass(cr,
 109                                             static_cast<size_t>(mdatoms.homenr),
 110                                             x,
 111                                             t);
 112         // dr = minimum_image_vector(b - a)
 113         pbc_dx(&pbc,
 114                b,
 115                a,
 116                dr);
 117         r2[0] += dr[0];
 118         r2[1] += dr[1];
 119         r2[2] += dr[2];
 120     }
 121     // In the case of a single-atom site, r1 and r2 are now correct if local or [0,0,0] if not local.
 122
 123
 124     // Master rank update call-back. This needs to be moved to a discrete place in the
 125     // time step to avoid extraneous barriers. The code would be prettier with "futures"...
 126     if ((cr.dd == nullptr) || MASTER(&cr))
 127     {
 128         restraint_->update(RVec(r1),
 129                            r2,
 130                            t);
 131     }
 132     // All ranks wait for the update to finish.
 133     // tMPI ranks are depending on structures that may have just been updated.
 134     if (DOMAINDECOMP(&cr))
 135     {
 136         // Note: this assumes that all ranks are hitting this line, which is not generally true.
 137         // I need to find the right subcommunicator. What I really want is a _scoped_ communicator...
 138         gmx_barrier(&cr);
 139     }
 140
 141     // Apply restraint on all thread ranks only after any updates have been made.
 142     auto result = restraint_->evaluate(RVec(r1),
 143                                        r2,
 144                                        t);
 145
 146     // This can easily be generalized for pair restraints that apply to selections instead of
 147     // individual indices, or to restraints that aren't pair restraints.
 148     const int  site1  = static_cast<int>(sites_.front().index());
 149     const int* aLocal = &site1;
 150     // Set forces using index `site1` if no domain decomposition, otherwise set with local index if available.
 151     auto      &force = forceProviderOutput->forceWithVirial_.force_;
 152     if ((cr.dd == nullptr) || (aLocal = cr.dd->ga2la->findHome(site1)))
 153     {
 154         force[static_cast<size_t>(*aLocal)] += result.force;
 155     }
 156
 157     // Note: Currently calculateForces is called once per restraint and each restraint
 158     // applies to a pair of atoms. Future optimizations may consolidate multiple restraints
 159     // with possibly duplicated sites, in which case we may prefer to iterate over non-frozen
 160     // sites to apply forces without explicitly expressing pairwise symmetry as in the
 161     // following logic.
 162     const int  site2  = static_cast<int>(sites_.back().index());
 163     const int* bLocal = &site2;
 164     if ((cr.dd == nullptr) || (bLocal = cr.dd->ga2la->findHome(site2)))
 165     {
 166         force[static_cast<size_t>(*bLocal)] -= result.force;
 167     }
 168 }
 169
 170 RestraintMDModuleImpl::~RestraintMDModuleImpl() = default;
 171
 172 RestraintMDModuleImpl::RestraintMDModuleImpl(std::shared_ptr<IRestraintPotential> restraint,
 173                                              const std::vector<int>              &sites) :
 174     forceProvider_(std::make_unique<RestraintForceProvider>(restraint,
 175                                                             sites))
 176 {
 177     GMX_ASSERT(forceProvider_, "Class invariant implies non-null ForceProvider.");
 178 }
 179
 180 IMdpOptionProvider* RestraintMDModuleImpl::mdpOptionProvider()
 181 {
 182     return nullptr;
 183 }
 184
 185 IMDOutputProvider* RestraintMDModuleImpl::outputProvider()
 186 {
 187     return nullptr;
 188 }
 189
 190 void RestraintMDModuleImpl::initForceProviders(ForceProviders* forceProviders)
 191 {
 192     GMX_ASSERT(forceProvider_, "Class invariant implies non-null ForceProvider member.");
 193     GMX_ASSERT(forceProviders, "Provided ForceProviders* assumed to be non-null.");
 194     forceProviders->addForceProvider(forceProvider_.get());
 195 }
 196
 197
 198 // Needs to be defined after implementation type is complete in order to have unique_ptr member.
 199 RestraintMDModule::~RestraintMDModule() = default;
 200
 201
 202 IMdpOptionProvider* RestraintMDModule::mdpOptionProvider()
 203 {
 204     return nullptr;
 205 }
 206
 207 IMDOutputProvider* RestraintMDModule::outputProvider()
 208 {
 209     return nullptr;
 210 }
 211
 212 void RestraintMDModule::initForceProviders(ForceProviders* forceProviders)
 213 {
 214     GMX_ASSERT(impl_, "Class invariant implies non-null implementation member.");
 215     impl_->initForceProviders(forceProviders);
 216 }
 217
 218 std::unique_ptr<RestraintMDModule>
 219 RestraintMDModule::create(std::shared_ptr<IRestraintPotential>  restraint,
 220                           const std::vector<int>               &sites)
 221 {
 222     auto implementation = std::make_unique<RestraintMDModuleImpl>(std::move(restraint),
 223                                                                   sites);
 224     auto newModule = std::make_unique<RestraintMDModule>(std::move(implementation));
 225     return newModule;
 226 }
 227
 228 // private constructor to implement static create() method.
 229 RestraintMDModule::RestraintMDModule(std::unique_ptr<RestraintMDModuleImpl> restraint) :
 230     impl_ {std::move(restraint)}
 231 {
 232 }
 233
 234 } // end namespace gmx