src/gromacs/restraint/restraintmdmodule.cpp

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