1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2014,2015,2017,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 #include "gmxpre.h"
36
37 #include "gromacs/mdlib/shake.h"
38
39 #include <assert.h>
40
41 #include <cmath>
42
43 #include <algorithm>
44 #include <vector>
45
46 #include <gtest/gtest.h>
47
48 #include "gromacs/utility/arrayref.h"
49
50 #include "testutils/refdata.h"
51 #include "testutils/testasserts.h"
52
53 namespace gmx
54 {
55 namespace
56 {
57
58 /*! \brief Stride of the vector of int used to describe each SHAKE
59 * constraint
60 *
61 * Like other such code, SHAKE is hard-wired to use t_ilist.iatoms as
62 * a flat vector of tuples of general data. Here, they are triples
63 * containing the index of the constraint type, and then the indices
64 * of the two atoms involved. So for each constraint, we must stride
65 * this vector by three to get access to its information. */
66 const int constraintStride = 3;
67
68 /*! \brief Compute the displacements between pairs of constrained
69 * atoms described in the iatom "topology". */
computeDisplacements(ArrayRef<const int> iatom,const std::vector<RVec> & positions)70 std::vector<RVec> computeDisplacements(ArrayRef<const int> iatom, const std::vector<RVec>& positions)
71 {
72 assert(0 == iatom.size() % constraintStride);
73 int numConstraints = iatom.size() / constraintStride;
74 std::vector<RVec> displacements;
75
76 for (int ll = 0; ll != numConstraints; ++ll)
77 {
78 int atom_i = iatom[ll * constraintStride + 1];
79 int atom_j = iatom[ll * constraintStride + 2];
80
81 displacements.push_back(positions[atom_i] - positions[atom_j]);
82 }
83
84 return displacements;
85 }
86
87 /*! \brief Compute half of the reduced mass of each pair of constrained
88 * atoms in the iatom "topology".
89 *
90 * The reduced mass is m = 1/(1/m_i + 1/m_j)) */
computeHalfOfReducedMasses(const std::vector<int> & iatom,const std::vector<real> & inverseMasses)91 std::vector<real> computeHalfOfReducedMasses(const std::vector<int>& iatom,
92 const std::vector<real>& inverseMasses)
93 {
94 int numConstraints = iatom.size() / constraintStride;
95 std::vector<real> halfOfReducedMasses;
96
97 for (int ll = 0; ll != numConstraints; ++ll)
98 {
99 int atom_i = iatom[ll * constraintStride + 1];
100 int atom_j = iatom[ll * constraintStride + 2];
101
102 halfOfReducedMasses.push_back(0.5 / (inverseMasses[atom_i] + inverseMasses[atom_j]));
103 }
104
105 return halfOfReducedMasses;
106 }
107
108 /*! \brief Compute the distances corresponding to the vector of displacements components */
computeDistancesSquared(ArrayRef<const RVec> displacements)109 std::vector<real> computeDistancesSquared(ArrayRef<const RVec> displacements)
110 {
111 int numDistancesSquared = displacements.size();
112 std::vector<real> distanceSquared;
113
114 for (int i = 0; i != numDistancesSquared; ++i)
115 {
116 distanceSquared.push_back(0.0);
117 for (int d = 0; d != DIM; ++d)
118 {
119 real displacement = displacements[i][d];
120 distanceSquared.back() += displacement * displacement;
121 }
122 }
123
124 return distanceSquared;
125 }
126
127 /*! \brief Test fixture for testing SHAKE */
128 class ShakeTest : public ::testing::Test
129 {
130 public:
131 /*! \brief Set up data for test cases to use when constructing
132 their input */
SetUp()133 void SetUp() override
134 {
135 inverseMassesDatabase_.push_back(2.0);
136 inverseMassesDatabase_.push_back(3.0);
137 inverseMassesDatabase_.push_back(4.0);
138 inverseMassesDatabase_.push_back(1.0);
139
140 positionsDatabase_.emplace_back(2.5, -3.1, 15.7);
141
142 positionsDatabase_.emplace_back(0.51, -3.02, 15.55);
143
144 positionsDatabase_.emplace_back(-0.5, -3.0, 15.2);
145
146 positionsDatabase_.emplace_back(-1.51, -2.95, 15.05);
147 }
148
149 //! Run the test
runTest(size_t gmx_unused numAtoms,size_t numConstraints,const std::vector<int> & iatom,const std::vector<real> & constrainedDistances,const std::vector<real> & inverseMasses,const std::vector<RVec> & positions)150 static void runTest(size_t gmx_unused numAtoms,
151 size_t numConstraints,
152 const std::vector<int>& iatom,
153 const std::vector<real>& constrainedDistances,
154 const std::vector<real>& inverseMasses,
155 const std::vector<RVec>& positions)
156 {
157 // Check the test input is consistent
158 assert(numConstraints * constraintStride == iatom.size());
159 assert(numConstraints == constrainedDistances.size());
160 assert(numAtoms == inverseMasses.size());
161 assert(numAtoms == positions.size());
162 for (size_t i = 0; i != numConstraints; ++i)
163 {
164 for (size_t j = 1; j < 3; j++)
165 {
166 // Check that the topology refers to atoms that have masses and positions
167 assert(iatom[i * constraintStride + j] >= 0);
168 assert(iatom[i * constraintStride + j] < static_cast<int>(numAtoms));
169 }
170 }
171 std::vector<real> distanceSquaredTolerances;
172 std::vector<real> lagrangianValues;
173 std::vector<real> constrainedDistancesSquared;
174
175 real coordMax = 0;
176 for (size_t i = 0; i != numConstraints; ++i)
177 {
178 constrainedDistancesSquared.push_back(constrainedDistances[i] * constrainedDistances[i]);
179 distanceSquaredTolerances.push_back(
180 0.5 / (constrainedDistancesSquared.back() * ShakeTest::tolerance_));
181 lagrangianValues.push_back(0.0);
182
183 for (size_t j = 1; j < constraintStride; j++)
184 {
185 for (int d = 0; d < DIM; d++)
186 {
187 coordMax = std::max(coordMax, std::abs(positions[iatom[i * constraintStride + j]][d]));
188 }
189 }
190 }
191 std::vector<real> halfOfReducedMasses = computeHalfOfReducedMasses(iatom, inverseMasses);
192 std::vector<RVec> initialDisplacements = computeDisplacements(iatom, positions);
193
194 std::vector<RVec> finalPositions = positions;
195 int numIterations = 0;
196 int numErrors = 0;
197
198 cshake(iatom.data(), numConstraints, &numIterations, ShakeTest::maxNumIterations_,
199 constrainedDistancesSquared, finalPositions, nullptr, initialDisplacements,
200 halfOfReducedMasses, omega_, inverseMasses.data(), distanceSquaredTolerances,
201 lagrangianValues, &numErrors);
202
203 std::vector<RVec> finalDisplacements = computeDisplacements(iatom, finalPositions);
204 std::vector<real> finalDistancesSquared = computeDistancesSquared(finalDisplacements);
205 assert(numConstraints == finalDistancesSquared.size());
206
207 EXPECT_EQ(0, numErrors);
208 EXPECT_GT(numIterations, 1);
209 EXPECT_LT(numIterations, ShakeTest::maxNumIterations_);
210 // TODO wrap this in a Google Mock matcher if there's
211 // other tests like it some time?
212 for (size_t i = 0; i != numConstraints; ++i)
213 {
214 // We need to allow for the requested tolerance plus rounding
215 // errors due to the absolute size of the coordinate values
216 test::FloatingPointTolerance constraintTolerance =
217 test::absoluteTolerance(std::sqrt(constrainedDistancesSquared[i]) * ShakeTest::tolerance_
218 + coordMax * GMX_REAL_EPS);
219 // Assert that the constrained distances are within the required tolerance
220 EXPECT_FLOAT_EQ_TOL(std::sqrt(constrainedDistancesSquared[i]),
221 std::sqrt(finalDistancesSquared[i]), constraintTolerance);
222 }
223 }
224
225 //! Tolerance for SHAKE conversion (ie. shake-tol .mdp setting)
226 static const real tolerance_;
227 //! Maximum number of iterations permitted in these tests
228 static const int maxNumIterations_;
229 //! SHAKE over-relaxation (SOR) factor
230 static const real omega_;
231 //! Database of inverse masses of atoms in the topology
232 std::vector<real> inverseMassesDatabase_;
233 //! Database of atom positions (three reals per atom)
234 std::vector<RVec> positionsDatabase_;
235 };
236
237 const real ShakeTest::tolerance_ = 1e-5;
238 const int ShakeTest::maxNumIterations_ = 30;
239 const real ShakeTest::omega_ = 1.0;
240
TEST_F(ShakeTest,ConstrainsOneBond)241 TEST_F(ShakeTest, ConstrainsOneBond)
242 {
243 int numAtoms = 2;
244 int numConstraints = 1;
245
246 std::vector<int> iatom;
247 iatom.push_back(-1); // unused
248 iatom.push_back(0); // i atom index
249 iatom.push_back(1); // j atom index
250
251 std::vector<real> constrainedDistances;
252 constrainedDistances.push_back(2.0);
253
254 std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
255 inverseMassesDatabase_.begin() + numAtoms);
256 std::vector<RVec> positions(positionsDatabase_.begin(), positionsDatabase_.begin() + numAtoms);
257
258 runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
259 }
260
TEST_F(ShakeTest,ConstrainsTwoDisjointBonds)261 TEST_F(ShakeTest, ConstrainsTwoDisjointBonds)
262 {
263 int numAtoms = 4;
264 int numConstraints = 2;
265
266 std::vector<int> iatom;
267 iatom.push_back(-1); // unused
268 iatom.push_back(0); // i atom index
269 iatom.push_back(1); // j atom index
270
271 iatom.push_back(-1); // unused
272 iatom.push_back(2); // i atom index
273 iatom.push_back(3); // j atom index
274
275 std::vector<real> constrainedDistances;
276 constrainedDistances.push_back(2.0);
277 constrainedDistances.push_back(1.0);
278
279 std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
280 inverseMassesDatabase_.begin() + numAtoms);
281 std::vector<RVec> positions(positionsDatabase_.begin(), positionsDatabase_.begin() + numAtoms);
282
283 runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
284 }
285
TEST_F(ShakeTest,ConstrainsTwoBondsWithACommonAtom)286 TEST_F(ShakeTest, ConstrainsTwoBondsWithACommonAtom)
287 {
288 int numAtoms = 3;
289 int numConstraints = 2;
290
291 std::vector<int> iatom;
292 iatom.push_back(-1); // unused
293 iatom.push_back(0); // i atom index
294 iatom.push_back(1); // j atom index
295
296 iatom.push_back(-1); // unused
297 iatom.push_back(1); // i atom index
298 iatom.push_back(2); // j atom index
299
300 std::vector<real> constrainedDistances;
301 constrainedDistances.push_back(2.0);
302 constrainedDistances.push_back(1.0);
303
304 std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
305 inverseMassesDatabase_.begin() + numAtoms);
306 std::vector<RVec> positions(positionsDatabase_.begin(), positionsDatabase_.begin() + numAtoms);
307
308 runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
309 }
310
TEST_F(ShakeTest,ConstrainsThreeBondsWithCommonAtoms)311 TEST_F(ShakeTest, ConstrainsThreeBondsWithCommonAtoms)
312 {
313 int numAtoms = 4;
314 int numConstraints = 3;
315
316 std::vector<int> iatom;
317 iatom.push_back(-1); // unused
318 iatom.push_back(0); // i atom index
319 iatom.push_back(1); // j atom index
320
321 iatom.push_back(-1); // unused
322 iatom.push_back(1); // i atom index
323 iatom.push_back(2); // j atom index
324
325 iatom.push_back(-1); // unused
326 iatom.push_back(2); // i atom index
327 iatom.push_back(3); // j atom index
328
329 std::vector<real> constrainedDistances;
330 constrainedDistances.push_back(2.0);
331 constrainedDistances.push_back(1.0);
332 constrainedDistances.push_back(1.0);
333
334 std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
335 inverseMassesDatabase_.begin() + numAtoms);
336 std::vector<RVec> positions(positionsDatabase_.begin(), positionsDatabase_.begin() + numAtoms);
337
338 runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
339 }
340
341 } // namespace
342 } // namespace gmx
343