?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Remove PImpl scaffolding from CUDA version of LINCS
[gromacs.git] / src / gromacs / mdlib / tests / constr.cpp
blob2f02a2ca4c9cf3885d81d3fe5c368e45fffde311
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2018,2019, 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 /*! \internal \file
36 * \brief SHAKE and LINCS tests.
37 *
38 * \todo Better tests for virial are needed.
39 * \todo Tests for bigger systems to test threads synchronization,
40 * reduction, etc. on the GPU.
41 * \todo Tests for algorithms for derivatives.
42 * \todo Free-energy perturbation tests
43 *
44 * \author Artem Zhmurov <zhmurov@gmail.com>
45 * \ingroup module_mdlib
46 */
47
48 #include "gmxpre.h"
49
50 #include "gromacs/mdlib/constr.h"
51
52 #include "config.h"
53
54 #include <assert.h>
55
56 #include <cmath>
57
58 #include <algorithm>
59 #include <unordered_map>
60 #include <vector>
61
62 #include <gtest/gtest.h>
63
64 #include "gromacs/fileio/gmxfio.h"
65 #include "gromacs/gmxlib/nrnb.h"
66 #include "gromacs/math/paddedvector.h"
67 #include "gromacs/math/vec.h"
68 #include "gromacs/math/vectypes.h"
69 #include "gromacs/mdlib/lincs.h"
70 #include "gromacs/mdlib/shake.h"
71 #include "gromacs/mdrunutility/multisim.h"
72 #include "gromacs/mdtypes/commrec.h"
73 #include "gromacs/mdtypes/inputrec.h"
74 #include "gromacs/mdtypes/mdatom.h"
75 #include "gromacs/pbcutil/pbc.h"
76 #include "gromacs/topology/block.h"
77 #include "gromacs/topology/idef.h"
78 #include "gromacs/topology/ifunc.h"
79 #include "gromacs/topology/topology.h"
80 #include "gromacs/utility/smalloc.h"
81 #include "gromacs/utility/stringutil.h"
82 #include "gromacs/utility/unique_cptr.h"
83
84 #include "testutils/refdata.h"
85 #include "testutils/testasserts.h"
86
87 #include "constr_impl.h"
88
89 namespace gmx
90 {
91 namespace test
92 {
93
94 ConstraintsTestData::ConstraintsTestData(const std::string &title,
95 int numAtoms, std::vector<real> masses,
96 std::vector<int> constraints, std::vector<real> constraintsR0,
97 bool computeVirial, tensor virialScaledRef,
98 bool compute_dHdLambda, float dHdLambdaRef,
99 real initialTime, real timestep,
100 const std::vector<RVec> &x, const std::vector<RVec> &xPrime, const std::vector<RVec> &v,
101 real shakeTolerance, gmx_bool shakeUseSOR,
102 int lincsNumIterations, int lincsExpansionOrder, real lincsWarnAngle)
103 {
104 // This is to trick Gerrit
105 {
106 {
107 title_ = title; // Human-friendly name of the system
108 numAtoms_ = numAtoms; // Number of atoms
109
110 // Masses of atoms
111 masses_ = masses;
112 invmass_.resize(numAtoms); // Vector of inverse masses
113
114 for (int i = 0; i < numAtoms; i++)
115 {
116 invmass_[i] = 1.0/masses.at(i);
117 }
118
119 // Saving constraints to check if they are satisfied after algorithm was applied
120 constraints_ = constraints; // Constraints indices (in type-i-j format)
121 constraintsR0_ = constraintsR0; // Equilibrium distances for each type of constraint
122
123 invdt_ = 1.0/timestep; // Inverse timestep
124
125 // Communication record
126 cr_.nnodes = 1;
127 cr_.dd = nullptr;
128
129 // Multisim data
130 ms_.sim = 0;
131 ms_.nsim = 1;
132
133 // Input record - data that usually comes from configuration file (.mdp)
134 ir_.efep = 0;
135 ir_.init_t = initialTime;
136 ir_.delta_t = timestep;
137 ir_.eI = 0;
138
139 // MD atoms data
140 md_.nMassPerturbed = 0;
141 md_.lambda = 0.0;
142 md_.invmass = invmass_.data();
143 md_.nr = numAtoms;
144 md_.homenr = numAtoms;
145
146 // Virial evaluation
147 computeVirial_ = computeVirial;
148 if (computeVirial)
149 {
150 for (int i = 0; i < DIM; i++)
151 {
152 for (int j = 0; j < DIM; j++)
153 {
154 virialScaled_[i][j] = 0;
155 virialScaledRef_[i][j] = virialScaledRef[i][j];
156 }
157 }
158 }
159
160
161 // Free energy evaluation
162 compute_dHdLambda_ = compute_dHdLambda;
163 dHdLambda_ = 0;
164 if (compute_dHdLambda_)
165 {
166 ir_.efep = efepYES;
167 dHdLambdaRef_ = dHdLambdaRef;
168 }
169 else
170 {
171 ir_.efep = efepNO;
172 dHdLambdaRef_ = 0;
173 }
174
175 // Constraints and their parameters (local topology)
176 for (int i = 0; i < F_NRE; i++)
177 {
178 idef_.il[i].nr = 0;
179 }
180 idef_.il[F_CONSTR].nr = constraints.size();
181
182 snew(idef_.il[F_CONSTR].iatoms, constraints.size());
183 int maxType = 0;
184 for (unsigned i = 0; i < constraints.size(); i++)
185 {
186 if (i % 3 == 0)
187 {
188 if (maxType < constraints.at(i))
189 {
190 maxType = constraints.at(i);
191 }
192 }
193 idef_.il[F_CONSTR].iatoms[i] = constraints.at(i);
194 }
195 snew(idef_.iparams, maxType + 1);
196 for (unsigned i = 0; i < constraints.size()/3; i++)
197 {
198 idef_.iparams[constraints.at(3*i)].constr.dA = constraintsR0.at(constraints.at(3*i));
199 idef_.iparams[constraints.at(3*i)].constr.dB = constraintsR0.at(constraints.at(3*i));
200 }
201
202 // Constraints and their parameters (global topology)
203 InteractionList interactionList;
204 interactionList.iatoms.resize(constraints.size());
205 for (unsigned i = 0; i < constraints.size(); i++)
206 {
207 interactionList.iatoms.at(i) = constraints.at(i);
208 }
209 InteractionList interactionListEmpty;
210 interactionListEmpty.iatoms.resize(0);
211
212 gmx_moltype_t molType;
213 molType.atoms.nr = numAtoms;
214 molType.ilist.at(F_CONSTR) = interactionList;
215 molType.ilist.at(F_CONSTRNC) = interactionListEmpty;
216 mtop_.moltype.push_back(molType);
217
218 gmx_molblock_t molBlock;
219 molBlock.type = 0;
220 molBlock.nmol = 1;
221 mtop_.molblock.push_back(molBlock);
222
223 mtop_.natoms = numAtoms;
224 mtop_.ffparams.iparams.resize(maxType + 1);
225 for (int i = 0; i <= maxType; i++)
226 {
227 mtop_.ffparams.iparams.at(i) = idef_.iparams[i];
228 }
229 mtop_.bIntermolecularInteractions = false;
230
231 // Coordinates and velocities
232 x_.resizeWithPadding(numAtoms);
233 xPrime_.resizeWithPadding(numAtoms);
234 xPrime0_.resizeWithPadding(numAtoms);
235 xPrime2_.resizeWithPadding(numAtoms);
236
237 v_.resizeWithPadding(numAtoms);
238 v0_.resizeWithPadding(numAtoms);
239
240 std::copy(x.begin(), x.end(), x_.begin());
241 std::copy(xPrime.begin(), xPrime.end(), xPrime_.begin());
242 std::copy(xPrime.begin(), xPrime.end(), xPrime0_.begin());
243 std::copy(xPrime.begin(), xPrime.end(), xPrime2_.begin());
244
245 std::copy(v.begin(), v.end(), v_.begin());
246 std::copy(v.begin(), v.end(), v0_.begin());
247
248 // SHAKE-specific parameters
249 ir_.shake_tol = shakeTolerance;
250 ir_.bShakeSOR = shakeUseSOR;
251
252 // LINCS-specific parameters
253 ir_.nLincsIter = lincsNumIterations;
254 ir_.nProjOrder = lincsExpansionOrder;
255 ir_.LincsWarnAngle = lincsWarnAngle;
256 }
257 }
258 }
259
260 /*! \brief
261 * Reset the data structure so it can be reused.
262 *
263 * Set the coordinates and velocities back to their values before
264 * constraining. The scaled virial tensor and dHdLambda are zeroed.
265 *
266 */
267 void ConstraintsTestData::reset()
268 {
269 xPrime_ = xPrime0_;
270 xPrime2_ = xPrime0_;
271 v_ = v0_;
272
273 if (computeVirial_)
274 {
275 for (int i = 0; i < DIM; i++)
276 {
277 for (int j = 0; j < DIM; j++)
278 {
279 virialScaled_[i][j] = 0;
280 }
281 }
282 }
283 dHdLambda_ = 0;
284 }
285
286 /*! \brief
287 * Cleaning up the memory.
288 */
289 ConstraintsTestData::~ConstraintsTestData()
290 {
291 sfree(idef_.il[F_CONSTR].iatoms);
292 sfree(idef_.iparams);
293 }
294
295 namespace
296 {
297
298 /*! \brief The two-dimensional parameter space for test.
299 *
300 * The test will run for all possible combinations of accessible
301 * values of the:
302 * 1. PBC setup ("PBCNONE" or "PBCXYZ")
303 * 2. The algorithm ("SHAKE", "LINCS" or "LINCS_GPU").
304 */
305 typedef std::tuple<std::string, std::string> ConstraintsTestParameters;
306
307 //! Names of all availible algorithms
308 std::vector<std::string> algorithmsNames;
309
310 //! Method that fills and returns algorithmNames to the test macros.
311 std::vector<std::string> getAlgorithmsNames()
312 {
313 algorithmsNames.emplace_back("SHAKE");
314 algorithmsNames.emplace_back("LINCS");
315 std::string errorMessage;
316 if (GMX_GPU == GMX_GPU_CUDA && canDetectGpus(&errorMessage))
317 {
318 algorithmsNames.emplace_back("LINCS_CUDA");
319 }
320 return algorithmsNames;
321 }
322
323 /*! \brief Test fixture for constraints.
324 *
325 * The fixture uses following test systems:
326 * 1. Two atoms, connected with one constraint (e.g. NH).
327 * 2. Three atoms, connected consequently with two constraints (e.g. CH2).
328 * 3. Three atoms, constrained to the fourth atom (e.g. CH3).
329 * 4. Four atoms, connected by two independent constraints.
330 * 5. Three atoms, connected by three constraints in a triangle
331 * (e.g. H2O with constrained H-O-H angle).
332 * 6. Four atoms, connected by three consequential constraints.
333 *
334 * For all systems, the final lengths of the constraints are tested against the
335 * reference values, the direction of each constraint is checked.
336 * Test also verifies that the center of mass has not been
337 * shifted by the constraints and that its velocity has not changed.
338 * For some systems, the value for scaled virial tensor is checked against
339 * pre-computed data.
340 */
341 class ConstraintsTest : public ::testing::TestWithParam<ConstraintsTestParameters>
342 {
343 public:
344 //! PBC setups
345 std::unordered_map <std::string, t_pbc> pbcs_;
346 //! Algorithms (SHAKE and LINCS)
347 std::unordered_map <std::string, void(*)(ConstraintsTestData *testData, t_pbc pbc)> algorithms_;
348
349 /*! \brief Test setup function.
350 *
351 * Setting up the pbcs and algorithms. Note, that corresponding string keywords
352 * have to be explicitly added at the end of this file when the tests are called.
353 *
354 */
355 void SetUp() override
356 {
357
358 //
359 // PBC initialization
360 //
361 t_pbc pbc;
362
363 // Infinitely small box
364 matrix boxNone = { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} };
365 set_pbc(&pbc, epbcNONE, boxNone);
366 pbcs_["PBCNone"] = pbc;
367
368 // Rectangular box
369 matrix boxXyz = { {10.0, 0.0, 0.0}, {0.0, 20.0, 0.0}, {0.0, 0.0, 15.0} };
370 set_pbc(&pbc, epbcXYZ, boxXyz);
371 pbcs_["PBCXYZ"] = pbc;
372
373 //
374 // Algorithms
375 //
376 // SHAKE
377 algorithms_["SHAKE"] = applyShake;
378 // LINCS
379 algorithms_["LINCS"] = applyLincs;
380 // LINCS using CUDA (will only be called if CUDA is available)
381 algorithms_["LINCS_CUDA"] = applyLincsCuda;
382 }
383
384 /*! \brief
385 * The test on the final length of constrained bonds.
386 *
387 * Goes through all the constraints and checks if the final length of all the constraints is equal
388 * to the target length with provided tolerance.
389 *
390 * \param[in] tolerance Allowed tolerance in final lengths.
391 * \param[in] testData Test data structure.
392 * \param[in] pbc Periodic boundary data.
393 */
394 void checkConstrainsLength(FloatingPointTolerance tolerance, const ConstraintsTestData &testData, t_pbc pbc)
395 {
396
397 // Test if all the constraints are satisfied
398 for (unsigned c = 0; c < testData.constraints_.size()/3; c++)
399 {
400 real r0 = testData.constraintsR0_.at(testData.constraints_.at(3*c));
401 int i = testData.constraints_.at(3*c + 1);
402 int j = testData.constraints_.at(3*c + 2);
403 RVec xij0, xij1;
404 real d0, d1;
405 if (pbc.ePBC == epbcXYZ)
406 {
407 pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
408 pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
409 }
410 else
411 {
412 rvec_sub(testData.x_[i], testData.x_[j], xij0);
413 rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
414 }
415 d0 = norm(xij0);
416 d1 = norm(xij1);
417 EXPECT_REAL_EQ_TOL(r0, d1, tolerance) << gmx::formatString(
418 "rij = %f, which is not equal to r0 = %f for constraint #%u, between atoms %d and %d"
419 " (before constraining rij was %f).", d1, r0, c, i, j, d0);
420 }
421 }
422
423 /*! \brief
424 * The test on the final length of constrained bonds.
425 *
426 * Goes through all the constraints and checks if the direction of constraint has not changed
427 * by the algorithm (i.e. the constraints algorithm arrived to the solution that is closest
428 * to the initial system conformation).
429 *
430 * \param[in] testData Test data structure.
431 * \param[in] pbc Periodic boundary data.
432 */
433 void checkConstrainsDirection(const ConstraintsTestData &testData, t_pbc pbc)
434 {
435
436 for (unsigned c = 0; c < testData.constraints_.size()/3; c++)
437 {
438 int i = testData.constraints_.at(3*c + 1);
439 int j = testData.constraints_.at(3*c + 2);
440 RVec xij0, xij1;
441 if (pbc.ePBC == epbcXYZ)
442 {
443 pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
444 pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
445 }
446 else
447 {
448 rvec_sub(testData.x_[i], testData.x_[j], xij0);
449 rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
450 }
451
452 real dot = xij0.dot(xij1);
453
454 EXPECT_GE(dot, 0.0) << gmx::formatString(
455 "The constraint %u changed direction. Constraining algorithm might have returned the wrong root "
456 "of the constraints equation.", c);
457
458 }
459 }
460
461 /*! \brief
462 * The test on the coordinates of the center of the mass (COM) of the system.
463 *
464 * Checks if the center of mass has not been shifted by the constraints. Note,
465 * that this test does not take into account the periodic boundary conditions.
466 * Hence it will not work should the constraints decide to move atoms across
467 * PBC borders.
468 *
469 * \param[in] tolerance Allowed tolerance in COM coordinates.
470 * \param[in] testData Test data structure.
471 */
472 void checkCOMCoordinates(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
473 {
474
475 RVec comPrime0({0.0, 0.0, 0.0});
476 RVec comPrime({0.0, 0.0, 0.0});
477 for (int i = 0; i < testData.numAtoms_; i++)
478 {
479 comPrime0 += testData.masses_[i]*testData.xPrime0_[i];
480 comPrime += testData.masses_[i]*testData.xPrime_[i];
481 }
482
483 comPrime0 /= testData.numAtoms_;
484 comPrime /= testData.numAtoms_;
485
486 EXPECT_REAL_EQ_TOL(comPrime[XX], comPrime0[XX], tolerance)
487 << "Center of mass was shifted by constraints in x-direction.";
488 EXPECT_REAL_EQ_TOL(comPrime[YY], comPrime0[YY], tolerance)
489 << "Center of mass was shifted by constraints in y-direction.";
490 EXPECT_REAL_EQ_TOL(comPrime[ZZ], comPrime0[ZZ], tolerance)
491 << "Center of mass was shifted by constraints in z-direction.";
492
493 }
494
495 /*! \brief
496 * The test on the velocity of the center of the mass (COM) of the system.
497 *
498 * Checks if the velocity of the center of mass has not changed.
499 *
500 * \param[in] tolerance Allowed tolerance in COM velocity components.
501 * \param[in] testData Test data structure.
502 */
503 void checkCOMVelocity(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
504 {
505
506 RVec comV0({0.0, 0.0, 0.0});
507 RVec comV({0.0, 0.0, 0.0});
508 for (int i = 0; i < testData.numAtoms_; i++)
509 {
510 comV0 += testData.masses_[i]*testData.v0_[i];
511 comV += testData.masses_[i]*testData.v_[i];
512 }
513 comV0 /= testData.numAtoms_;
514 comV /= testData.numAtoms_;
515
516 EXPECT_REAL_EQ_TOL(comV[XX], comV0[XX], tolerance)
517 << "Velocity of the center of mass in x-direction has been changed by constraints.";
518 EXPECT_REAL_EQ_TOL(comV[YY], comV0[YY], tolerance)
519 << "Velocity of the center of mass in y-direction has been changed by constraints.";
520 EXPECT_REAL_EQ_TOL(comV[ZZ], comV0[ZZ], tolerance)
521 << "Velocity of the center of mass in z-direction has been changed by constraints.";
522 }
523
524 /*! \brief
525 * The test on the final coordinates (not used).
526 *
527 * Goes through all atoms and checks if the final positions correspond to the
528 * provided reference set of coordinates.
529 *
530 * \param[in] xPrimeRef The reference set of coordinates.
531 * \param[in] tolerance Tolerance for the coordinates test.
532 * \param[in] testData Test data structure.
533 */
534 void checkFinalCoordinates(std::vector<RVec> xPrimeRef, FloatingPointTolerance tolerance,
535 const ConstraintsTestData &testData)
536 {
537 for (int i = 0; i < testData.numAtoms_; i++)
538 {
539 for (int d = 0; d < DIM; d++)
540 {
541 EXPECT_REAL_EQ_TOL(xPrimeRef.at(i)[d], testData.xPrime_[i][d], tolerance) << gmx::formatString(
542 "Coordinates after constrains were applied differ from these in the reference set for atom #%d.", i);
543 }
544 }
545 }
546
547 /*! \brief
548 * The test of virial tensor.
549 *
550 * Checks if the values in the scaled virial tensor are equal to pre-computed values.
551 *
552 * \param[in] tolerance Tolerance for the tensor values.
553 * \param[in] testData Test data structure.
554 */
555 void checkVirialTensor(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
556 {
557 for (int i = 0; i < DIM; i++)
558 {
559 for (int j = 0; j < DIM; j++)
560 {
561 EXPECT_REAL_EQ_TOL(testData.virialScaledRef_[i][j], testData.virialScaled_[i][j],
562 tolerance) << gmx::formatString(
563 "Values in virial tensor at [%d][%d] are not within the tolerance from reference value.", i, j);
564 }
565 }
566 }
567
568 /*! \brief
569 * The test for FEP (not used).
570 *
571 * Checks if the value of dH/dLambda is equal to the reference value.
572 * \todo Add tests for dHdLambda values.
573 *
574 * \param[in] dHdLambdaRef Reference value.
575 * \param[in] tolerance Tolerance.
576 * \param[in] testData Test data structure.
577 */
578 void checkFEP(const real dHdLambdaRef, FloatingPointTolerance tolerance,
579 const ConstraintsTestData &testData)
580 {
581 EXPECT_REAL_EQ_TOL(dHdLambdaRef, testData.dHdLambda_, tolerance) <<
582 "Computed value for dV/dLambda is not equal to the reference value. ";
583 }
584
585 };
586
587 TEST_P(ConstraintsTest, SingleConstraint){
588 std::string title = "one constraint (e.g. OH)";
589 int numAtoms = 2;
590
591 std::vector<real> masses = {1.0, 12.0};
592 std::vector<int> constraints = {0, 0, 1};
593 std::vector<real> constraintsR0 = {0.1};
594
595 real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
596
597 std::vector<RVec> x = {{ 0.0, oneTenthOverSqrtTwo, 0.0 },
598 { oneTenthOverSqrtTwo, 0.0, 0.0 }};
599 std::vector<RVec> xPrime = {{ 0.01, 0.08, 0.01 },
600 { 0.06, 0.01, -0.01 }};
601 std::vector<RVec> v = {{ 1.0, 2.0, 3.0 },
602 { 3.0, 2.0, 1.0 }};
603
604 tensor virialScaledRef = {{-5.58e-04, 5.58e-04, 0.00e+00 },
605 { 5.58e-04, -5.58e-04, 0.00e+00 },
606 { 0.00e+00, 0.00e+00, 0.00e+00 }};
607
608 real shakeTolerance = 0.0001;
609 gmx_bool shakeUseSOR = false;
610
611 int lincsNIter = 1;
612 int lincslincsExpansionOrder = 4;
613 real lincsWarnAngle = 30.0;
614
615 std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
616 (title, numAtoms, masses,
617 constraints, constraintsR0,
618 true, virialScaledRef,
619 false, 0,
620 real(0.0), real(0.001),
621 x, xPrime, v,
622 shakeTolerance, shakeUseSOR,
623 lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
624 std::string pbcName;
625 std::string algorithmName;
626 std::tie(pbcName, algorithmName) = GetParam();
627 t_pbc pbc = pbcs_.at(pbcName);
628
629 // Apply constraints
630 algorithms_.at(algorithmName)(testData.get(), pbc);
631
632 checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
633 checkConstrainsDirection(*testData, pbc);
634 checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
635 checkCOMVelocity(absoluteTolerance(0.0001), *testData);
636
637 checkVirialTensor(absoluteTolerance(0.0001), *testData);
638
639 }
640
641 TEST_P(ConstraintsTest, TwoDisjointConstraints){
642
643 std::string title = "two disjoint constraints";
644 int numAtoms = 4;
645 std::vector<real> masses = {0.5, 1.0/3.0, 0.25, 1.0};
646 std::vector<int> constraints = {0, 0, 1, 1, 2, 3};
647 std::vector<real> constraintsR0 = {2.0, 1.0};
648
649
650 std::vector<RVec> x = {{ 2.50, -3.10, 15.70 },
651 { 0.51, -3.02, 15.55 },
652 { -0.50, -3.00, 15.20 },
653 { -1.51, -2.95, 15.05 }};
654
655 std::vector<RVec> xPrime = {{ 2.50, -3.10, 15.70 },
656 { 0.51, -3.02, 15.55 },
657 { -0.50, -3.00, 15.20 },
658 { -1.51, -2.95, 15.05 }};
659
660 std::vector<RVec> v = {{ 0.0, 1.0, 0.0 },
661 { 1.0, 0.0, 0.0 },
662 { 0.0, 0.0, 1.0 },
663 { 0.0, 0.0, 0.0 }};
664
665 tensor virialScaledRef = {{ 3.3e-03, -1.7e-04, 5.6e-04 },
666 {-1.7e-04, 8.9e-06, -2.8e-05 },
667 { 5.6e-04, -2.8e-05, 8.9e-05 }};
668
669 real shakeTolerance = 0.0001;
670 gmx_bool shakeUseSOR = false;
671
672 int lincsNIter = 1;
673 int lincslincsExpansionOrder = 4;
674 real lincsWarnAngle = 30.0;
675
676 std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
677 (title, numAtoms, masses,
678 constraints, constraintsR0,
679 true, virialScaledRef,
680 false, 0,
681 real(0.0), real(0.001),
682 x, xPrime, v,
683 shakeTolerance, shakeUseSOR,
684 lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
685
686 std::string pbcName;
687 std::string algorithmName;
688 std::tie(pbcName, algorithmName) = GetParam();
689 t_pbc pbc = pbcs_.at(pbcName);
690
691 // Apply constraints
692 algorithms_.at(algorithmName)(testData.get(), pbc);
693
694 checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
695 checkConstrainsDirection(*testData, pbc);
696 checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
697 checkCOMVelocity(absoluteTolerance(0.0001), *testData);
698
699 checkVirialTensor(absoluteTolerance(0.0001), *testData);
700
701 }
702
703 TEST_P(ConstraintsTest, ThreeSequentialConstraints){
704
705 std::string title = "three atoms, connected longitudinally (e.g. CH2)";
706 int numAtoms = 3;
707 std::vector<real> masses = {1.0, 12.0, 16.0 };
708 std::vector<int> constraints = {0, 0, 1, 1, 1, 2};
709 std::vector<real> constraintsR0 = {0.1, 0.2};
710
711 real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
712 real twoTenthsOverSqrtThree = 0.2_real / std::sqrt(3.0_real);
713
714 std::vector<RVec> x = {{ oneTenthOverSqrtTwo, oneTenthOverSqrtTwo, 0.0 },
715 { 0.0, 0.0, 0.0 },
716 { twoTenthsOverSqrtThree, twoTenthsOverSqrtThree, twoTenthsOverSqrtThree }};
717
718 std::vector<RVec> xPrime = {{ 0.08, 0.07, 0.01 },
719 { -0.02, 0.01, -0.02 },
720 { 0.10, 0.12, 0.11 }};
721
722 std::vector<RVec> v = {{ 1.0, 0.0, 0.0 },
723 { 0.0, 1.0, 0.0 },
724 { 0.0, 0.0, 1.0 }};
725
726 tensor virialScaledRef = {{ 4.14e-03, 4.14e-03, 3.31e-03},
727 { 4.14e-03, 4.14e-03, 3.31e-03},
728 { 3.31e-03, 3.31e-03, 3.31e-03}};
729
730 real shakeTolerance = 0.0001;
731 gmx_bool shakeUseSOR = false;
732
733 int lincsNIter = 1;
734 int lincslincsExpansionOrder = 4;
735 real lincsWarnAngle = 30.0;
736
737 std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
738 (title, numAtoms, masses,
739 constraints, constraintsR0,
740 true, virialScaledRef,
741 false, 0,
742 real(0.0), real(0.001),
743 x, xPrime, v,
744 shakeTolerance, shakeUseSOR,
745 lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
746
747 std::string pbcName;
748 std::string algorithmName;
749 std::tie(pbcName, algorithmName) = GetParam();
750 t_pbc pbc = pbcs_.at(pbcName);
751
752 // Apply constraints
753 algorithms_.at(algorithmName)(testData.get(), pbc);
754
755 checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
756 checkConstrainsDirection(*testData, pbc);
757 checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
758 checkCOMVelocity(absoluteTolerance(0.0001), *testData);
759
760 checkVirialTensor(absoluteTolerance(0.0001), *testData);
761
762 }
763
764 TEST_P(ConstraintsTest, ThreeConstraintsWithCentralAtom){
765
766 std::string title = "three atoms, connected to the central atom (e.g. CH3)";
767 int numAtoms = 4;
768 std::vector<real> masses = {12.0, 1.0, 1.0, 1.0 };
769 std::vector<int> constraints = {0, 0, 1, 0, 0, 2, 0, 0, 3};
770 std::vector<real> constraintsR0 = {0.1};
771
772
773 std::vector<RVec> x = {{ 0.00, 0.00, 0.00 },
774 { 0.10, 0.00, 0.00 },
775 { 0.00, -0.10, 0.00 },
776 { 0.00, 0.00, 0.10 }};
777
778 std::vector<RVec> xPrime = {{ 0.004, 0.009, -0.010 },
779 { 0.110, -0.006, 0.003 },
780 {-0.007, -0.102, -0.007 },
781 {-0.005, 0.011, 0.102 }};
782
783 std::vector<RVec> v = {{ 1.0, 0.0, 0.0 },
784 { 1.0, 0.0, 0.0 },
785 { 1.0, 0.0, 0.0 },
786 { 1.0, 0.0, 0.0 }};
787
788 tensor virialScaledRef = {{7.14e-04, 0.00e+00, 0.00e+00},
789 {0.00e+00, 1.08e-03, 0.00e+00},
790 {0.00e+00, 0.00e+00, 1.15e-03}};
791
792 real shakeTolerance = 0.0001;
793 gmx_bool shakeUseSOR = false;
794
795 int lincsNIter = 1;
796 int lincslincsExpansionOrder = 4;
797 real lincsWarnAngle = 30.0;
798
799 std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
800 (title, numAtoms, masses,
801 constraints, constraintsR0,
802 true, virialScaledRef,
803 false, 0,
804 real(0.0), real(0.001),
805 x, xPrime, v,
806 shakeTolerance, shakeUseSOR,
807 lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
808
809 std::string pbcName;
810 std::string algorithmName;
811 std::tie(pbcName, algorithmName) = GetParam();
812 t_pbc pbc = pbcs_.at(pbcName);
813
814 // Apply constraints
815 algorithms_.at(algorithmName)(testData.get(), pbc);
816
817 checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
818 checkConstrainsDirection(*testData, pbc);
819 checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
820 checkCOMVelocity(absoluteTolerance(0.0001), *testData);
821
822 checkVirialTensor(absoluteTolerance(0.0001), *testData);
823 }
824
825 TEST_P(ConstraintsTest, FourSequentialConstraints){
826
827 std::string title = "four atoms, connected longitudinally";
828 int numAtoms = 4;
829 std::vector<real> masses = {0.5, 1.0/3.0, 0.25, 1.0};
830 std::vector<int> constraints = {0, 0, 1, 1, 1, 2, 2, 2, 3};
831 std::vector<real> constraintsR0 = {2.0, 1.0, 1.0};
832
833
834 std::vector<RVec> x = {{ 2.50, -3.10, 15.70 },
835 { 0.51, -3.02, 15.55 },
836 { -0.50, -3.00, 15.20 },
837 { -1.51, -2.95, 15.05 }};
838
839 std::vector<RVec> xPrime = {{ 2.50, -3.10, 15.70 },
840 { 0.51, -3.02, 15.55 },
841 { -0.50, -3.00, 15.20 },
842 { -1.51, -2.95, 15.05 }};
843
844 std::vector<RVec> v = {{ 0.0, 0.0, 2.0 },
845 { 0.0, 0.0, 3.0 },
846 { 0.0, 0.0, -4.0 },
847 { 0.0, 0.0, -1.0 }};
848
849 tensor virialScaledRef = {{ 1.15e-01, -4.20e-03, 2.12e-02},
850 {-4.20e-03, 1.70e-04, -6.41e-04},
851 { 2.12e-02, -6.41e-04, 5.45e-03}};
852
853 real shakeTolerance = 0.0001;
854 gmx_bool shakeUseSOR = false;
855
856 int lincsNIter = 4;
857 int lincslincsExpansionOrder = 8;
858 real lincsWarnAngle = 30.0;
859
860 std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
861 (title, numAtoms, masses,
862 constraints, constraintsR0,
863 true, virialScaledRef,
864 false, 0,
865 real(0.0), real(0.001),
866 x, xPrime, v,
867 shakeTolerance, shakeUseSOR,
868 lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
869
870 std::string pbcName;
871 std::string algorithmName;
872 std::tie(pbcName, algorithmName) = GetParam();
873 t_pbc pbc = pbcs_.at(pbcName);
874
875 // Apply constraints
876 algorithms_.at(algorithmName)(testData.get(), pbc);
877
878 checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
879 checkConstrainsDirection(*testData, pbc);
880 checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
881 checkCOMVelocity(absoluteTolerance(0.0001), *testData);
882
883 checkVirialTensor(absoluteTolerance(0.01), *testData);
884
885 }
886
887 TEST_P(ConstraintsTest, TriangleOfConstraints){
888
889 std::string title = "basic triangle (tree atoms, connected to each other)";
890 int numAtoms = 3;
891 std::vector<real> masses = {1.0, 1.0, 1.0};
892 std::vector<int> constraints = {0, 0, 1, 2, 0, 2, 1, 1, 2};
893 std::vector<real> constraintsR0 = {0.1, 0.1, 0.1};
894
895 real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
896
897 std::vector<RVec> x = {{ oneTenthOverSqrtTwo, 0.0, 0.0 },
898 { 0.0, oneTenthOverSqrtTwo, 0.0 },
899 { 0.0, 0.0, oneTenthOverSqrtTwo }};
900
901 std::vector<RVec> xPrime = {{ 0.09, -0.02, 0.01 },
902 { -0.02, 0.10, -0.02 },
903 { 0.03, -0.01, 0.07 }};
904
905 std::vector<RVec> v = {{ 1.0, 1.0, 1.0 },
906 { -2.0, -2.0, -2.0 },
907 { 1.0, 1.0, 1.0 }};
908
909 tensor virialScaledRef = {{ 6.00e-04, -1.61e-03, 1.01e-03},
910 {-1.61e-03, 2.53e-03, -9.25e-04},
911 { 1.01e-03, -9.25e-04, -8.05e-05}};
912
913 real shakeTolerance = 0.0001;
914 gmx_bool shakeUseSOR = false;
915
916 int lincsNIter = 1;
917 int lincslincsExpansionOrder = 4;
918 real lincsWarnAngle = 30.0;
919
920 std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
921 (title, numAtoms, masses,
922 constraints, constraintsR0,
923 true, virialScaledRef,
924 false, 0,
925 real(0.0), real(0.001),
926 x, xPrime, v,
927 shakeTolerance, shakeUseSOR,
928 lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
929
930 std::string pbcName;
931 std::string algorithmName;
932 std::tie(pbcName, algorithmName) = GetParam();
933 t_pbc pbc = pbcs_.at(pbcName);
934
935 // Apply constraints
936 algorithms_.at(algorithmName)(testData.get(), pbc);
937
938 checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
939 checkConstrainsDirection(*testData, pbc);
940 checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
941 checkCOMVelocity(absoluteTolerance(0.0001), *testData);
942
943 checkVirialTensor(absoluteTolerance(0.00001), *testData);
944
945 }
946
947
948 INSTANTIATE_TEST_CASE_P(WithParameters, ConstraintsTest,
949 ::testing::Combine(::testing::Values("PBCNone", "PBCXYZ"),
950 ::testing::ValuesIn(getAlgorithmsNames())));
951
952 } // namespace
953 } // namespace test
954 } // namespace gmx
The ultimate molecular dynamics simulation package