2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016, 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.
40 #include <gtest/gtest.h>
42 #include "gromacs/math/vec.h"
43 #include "gromacs/math/vectypes.h"
44 #include "gromacs/mdlib/constr.h"
45 #include "gromacs/mdtypes/mdatom.h"
46 #include "gromacs/pbcutil/pbc.h"
47 #include "gromacs/topology/idef.h"
48 #include "gromacs/topology/topology.h"
49 #include "gromacs/utility/smalloc.h"
50 #include "gromacs/utility/stringutil.h"
51 #include "gromacs/utility/unique_cptr.h"
53 #include "testutils/testasserts.h"
61 //! Database of 51 water atom input positions (DIM reals per atom, taken from spc216.gro) for use as test inputs.
62 const double g_positions[] = {
116 //! Simple cubic simulation box to use in tests
117 matrix g_box = {{real(1.86206), 0, 0}, {0, real(1.86206), 0}, {0, 0, real(1.86206)}};
119 //! Convenience typedef
120 typedef std::tuple<int, bool, bool, bool> SettleTestParameters;
122 /*! \brief Test fixture for testing SETTLE position updates
124 * \todo This also tests that if the calling code requires velocities
125 * and virial updates, that those outputs do change, but does not test
126 * that those changes are correct.
128 * \todo Only no-PBC and cubic-PBC are tested here, but the correct
129 * function of the SIMD version of set_pbx_auic in all cases should be
132 class SettleTest : public ::testing::TestWithParam<SettleTestParameters>
135 //! Updated water atom positions to constrain (DIM reals per atom)
136 std::vector<real> updatedPositions_;
137 //! Water atom velocities to constrain (DIM reals per atom)
138 std::vector<real> velocities_;
139 //! PBC option to test
141 //! PBC option to test
146 updatedPositions_(std::begin(g_positions), std::end(g_positions)),
147 velocities_(updatedPositions_.size(), 0)
149 set_pbc(&pbcNone_, epbcNONE, g_box);
150 set_pbc(&pbcXYZ_, epbcXYZ, g_box);
152 // Perturb the atom positions, to appear like an
153 // "update," and where there is definitely constraining
155 for (size_t i = 0; i != updatedPositions_.size(); ++i)
159 updatedPositions_[i] += 0.01;
163 updatedPositions_[i] -= 0.01;
167 updatedPositions_[i] += 0.02;
171 updatedPositions_[i] -= 0.02;
177 TEST_P(SettleTest, SatisfiesConstraints)
180 bool usePbc, useVelocities, calcVirial;
181 // Make some symbolic names for the parameter combination under
183 std::tie(numSettles, usePbc, useVelocities, calcVirial) = GetParam();
185 // Make a string that describes which parameter combination is
186 // being tested, to help make failing tests comprehensible.
187 std::string testDescription = formatString("while testing %d SETTLEs, %sPBC, %svelocities and %scalculating the virial",
189 usePbc ? "with " : "without ",
190 useVelocities ? "with " : "without ",
191 calcVirial ? "" : "not ");
193 const int settleType = 0;
194 const int atomsPerSettle = NRAL(F_SETTLE);
195 ASSERT_LE(numSettles, updatedPositions_.size() / (atomsPerSettle * DIM)) << "cannot test that many SETTLEs " << testDescription;
197 // Set up the topology. We still have to make some raw pointers,
198 // but they are put into scope guards for automatic cleanup.
201 const unique_cptr<gmx_mtop_t> mtopGuard(mtop);
204 snew(mtop->moltype, mtop->nmoltype);
205 const unique_cptr<gmx_moltype_t> moltypeGuard(mtop->moltype);
207 snew(mtop->molblock, mtop->nmolblock);
208 const unique_cptr<gmx_molblock_t> molblockGuard(mtop->molblock);
209 mtop->molblock[0].type = 0;
210 std::vector<int> iatoms;
211 for (int i = 0; i < numSettles; ++i)
213 iatoms.push_back(settleType);
214 iatoms.push_back(i*atomsPerSettle+0);
215 iatoms.push_back(i*atomsPerSettle+1);
216 iatoms.push_back(i*atomsPerSettle+2);
218 mtop->moltype[0].ilist[F_SETTLE].iatoms = iatoms.data();
219 mtop->moltype[0].ilist[F_SETTLE].nr = iatoms.size();
221 // Set up the SETTLE parameters.
222 mtop->ffparams.ntypes = 1;
223 snew(mtop->ffparams.iparams, mtop->ffparams.ntypes);
224 const unique_cptr<t_iparams> iparamsGuard(mtop->ffparams.iparams);
225 const real dOH = 0.09572;
226 const real dHH = 0.15139;
227 mtop->ffparams.iparams[settleType].settle.doh = dOH;
228 mtop->ffparams.iparams[settleType].settle.dhh = dHH;
230 // Set up the masses.
232 std::vector<real> mass, massReciprocal;
233 const real oxygenMass = 15.9994, hydrogenMass = 1.008;
234 for (int i = 0; i < numSettles; ++i)
236 mass.push_back(oxygenMass);
237 mass.push_back(hydrogenMass);
238 mass.push_back(hydrogenMass);
239 massReciprocal.push_back(1./oxygenMass);
240 massReciprocal.push_back(1./hydrogenMass);
241 massReciprocal.push_back(1./hydrogenMass);
243 mdatoms.massT = mass.data();
244 mdatoms.invmass = massReciprocal.data();
245 mdatoms.homenr = numSettles * atomsPerSettle;
247 // Finally make the settle data structures
248 gmx_settledata_t settled = settle_init(mtop);
249 settle_set_constraints(settled, &mtop->moltype[0].ilist[F_SETTLE], &mdatoms);
251 // Copy the original positions from the array of doubles to a vector of reals
252 std::vector<real> startingPositions(std::begin(g_positions), std::end(g_positions));
256 tensor virial = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
257 int numThreads = 1, threadIndex = 0;
258 const real reciprocalTimeStep = 1.0/0.002;
259 csettle(settled, numThreads, threadIndex,
260 usePbc ? &pbcXYZ_ : &pbcNone_,
261 startingPositions.data(), updatedPositions_.data(), reciprocalTimeStep,
262 useVelocities ? velocities_.data() : nullptr,
263 calcVirial, virial, &errorOccured);
264 settle_free(settled);
265 EXPECT_FALSE(errorOccured) << testDescription;
267 // The necessary tolerances for the test to pass were determined
268 // empirically. This isn't nice, but the required behaviour that
269 // SETTLE produces constrained coordinates consistent with
270 // sensible sampling needs to be tested at a much higher level.
271 FloatingPointTolerance tolerance =
272 relativeToleranceAsPrecisionDependentUlp(dOH*dOH, 80, 380);
274 // Verify the updated coordinates match the requirements
275 int positionIndex = 0, velocityIndex = 0;
276 for (int i = 0; i < numSettles; ++i)
279 updatedPositions_[positionIndex++], updatedPositions_[positionIndex++], updatedPositions_[positionIndex++]
282 updatedPositions_[positionIndex++], updatedPositions_[positionIndex++], updatedPositions_[positionIndex++]
285 updatedPositions_[positionIndex++], updatedPositions_[positionIndex++], updatedPositions_[positionIndex++]
288 EXPECT_REAL_EQ_TOL(dOH*dOH, distance2(positionO, positionH1), tolerance) << formatString("for water %d ", i) << testDescription;
289 EXPECT_REAL_EQ_TOL(dOH*dOH, distance2(positionO, positionH2), tolerance) << formatString("for water %d ", i) << testDescription;
290 EXPECT_REAL_EQ_TOL(dHH*dHH, distance2(positionH1, positionH2), tolerance) << formatString("for water %d ", i) << testDescription;
292 // This merely tests whether the velocities were
293 // updated from the starting values of zero (or not),
294 // but not whether the update was correct.
295 for (int j = 0; j < atomsPerSettle * DIM; ++j, ++velocityIndex)
297 EXPECT_TRUE(useVelocities == (0. != velocities_[velocityIndex])) << formatString("for water %d velocity coordinate %d ", i, j) << testDescription;
301 // This merely tests whether the viral was updated from
302 // the starting values of zero (or not), but not whether
303 // any update was correct.
304 for (int d = 0; d < DIM; ++d)
306 for (int dd = 0; dd < DIM; ++dd)
308 EXPECT_TRUE(calcVirial == (0. != virial[d][dd])) << formatString("for virial component[%d][%d] ", d, dd) << testDescription;
313 // Scan the full Cartesian product of numbers of SETTLE interactions
314 // (4 and 17 are chosen to test cases that do and do not match
315 // hardware SIMD widths), and whether or not we use PBC, velocities or
316 // calculate the virial contribution.
317 INSTANTIATE_TEST_CASE_P(WithParameters, SettleTest,
318 ::testing::Combine(::testing::Values(1, 4, 7),