2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2018,2019,2020, 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/paddedvector.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/math/vectypes.h"
58 #include "gromacs/mdtypes/mdatom.h"
59 #include "gromacs/pbcutil/ishift.h"
60 #include "gromacs/pbcutil/pbc.h"
61 #include "gromacs/topology/idef.h"
62 #include "gromacs/utility/strconvert.h"
63 #include "gromacs/utility/stringstream.h"
64 #include "gromacs/utility/textwriter.h"
66 #include "testutils/refdata.h"
67 #include "testutils/testasserts.h"
74 //! Number of atoms used in these tests.
75 constexpr int c_numAtoms = 4;
77 /*! \brief Output from bonded kernels
79 * \todo Later this might turn into the actual output struct. */
80 struct OutputQuantities
82 //! Energy of this interaction
84 //! Derivative with respect to lambda
87 rvec fshift[N_IVEC] = { { 0 } };
89 alignas(GMX_REAL_MAX_SIMD_WIDTH * sizeof(real)) rvec4 f[c_numAtoms] = { { 0 } };
92 /*! \brief Utility to check the output from bonded tests
94 * \param[in] checker Reference checker
95 * \param[in] output The output from the test to check
96 * \param[in] bondedKernelFlavor Flavor for determining what output to check
98 void checkOutput(test::TestReferenceChecker* checker,
99 const OutputQuantities& output,
100 const BondedKernelFlavor bondedKernelFlavor)
102 if (computeEnergy(bondedKernelFlavor))
104 checker->checkReal(output.energy, "Epot ");
105 // Should still be zero when not doing FEP, so may as well test it.
106 checker->checkReal(output.dvdlambda, "dVdlambda ");
108 if (computeVirial(bondedKernelFlavor))
110 checker->checkVector(output.fshift[CENTRAL], "Central shift forces");
112 checker->checkSequence(std::begin(output.f), std::end(output.f), "Forces");
115 /*! \brief Input structure for listed forces tests
122 //! Tolerance for float evaluation
124 //! Tolerance for double evaluation
125 double dtoler = 1e-8;
126 //! Do free energy perturbation?
128 //! Interaction parameters
129 t_iparams iparams = { { 0 } };
131 friend std::ostream& operator<<(std::ostream& out, const iListInput& input);
136 /*! \brief Constructor with tolerance
138 * \param[in] ftol Single precision tolerance
139 * \param[in] dtol Double precision tolerance
141 iListInput(float ftol, double dtol)
146 /*! \brief Set parameters for harmonic potential
148 * Free energy perturbation is turned on when A
149 * and B parameters are different.
150 * \param[in] ft Function type
151 * \param[in] rA Equilibrium value A
152 * \param[in] krA Force constant A
153 * \param[in] rB Equilibrium value B
154 * \param[in] krB Force constant B
155 * \return The structure itself.
157 iListInput setHarmonic(int ft, real rA, real krA, real rB, real krB)
159 iparams.harmonic.rA = rA;
160 iparams.harmonic.rB = rB;
161 iparams.harmonic.krA = krA;
162 iparams.harmonic.krB = krB;
164 fep = (rA != rB || krA != krB);
167 /*! \brief Set parameters for harmonic potential
169 * \param[in] ft Function type
170 * \param[in] rA Equilibrium value
171 * \param[in] krA Force constant
172 * \return The structure itself.
174 iListInput setHarmonic(int ft, real rA, real krA) { return setHarmonic(ft, rA, krA, rA, krA); }
175 /*! \brief Set parameters for cubic potential
177 * \param[in] b0 Equilibrium bond length
178 * \param[in] kb Harmonic force constant
179 * \param[in] kcub Cubic force constant
180 * \return The structure itself.
182 iListInput setCubic(real b0, real kb, real kcub)
184 ftype = F_CUBICBONDS;
185 iparams.cubic.b0 = b0;
186 iparams.cubic.kb = kb;
187 iparams.cubic.kcub = kcub;
190 /*! \brief Set parameters for morse potential
192 * Free energy perturbation is turned on when A
193 * and B parameters are different.
194 * \param[in] b0A Equilibrium value A
195 * \param[in] cbA Force constant A
196 * \param[in] betaA Steepness parameter A
197 * \param[in] b0B Equilibrium value B
198 * \param[in] cbB Force constant B
199 * \param[in] betaB Steepness parameter B
200 * \return The structure itself.
202 iListInput setMorse(real b0A, real cbA, real betaA, real b0B, real cbB, real betaB)
205 iparams.morse.b0A = b0A;
206 iparams.morse.cbA = cbA;
207 iparams.morse.betaA = betaA;
208 iparams.morse.b0B = b0B;
209 iparams.morse.cbB = cbB;
210 iparams.morse.betaB = betaB;
211 fep = (b0A != b0B || cbA != cbB || betaA != betaB);
214 /*! \brief Set parameters for morse potential
216 * \param[in] b0A Equilibrium value
217 * \param[in] cbA Force constant
218 * \param[in] betaA Steepness parameter
219 * \return The structure itself.
221 iListInput setMorse(real b0A, real cbA, real betaA)
223 return setMorse(b0A, cbA, betaA, b0A, cbA, betaA);
225 /*! \brief Set parameters for fene potential
227 * \param[in] bm Equilibrium bond length
228 * \param[in] kb Force constant
229 * \return The structure itself.
231 iListInput setFene(real bm, real kb)
234 iparams.fene.bm = bm;
235 iparams.fene.kb = kb;
238 /*! \brief Set parameters for linear angle potential
240 * Free energy perturbation is turned on when A
241 * and B parameters are different.
242 * \param[in] klinA Force constant A
243 * \param[in] aA The position of the central atom A
244 * \param[in] klinB Force constant B
245 * \param[in] aB The position of the central atom B
246 * \return The structure itself.
248 iListInput setLinearAngle(real klinA, real aA, real klinB, real aB)
250 ftype = F_LINEAR_ANGLES;
251 iparams.linangle.klinA = klinA;
252 iparams.linangle.aA = aA;
253 iparams.linangle.klinB = klinB;
254 iparams.linangle.aB = aB;
255 fep = (klinA != klinB || aA != aB);
258 /*! \brief Set parameters for linear angle potential
260 * \param[in] klinA Force constant
261 * \param[in] aA The position of the central atom
262 * \return The structure itself.
264 iListInput setLinearAngle(real klinA, real aA) { return setLinearAngle(klinA, aA, klinA, aA); }
265 /*! \brief Set parameters for Urey Bradley potential
267 * Free energy perturbation is turned on when A
268 * and B parameters are different.
269 * \param[in] thetaA Equilibrium angle A
270 * \param[in] kthetaA Force constant A
271 * \param[in] r13A The distance between i and k atoms A
272 * \param[in] kUBA The force constant for 1-3 distance A
273 * \param[in] thetaB Equilibrium angle B
274 * \param[in] kthetaB Force constant B
275 * \param[in] r13B The distance between i and k atoms B
276 * \param[in] kUBB The force constant for 1-3 distance B
277 * \return The structure itself.
280 setUreyBradley(real thetaA, real kthetaA, real r13A, real kUBA, real thetaB, real kthetaB, real r13B, real kUBB)
282 ftype = F_UREY_BRADLEY;
283 iparams.u_b.thetaA = thetaA;
284 iparams.u_b.kthetaA = kthetaA;
285 iparams.u_b.r13A = r13A;
286 iparams.u_b.kUBA = kUBA;
287 iparams.u_b.thetaB = thetaB;
288 iparams.u_b.kthetaB = kthetaB;
289 iparams.u_b.r13B = r13B;
290 iparams.u_b.kUBB = kUBB;
291 fep = (thetaA != thetaB || kthetaA != kthetaB || r13A != r13B || kUBA != kUBB);
294 /*! \brief Set parameters for Urey Bradley potential
296 * \param[in] thetaA Equilibrium angle
297 * \param[in] kthetaA Force constant
298 * \param[in] r13A The distance between i and k atoms
299 * \param[in] kUBA The force constant for 1-3 distance
300 * \return The structure itself.
302 iListInput setUreyBradley(real thetaA, real kthetaA, real r13A, real kUBA)
304 return setUreyBradley(thetaA, kthetaA, r13A, kUBA, thetaA, kthetaA, r13A, kUBA);
306 /*! \brief Set parameters for Cross Bond Bonds potential
308 * \param[in] r1e First bond length i-j
309 * \param[in] r2e Second bond length i-k
310 * \param[in] krr The force constant
311 * \return The structure itself.
313 iListInput setCrossBondBonds(real r1e, real r2e, real krr)
315 ftype = F_CROSS_BOND_BONDS;
316 iparams.cross_bb.r1e = r1e;
317 iparams.cross_bb.r2e = r2e;
318 iparams.cross_bb.krr = krr;
321 /*! \brief Set parameters for Cross Bond Angles potential
323 * \param[in] r1e First bond length i-j
324 * \param[in] r2e Second bond length j-k
325 * \param[in] r3e Third bond length i-k
326 * \param[in] krt The force constant
327 * \return The structure itself.
329 iListInput setCrossBondAngles(real r1e, real r2e, real r3e, real krt)
331 ftype = F_CROSS_BOND_ANGLES;
332 iparams.cross_ba.r1e = r1e;
333 iparams.cross_ba.r2e = r2e;
334 iparams.cross_ba.r3e = r3e;
335 iparams.cross_ba.krt = krt;
338 /*! \brief Set parameters for Quartic Angles potential
340 * \param[in] theta Angle
341 * \param[in] c Array of parameters
342 * \return The structure itself.
344 iListInput setQuarticAngles(real theta, const real c[5])
346 ftype = F_QUARTIC_ANGLES;
347 iparams.qangle.theta = theta;
348 iparams.qangle.c[0] = c[0];
349 iparams.qangle.c[1] = c[1];
350 iparams.qangle.c[2] = c[2];
351 iparams.qangle.c[3] = c[3];
352 iparams.qangle.c[4] = c[4];
355 /*! \brief Set parameters for proper dihedrals potential
357 * Free energy perturbation is turned on when A
358 * and B parameters are different.
359 * \param[in] ft Function type
360 * \param[in] phiA Dihedral angle A
361 * \param[in] cpA Force constant A
362 * \param[in] mult Multiplicity of the angle
363 * \param[in] phiB Dihedral angle B
364 * \param[in] cpB Force constant B
365 * \return The structure itself.
367 iListInput setPDihedrals(int ft, real phiA, real cpA, int mult, real phiB, real cpB)
370 iparams.pdihs.phiA = phiA;
371 iparams.pdihs.cpA = cpA;
372 iparams.pdihs.phiB = phiB;
373 iparams.pdihs.cpB = cpB;
374 iparams.pdihs.mult = mult;
375 fep = (phiA != phiB || cpA != cpB);
378 /*! \brief Set parameters for proper dihedrals potential
380 * \param[in] ft Function type
381 * \param[in] phiA Dihedral angle
382 * \param[in] cpA Force constant
383 * \param[in] mult Multiplicity of the angle
384 * \return The structure itself.
386 iListInput setPDihedrals(int ft, real phiA, real cpA, int mult)
388 return setPDihedrals(ft, phiA, cpA, mult, phiA, cpA);
390 /*! \brief Set parameters for Ryckaert-Bellemans dihedrals potential
392 * Free energy perturbation is turned on when A
393 * and B parameters are different.
394 * \param[in] rbcA Force constants A
395 * \param[in] rbcB Force constants B
396 * \return The structure itself.
398 iListInput setRbDihedrals(const real rbcA[NR_RBDIHS], const real rbcB[NR_RBDIHS])
402 for (int i = 0; i < NR_RBDIHS; i++)
404 iparams.rbdihs.rbcA[i] = rbcA[i];
405 iparams.rbdihs.rbcB[i] = rbcB[i];
406 fep = fep || (rbcA[i] != rbcB[i]);
410 /*! \brief Set parameters for Ryckaert-Bellemans dihedrals potential
412 * \param[in] rbc Force constants
413 * \return The structure itself.
415 iListInput setRbDihedrals(const real rbc[NR_RBDIHS]) { return setRbDihedrals(rbc, rbc); }
416 /*! \brief Set parameters for Polarization
418 * \param[in] alpha Polarizability
419 * \return The structure itself.
421 iListInput setPolarization(real alpha)
423 ftype = F_POLARIZATION;
425 iparams.polarize.alpha = alpha;
428 /*! \brief Set parameters for Anharmonic Polarization
430 * \param[in] alpha Polarizability (nm^3)
431 * \param[in] drcut The cut-off distance (nm) after which the
432 * fourth power kicks in
433 * \param[in] khyp The force constant for the fourth power
434 * \return The structure itself.
436 iListInput setAnharmPolarization(real alpha, real drcut, real khyp)
438 ftype = F_ANHARM_POL;
440 iparams.anharm_polarize.alpha = alpha;
441 iparams.anharm_polarize.drcut = drcut;
442 iparams.anharm_polarize.khyp = khyp;
445 /*! \brief Set parameters for Thole Polarization
447 * \param[in] a Thole factor
448 * \param[in] alpha1 Polarizability 1 (nm^3)
449 * \param[in] alpha2 Polarizability 2 (nm^3)
450 * \param[in] rfac Distance factor
451 * \return The structure itself.
453 iListInput setTholePolarization(real a, real alpha1, real alpha2, real rfac)
458 iparams.thole.alpha1 = alpha1;
459 iparams.thole.alpha2 = alpha2;
460 iparams.thole.rfac = rfac;
463 /*! \brief Set parameters for Water Polarization
465 * \param[in] alpha_x Polarizability X (nm^3)
466 * \param[in] alpha_y Polarizability Y (nm^3)
467 * \param[in] alpha_z Polarizability Z (nm^3)
468 * \param[in] rOH Oxygen-Hydrogen distance
469 * \param[in] rHH Hydrogen-Hydrogen distance
470 * \param[in] rOD Oxygen-Dummy distance
471 * \return The structure itself.
473 iListInput setWaterPolarization(real alpha_x, real alpha_y, real alpha_z, real rOH, real rHH, real rOD)
477 iparams.wpol.al_x = alpha_x;
478 iparams.wpol.al_y = alpha_y;
479 iparams.wpol.al_z = alpha_z;
480 iparams.wpol.rOH = rOH;
481 iparams.wpol.rHH = rHH;
482 iparams.wpol.rOD = rOD;
487 //! Prints the interaction and parameters to a stream
488 std::ostream& operator<<(std::ostream& out, const iListInput& input)
491 out << "Function type " << input.ftype << " called " << interaction_function[input.ftype].name
492 << " ie. labelled '" << interaction_function[input.ftype].longname << "' in an energy file"
495 // Organize to print the legacy C union t_iparams, whose
496 // relevant contents vary with ftype.
497 StringOutputStream stream;
499 TextWriter writer(&stream);
500 printInteractionParameters(&writer, input.ftype, input.iparams);
502 out << "Function parameters " << stream.toString();
503 out << "Parameters trigger FEP? " << (input.fep ? "true" : "false") << endl;
507 /*! \brief Utility to fill iatoms struct
509 * \param[in] ftype Function type
510 * \param[out] iatoms Pointer to iatoms struct
512 void fillIatoms(int ftype, std::vector<t_iatom>* iatoms)
514 std::unordered_map<int, std::vector<int>> ia = { { 2, { 0, 0, 1, 0, 1, 2, 0, 2, 3 } },
515 { 3, { 0, 0, 1, 2, 0, 1, 2, 3 } },
516 { 4, { 0, 0, 1, 2, 3 } },
517 { 5, { 0, 0, 1, 2, 3, 0 } } };
518 EXPECT_TRUE(ftype >= 0 && ftype < F_NRE);
519 int nral = interaction_function[ftype].nratoms;
520 for (auto& i : ia[nral])
522 iatoms->push_back(i);
526 class ListedForcesTest :
527 public ::testing::TestWithParam<std::tuple<iListInput, PaddedVector<RVec>, PbcType>>
532 PaddedVector<RVec> x_;
535 test::TestReferenceData refData_;
536 test::TestReferenceChecker checker_;
537 ListedForcesTest() : checker_(refData_.rootChecker())
539 input_ = std::get<0>(GetParam());
540 x_ = std::get<1>(GetParam());
541 pbcType_ = std::get<2>(GetParam());
543 box_[0][0] = box_[1][1] = box_[2][2] = 1.5;
544 set_pbc(&pbc_, pbcType_, box_);
545 // We need quite specific tolerances here since angle functions
546 // etc. are not very precise and reproducible.
547 test::FloatingPointTolerance tolerance(test::FloatingPointTolerance(
548 input_.ftoler, 1.0e-6, input_.dtoler, 1.0e-12, 10000, 100, false));
549 checker_.setDefaultTolerance(tolerance);
551 void testOneIfunc(test::TestReferenceChecker* checker, const std::vector<t_iatom>& iatoms, const real lambda)
553 SCOPED_TRACE(std::string("Testing PBC ") + c_pbcTypeNames[pbcType_]);
554 std::vector<int> ddgatindex = { 0, 1, 2, 3 };
555 std::vector<real> chargeA = { 1.5, -2.0, 1.5, -1.0 };
556 t_mdatoms mdatoms = { 0 };
557 mdatoms.chargeA = chargeA.data();
558 /* Here we run both the standard, plain-C force+shift-forcs+energy+free-energy
559 * kernel flavor and the potentially optimized, with SIMD and less output,
560 * force only kernels. Note that we also run the optimized kernel for free-energy
561 * input when lambda=0, as the force output should match the non free-energy case.
563 std::vector<BondedKernelFlavor> flavors = { BondedKernelFlavor::ForcesAndVirialAndEnergy };
564 if (!input_.fep || lambda == 0)
566 flavors.push_back(BondedKernelFlavor::ForcesSimdWhenAvailable);
568 for (const auto flavor : flavors)
570 OutputQuantities output;
572 calculateSimpleBond(input_.ftype, iatoms.size(), iatoms.data(), &input_.iparams,
573 as_rvec_array(x_.data()), output.f, output.fshift, &pbc_,
574 lambda, &output.dvdlambda, &mdatoms,
575 /* struct t_fcdata * */ nullptr, ddgatindex.data(), flavor);
576 // Internal consistency test of both test input
577 // and bonded functions.
578 EXPECT_TRUE((input_.fep || (output.dvdlambda == 0.0))) << "dvdlambda was " << output.dvdlambda;
579 checkOutput(checker, output, flavor);
584 test::TestReferenceChecker thisChecker =
585 checker_.checkCompound("FunctionType", interaction_function[input_.ftype].name)
586 .checkCompound("FEP", (input_.fep ? "Yes" : "No"));
587 std::vector<t_iatom> iatoms;
588 fillIatoms(input_.ftype, &iatoms);
591 const int numLambdas = 3;
592 for (int i = 0; i < numLambdas; ++i)
594 const real lambda = i / (numLambdas - 1.0);
595 auto valueChecker = thisChecker.checkCompound("Lambda", toString(lambda));
596 testOneIfunc(&valueChecker, iatoms, lambda);
601 testOneIfunc(&thisChecker, iatoms, 0.0);
606 TEST_P(ListedForcesTest, Ifunc)
611 //! Function types for testing bonds. Add new terms at the end.
612 std::vector<iListInput> c_InputBonds = {
613 { iListInput().setHarmonic(F_BONDS, 0.15, 500.0) },
614 { iListInput(2e-6F, 1e-8).setHarmonic(F_BONDS, 0.15, 500.0, 0.17, 400.0) },
615 { iListInput(1e-4F, 1e-8).setHarmonic(F_G96BONDS, 0.15, 50.0) },
616 { iListInput().setHarmonic(F_G96BONDS, 0.15, 50.0, 0.17, 40.0) },
617 { iListInput().setCubic(0.16, 50.0, 2.0) },
618 { iListInput(2e-6F, 1e-8).setMorse(0.15, 50.0, 2.0, 0.17, 40.0, 1.6) },
619 { iListInput(2e-6F, 1e-8).setMorse(0.15, 30.0, 2.7) },
620 { iListInput().setFene(0.4, 5.0) }
623 //! Constants for Quartic Angles
624 const real cQuarticAngles[5] = { 1.1, 2.3, 4.6, 7.8, 9.2 };
626 //! Function types for testing angles. Add new terms at the end.
627 std::vector<iListInput> c_InputAngles = {
628 { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES, 100.0, 50.0) },
629 { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES, 100.15, 50.0, 95.0, 30.0) },
630 { iListInput(8e-3, 1e-8).setHarmonic(F_G96ANGLES, 100.0, 50.0) },
631 { iListInput(8e-3, 1e-8).setHarmonic(F_G96ANGLES, 100.0, 50.0, 95.0, 30.0) },
632 { iListInput().setLinearAngle(50.0, 0.4) },
633 { iListInput().setLinearAngle(50.0, 0.4, 40.0, 0.6) },
634 { iListInput(2e-6, 1e-8).setCrossBondBonds(0.8, 0.7, 45.0) },
635 { iListInput(3e-6, 1e-8).setCrossBondAngles(0.8, 0.7, 0.3, 45.0) },
636 { iListInput(2e-2, 1e-8).setUreyBradley(950.0, 46.0, 0.3, 5.0) },
637 { iListInput(2e-2, 1e-8).setUreyBradley(100.0, 45.0, 0.3, 5.0, 90.0, 47.0, 0.32, 7.0) },
638 { iListInput(2e-3, 1e-8).setQuarticAngles(87.0, cQuarticAngles) }
641 //! Constants for Ryckaert-Bellemans A
642 const real rbcA[NR_RBDIHS] = { -5.35, 13.6, 8.4, -16.7, 0.3, 12.4 };
644 //! Constants for Ryckaert-Bellemans B
645 const real rbcB[NR_RBDIHS] = { -6.35, 12.6, 8.1, -10.7, 0.9, 15.4 };
647 //! Constants for Ryckaert-Bellemans without FEP
648 const real rbc[NR_RBDIHS] = { -7.35, 13.6, 8.4, -16.7, 1.3, 12.4 };
650 //! Function types for testing dihedrals. Add new terms at the end.
651 std::vector<iListInput> c_InputDihs = {
652 { iListInput(5e-4, 1e-8).setPDihedrals(F_PDIHS, -100.0, 10.0, 2, -80.0, 20.0) },
653 { iListInput(1e-4, 1e-8).setPDihedrals(F_PDIHS, -105.0, 15.0, 2) },
654 { iListInput(2e-4, 1e-8).setHarmonic(F_IDIHS, 100.0, 50.0) },
655 { iListInput(2e-4, 1e-8).setHarmonic(F_IDIHS, 100.15, 50.0, 95.0, 30.0) },
656 { iListInput(4e-4, 1e-8).setRbDihedrals(rbcA, rbcB) },
657 { iListInput(4e-4, 1e-8).setRbDihedrals(rbc) }
660 //! Function types for testing polarization. Add new terms at the end.
661 std::vector<iListInput> c_InputPols = {
662 { iListInput(2e-5, 1e-8).setPolarization(0.12) },
663 { iListInput(2e-3, 1e-8).setAnharmPolarization(0.0013, 0.02, 1235.6) },
664 { iListInput(1.4e-3, 1e-8).setTholePolarization(0.26, 0.07, 0.09, 1.6) },
665 { iListInput(2e-3, 1e-8).setWaterPolarization(0.001, 0.0012, 0.0016, 0.095, 0.15, 0.02) },
668 //! Function types for testing polarization. Add new terms at the end.
669 std::vector<iListInput> c_InputRestraints = {
670 { iListInput(1e-4, 1e-8).setPDihedrals(F_ANGRES, -100.0, 10.0, 2, -80.0, 20.0) },
671 { iListInput(1e-4, 1e-8).setPDihedrals(F_ANGRES, -105.0, 15.0, 2) },
672 { iListInput(1e-4, 1e-8).setPDihedrals(F_ANGRESZ, -100.0, 10.0, 2, -80.0, 20.0) },
673 { iListInput(1e-4, 1e-8).setPDihedrals(F_ANGRESZ, -105.0, 15.0, 2) },
674 { iListInput(2e-3, 1e-8).setHarmonic(F_RESTRANGLES, 100.0, 50.0) },
675 { iListInput(2e-3, 1e-8).setHarmonic(F_RESTRANGLES, 100.0, 50.0, 110.0, 45.0) }
678 //! Coordinates for testing
679 std::vector<PaddedVector<RVec>> c_coordinatesForTests = {
680 { { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.2 }, { 0.005, 0.0, 0.1 }, { -0.001, 0.1, 0.0 } },
681 { { 0.5, 0.0, 0.0 }, { 0.5, 0.0, 0.15 }, { 0.5, 0.07, 0.22 }, { 0.5, 0.18, 0.22 } },
682 { { -0.1143, -0.0282, 0.0 }, { 0.0, 0.0434, 0.0 }, { 0.1185, -0.0138, 0.0 }, { -0.0195, 0.1498, 0.0 } }
685 //! PBC values for testing
686 std::vector<PbcType> c_pbcForTests = { PbcType::No, PbcType::XY, PbcType::Xyz };
688 // Those tests give errors with the intel compiler and nothing else, so we disable them only there.
689 #ifndef __INTEL_COMPILER
690 INSTANTIATE_TEST_CASE_P(Bond,
692 ::testing::Combine(::testing::ValuesIn(c_InputBonds),
693 ::testing::ValuesIn(c_coordinatesForTests),
694 ::testing::ValuesIn(c_pbcForTests)));
696 INSTANTIATE_TEST_CASE_P(Angle,
698 ::testing::Combine(::testing::ValuesIn(c_InputAngles),
699 ::testing::ValuesIn(c_coordinatesForTests),
700 ::testing::ValuesIn(c_pbcForTests)));
702 INSTANTIATE_TEST_CASE_P(Dihedral,
704 ::testing::Combine(::testing::ValuesIn(c_InputDihs),
705 ::testing::ValuesIn(c_coordinatesForTests),
706 ::testing::ValuesIn(c_pbcForTests)));
708 INSTANTIATE_TEST_CASE_P(Polarize,
710 ::testing::Combine(::testing::ValuesIn(c_InputPols),
711 ::testing::ValuesIn(c_coordinatesForTests),
712 ::testing::ValuesIn(c_pbcForTests)));
714 INSTANTIATE_TEST_CASE_P(Restraints,
716 ::testing::Combine(::testing::ValuesIn(c_InputRestraints),
717 ::testing::ValuesIn(c_coordinatesForTests),
718 ::testing::ValuesIn(c_pbcForTests)));