2 * This file is part of the GROMACS molecular simulation package.
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.
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.
36 * \brief SHAKE and LINCS tests.
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
44 * \author Artem Zhmurov <zhmurov@gmail.com>
45 * \ingroup module_mdlib
50 #include "gromacs/mdlib/constr.h"
59 #include <unordered_map>
62 #include <gtest/gtest.h>
64 #include "gromacs/fileio/gmxfio.h"
65 #include "gromacs/gmxlib/nrnb.h"
66 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
67 #include "gromacs/gpu_utils/gpu_utils.h"
68 #include "gromacs/math/paddedvector.h"
69 #include "gromacs/math/vec.h"
70 #include "gromacs/math/vectypes.h"
71 #include "gromacs/mdlib/gmx_omp_nthreads.h"
72 #include "gromacs/mdlib/lincs.h"
73 #include "gromacs/mdlib/lincs_cuda.h"
74 #include "gromacs/mdlib/shake.h"
75 #include "gromacs/mdtypes/commrec.h"
76 #include "gromacs/mdtypes/inputrec.h"
77 #include "gromacs/mdtypes/mdatom.h"
78 #include "gromacs/pbcutil/pbc.h"
79 #include "gromacs/topology/block.h"
80 #include "gromacs/topology/idef.h"
81 #include "gromacs/topology/ifunc.h"
82 #include "gromacs/topology/topology.h"
83 #include "gromacs/utility/smalloc.h"
84 #include "gromacs/utility/stringutil.h"
85 #include "gromacs/utility/unique_cptr.h"
87 #include "testutils/refdata.h"
88 #include "testutils/testasserts.h"
96 * Constraints test data structure.
98 * Structure to collect all the necessary data, including system coordinates and topology,
99 * constraints information, etc. The structure can be reset and reused.
101 struct ConstraintsTestData
104 //! Human-friendly name for a system
111 std::vector
<real
> masses_
;
113 std::vector
<real
> invmass_
;
114 //! Communication record
116 //! Input record (info that usually in .mdp file)
124 //! Computational time array (normally used to benchmark performance)
129 //! Number of flexible constraints
131 //! Whether the virial should be computed
134 tensor virialScaled_
;
135 //! Scaled virial (reference values)
136 tensor virialScaledRef_
;
137 //! If the free energy is computed
138 bool compute_dHdLambda_
;
139 //! For free energy computation
141 //! For free energy computation (reference value)
144 //! Coordinates before the timestep
145 PaddedVector
<RVec
> x_
;
146 //! Coordinates after timestep, output for the constraints
147 PaddedVector
<RVec
> xPrime_
;
148 //! Backup for coordinates (for reset)
149 PaddedVector
<RVec
> xPrime0_
;
150 //! Intermediate set of coordinates (normally used for projection correction)
151 PaddedVector
<RVec
> xPrime2_
;
153 PaddedVector
<RVec
> v_
;
154 //! Backup for velocities (for reset)
155 PaddedVector
<RVec
> v0_
;
157 //! Constraints data (type1-i1-j1-type2-i2-j2-...)
158 std::vector
<int> constraints_
;
159 //! Target lengths for all constraint types
160 std::vector
<real
> constraintsR0_
;
163 * Constructor for the object with all parameters and variables needed by constraints algorithms.
165 * This constructor assembles stubs for all the data structures, required to initialize
166 * and apply LINCS and SHAKE constraints. The coordinates and velocities before constraining
167 * are saved to allow for reset. The constraints data are stored for testing after constraints
170 * \param[in] title Human-friendly name of the system.
171 * \param[in] numAtoms Number of atoms in the system.
172 * \param[in] masses Atom masses. Size of this vector should be equal to numAtoms.
173 * \param[in] constraints List of constraints, organized in triples of integers.
174 * First integer is the index of type for a constraint, second
175 * and third are the indices of constrained atoms. The types
176 * of constraints should be sequential but not necessarily
177 * start from zero (which is the way they normally are in
179 * \param[in] constraintsR0 Target values for bond lengths for bonds of each type. The
180 * size of this vector should be equal to the total number of
181 * unique types in constraints vector.
182 * \param[in] computeVirial Whether the virial should be computed.
183 * \param[in] virialScaledRef Reference values for scaled virial tensor.
184 * \param[in] compute_dHdLambda Whether free energy should be computed.
185 * \param[in] dHdLambdaRef Reference value for dHdLambda.
186 * \param[in] initialTime Initial time.
187 * \param[in] timestep Timestep.
188 * \param[in] x Coordinates before integration step.
189 * \param[in] xPrime Coordinates after integration step, but before constraining.
190 * \param[in] v Velocities before constraining.
191 * \param[in] shakeTolerance Target tolerance for SHAKE.
192 * \param[in] shakeUseSOR Use successive over-relaxation method for SHAKE iterations.
193 * The general formula is:
194 * x_n+1 = (1-omega)*x_n + omega*f(x_n),
195 * where omega = 1 if SOR is off and may be < 1 if SOR is on.
196 * \param[in] lincsNumIterations Number of iterations used to compute the inverse matrix.
197 * \param[in] lincsExpansionOrder The order for algorithm that adjusts the direction of the
198 * bond after constraints are applied.
199 * \param[in] lincsWarnAngle The threshold value for the change in bond angle. When
200 * exceeded the program will issue a warning.
203 ConstraintsTestData(const std::string
&title
,
204 int numAtoms
, std::vector
<real
> masses
,
205 std::vector
<int> constraints
, std::vector
<real
> constraintsR0
,
206 bool computeVirial
, tensor virialScaledRef
,
207 bool compute_dHdLambda
, float dHdLambdaRef
,
208 real initialTime
, real timestep
,
209 const std::vector
<RVec
> &x
, const std::vector
<RVec
> &xPrime
, const std::vector
<RVec
> &v
,
210 real shakeTolerance
, gmx_bool shakeUseSOR
,
211 int lincsNumIterations
, int lincsExpansionOrder
, real lincsWarnAngle
)
213 title_
= title
; // Human-friendly name of the system
214 numAtoms_
= numAtoms
; // Number of atoms
218 invmass_
.resize(numAtoms
); // Vector of inverse masses
220 for (int i
= 0; i
< numAtoms
; i
++)
222 invmass_
[i
] = 1.0/masses
.at(i
);
225 // Saving constraints to check if they are satisfied after algorithm was applied
226 constraints_
= constraints
; // Constraints indices (in type-i-j format)
227 constraintsR0_
= constraintsR0
; // Equilibrium distances for each type of constraint
229 invdt_
= 1.0/timestep
; // Inverse timestep
231 // Communication record
239 // Input record - data that usually comes from configuration file (.mdp)
241 ir_
.init_t
= initialTime
;
242 ir_
.delta_t
= timestep
;
246 md_
.nMassPerturbed
= 0;
248 md_
.invmass
= invmass_
.data();
250 md_
.homenr
= numAtoms
;
253 computeVirial_
= computeVirial
;
256 for (int i
= 0; i
< DIM
; i
++)
258 for (int j
= 0; j
< DIM
; j
++)
260 virialScaled_
[i
][j
] = 0;
261 virialScaledRef_
[i
][j
] = virialScaledRef
[i
][j
];
267 // Free energy evaluation
268 compute_dHdLambda_
= compute_dHdLambda
;
270 if (compute_dHdLambda_
)
273 dHdLambdaRef_
= dHdLambdaRef
;
281 // Constraints and their parameters (local topology)
282 for (int i
= 0; i
< F_NRE
; i
++)
286 idef_
.il
[F_CONSTR
].nr
= constraints
.size();
288 snew(idef_
.il
[F_CONSTR
].iatoms
, constraints
.size());
290 for (unsigned i
= 0; i
< constraints
.size(); i
++)
294 if (maxType
< constraints
.at(i
))
296 maxType
= constraints
.at(i
);
299 idef_
.il
[F_CONSTR
].iatoms
[i
] = constraints
.at(i
);
301 snew(idef_
.iparams
, maxType
+ 1);
302 for (unsigned i
= 0; i
< constraints
.size()/3; i
++)
304 idef_
.iparams
[constraints
.at(3*i
)].constr
.dA
= constraintsR0
.at(constraints
.at(3*i
));
305 idef_
.iparams
[constraints
.at(3*i
)].constr
.dB
= constraintsR0
.at(constraints
.at(3*i
));
308 // Constraints and their parameters (global topology)
309 InteractionList interactionList
;
310 interactionList
.iatoms
.resize(constraints
.size());
311 for (unsigned i
= 0; i
< constraints
.size(); i
++)
313 interactionList
.iatoms
.at(i
) = constraints
.at(i
);
315 InteractionList interactionListEmpty
;
316 interactionListEmpty
.iatoms
.resize(0);
318 gmx_moltype_t molType
;
319 molType
.atoms
.nr
= numAtoms
;
320 molType
.ilist
.at(F_CONSTR
) = interactionList
;
321 molType
.ilist
.at(F_CONSTRNC
) = interactionListEmpty
;
322 mtop_
.moltype
.push_back(molType
);
324 gmx_molblock_t molBlock
;
327 mtop_
.molblock
.push_back(molBlock
);
329 mtop_
.natoms
= numAtoms
;
330 mtop_
.ffparams
.iparams
.resize(maxType
+ 1);
331 for (int i
= 0; i
<= maxType
; i
++)
333 mtop_
.ffparams
.iparams
.at(i
) = idef_
.iparams
[i
];
335 mtop_
.bIntermolecularInteractions
= false;
337 // Coordinates and velocities
338 x_
.resizeWithPadding(numAtoms
);
339 xPrime_
.resizeWithPadding(numAtoms
);
340 xPrime0_
.resizeWithPadding(numAtoms
);
341 xPrime2_
.resizeWithPadding(numAtoms
);
343 v_
.resizeWithPadding(numAtoms
);
344 v0_
.resizeWithPadding(numAtoms
);
346 std::copy(x
.begin(), x
.end(), x_
.begin());
347 std::copy(xPrime
.begin(), xPrime
.end(), xPrime_
.begin());
348 std::copy(xPrime
.begin(), xPrime
.end(), xPrime0_
.begin());
349 std::copy(xPrime
.begin(), xPrime
.end(), xPrime2_
.begin());
351 std::copy(v
.begin(), v
.end(), v_
.begin());
352 std::copy(v
.begin(), v
.end(), v0_
.begin());
354 // SHAKE-specific parameters
355 ir_
.shake_tol
= shakeTolerance
;
356 ir_
.bShakeSOR
= shakeUseSOR
;
358 // LINCS-specific parameters
359 ir_
.nLincsIter
= lincsNumIterations
;
360 ir_
.nProjOrder
= lincsExpansionOrder
;
361 ir_
.LincsWarnAngle
= lincsWarnAngle
;
365 * Reset the data structure so it can be reused.
367 * Set the coordinates and velocities back to their values before
368 * constraining. The scaled virial tensor and dHdLambda are zeroed.
379 for (int i
= 0; i
< DIM
; i
++)
381 for (int j
= 0; j
< DIM
; j
++)
383 virialScaled_
[i
][j
] = 0;
391 * Cleaning up the memory.
393 ~ConstraintsTestData()
395 sfree(idef_
.il
[F_CONSTR
].iatoms
);
396 sfree(idef_
.iparams
);
401 /*! \brief The two-dimensional parameter space for test.
403 * The test will run for all possible combinations of accessible
405 * 1. PBC setup ("PBCNONE" or "PBCXYZ")
406 * 2. The algorithm ("SHAKE", "LINCS" or "LINCS_GPU").
408 typedef std::tuple
<std::string
, std::string
> ConstraintsTestParameters
;
410 /*! \brief Names of all availible algorithms
412 * Constructed from the algorithms_ field of the test class.
413 * Used as the list of values of second parameter in ConstraintsTestParameters.
415 static std::vector
<std::string
> algorithmsNames
;
418 /*! \brief Test fixture for constraints.
420 * The fixture uses following test systems:
421 * 1. Two atoms, connected with one constraint (e.g. NH).
422 * 2. Three atoms, connected consequently with two constraints (e.g. CH2).
423 * 3. Three atoms, constrained to the fourth atom (e.g. CH3).
424 * 4. Four atoms, connected by two independent constraints.
425 * 5. Three atoms, connected by three constraints in a triangle
426 * (e.g. H2O with constrained H-O-H angle).
427 * 6. Four atoms, connected by three consequential constraints.
429 * For all systems, the final lengths of the constraints are tested against the
430 * reference values, the direction of each constraint is checked.
431 * Test also verifies that the center of mass has not been
432 * shifted by the constraints and that its velocity has not changed.
433 * For some systems, the value for scaled virial tensor is checked against
436 class ConstraintsTest
: public ::testing::TestWithParam
<ConstraintsTestParameters
>
440 std::unordered_map
<std::string
, t_pbc
> pbcs_
;
441 //! Algorithms (SHAKE and LINCS)
442 std::unordered_map
<std::string
, void(*)(ConstraintsTestData
*testData
, t_pbc pbc
)> algorithms_
;
444 /*! \brief Test setup function.
446 * Setting up the pbcs and algorithms. Note, that corresponding string keywords
447 * have to be explicitly added at the end of this file when the tests are called.
450 void SetUp() override
454 // PBC initialization
458 // Infinitely small box
459 matrix boxNone
= { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} };
460 set_pbc(&pbc
, epbcNONE
, boxNone
);
461 pbcs_
["PBCNone"] = pbc
;
464 matrix boxXyz
= { {10.0, 0.0, 0.0}, {0.0, 20.0, 0.0}, {0.0, 0.0, 15.0} };
465 set_pbc(&pbc
, epbcXYZ
, boxXyz
);
466 pbcs_
["PBCXYZ"] = pbc
;
472 algorithms_
["SHAKE"] = applyShake
;
474 algorithms_
["LINCS"] = applyLincs
;
475 // LINCS using CUDA (if CUDA is available)
476 #if GMX_GPU == GMX_GPU_CUDA
477 std::string errorMessage
;
478 if (canDetectGpus(&errorMessage
))
480 algorithms_
["LINCS_CUDA"] = applyLincsCuda
;
483 for (const auto &algorithm
: algorithms_
)
485 algorithmsNames
.push_back(algorithm
.first
);
490 * Initialize and apply SHAKE constraints.
492 * \param[in] testData Test data structure.
493 * \param[in] pbc Periodic boundary data (not used in SHAKE).
495 static void applyShake(ConstraintsTestData
*testData
, t_pbc gmx_unused pbc
)
497 shakedata
* shaked
= shake_init();
498 make_shake_sblock_serial(shaked
, &testData
->idef_
, testData
->md_
);
499 bool success
= constrain_shake(
502 testData
->invmass_
.data(),
505 as_rvec_array(testData
->x_
.data()),
506 as_rvec_array(testData
->xPrime_
.data()),
507 as_rvec_array(testData
->xPrime2_
.data()),
509 testData
->md_
.lambda
,
510 &testData
->dHdLambda_
,
512 as_rvec_array(testData
->v_
.data()),
513 testData
->computeVirial_
,
514 testData
->virialScaled_
,
516 gmx::ConstraintVariable::Positions
);
517 EXPECT_TRUE(success
) << "Test failed with a false return value in SHAKE.";
522 * Initialize and apply LINCS constraints.
524 * \param[in] testData Test data structure.
525 * \param[in] pbc Periodic boundary data.
527 static void applyLincs(ConstraintsTestData
*testData
, t_pbc pbc
)
532 int warncount_lincs
= 0;
533 gmx_omp_nthreads_set(emntLINCS
, 1);
535 // Make blocka structure for faster LINCS setup
536 std::vector
<t_blocka
> at2con_mt
;
537 at2con_mt
.reserve(testData
->mtop_
.moltype
.size());
538 for (const gmx_moltype_t
&moltype
: testData
->mtop_
.moltype
)
540 // This function is in constr.cpp
541 at2con_mt
.push_back(make_at2con(moltype
,
542 testData
->mtop_
.ffparams
.iparams
,
543 flexibleConstraintTreatment(EI_DYNAMICS(testData
->ir_
.eI
))));
546 lincsd
= init_lincs(nullptr, testData
->mtop_
,
547 testData
->nflexcon_
, at2con_mt
,
549 testData
->ir_
.nLincsIter
, testData
->ir_
.nProjOrder
);
550 set_lincs(testData
->idef_
, testData
->md_
, EI_DYNAMICS(testData
->ir_
.eI
), &testData
->cr_
, lincsd
);
552 // Evaluate constraints
553 bool sucess
= constrain_lincs(false, testData
->ir_
, 0, lincsd
, testData
->md_
,
556 as_rvec_array(testData
->x_
.data()),
557 as_rvec_array(testData
->xPrime_
.data()),
558 as_rvec_array(testData
->xPrime2_
.data()),
560 testData
->md_
.lambda
, &testData
->dHdLambda_
,
562 as_rvec_array(testData
->v_
.data()),
563 testData
->computeVirial_
, testData
->virialScaled_
,
564 gmx::ConstraintVariable::Positions
,
566 maxwarn
, &warncount_lincs
);
567 EXPECT_TRUE(sucess
) << "Test failed with a false return value in LINCS.";
568 EXPECT_EQ(warncount_lincs
, 0) << "There were warnings in LINCS.";
569 for (unsigned int i
= 0; i
< testData
->mtop_
.moltype
.size(); i
++)
571 sfree(at2con_mt
.at(i
).index
);
572 sfree(at2con_mt
.at(i
).a
);
578 * Initialize and apply LINCS constraints on CUDA-enabled GPU.
580 * \param[in] testData Test data structure.
581 * \param[in] pbc Periodic boundary data.
583 static void applyLincsCuda(ConstraintsTestData
*testData
, t_pbc pbc
)
585 auto lincsCuda
= std::make_unique
<LincsCuda
>(testData
->numAtoms_
,
586 testData
->ir_
.nLincsIter
,
587 testData
->ir_
.nProjOrder
);
588 lincsCuda
->set(testData
->idef_
, testData
->md_
);
589 lincsCuda
->setPbc(&pbc
);
590 lincsCuda
->copyCoordinatesToGpu(as_rvec_array(testData
->x_
.data()),
591 as_rvec_array(testData
->xPrime_
.data()));
592 lincsCuda
->copyVelocitiesToGpu(as_rvec_array(testData
->v_
.data()));
593 lincsCuda
->apply(true, testData
->invdt_
, testData
->computeVirial_
, testData
->virialScaled_
);
594 lincsCuda
->copyCoordinatesFromGpu(as_rvec_array(testData
->xPrime_
.data()));
595 lincsCuda
->copyVelocitiesFromGpu(as_rvec_array(testData
->v_
.data()));
599 * The test on the final length of constrained bonds.
601 * Goes through all the constraints and checks if the final length of all the constraints is equal
602 * to the target length with provided tolerance.
604 * \param[in] tolerance Allowed tolerance in final lengths.
605 * \param[in] testData Test data structure.
606 * \param[in] pbc Periodic boundary data.
608 void checkConstrainsLength(FloatingPointTolerance tolerance
, const ConstraintsTestData
&testData
, t_pbc pbc
)
611 // Test if all the constraints are satisfied
612 for (unsigned c
= 0; c
< testData
.constraints_
.size()/3; c
++)
614 real r0
= testData
.constraintsR0_
.at(testData
.constraints_
.at(3*c
));
615 int i
= testData
.constraints_
.at(3*c
+ 1);
616 int j
= testData
.constraints_
.at(3*c
+ 2);
619 if (pbc
.ePBC
== epbcXYZ
)
621 pbc_dx_aiuc(&pbc
, testData
.x_
[i
], testData
.x_
[j
], xij0
);
622 pbc_dx_aiuc(&pbc
, testData
.xPrime_
[i
], testData
.xPrime_
[j
], xij1
);
626 rvec_sub(testData
.x_
[i
], testData
.x_
[j
], xij0
);
627 rvec_sub(testData
.xPrime_
[i
], testData
.xPrime_
[j
], xij1
);
631 EXPECT_REAL_EQ_TOL(r0
, d1
, tolerance
) << gmx::formatString(
632 "rij = %f, which is not equal to r0 = %f for constraint #%u, between atoms %d and %d"
633 " (before constraining rij was %f).", d1
, r0
, c
, i
, j
, d0
);
638 * The test on the final length of constrained bonds.
640 * Goes through all the constraints and checks if the direction of constraint has not changed
641 * by the algorithm (i.e. the constraints algorithm arrived to the solution that is closest
642 * to the initial system conformation).
644 * \param[in] testData Test data structure.
645 * \param[in] pbc Periodic boundary data.
647 void checkConstrainsDirection(const ConstraintsTestData
&testData
, t_pbc pbc
)
650 for (unsigned c
= 0; c
< testData
.constraints_
.size()/3; c
++)
652 int i
= testData
.constraints_
.at(3*c
+ 1);
653 int j
= testData
.constraints_
.at(3*c
+ 2);
655 if (pbc
.ePBC
== epbcXYZ
)
657 pbc_dx_aiuc(&pbc
, testData
.x_
[i
], testData
.x_
[j
], xij0
);
658 pbc_dx_aiuc(&pbc
, testData
.xPrime_
[i
], testData
.xPrime_
[j
], xij1
);
662 rvec_sub(testData
.x_
[i
], testData
.x_
[j
], xij0
);
663 rvec_sub(testData
.xPrime_
[i
], testData
.xPrime_
[j
], xij1
);
666 real dot
= xij0
.dot(xij1
);
668 EXPECT_GE(dot
, 0.0) << gmx::formatString(
669 "The constraint %u changed direction. Constraining algorithm might have returned the wrong root "
670 "of the constraints equation.", c
);
676 * The test on the coordinates of the center of the mass (COM) of the system.
678 * Checks if the center of mass has not been shifted by the constraints. Note,
679 * that this test does not take into account the periodic boundary conditions.
680 * Hence it will not work should the constraints decide to move atoms across
683 * \param[in] tolerance Allowed tolerance in COM coordinates.
684 * \param[in] testData Test data structure.
686 void checkCOMCoordinates(FloatingPointTolerance tolerance
, const ConstraintsTestData
&testData
)
689 RVec
comPrime0({0.0, 0.0, 0.0});
690 RVec
comPrime({0.0, 0.0, 0.0});
691 for (int i
= 0; i
< testData
.numAtoms_
; i
++)
693 comPrime0
+= testData
.masses_
[i
]*testData
.xPrime0_
[i
];
694 comPrime
+= testData
.masses_
[i
]*testData
.xPrime_
[i
];
697 comPrime0
/= testData
.numAtoms_
;
698 comPrime
/= testData
.numAtoms_
;
700 EXPECT_REAL_EQ_TOL(comPrime
[XX
], comPrime0
[XX
], tolerance
)
701 << "Center of mass was shifted by constraints in x-direction.";
702 EXPECT_REAL_EQ_TOL(comPrime
[YY
], comPrime0
[YY
], tolerance
)
703 << "Center of mass was shifted by constraints in y-direction.";
704 EXPECT_REAL_EQ_TOL(comPrime
[ZZ
], comPrime0
[ZZ
], tolerance
)
705 << "Center of mass was shifted by constraints in z-direction.";
710 * The test on the velocity of the center of the mass (COM) of the system.
712 * Checks if the velocity of the center of mass has not changed.
714 * \param[in] tolerance Allowed tolerance in COM velocity components.
715 * \param[in] testData Test data structure.
717 void checkCOMVelocity(FloatingPointTolerance tolerance
, const ConstraintsTestData
&testData
)
720 RVec
comV0({0.0, 0.0, 0.0});
721 RVec
comV({0.0, 0.0, 0.0});
722 for (int i
= 0; i
< testData
.numAtoms_
; i
++)
724 comV0
+= testData
.masses_
[i
]*testData
.v0_
[i
];
725 comV
+= testData
.masses_
[i
]*testData
.v_
[i
];
727 comV0
/= testData
.numAtoms_
;
728 comV
/= testData
.numAtoms_
;
730 EXPECT_REAL_EQ_TOL(comV
[XX
], comV0
[XX
], tolerance
)
731 << "Velocity of the center of mass in x-direction has been changed by constraints.";
732 EXPECT_REAL_EQ_TOL(comV
[YY
], comV0
[YY
], tolerance
)
733 << "Velocity of the center of mass in y-direction has been changed by constraints.";
734 EXPECT_REAL_EQ_TOL(comV
[ZZ
], comV0
[ZZ
], tolerance
)
735 << "Velocity of the center of mass in z-direction has been changed by constraints.";
739 * The test on the final coordinates (not used).
741 * Goes through all atoms and checks if the final positions correspond to the
742 * provided reference set of coordinates.
744 * \param[in] xPrimeRef The reference set of coordinates.
745 * \param[in] tolerance Tolerance for the coordinates test.
746 * \param[in] testData Test data structure.
748 void checkFinalCoordinates(std::vector
<RVec
> xPrimeRef
, FloatingPointTolerance tolerance
,
749 const ConstraintsTestData
&testData
)
751 for (int i
= 0; i
< testData
.numAtoms_
; i
++)
753 for (int d
= 0; d
< DIM
; d
++)
755 EXPECT_REAL_EQ_TOL(xPrimeRef
.at(i
)[d
], testData
.xPrime_
[i
][d
], tolerance
) << gmx::formatString(
756 "Coordinates after constrains were applied differ from these in the reference set for atom #%d.", i
);
762 * The test of virial tensor.
764 * Checks if the values in the scaled virial tensor are equal to pre-computed values.
766 * \param[in] tolerance Tolerance for the tensor values.
767 * \param[in] testData Test data structure.
769 void checkVirialTensor(FloatingPointTolerance tolerance
, const ConstraintsTestData
&testData
)
771 for (int i
= 0; i
< DIM
; i
++)
773 for (int j
= 0; j
< DIM
; j
++)
775 EXPECT_REAL_EQ_TOL(testData
.virialScaledRef_
[i
][j
], testData
.virialScaled_
[i
][j
],
776 tolerance
) << gmx::formatString(
777 "Values in virial tensor at [%d][%d] are not within the tolerance from reference value.", i
, j
);
783 * The test for FEP (not used).
785 * Checks if the value of dH/dLambda is equal to the reference value.
786 * \todo Add tests for dHdLambda values.
788 * \param[in] dHdLambdaRef Reference value.
789 * \param[in] tolerance Tolerance.
790 * \param[in] testData Test data structure.
792 void checkFEP(const real dHdLambdaRef
, FloatingPointTolerance tolerance
,
793 const ConstraintsTestData
&testData
)
795 EXPECT_REAL_EQ_TOL(dHdLambdaRef
, testData
.dHdLambda_
, tolerance
) <<
796 "Computed value for dV/dLambda is not equal to the reference value. ";
801 TEST_P(ConstraintsTest
, SingleConstraint
){
802 std::string title
= "one constraint (e.g. OH)";
805 std::vector
<real
> masses
= {1.0, 12.0};
806 std::vector
<int> constraints
= {0, 0, 1};
807 std::vector
<real
> constraintsR0
= {0.1};
809 real oneTenthOverSqrtTwo
= 0.1_real
/ std::sqrt(2.0_real
);
811 std::vector
<RVec
> x
= {{ 0.0, oneTenthOverSqrtTwo
, 0.0 },
812 { oneTenthOverSqrtTwo
, 0.0, 0.0 }};
813 std::vector
<RVec
> xPrime
= {{ 0.01, 0.08, 0.01 },
814 { 0.06, 0.01, -0.01 }};
815 std::vector
<RVec
> v
= {{ 1.0, 2.0, 3.0 },
818 tensor virialScaledRef
= {{-5.58e-04, 5.58e-04, 0.00e+00 },
819 { 5.58e-04, -5.58e-04, 0.00e+00 },
820 { 0.00e+00, 0.00e+00, 0.00e+00 }};
822 real shakeTolerance
= 0.0001;
823 gmx_bool shakeUseSOR
= false;
826 int lincslincsExpansionOrder
= 4;
827 real lincsWarnAngle
= 30.0;
829 std::unique_ptr
<ConstraintsTestData
> testData
= std::make_unique
<ConstraintsTestData
>
830 (title
, numAtoms
, masses
,
831 constraints
, constraintsR0
,
832 true, virialScaledRef
,
834 real(0.0), real(0.001),
836 shakeTolerance
, shakeUseSOR
,
837 lincsNIter
, lincslincsExpansionOrder
, lincsWarnAngle
);
839 std::string algorithmName
;
840 std::tie(pbcName
, algorithmName
) = GetParam();
841 t_pbc pbc
= pbcs_
.at(pbcName
);
844 algorithms_
.at(algorithmName
)(testData
.get(), pbc
);
846 checkConstrainsLength(absoluteTolerance(0.0002), *testData
, pbc
);
847 checkConstrainsDirection(*testData
, pbc
);
848 checkCOMCoordinates(absoluteTolerance(0.0001), *testData
);
849 checkCOMVelocity(absoluteTolerance(0.0001), *testData
);
851 checkVirialTensor(absoluteTolerance(0.0001), *testData
);
855 TEST_P(ConstraintsTest
, TwoDisjointConstraints
){
857 std::string title
= "two disjoint constraints";
859 std::vector
<real
> masses
= {0.5, 1.0/3.0, 0.25, 1.0};
860 std::vector
<int> constraints
= {0, 0, 1, 1, 2, 3};
861 std::vector
<real
> constraintsR0
= {2.0, 1.0};
864 std::vector
<RVec
> x
= {{ 2.50, -3.10, 15.70 },
865 { 0.51, -3.02, 15.55 },
866 { -0.50, -3.00, 15.20 },
867 { -1.51, -2.95, 15.05 }};
869 std::vector
<RVec
> xPrime
= {{ 2.50, -3.10, 15.70 },
870 { 0.51, -3.02, 15.55 },
871 { -0.50, -3.00, 15.20 },
872 { -1.51, -2.95, 15.05 }};
874 std::vector
<RVec
> v
= {{ 0.0, 1.0, 0.0 },
879 tensor virialScaledRef
= {{ 3.3e-03, -1.7e-04, 5.6e-04 },
880 {-1.7e-04, 8.9e-06, -2.8e-05 },
881 { 5.6e-04, -2.8e-05, 8.9e-05 }};
883 real shakeTolerance
= 0.0001;
884 gmx_bool shakeUseSOR
= false;
887 int lincslincsExpansionOrder
= 4;
888 real lincsWarnAngle
= 30.0;
890 std::unique_ptr
<ConstraintsTestData
> testData
= std::make_unique
<ConstraintsTestData
>
891 (title
, numAtoms
, masses
,
892 constraints
, constraintsR0
,
893 true, virialScaledRef
,
895 real(0.0), real(0.001),
897 shakeTolerance
, shakeUseSOR
,
898 lincsNIter
, lincslincsExpansionOrder
, lincsWarnAngle
);
901 std::string algorithmName
;
902 std::tie(pbcName
, algorithmName
) = GetParam();
903 t_pbc pbc
= pbcs_
.at(pbcName
);
906 algorithms_
.at(algorithmName
)(testData
.get(), pbc
);
908 checkConstrainsLength(absoluteTolerance(0.0002), *testData
, pbc
);
909 checkConstrainsDirection(*testData
, pbc
);
910 checkCOMCoordinates(absoluteTolerance(0.0001), *testData
);
911 checkCOMVelocity(absoluteTolerance(0.0001), *testData
);
913 checkVirialTensor(absoluteTolerance(0.0001), *testData
);
917 TEST_P(ConstraintsTest
, ThreeSequentialConstraints
){
919 std::string title
= "three atoms, connected longitudinally (e.g. CH2)";
921 std::vector
<real
> masses
= {1.0, 12.0, 16.0 };
922 std::vector
<int> constraints
= {0, 0, 1, 1, 1, 2};
923 std::vector
<real
> constraintsR0
= {0.1, 0.2};
925 real oneTenthOverSqrtTwo
= 0.1_real
/ std::sqrt(2.0_real
);
926 real twoTenthsOverSqrtThree
= 0.2_real
/ std::sqrt(3.0_real
);
928 std::vector
<RVec
> x
= {{ oneTenthOverSqrtTwo
, oneTenthOverSqrtTwo
, 0.0 },
930 { twoTenthsOverSqrtThree
, twoTenthsOverSqrtThree
, twoTenthsOverSqrtThree
}};
932 std::vector
<RVec
> xPrime
= {{ 0.08, 0.07, 0.01 },
933 { -0.02, 0.01, -0.02 },
934 { 0.10, 0.12, 0.11 }};
936 std::vector
<RVec
> v
= {{ 1.0, 0.0, 0.0 },
940 tensor virialScaledRef
= {{ 4.14e-03, 4.14e-03, 3.31e-03},
941 { 4.14e-03, 4.14e-03, 3.31e-03},
942 { 3.31e-03, 3.31e-03, 3.31e-03}};
944 real shakeTolerance
= 0.0001;
945 gmx_bool shakeUseSOR
= false;
948 int lincslincsExpansionOrder
= 4;
949 real lincsWarnAngle
= 30.0;
951 std::unique_ptr
<ConstraintsTestData
> testData
= std::make_unique
<ConstraintsTestData
>
952 (title
, numAtoms
, masses
,
953 constraints
, constraintsR0
,
954 true, virialScaledRef
,
956 real(0.0), real(0.001),
958 shakeTolerance
, shakeUseSOR
,
959 lincsNIter
, lincslincsExpansionOrder
, lincsWarnAngle
);
962 std::string algorithmName
;
963 std::tie(pbcName
, algorithmName
) = GetParam();
964 t_pbc pbc
= pbcs_
.at(pbcName
);
967 algorithms_
.at(algorithmName
)(testData
.get(), pbc
);
969 checkConstrainsLength(absoluteTolerance(0.0002), *testData
, pbc
);
970 checkConstrainsDirection(*testData
, pbc
);
971 checkCOMCoordinates(absoluteTolerance(0.0001), *testData
);
972 checkCOMVelocity(absoluteTolerance(0.0001), *testData
);
974 checkVirialTensor(absoluteTolerance(0.0001), *testData
);
978 TEST_P(ConstraintsTest
, ThreeConstraintsWithCentralAtom
){
980 std::string title
= "three atoms, connected to the central atom (e.g. CH3)";
982 std::vector
<real
> masses
= {12.0, 1.0, 1.0, 1.0 };
983 std::vector
<int> constraints
= {0, 0, 1, 0, 0, 2, 0, 0, 3};
984 std::vector
<real
> constraintsR0
= {0.1};
987 std::vector
<RVec
> x
= {{ 0.00, 0.00, 0.00 },
988 { 0.10, 0.00, 0.00 },
989 { 0.00, -0.10, 0.00 },
990 { 0.00, 0.00, 0.10 }};
992 std::vector
<RVec
> xPrime
= {{ 0.004, 0.009, -0.010 },
993 { 0.110, -0.006, 0.003 },
994 {-0.007, -0.102, -0.007 },
995 {-0.005, 0.011, 0.102 }};
997 std::vector
<RVec
> v
= {{ 1.0, 0.0, 0.0 },
1002 tensor virialScaledRef
= {{7.14e-04, 0.00e+00, 0.00e+00},
1003 {0.00e+00, 1.08e-03, 0.00e+00},
1004 {0.00e+00, 0.00e+00, 1.15e-03}};
1006 real shakeTolerance
= 0.0001;
1007 gmx_bool shakeUseSOR
= false;
1010 int lincslincsExpansionOrder
= 4;
1011 real lincsWarnAngle
= 30.0;
1013 std::unique_ptr
<ConstraintsTestData
> testData
= std::make_unique
<ConstraintsTestData
>
1014 (title
, numAtoms
, masses
,
1015 constraints
, constraintsR0
,
1016 true, virialScaledRef
,
1018 real(0.0), real(0.001),
1020 shakeTolerance
, shakeUseSOR
,
1021 lincsNIter
, lincslincsExpansionOrder
, lincsWarnAngle
);
1023 std::string pbcName
;
1024 std::string algorithmName
;
1025 std::tie(pbcName
, algorithmName
) = GetParam();
1026 t_pbc pbc
= pbcs_
.at(pbcName
);
1028 // Apply constraints
1029 algorithms_
.at(algorithmName
)(testData
.get(), pbc
);
1031 checkConstrainsLength(absoluteTolerance(0.0002), *testData
, pbc
);
1032 checkConstrainsDirection(*testData
, pbc
);
1033 checkCOMCoordinates(absoluteTolerance(0.0001), *testData
);
1034 checkCOMVelocity(absoluteTolerance(0.0001), *testData
);
1036 checkVirialTensor(absoluteTolerance(0.0001), *testData
);
1039 TEST_P(ConstraintsTest
, FourSequentialConstraints
){
1041 std::string title
= "four atoms, connected longitudinally";
1043 std::vector
<real
> masses
= {0.5, 1.0/3.0, 0.25, 1.0};
1044 std::vector
<int> constraints
= {0, 0, 1, 1, 1, 2, 2, 2, 3};
1045 std::vector
<real
> constraintsR0
= {2.0, 1.0, 1.0};
1048 std::vector
<RVec
> x
= {{ 2.50, -3.10, 15.70 },
1049 { 0.51, -3.02, 15.55 },
1050 { -0.50, -3.00, 15.20 },
1051 { -1.51, -2.95, 15.05 }};
1053 std::vector
<RVec
> xPrime
= {{ 2.50, -3.10, 15.70 },
1054 { 0.51, -3.02, 15.55 },
1055 { -0.50, -3.00, 15.20 },
1056 { -1.51, -2.95, 15.05 }};
1058 std::vector
<RVec
> v
= {{ 0.0, 0.0, 2.0 },
1061 { 0.0, 0.0, -1.0 }};
1063 tensor virialScaledRef
= {{ 1.15e-01, -4.20e-03, 2.12e-02},
1064 {-4.20e-03, 1.70e-04, -6.41e-04},
1065 { 2.12e-02, -6.41e-04, 5.45e-03}};
1067 real shakeTolerance
= 0.0001;
1068 gmx_bool shakeUseSOR
= false;
1071 int lincslincsExpansionOrder
= 8;
1072 real lincsWarnAngle
= 30.0;
1074 std::unique_ptr
<ConstraintsTestData
> testData
= std::make_unique
<ConstraintsTestData
>
1075 (title
, numAtoms
, masses
,
1076 constraints
, constraintsR0
,
1077 true, virialScaledRef
,
1079 real(0.0), real(0.001),
1081 shakeTolerance
, shakeUseSOR
,
1082 lincsNIter
, lincslincsExpansionOrder
, lincsWarnAngle
);
1084 std::string pbcName
;
1085 std::string algorithmName
;
1086 std::tie(pbcName
, algorithmName
) = GetParam();
1087 t_pbc pbc
= pbcs_
.at(pbcName
);
1089 // Apply constraints
1090 algorithms_
.at(algorithmName
)(testData
.get(), pbc
);
1092 checkConstrainsLength(absoluteTolerance(0.0002), *testData
, pbc
);
1093 checkConstrainsDirection(*testData
, pbc
);
1094 checkCOMCoordinates(absoluteTolerance(0.0001), *testData
);
1095 checkCOMVelocity(absoluteTolerance(0.0001), *testData
);
1097 checkVirialTensor(absoluteTolerance(0.01), *testData
);
1101 TEST_P(ConstraintsTest
, TriangleOfConstraints
){
1103 std::string title
= "basic triangle (tree atoms, connected to each other)";
1105 std::vector
<real
> masses
= {1.0, 1.0, 1.0};
1106 std::vector
<int> constraints
= {0, 0, 1, 2, 0, 2, 1, 1, 2};
1107 std::vector
<real
> constraintsR0
= {0.1, 0.1, 0.1};
1109 real oneTenthOverSqrtTwo
= 0.1_real
/ std::sqrt(2.0_real
);
1111 std::vector
<RVec
> x
= {{ oneTenthOverSqrtTwo
, 0.0, 0.0 },
1112 { 0.0, oneTenthOverSqrtTwo
, 0.0 },
1113 { 0.0, 0.0, oneTenthOverSqrtTwo
}};
1115 std::vector
<RVec
> xPrime
= {{ 0.09, -0.02, 0.01 },
1116 { -0.02, 0.10, -0.02 },
1117 { 0.03, -0.01, 0.07 }};
1119 std::vector
<RVec
> v
= {{ 1.0, 1.0, 1.0 },
1120 { -2.0, -2.0, -2.0 },
1123 tensor virialScaledRef
= {{ 6.00e-04, -1.61e-03, 1.01e-03},
1124 {-1.61e-03, 2.53e-03, -9.25e-04},
1125 { 1.01e-03, -9.25e-04, -8.05e-05}};
1127 real shakeTolerance
= 0.0001;
1128 gmx_bool shakeUseSOR
= false;
1131 int lincslincsExpansionOrder
= 4;
1132 real lincsWarnAngle
= 30.0;
1134 std::unique_ptr
<ConstraintsTestData
> testData
= std::make_unique
<ConstraintsTestData
>
1135 (title
, numAtoms
, masses
,
1136 constraints
, constraintsR0
,
1137 true, virialScaledRef
,
1139 real(0.0), real(0.001),
1141 shakeTolerance
, shakeUseSOR
,
1142 lincsNIter
, lincslincsExpansionOrder
, lincsWarnAngle
);
1144 std::string pbcName
;
1145 std::string algorithmName
;
1146 std::tie(pbcName
, algorithmName
) = GetParam();
1147 t_pbc pbc
= pbcs_
.at(pbcName
);
1149 // Apply constraints
1150 algorithms_
.at(algorithmName
)(testData
.get(), pbc
);
1152 checkConstrainsLength(absoluteTolerance(0.0002), *testData
, pbc
);
1153 checkConstrainsDirection(*testData
, pbc
);
1154 checkCOMCoordinates(absoluteTolerance(0.0001), *testData
);
1155 checkCOMVelocity(absoluteTolerance(0.0001), *testData
);
1157 checkVirialTensor(absoluteTolerance(0.00001), *testData
);
1162 #if GMX_GPU == GMX_GPU_CUDA
1163 INSTANTIATE_TEST_CASE_P(WithParameters
, ConstraintsTest
,
1164 ::testing::Combine(::testing::Values("PBCNone", "PBCXYZ"),
1165 ::testing::ValuesIn(algorithmsNames
)));
1168 #if GMX_GPU != GMX_GPU_CUDA
1169 INSTANTIATE_TEST_CASE_P(WithParameters
, ConstraintsTest
,
1170 ::testing::Combine(::testing::Values("PBCNone", "PBCXYZ"),
1171 ::testing::Values("SHAKE", "LINCS")));