2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2019,2020,2021, 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.
37 * \brief Implements SETTLE using GPU
39 * This file contains implementation for the data management of GPU version of SETTLE constraints algorithm.
41 * \author Artem Zhmurov <zhmurov@gmail.com>
43 * \ingroup module_mdlib
47 #include "settle_gpu.h"
56 #include "gromacs/gpu_utils/devicebuffer.h"
57 #include "gromacs/gpu_utils/gputraits.h"
58 #include "gromacs/math/functions.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/mdlib/settle_gpu_internal.h"
61 #include "gromacs/pbcutil/pbc.h"
66 void SettleGpu::apply(const DeviceBuffer<Float3>& d_x,
67 DeviceBuffer<Float3> d_xp,
68 const bool updateVelocities,
69 DeviceBuffer<Float3> d_v,
71 const bool computeVirial,
73 const PbcAiuc& pbcAiuc)
76 // Early exit if no settles
84 // Fill with zeros so the values can be reduced to it
85 // Only 6 values are needed because virial is symmetrical
86 clearDeviceBufferAsync(&d_virialScaled_, 0, 6, deviceStream_);
89 launchSettleGpuKernel(numSettles_,
105 copyFromDeviceBuffer(
106 h_virialScaled_.data(), &d_virialScaled_, 0, 6, deviceStream_, GpuApiCallBehavior::Sync, nullptr);
108 // Mapping [XX, XY, XZ, YY, YZ, ZZ] internal format to a tensor object
109 virialScaled[XX][XX] += h_virialScaled_[0];
110 virialScaled[XX][YY] += h_virialScaled_[1];
111 virialScaled[XX][ZZ] += h_virialScaled_[2];
113 virialScaled[YY][XX] += h_virialScaled_[1];
114 virialScaled[YY][YY] += h_virialScaled_[3];
115 virialScaled[YY][ZZ] += h_virialScaled_[4];
117 virialScaled[ZZ][XX] += h_virialScaled_[2];
118 virialScaled[ZZ][YY] += h_virialScaled_[4];
119 virialScaled[ZZ][ZZ] += h_virialScaled_[5];
123 SettleGpu::SettleGpu(const gmx_mtop_t& mtop, const DeviceContext& deviceContext, const DeviceStream& deviceStream) :
124 deviceContext_(deviceContext), deviceStream_(deviceStream)
126 static_assert(sizeof(real) == sizeof(float),
127 "Real numbers should be in single precision in GPU code.");
129 // This is to prevent the assertion failure for the systems without water
130 int totalSettles = 0;
131 for (unsigned mt = 0; mt < mtop.moltype.size(); mt++)
133 const int nral1 = 1 + NRAL(F_SETTLE);
134 InteractionList interactionList = mtop.moltype[mt].ilist[F_SETTLE];
135 std::vector<int> iatoms = interactionList.iatoms;
136 totalSettles += iatoms.size() / nral1;
138 if (totalSettles == 0)
143 // TODO This should be lifted to a separate subroutine that gets the values of Oxygen and
144 // Hydrogen masses, checks if they are consistent across the topology and if there is no more
145 // than two values for each mass if the free energy perturbation is enabled. In later case,
146 // masses may need to be updated on a regular basis (i.e. in set(...) method).
147 // TODO Do the checks for FEP
151 for (unsigned mt = 0; mt < mtop.moltype.size(); mt++)
153 const int nral1 = 1 + NRAL(F_SETTLE);
154 InteractionList interactionList = mtop.moltype[mt].ilist[F_SETTLE];
155 std::vector<int> iatoms = interactionList.iatoms;
156 for (unsigned i = 0; i < iatoms.size() / nral1; i++)
158 WaterMolecule settler;
159 settler.ow1 = iatoms[i * nral1 + 1]; // Oxygen index
160 settler.hw2 = iatoms[i * nral1 + 2]; // First hydrogen index
161 settler.hw3 = iatoms[i * nral1 + 3]; // Second hydrogen index
162 t_atom ow1 = mtop.moltype[mt].atoms.atom[settler.ow1];
163 t_atom hw2 = mtop.moltype[mt].atoms.atom[settler.hw2];
164 t_atom hw3 = mtop.moltype[mt].atoms.atom[settler.hw3];
174 GMX_RELEASE_ASSERT(mO == ow1.m,
175 "Topology has different values for oxygen mass. Should be identical "
176 "in order to use SETTLE.");
177 GMX_RELEASE_ASSERT(hw2.m == hw3.m && hw2.m == mH,
178 "Topology has different values for hydrogen mass. Should be "
179 "identical in order to use SETTLE.");
183 GMX_RELEASE_ASSERT(mO > 0, "Could not find oxygen mass in the topology. Needed in SETTLE.");
184 GMX_RELEASE_ASSERT(mH > 0, "Could not find hydrogen mass in the topology. Needed in SETTLE.");
186 // TODO Very similar to SETTLE initialization on CPU. Should be lifted to a separate method
187 // (one that gets dOH and dHH values and checks them for consistency)
188 int settle_type = -1;
189 for (unsigned mt = 0; mt < mtop.moltype.size(); mt++)
191 const int nral1 = 1 + NRAL(F_SETTLE);
192 InteractionList interactionList = mtop.moltype[mt].ilist[F_SETTLE];
193 for (int i = 0; i < interactionList.size(); i += nral1)
195 if (settle_type == -1)
197 settle_type = interactionList.iatoms[i];
199 else if (interactionList.iatoms[i] != settle_type)
202 "The [molecules] section of your topology specifies more than one block "
204 "a [moleculetype] with a [settles] block. Only one such is allowed.\n"
205 "If you are trying to partition your solvent into different *groups*\n"
206 "(e.g. for freezing, T-coupling, etc.), you are using the wrong "
208 "files specify groups. Otherwise, you may wish to change the least-used\n"
209 "block of molecules with SETTLE constraints into 3 normal constraints.");
214 GMX_RELEASE_ASSERT(settle_type >= 0, "settle_init called without settles");
216 real dOH = mtop.ffparams.iparams[settle_type].settle.doh;
217 real dHH = mtop.ffparams.iparams[settle_type].settle.dhh;
219 settleParameters_ = settleParameters(mO, mH, 1.0 / mO, 1.0 / mH, dOH, dHH);
221 allocateDeviceBuffer(&d_virialScaled_, 6, deviceContext_);
222 h_virialScaled_.resize(6);
225 SettleGpu::~SettleGpu()
227 // Early exit if there is no settles
228 if (numSettles_ == 0)
232 freeDeviceBuffer(&d_virialScaled_);
233 if (numAtomIdsAlloc_ > 0)
235 freeDeviceBuffer(&d_atomIds_);
239 void SettleGpu::set(const InteractionDefinitions& idef)
241 const int nral1 = 1 + NRAL(F_SETTLE);
242 const InteractionList& il_settle = idef.il[F_SETTLE];
243 ArrayRef<const int> iatoms = il_settle.iatoms;
244 numSettles_ = il_settle.size() / nral1;
246 reallocateDeviceBuffer(&d_atomIds_, numSettles_, &numAtomIds_, &numAtomIdsAlloc_, deviceContext_);
247 h_atomIds_.resize(numSettles_);
248 for (int i = 0; i < numSettles_; i++)
250 WaterMolecule settler;
251 settler.ow1 = iatoms[i * nral1 + 1]; // Oxygen index
252 settler.hw2 = iatoms[i * nral1 + 2]; // First hydrogen index
253 settler.hw3 = iatoms[i * nral1 + 3]; // Second hydrogen index
254 h_atomIds_[i] = settler;
257 &d_atomIds_, h_atomIds_.data(), 0, numSettles_, deviceStream_, GpuApiCallBehavior::Sync, nullptr);