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.
37 * Implements test of bonded force routines
39 * \author David van der Spoel <david.vanderspoel@icm.uu.se>
40 * \ingroup module_listed_forces
44 #include "gromacs/listed_forces/bonded.h"
49 #include <unordered_map>
51 #include <gtest/gtest.h>
53 #include "gromacs/listed_forces/listed_forces.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/math/vectypes.h"
57 #include "gromacs/pbcutil/ishift.h"
58 #include "gromacs/pbcutil/pbc.h"
59 #include "gromacs/topology/idef.h"
60 #include "gromacs/utility/strconvert.h"
62 #include "testutils/refdata.h"
63 #include "testutils/testasserts.h"
70 //! Number of atoms used in these tests.
71 constexpr int c_numAtoms
= 4;
73 /*! \brief Output from bonded kernels
75 * \todo Later this might turn into the actual output struct. */
76 struct OutputQuantities
78 //! Energy of this interaction
80 //! Derivative with respect to lambda
83 rvec fshift
[N_IVEC
] = {{0}};
85 rvec4 f
[c_numAtoms
] = {{0}};
88 /*! \brief Utility to check the output from bonded tests
90 * \param[in] checker Reference checker
91 * \param[in] output The output from the test to check
93 void checkOutput(test::TestReferenceChecker
*checker
,
94 const OutputQuantities
&output
)
96 checker
->checkReal(output
.energy
, "Epot ");
97 // Should still be zero if not doing FEP, so may as well test it.
98 checker
->checkReal(output
.dvdlambda
, "dVdlambda ");
99 checker
->checkVector(output
.fshift
[CENTRAL
], "Central shift forces");
100 checker
->checkSequence(std::begin(output
.f
), std::end(output
.f
), "Forces");
103 /*! \brief Input structure for listed forces tests
110 //! Tolerance for float evaluation
112 //! Tolerance for double evaluation
113 double dtoler
= 1e-8;
114 //! Do free energy perturbation?
116 //! Interaction parameters
117 t_iparams iparams
= {{ 0 }};
122 /*! \brief Constructor with tolerance
124 * \param[in] ftol Single precision tolerance
125 * \param[in] dtol Double precision tolerance
127 iListInput(float ftol
, double dtol
)
132 /*! \brief Set parameters for harmonic potential
134 * Free energy perturbation is turned on when A
135 * and B parameters are different.
136 * \param[in] ft Function type
137 * \param[in] rA Equilibrium value A
138 * \param[in] krA Force constant A
139 * \param[in] rB Equilibrium value B
140 * \param[in] krB Force constant B
141 * \return The structure itself.
143 iListInput
setHarmonic(int ft
, real rA
, real krA
, real rB
, real krB
)
145 iparams
.harmonic
.rA
= rA
;
146 iparams
.harmonic
.rB
= rB
;
147 iparams
.harmonic
.krA
= krA
;
148 iparams
.harmonic
.krB
= krB
;
150 fep
= (rA
!= rB
|| krA
!= krB
);
153 /*! \brief Set parameters for harmonic potential
155 * \param[in] ft Function type
156 * \param[in] rA Equilibrium value
157 * \param[in] krA Force constant
158 * \return The structure itself.
160 iListInput
setHarmonic(int ft
, real rA
, real krA
)
162 return setHarmonic(ft
, rA
, krA
, rA
, krA
);
164 /*! \brief Set parameters for cubic potential
166 * \param[in] b0 Equilibrium bond length
167 * \param[in] kb Harmonic force constant
168 * \param[in] kcub Cubic force constant
169 * \return The structure itself.
171 iListInput
setCubic(real b0
, real kb
, real kcub
)
173 ftype
= F_CUBICBONDS
;
174 iparams
.cubic
.b0
= b0
;
175 iparams
.cubic
.kb
= kb
;
176 iparams
.cubic
.kcub
= kcub
;
179 /*! \brief Set parameters for morse potential
181 * Free energy perturbation is turned on when A
182 * and B parameters are different.
183 * \param[in] b0A Equilibrium value A
184 * \param[in] cbA Force constant A
185 * \param[in] betaA Steepness parameter A
186 * \param[in] b0B Equilibrium value B
187 * \param[in] cbB Force constant B
188 * \param[in] betaB Steepness parameter B
189 * \return The structure itself.
191 iListInput
setMorse(real b0A
, real cbA
, real betaA
, real b0B
, real cbB
, real betaB
)
194 iparams
.morse
.b0A
= b0A
;
195 iparams
.morse
.cbA
= cbA
;
196 iparams
.morse
.betaA
= betaA
;
197 iparams
.morse
.b0B
= b0B
;
198 iparams
.morse
.cbB
= cbB
;
199 iparams
.morse
.betaB
= betaB
;
200 fep
= (b0A
!= b0B
|| cbA
!= cbB
|| betaA
!= betaB
);
203 /*! \brief Set parameters for morse potential
205 * \param[in] b0A Equilibrium value
206 * \param[in] cbA Force constant
207 * \param[in] betaA Steepness parameter
208 * \return The structure itself.
210 iListInput
setMorse(real b0A
, real cbA
, real betaA
)
212 return setMorse(b0A
, cbA
, betaA
, b0A
, cbA
, betaA
);
214 /*! \brief Set parameters for fene potential
216 * \param[in] bm Equilibrium bond length
217 * \param[in] kb Force constant
218 * \return The structure itself.
220 iListInput
setFene(real bm
, real kb
)
223 iparams
.fene
.bm
= bm
;
224 iparams
.fene
.kb
= kb
;
227 /*! \brief Set parameters for linear angle potential
229 * Free energy perturbation is turned on when A
230 * and B parameters are different.
231 * \param[in] klinA Force constant A
232 * \param[in] aA The position of the central atom A
233 * \param[in] klinB Force constant B
234 * \param[in] aB The position of the central atom B
235 * \return The structure itself.
237 iListInput
setLinearAngle(real klinA
, real aA
, real klinB
, real aB
)
239 ftype
= F_LINEAR_ANGLES
;
240 iparams
.linangle
.klinA
= klinA
;
241 iparams
.linangle
.aA
= aA
;
242 iparams
.linangle
.klinB
= klinB
;
243 iparams
.linangle
.aB
= aB
;
244 fep
= (klinA
!= klinB
|| aA
!= aB
);
247 /*! \brief Set parameters for linear angle potential
249 * \param[in] klinA Force constant
250 * \param[in] aA The position of the central atom
251 * \return The structure itself.
253 iListInput
setLinearAngle(real klinA
, real aA
)
255 return setLinearAngle(klinA
, aA
, klinA
, aA
);
257 /*! \brief Set parameters for Urey Bradley potential
259 * Free energy perturbation is turned on when A
260 * and B parameters are different.
261 * \param[in] thetaA Equilibrium angle A
262 * \param[in] kthetaA Force constant A
263 * \param[in] r13A The distance between i and k atoms A
264 * \param[in] kUBA The force constant for 1-3 distance A
265 * \param[in] thetaB Equilibrium angle B
266 * \param[in] kthetaB Force constant B
267 * \param[in] r13B The distance between i and k atoms B
268 * \param[in] kUBB The force constant for 1-3 distance B
269 * \return The structure itself.
271 iListInput
setUreyBradley(real thetaA
, real kthetaA
, real r13A
, real kUBA
,
272 real thetaB
, real kthetaB
, real r13B
, real kUBB
)
274 ftype
= F_UREY_BRADLEY
;
275 iparams
.u_b
.thetaA
= thetaA
;
276 iparams
.u_b
.kthetaA
= kthetaA
;
277 iparams
.u_b
.r13A
= r13A
;
278 iparams
.u_b
.kUBA
= kUBA
;
279 iparams
.u_b
.thetaB
= thetaB
;
280 iparams
.u_b
.kthetaB
= kthetaB
;
281 iparams
.u_b
.r13B
= r13B
;
282 iparams
.u_b
.kUBB
= kUBB
;
283 fep
= (thetaA
!= thetaB
|| kthetaA
!= kthetaB
|| r13A
!= r13B
|| kUBA
!= kUBB
);
286 /*! \brief Set parameters for Urey Bradley potential
288 * \param[in] thetaA Equilibrium angle
289 * \param[in] kthetaA Force constant
290 * \param[in] r13A The distance between i and k atoms
291 * \param[in] kUBA The force constant for 1-3 distance
292 * \return The structure itself.
294 iListInput
setUreyBradley(real thetaA
, real kthetaA
, real r13A
, real kUBA
)
296 return setUreyBradley(thetaA
, kthetaA
, r13A
, kUBA
,
297 thetaA
, kthetaA
, r13A
, kUBA
);
299 /*! \brief Set parameters for Cross Bond Bonds potential
301 * \param[in] r1e First bond length i-j
302 * \param[in] r2e Second bond length i-k
303 * \param[in] krr The force constant
304 * \return The structure itself.
306 iListInput
setCrossBondBonds(real r1e
, real r2e
, real krr
)
308 ftype
= F_CROSS_BOND_BONDS
;
309 iparams
.cross_bb
.r1e
= r1e
;
310 iparams
.cross_bb
.r2e
= r2e
;
311 iparams
.cross_bb
.krr
= krr
;
314 /*! \brief Set parameters for Cross Bond Angles potential
316 * \param[in] r1e First bond length i-j
317 * \param[in] r2e Second bond length j-k
318 * \param[in] r3e Third bond length i-k
319 * \param[in] krt The force constant
320 * \return The structure itself.
322 iListInput
setCrossBondAngles(real r1e
, real r2e
, real r3e
, real krt
)
324 ftype
= F_CROSS_BOND_ANGLES
;
325 iparams
.cross_ba
.r1e
= r1e
;
326 iparams
.cross_ba
.r2e
= r2e
;
327 iparams
.cross_ba
.r3e
= r3e
;
328 iparams
.cross_ba
.krt
= krt
;
331 /*! \brief Set parameters for Quartic Angles potential
333 * \param[in] theta Angle
334 * \param[in] c Array of parameters
335 * \return The structure itself.
337 iListInput
setQuarticAngles(real theta
, const real c
[5])
339 ftype
= F_QUARTIC_ANGLES
;
340 iparams
.qangle
.theta
= theta
;
341 iparams
.qangle
.c
[0] = c
[0];
342 iparams
.qangle
.c
[1] = c
[1];
343 iparams
.qangle
.c
[2] = c
[2];
344 iparams
.qangle
.c
[3] = c
[3];
345 iparams
.qangle
.c
[4] = c
[4];
348 /*! \brief Set parameters for proper dihedrals potential
350 * Free energy perturbation is turned on when A
351 * and B parameters are different.
352 * \param[in] phiA Dihedral angle A
353 * \param[in] cpA Force constant A
354 * \param[in] mult Multiplicity of the angle
355 * \param[in] phiB Dihedral angle B
356 * \param[in] cpB Force constant B
357 * \return The structure itself.
359 iListInput
setProperDihedrals(real phiA
, real cpA
, int mult
, real phiB
, real cpB
)
362 iparams
.pdihs
.phiA
= phiA
;
363 iparams
.pdihs
.cpA
= cpA
;
364 iparams
.pdihs
.phiB
= phiB
;
365 iparams
.pdihs
.cpB
= cpB
;
366 iparams
.pdihs
.mult
= mult
;
367 fep
= (phiA
!= phiB
|| cpA
!= cpB
);
370 /*! \brief Set parameters for proper dihedrals potential
372 * \param[in] phiA Dihedral angle
373 * \param[in] cpA Force constant
374 * \param[in] mult Multiplicity of the angle
375 * \return The structure itself.
377 iListInput
setProperDihedrals(real phiA
, real cpA
, int mult
)
379 return setProperDihedrals(phiA
, cpA
, mult
, phiA
, cpA
);
381 /*! \brief Set parameters for Ryckaert-Bellemans dihedrals potential
383 * Free energy perturbation is turned on when A
384 * and B parameters are different.
385 * \param[in] rbcA Force constants A
386 * \param[in] rbcB Force constants B
387 * \return The structure itself.
389 iListInput
setRbDihedrals(const real rbcA
[NR_RBDIHS
], const real rbcB
[NR_RBDIHS
])
393 for (int i
= 0; i
< NR_RBDIHS
; i
++)
395 iparams
.rbdihs
.rbcA
[i
] = rbcA
[i
];
396 iparams
.rbdihs
.rbcB
[i
] = rbcB
[i
];
397 fep
= fep
|| (rbcA
[i
] != rbcB
[i
]);
401 /*! \brief Set parameters for Ryckaert-Bellemans dihedrals potential
403 * \param[in] rbc Force constants
404 * \return The structure itself.
406 iListInput
setRbDihedrals(const real rbc
[NR_RBDIHS
])
408 return setRbDihedrals(rbc
, rbc
);
412 /*! \brief Utility to fill iatoms struct
414 * \param[in] ftype Function type
415 * \param[out] iatoms Pointer to iatoms struct
417 void fillIatoms(int ftype
, std::vector
<t_iatom
> *iatoms
)
419 std::unordered_map
<int, std::vector
<int> > ia
=
420 { { 2, { 0, 0, 1, 0, 1, 2, 0, 2, 3 } },
421 { 3, { 0, 0, 1, 2, 0, 1, 2, 3 } },
422 { 4, { 0, 0, 1, 2, 3 } }};
423 EXPECT_TRUE(ftype
>= 0 && ftype
< F_NRE
);
424 int nral
= interaction_function
[ftype
].nratoms
;
425 for (auto &i
: ia
[nral
])
427 iatoms
->push_back(i
);
431 class ListedForcesTest
: public ::testing::TestWithParam
<std::tuple
<iListInput
, std::vector
<gmx::RVec
>, int> >
436 std::vector
<gmx::RVec
> x_
;
439 test::TestReferenceData refData_
;
440 test::TestReferenceChecker checker_
;
441 ListedForcesTest( ) :
442 checker_(refData_
.rootChecker())
444 input_
= std::get
<0>(GetParam());
445 x_
= std::get
<1>(GetParam());
446 epbc_
= std::get
<2>(GetParam());
448 box_
[0][0] = box_
[1][1] = box_
[2][2] = 1.5;
449 set_pbc(&pbc_
, epbc_
, box_
);
451 // We need quite specific tolerances here since angle functions
452 // etc. are not very precise and reproducible.
453 test::FloatingPointTolerance
tolerance(test::FloatingPointTolerance(input_
.ftoler
, 1.0e-6,
454 input_
.dtoler
, 1.0e-12,
456 checker_
.setDefaultTolerance(tolerance
);
458 void testOneIfunc(test::TestReferenceChecker
*checker
,
459 const std::vector
<t_iatom
> &iatoms
,
462 SCOPED_TRACE(std::string("Testing PBC ") + epbc_names
[epbc_
]);
463 std::vector
<int> ddgatindex
= { 0, 1, 2, 3 };
464 OutputQuantities output
;
465 output
.energy
= bondedFunction(input_
.ftype
) (iatoms
.size(),
468 as_rvec_array(x_
.data()),
469 output
.f
, output
.fshift
,
471 /* const struct t_graph *g */ nullptr,
472 lambda
, &output
.dvdlambda
,
473 /* const struct t_mdatoms *md */ nullptr,
474 /* struct t_fcdata *fcd */ nullptr,
476 // Internal consistency test of both test input
477 // and bonded functions.
478 EXPECT_TRUE((input_
.fep
|| (output
.dvdlambda
== 0.0)));
479 checkOutput(checker
, output
);
483 test::TestReferenceChecker thisChecker
=
484 checker_
.checkCompound("FunctionType",
485 interaction_function
[input_
.ftype
].name
).checkCompound("FEP", (input_
.fep
? "Yes" : "No"));
486 std::vector
<t_iatom
> iatoms
;
487 fillIatoms(input_
.ftype
, &iatoms
);
490 const int numLambdas
= 3;
491 for (int i
= 0; i
< numLambdas
; ++i
)
493 const real lambda
= i
/ (numLambdas
- 1.0);
494 auto valueChecker
= thisChecker
.checkCompound("Lambda", toString(lambda
));
495 testOneIfunc(&valueChecker
, iatoms
, lambda
);
500 testOneIfunc(&thisChecker
, iatoms
, 0.0);
505 TEST_P (ListedForcesTest
, Ifunc
)
510 //! Function types for testing bonds. Add new terms at the end.
511 std::vector
<iListInput
> c_InputBonds
=
513 { iListInput().setHarmonic(F_BONDS
, 0.15, 500.0) },
514 { iListInput(2e-6f
, 1e-8).setHarmonic(F_BONDS
, 0.15, 500.0, 0.17, 400.0) },
515 { iListInput(1e-4f
, 1e-8).setHarmonic(F_G96BONDS
, 0.15, 50.0) },
516 { iListInput().setHarmonic(F_G96BONDS
, 0.15, 50.0, 0.17, 40.0) },
517 { iListInput().setCubic(0.16, 50.0, 2.0) },
518 { iListInput(2e-6f
, 1e-8).setMorse(0.15, 50.0, 2.0, 0.17, 40.0, 1.6) },
519 { iListInput(2e-6f
, 1e-8).setMorse(0.15, 30.0, 2.7) },
520 { iListInput().setFene(0.4, 5.0) }
523 //! Constants for Quartic Angles
524 const real cQuarticAngles
[5] = { 1.1, 2.3, 4.6, 7.8, 9.2 };
526 //! Function types for testing angles. Add new terms at the end.
527 std::vector
<iListInput
> c_InputAngles
=
529 { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES
, 100.0, 50.0) },
530 { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES
, 100.15, 50.0, 95.0, 30.0) },
531 { iListInput(8e-3, 1e-8).setHarmonic(F_G96ANGLES
, 100.0, 50.0) },
532 { iListInput(8e-3, 1e-8).setHarmonic(F_G96ANGLES
, 100.0, 50.0, 95.0, 30.0) },
533 { iListInput().setLinearAngle(50.0, 0.4) },
534 { iListInput().setLinearAngle(50.0, 0.4, 40.0, 0.6) },
535 { iListInput(2e-6, 1e-8).setCrossBondBonds(0.8, 0.7, 45.0) },
536 { iListInput(2e-6, 1e-8).setCrossBondAngles(0.8, 0.7, 0.3, 45.0) },
537 { iListInput(8e-3, 1e-8).setUreyBradley(950.0, 46.0, 0.3, 5.0) },
538 { iListInput(8e-3, 1e-8).setUreyBradley(100.0, 45.0, 0.3, 5.0, 90.0, 47.0, 0.32, 7.0) },
539 { iListInput(7e-4, 1e-8).setQuarticAngles(87.0, cQuarticAngles
) }
542 //! Constants for Ryckaert-Bellemans A
543 const real rbcA
[NR_RBDIHS
] = { -5.35, 13.6, 8.4, -16.7, 0.3, 12.4 };
545 //! Constants for Ryckaert-Bellemans B
546 const real rbcB
[NR_RBDIHS
] = { -6.35, 12.6, 8.1, -10.7, 0.9, 15.4 };
548 //! Constants for Ryckaert-Bellemans without FEP
549 const real rbc
[NR_RBDIHS
] = { -7.35, 13.6, 8.4, -16.7, 1.3, 12.4 };
551 //! Function types for testing dihedrals. Add new terms at the end.
552 std::vector
<iListInput
> c_InputDihs
=
554 { iListInput(1e-4, 1e-8).setProperDihedrals(-100.0, 10.0, 2, -80.0, 20.0) },
555 { iListInput(1e-4, 1e-8).setProperDihedrals(-105.0, 15.0, 2) },
556 { iListInput(2e-4, 1e-8).setHarmonic(F_IDIHS
, 100.0, 50.0) },
557 { iListInput(2e-4, 1e-8).setHarmonic(F_IDIHS
, 100.15, 50.0, 95.0, 30.0) },
558 { iListInput(3e-4, 1e-8).setRbDihedrals(rbcA
, rbcB
) },
559 { iListInput(3e-4, 1e-8).setRbDihedrals(rbc
) }
562 //! Coordinates for testing
563 std::vector
<std::vector
<gmx::RVec
> > c_coordinatesForTests
=
565 {{ 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.2 }, { 0.005, 0.0, 0.1 }, { -0.001, 0.1, 0.0 }},
566 {{ 0.5, 0.0, 0.0 }, { 0.5, 0.0, 0.15 }, { 0.5, 0.07, 0.22 }, { 0.5, 0.18, 0.22 }},
567 {{ -0.1143, -0.0282, 0.0 }, { 0.0, 0.0434, 0.0 }, { 0.1185, -0.0138, 0.0 }, { -0.0195, 0.1498, 0.0 }}
570 //! PBC values for testing
571 std::vector
<int> c_pbcForTests
= { epbcNONE
, epbcXY
, epbcXYZ
};
573 INSTANTIATE_TEST_CASE_P(Bond
, ListedForcesTest
, ::testing::Combine(::testing::ValuesIn(c_InputBonds
), ::testing::ValuesIn(c_coordinatesForTests
), ::testing::ValuesIn(c_pbcForTests
)));
575 INSTANTIATE_TEST_CASE_P(Angle
, ListedForcesTest
, ::testing::Combine(::testing::ValuesIn(c_InputAngles
), ::testing::ValuesIn(c_coordinatesForTests
), ::testing::ValuesIn(c_pbcForTests
)));
577 INSTANTIATE_TEST_CASE_P(Dihedral
, ListedForcesTest
, ::testing::Combine(::testing::ValuesIn(c_InputDihs
), ::testing::ValuesIn(c_coordinatesForTests
), ::testing::ValuesIn(c_pbcForTests
)));