2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2017,2018, 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 #include "gromacs/mdlib/shake.h"
46 #include <gtest/gtest.h>
48 #include "testutils/refdata.h"
49 #include "testutils/testasserts.h"
56 /*! \brief Stride of the vector of int used to describe each SHAKE
59 * Like other such code, SHAKE is hard-wired to use t_ilist.iatoms as
60 * a flat vector of tuples of general data. Here, they are triples
61 * containing the index of the constraint type, and then the indices
62 * of the two atoms involved. So for each constraint, we must stride
63 * this vector by three to get access to its information. */
64 const int constraintStride = 3;
66 /*! \brief Compute the displacements between pairs of constrained
67 * atoms described in the iatom "topology". */
69 computeDisplacements(const std::vector<int> &iatom,
70 const std::vector<real> &positions)
72 assert(0 == iatom.size() % constraintStride);
73 int numConstraints = iatom.size() / constraintStride;
74 std::vector<real> displacements;
76 for (int ll = 0; ll != numConstraints; ++ll)
78 int atom_i = iatom[ll*constraintStride + 1];
79 int atom_j = iatom[ll*constraintStride + 2];
81 for (int d = 0; d != DIM; d++)
83 displacements.push_back(positions[atom_i*DIM + d] - positions[atom_j*DIM + d]);
90 /*! \brief Compute half of the reduced mass of each pair of constrained
91 * atoms in the iatom "topology".
93 * The reduced mass is m = 1/(1/m_i + 1/m_j)) */
95 computeHalfOfReducedMasses(const std::vector<int> &iatom,
96 const std::vector<real> &inverseMasses)
98 int numConstraints = iatom.size() / constraintStride;
99 std::vector<real> halfOfReducedMasses;
101 for (int ll = 0; ll != numConstraints; ++ll)
103 int atom_i = iatom[ll*constraintStride + 1];
104 int atom_j = iatom[ll*constraintStride + 2];
106 halfOfReducedMasses.push_back(0.5 / (inverseMasses[atom_i] + inverseMasses[atom_j]));
109 return halfOfReducedMasses;
112 /*! \brief Compute the distances corresponding to the vector of displacements components */
114 computeDistancesSquared(const std::vector<real> &displacements)
116 assert(0 == displacements.size() % DIM);
117 int numDistancesSquared = displacements.size() / DIM;
118 std::vector<real> distanceSquared;
120 for (int i = 0; i != numDistancesSquared; ++i)
122 distanceSquared.push_back(0.0);
123 for (int d = 0; d != DIM; ++d)
125 real displacement = displacements[i*DIM + d];
126 distanceSquared.back() += displacement * displacement;
130 return distanceSquared;
133 /*! \brief Test fixture for testing SHAKE */
134 class ShakeTest : public ::testing::Test
137 /*! \brief Set up data for test cases to use when constructing
139 void SetUp() override
141 inverseMassesDatabase_.push_back(2.0);
142 inverseMassesDatabase_.push_back(3.0);
143 inverseMassesDatabase_.push_back(4.0);
144 inverseMassesDatabase_.push_back(1.0);
146 positionsDatabase_.push_back(2.5);
147 positionsDatabase_.push_back(-3.1);
148 positionsDatabase_.push_back(15.7);
150 positionsDatabase_.push_back(0.51);
151 positionsDatabase_.push_back(-3.02);
152 positionsDatabase_.push_back(15.55);
154 positionsDatabase_.push_back(-0.5);
155 positionsDatabase_.push_back(-3.0);
156 positionsDatabase_.push_back(15.2);
158 positionsDatabase_.push_back(-1.51);
159 positionsDatabase_.push_back(-2.95);
160 positionsDatabase_.push_back(15.05);
164 void runTest(size_t gmx_unused numAtoms,
165 size_t numConstraints,
166 const std::vector<int> &iatom,
167 const std::vector<real> &constrainedDistances,
168 const std::vector<real> &inverseMasses,
169 const std::vector<real> &positions)
171 // Check the test input is consistent
172 assert(numConstraints * constraintStride == iatom.size());
173 assert(numConstraints == constrainedDistances.size());
174 assert(numAtoms == inverseMasses.size());
175 assert(numAtoms * DIM == positions.size());
176 for (size_t i = 0; i != numConstraints; ++i)
178 for (size_t j = 1; j < 3; j++)
180 // Check that the topology refers to atoms that have masses and positions
181 assert(iatom[i*constraintStride + j] >= 0);
182 assert(iatom[i*constraintStride + j] < static_cast<int>(numAtoms));
185 std::vector<real> distanceSquaredTolerances;
186 std::vector<real> lagrangianValues;
187 std::vector<real> constrainedDistancesSquared;
190 for (size_t i = 0; i != numConstraints; ++i)
192 constrainedDistancesSquared.push_back(constrainedDistances[i] * constrainedDistances[i]);
193 distanceSquaredTolerances.push_back(0.5 / (constrainedDistancesSquared.back() * ShakeTest::tolerance_));
194 lagrangianValues.push_back(0.0);
196 for (size_t j = 1; j < constraintStride; j++)
198 for (int d = 0; d < DIM; d++)
200 coordMax = std::max(coordMax, std::abs(positions[iatom[i*constraintStride + j]*DIM + d]));
204 std::vector<real> halfOfReducedMasses = computeHalfOfReducedMasses(iatom, inverseMasses);
205 std::vector<real> initialDisplacements = computeDisplacements(iatom, positions);
207 std::vector<real> finalPositions = positions;
208 int numIterations = 0;
211 cshake(iatom.data(), numConstraints, &numIterations,
212 ShakeTest::maxNumIterations_, constrainedDistancesSquared.data(),
213 finalPositions.data(), initialDisplacements.data(),
214 halfOfReducedMasses.data(), omega_, inverseMasses.data(),
215 distanceSquaredTolerances.data(),
216 lagrangianValues.data(),
219 std::vector<real> finalDisplacements = computeDisplacements(iatom, finalPositions);
220 std::vector<real> finalDistancesSquared = computeDistancesSquared(finalDisplacements);
221 assert(numConstraints == finalDistancesSquared.size());
223 EXPECT_EQ(0, numErrors);
224 EXPECT_GT(numIterations, 1);
225 EXPECT_LT(numIterations, ShakeTest::maxNumIterations_);
226 // TODO wrap this in a Google Mock matcher if there's
227 // other tests like it some time?
228 for (size_t i = 0; i != numConstraints; ++i)
230 // We need to allow for the requested tolerance plus rounding
231 // errors due to the absolute size of the coordinate values
232 test::FloatingPointTolerance constraintTolerance =
233 test::absoluteTolerance(std::sqrt(constrainedDistancesSquared[i])*ShakeTest::tolerance_ + coordMax*GMX_REAL_EPS);
234 // Assert that the constrained distances are within the required tolerance
235 EXPECT_FLOAT_EQ_TOL(std::sqrt(constrainedDistancesSquared[i]),
236 std::sqrt(finalDistancesSquared[i]),
237 constraintTolerance);
241 //! Tolerance for SHAKE conversion (ie. shake-tol .mdp setting)
242 static const real tolerance_;
243 //! Maximum number of iterations permitted in these tests
244 static const int maxNumIterations_;
245 //! SHAKE over-relaxation (SOR) factor
246 static const real omega_;
247 //! Database of inverse masses of atoms in the topology
248 std::vector<real> inverseMassesDatabase_;
249 //! Database of atom positions (three reals per atom)
250 std::vector<real> positionsDatabase_;
253 const real ShakeTest::tolerance_ = 1e-5;
254 const int ShakeTest::maxNumIterations_ = 30;
255 const real ShakeTest::omega_ = 1.0;
257 TEST_F(ShakeTest, ConstrainsOneBond)
260 int numConstraints = 1;
262 std::vector<int> iatom;
263 iatom.push_back(-1); // unused
264 iatom.push_back(0); // i atom index
265 iatom.push_back(1); // j atom index
267 std::vector<real> constrainedDistances;
268 constrainedDistances.push_back(2.0);
270 std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
271 inverseMassesDatabase_.begin() + numAtoms);
272 std::vector<real> positions(positionsDatabase_.begin(),
273 positionsDatabase_.begin() + numAtoms * DIM);
275 runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
278 TEST_F(ShakeTest, ConstrainsTwoDisjointBonds)
281 int numConstraints = 2;
283 std::vector<int> iatom;
284 iatom.push_back(-1); // unused
285 iatom.push_back(0); // i atom index
286 iatom.push_back(1); // j atom index
288 iatom.push_back(-1); // unused
289 iatom.push_back(2); // i atom index
290 iatom.push_back(3); // j atom index
292 std::vector<real> constrainedDistances;
293 constrainedDistances.push_back(2.0);
294 constrainedDistances.push_back(1.0);
296 std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
297 inverseMassesDatabase_.begin() + numAtoms);
298 std::vector<real> positions(positionsDatabase_.begin(),
299 positionsDatabase_.begin() + numAtoms * DIM);
301 runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
304 TEST_F(ShakeTest, ConstrainsTwoBondsWithACommonAtom)
307 int numConstraints = 2;
309 std::vector<int> iatom;
310 iatom.push_back(-1); // unused
311 iatom.push_back(0); // i atom index
312 iatom.push_back(1); // j atom index
314 iatom.push_back(-1); // unused
315 iatom.push_back(1); // i atom index
316 iatom.push_back(2); // j atom index
318 std::vector<real> constrainedDistances;
319 constrainedDistances.push_back(2.0);
320 constrainedDistances.push_back(1.0);
322 std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
323 inverseMassesDatabase_.begin() + numAtoms);
324 std::vector<real> positions(positionsDatabase_.begin(),
325 positionsDatabase_.begin() + numAtoms * DIM);
327 runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
330 TEST_F(ShakeTest, ConstrainsThreeBondsWithCommonAtoms)
333 int numConstraints = 3;
335 std::vector<int> iatom;
336 iatom.push_back(-1); // unused
337 iatom.push_back(0); // i atom index
338 iatom.push_back(1); // j atom index
340 iatom.push_back(-1); // unused
341 iatom.push_back(1); // i atom index
342 iatom.push_back(2); // j atom index
344 iatom.push_back(-1); // unused
345 iatom.push_back(2); // i atom index
346 iatom.push_back(3); // j atom index
348 std::vector<real> constrainedDistances;
349 constrainedDistances.push_back(2.0);
350 constrainedDistances.push_back(1.0);
351 constrainedDistances.push_back(1.0);
353 std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
354 inverseMassesDatabase_.begin() + numAtoms);
355 std::vector<real> positions(positionsDatabase_.begin(),
356 positionsDatabase_.begin() + numAtoms * DIM);
358 runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);