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