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/mdtypes/mdatom.h"
58 #include "gromacs/pbcutil/ishift.h"
59 #include "gromacs/pbcutil/pbc.h"
60 #include "gromacs/topology/idef.h"
61 #include "gromacs/utility/strconvert.h"
63 #include "testutils/refdata.h"
64 #include "testutils/testasserts.h"
71 //! Number of atoms used in these tests.
72 constexpr int c_numAtoms = 4;
74 /*! \brief Output from bonded kernels
76 * \todo Later this might turn into the actual output struct. */
77 struct OutputQuantities
79 //! Energy of this interaction
81 //! Derivative with respect to lambda
84 rvec fshift[N_IVEC] = {{0}};
86 rvec4 f[c_numAtoms] = {{0}};
89 /*! \brief Utility to check the output from bonded tests
91 * \param[in] checker Reference checker
92 * \param[in] output The output from the test to check
94 void checkOutput(test::TestReferenceChecker *checker,
95 const OutputQuantities &output)
97 checker->checkReal(output.energy, "Epot ");
98 // Should still be zero if not doing FEP, so may as well test it.
99 checker->checkReal(output.dvdlambda, "dVdlambda ");
100 checker->checkVector(output.fshift[CENTRAL], "Central shift forces");
101 checker->checkSequence(std::begin(output.f), std::end(output.f), "Forces");
104 /*! \brief Input structure for listed forces tests
111 //! Tolerance for float evaluation
113 //! Tolerance for double evaluation
114 double dtoler = 1e-8;
115 //! Do free energy perturbation?
117 //! Interaction parameters
118 t_iparams iparams = {{ 0 }};
123 /*! \brief Constructor with tolerance
125 * \param[in] ftol Single precision tolerance
126 * \param[in] dtol Double precision tolerance
128 iListInput(float ftol, double dtol)
133 /*! \brief Set parameters for harmonic potential
135 * Free energy perturbation is turned on when A
136 * and B parameters are different.
137 * \param[in] ft Function type
138 * \param[in] rA Equilibrium value A
139 * \param[in] krA Force constant A
140 * \param[in] rB Equilibrium value B
141 * \param[in] krB Force constant B
142 * \return The structure itself.
144 iListInput setHarmonic(int ft, real rA, real krA, real rB, real krB)
146 iparams.harmonic.rA = rA;
147 iparams.harmonic.rB = rB;
148 iparams.harmonic.krA = krA;
149 iparams.harmonic.krB = krB;
151 fep = (rA != rB || krA != krB);
154 /*! \brief Set parameters for harmonic potential
156 * \param[in] ft Function type
157 * \param[in] rA Equilibrium value
158 * \param[in] krA Force constant
159 * \return The structure itself.
161 iListInput setHarmonic(int ft, real rA, real krA)
163 return setHarmonic(ft, rA, krA, rA, krA);
165 /*! \brief Set parameters for cubic potential
167 * \param[in] b0 Equilibrium bond length
168 * \param[in] kb Harmonic force constant
169 * \param[in] kcub Cubic force constant
170 * \return The structure itself.
172 iListInput setCubic(real b0, real kb, real kcub)
174 ftype = F_CUBICBONDS;
175 iparams.cubic.b0 = b0;
176 iparams.cubic.kb = kb;
177 iparams.cubic.kcub = kcub;
180 /*! \brief Set parameters for morse potential
182 * Free energy perturbation is turned on when A
183 * and B parameters are different.
184 * \param[in] b0A Equilibrium value A
185 * \param[in] cbA Force constant A
186 * \param[in] betaA Steepness parameter A
187 * \param[in] b0B Equilibrium value B
188 * \param[in] cbB Force constant B
189 * \param[in] betaB Steepness parameter B
190 * \return The structure itself.
192 iListInput setMorse(real b0A, real cbA, real betaA, real b0B, real cbB, real betaB)
195 iparams.morse.b0A = b0A;
196 iparams.morse.cbA = cbA;
197 iparams.morse.betaA = betaA;
198 iparams.morse.b0B = b0B;
199 iparams.morse.cbB = cbB;
200 iparams.morse.betaB = betaB;
201 fep = (b0A != b0B || cbA != cbB || betaA != betaB);
204 /*! \brief Set parameters for morse potential
206 * \param[in] b0A Equilibrium value
207 * \param[in] cbA Force constant
208 * \param[in] betaA Steepness parameter
209 * \return The structure itself.
211 iListInput setMorse(real b0A, real cbA, real betaA)
213 return setMorse(b0A, cbA, betaA, b0A, cbA, betaA);
215 /*! \brief Set parameters for fene potential
217 * \param[in] bm Equilibrium bond length
218 * \param[in] kb Force constant
219 * \return The structure itself.
221 iListInput setFene(real bm, real kb)
224 iparams.fene.bm = bm;
225 iparams.fene.kb = kb;
228 /*! \brief Set parameters for linear angle potential
230 * Free energy perturbation is turned on when A
231 * and B parameters are different.
232 * \param[in] klinA Force constant A
233 * \param[in] aA The position of the central atom A
234 * \param[in] klinB Force constant B
235 * \param[in] aB The position of the central atom B
236 * \return The structure itself.
238 iListInput setLinearAngle(real klinA, real aA, real klinB, real aB)
240 ftype = F_LINEAR_ANGLES;
241 iparams.linangle.klinA = klinA;
242 iparams.linangle.aA = aA;
243 iparams.linangle.klinB = klinB;
244 iparams.linangle.aB = aB;
245 fep = (klinA != klinB || aA != aB);
248 /*! \brief Set parameters for linear angle potential
250 * \param[in] klinA Force constant
251 * \param[in] aA The position of the central atom
252 * \return The structure itself.
254 iListInput setLinearAngle(real klinA, real aA)
256 return setLinearAngle(klinA, aA, klinA, aA);
258 /*! \brief Set parameters for Urey Bradley potential
260 * Free energy perturbation is turned on when A
261 * and B parameters are different.
262 * \param[in] thetaA Equilibrium angle A
263 * \param[in] kthetaA Force constant A
264 * \param[in] r13A The distance between i and k atoms A
265 * \param[in] kUBA The force constant for 1-3 distance A
266 * \param[in] thetaB Equilibrium angle B
267 * \param[in] kthetaB Force constant B
268 * \param[in] r13B The distance between i and k atoms B
269 * \param[in] kUBB The force constant for 1-3 distance B
270 * \return The structure itself.
272 iListInput setUreyBradley(real thetaA, real kthetaA, real r13A, real kUBA,
273 real thetaB, real kthetaB, real r13B, real kUBB)
275 ftype = F_UREY_BRADLEY;
276 iparams.u_b.thetaA = thetaA;
277 iparams.u_b.kthetaA = kthetaA;
278 iparams.u_b.r13A = r13A;
279 iparams.u_b.kUBA = kUBA;
280 iparams.u_b.thetaB = thetaB;
281 iparams.u_b.kthetaB = kthetaB;
282 iparams.u_b.r13B = r13B;
283 iparams.u_b.kUBB = kUBB;
284 fep = (thetaA != thetaB || kthetaA != kthetaB || r13A != r13B || kUBA != kUBB);
287 /*! \brief Set parameters for Urey Bradley potential
289 * \param[in] thetaA Equilibrium angle
290 * \param[in] kthetaA Force constant
291 * \param[in] r13A The distance between i and k atoms
292 * \param[in] kUBA The force constant for 1-3 distance
293 * \return The structure itself.
295 iListInput setUreyBradley(real thetaA, real kthetaA, real r13A, real kUBA)
297 return setUreyBradley(thetaA, kthetaA, r13A, kUBA,
298 thetaA, kthetaA, r13A, kUBA);
300 /*! \brief Set parameters for Cross Bond Bonds potential
302 * \param[in] r1e First bond length i-j
303 * \param[in] r2e Second bond length i-k
304 * \param[in] krr The force constant
305 * \return The structure itself.
307 iListInput setCrossBondBonds(real r1e, real r2e, real krr)
309 ftype = F_CROSS_BOND_BONDS;
310 iparams.cross_bb.r1e = r1e;
311 iparams.cross_bb.r2e = r2e;
312 iparams.cross_bb.krr = krr;
315 /*! \brief Set parameters for Cross Bond Angles potential
317 * \param[in] r1e First bond length i-j
318 * \param[in] r2e Second bond length j-k
319 * \param[in] r3e Third bond length i-k
320 * \param[in] krt The force constant
321 * \return The structure itself.
323 iListInput setCrossBondAngles(real r1e, real r2e, real r3e, real krt)
325 ftype = F_CROSS_BOND_ANGLES;
326 iparams.cross_ba.r1e = r1e;
327 iparams.cross_ba.r2e = r2e;
328 iparams.cross_ba.r3e = r3e;
329 iparams.cross_ba.krt = krt;
332 /*! \brief Set parameters for Quartic Angles potential
334 * \param[in] theta Angle
335 * \param[in] c Array of parameters
336 * \return The structure itself.
338 iListInput setQuarticAngles(real theta, const real c[5])
340 ftype = F_QUARTIC_ANGLES;
341 iparams.qangle.theta = theta;
342 iparams.qangle.c[0] = c[0];
343 iparams.qangle.c[1] = c[1];
344 iparams.qangle.c[2] = c[2];
345 iparams.qangle.c[3] = c[3];
346 iparams.qangle.c[4] = c[4];
349 /*! \brief Set parameters for proper dihedrals potential
351 * Free energy perturbation is turned on when A
352 * and B parameters are different.
353 * \param[in] phiA Dihedral angle A
354 * \param[in] cpA Force constant A
355 * \param[in] mult Multiplicity of the angle
356 * \param[in] phiB Dihedral angle B
357 * \param[in] cpB Force constant B
358 * \return The structure itself.
360 iListInput setProperDihedrals(real phiA, real cpA, int mult, real phiB, real cpB)
363 iparams.pdihs.phiA = phiA;
364 iparams.pdihs.cpA = cpA;
365 iparams.pdihs.phiB = phiB;
366 iparams.pdihs.cpB = cpB;
367 iparams.pdihs.mult = mult;
368 fep = (phiA != phiB || cpA != cpB);
371 /*! \brief Set parameters for proper dihedrals potential
373 * \param[in] phiA Dihedral angle
374 * \param[in] cpA Force constant
375 * \param[in] mult Multiplicity of the angle
376 * \return The structure itself.
378 iListInput setProperDihedrals(real phiA, real cpA, int mult)
380 return setProperDihedrals(phiA, cpA, mult, phiA, cpA);
382 /*! \brief Set parameters for Ryckaert-Bellemans dihedrals potential
384 * Free energy perturbation is turned on when A
385 * and B parameters are different.
386 * \param[in] rbcA Force constants A
387 * \param[in] rbcB Force constants B
388 * \return The structure itself.
390 iListInput setRbDihedrals(const real rbcA[NR_RBDIHS], const real rbcB[NR_RBDIHS])
394 for (int i = 0; i < NR_RBDIHS; i++)
396 iparams.rbdihs.rbcA[i] = rbcA[i];
397 iparams.rbdihs.rbcB[i] = rbcB[i];
398 fep = fep || (rbcA[i] != rbcB[i]);
402 /*! \brief Set parameters for Ryckaert-Bellemans dihedrals potential
404 * \param[in] rbc Force constants
405 * \return The structure itself.
407 iListInput setRbDihedrals(const real rbc[NR_RBDIHS])
409 return setRbDihedrals(rbc, rbc);
411 /*! \brief Set parameters for Polarization
413 * \param[in] alpha Polarizability
414 * \return The structure itself.
416 iListInput setPolarization(real alpha)
418 ftype = F_POLARIZATION;
420 iparams.polarize.alpha = alpha;
423 /*! \brief Set parameters for Anharmonic Polarization
425 * \param[in] alpha Polarizability (nm^3)
426 * \param[in] drcut The cut-off distance (nm) after which the
427 * fourth power kicks in
428 * \param[in] khyp The force constant for the fourth power
429 * \return The structure itself.
431 iListInput setAnharmPolarization(real alpha,
435 ftype = F_ANHARM_POL;
437 iparams.anharm_polarize.alpha = alpha;
438 iparams.anharm_polarize.drcut = drcut;
439 iparams.anharm_polarize.khyp = khyp;
442 /*! \brief Set parameters for Thole Polarization
444 * \param[in] a Thole factor
445 * \param[in] alpha1 Polarizability 1 (nm^3)
446 * \param[in] alpha2 Polarizability 2 (nm^3)
447 * \param[in] rfac Distance factor
448 * \return The structure itself.
450 iListInput setTholePolarization(real a,
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,
482 iparams.wpol.al_x = alpha_x;
483 iparams.wpol.al_y = alpha_y;
484 iparams.wpol.al_z = alpha_z;
485 iparams.wpol.rOH = rOH;
486 iparams.wpol.rHH = rHH;
487 iparams.wpol.rOD = rOD;
492 /*! \brief Utility to fill iatoms struct
494 * \param[in] ftype Function type
495 * \param[out] iatoms Pointer to iatoms struct
497 void fillIatoms(int ftype, std::vector<t_iatom> *iatoms)
499 std::unordered_map<int, std::vector<int> > ia =
500 { { 2, { 0, 0, 1, 0, 1, 2, 0, 2, 3 } },
501 { 3, { 0, 0, 1, 2, 0, 1, 2, 3 } },
502 { 4, { 0, 0, 1, 2, 3 } },
503 { 5, { 0, 0, 1, 2, 3, 0 } } };
504 EXPECT_TRUE(ftype >= 0 && ftype < F_NRE);
505 int nral = interaction_function[ftype].nratoms;
506 for (auto &i : ia[nral])
508 iatoms->push_back(i);
512 class ListedForcesTest : public ::testing::TestWithParam<std::tuple<iListInput, std::vector<gmx::RVec>, int> >
517 std::vector<gmx::RVec> x_;
520 test::TestReferenceData refData_;
521 test::TestReferenceChecker checker_;
522 ListedForcesTest( ) :
523 checker_(refData_.rootChecker())
525 input_ = std::get<0>(GetParam());
526 x_ = std::get<1>(GetParam());
527 epbc_ = std::get<2>(GetParam());
529 box_[0][0] = box_[1][1] = box_[2][2] = 1.5;
530 set_pbc(&pbc_, epbc_, box_);
532 // We need quite specific tolerances here since angle functions
533 // etc. are not very precise and reproducible.
534 test::FloatingPointTolerance tolerance(test::FloatingPointTolerance(input_.ftoler, 1.0e-6,
535 input_.dtoler, 1.0e-12,
537 checker_.setDefaultTolerance(tolerance);
539 void testOneIfunc(test::TestReferenceChecker *checker,
540 const std::vector<t_iatom> &iatoms,
543 SCOPED_TRACE(std::string("Testing PBC ") + epbc_names[epbc_]);
544 std::vector<int> ddgatindex = { 0, 1, 2, 3 };
545 std::vector<real> chargeA = { 1.5, -2.0, 1.5, -1.0 };
546 t_mdatoms mdatoms = {0};
547 mdatoms.chargeA = chargeA.data();
548 OutputQuantities output;
549 output.energy = bondedFunction(input_.ftype) (iatoms.size(),
552 as_rvec_array(x_.data()),
553 output.f, output.fshift,
555 /* const struct t_graph *g */ nullptr,
556 lambda, &output.dvdlambda,
558 /* struct t_fcdata *fcd */ nullptr,
560 // Internal consistency test of both test input
561 // and bonded functions.
562 EXPECT_TRUE((input_.fep || (output.dvdlambda == 0.0)));
563 checkOutput(checker, output);
567 test::TestReferenceChecker thisChecker =
568 checker_.checkCompound("FunctionType",
569 interaction_function[input_.ftype].name).checkCompound("FEP", (input_.fep ? "Yes" : "No"));
570 std::vector<t_iatom> iatoms;
571 fillIatoms(input_.ftype, &iatoms);
574 const int numLambdas = 3;
575 for (int i = 0; i < numLambdas; ++i)
577 const real lambda = i / (numLambdas - 1.0);
578 auto valueChecker = thisChecker.checkCompound("Lambda", toString(lambda));
579 testOneIfunc(&valueChecker, iatoms, lambda);
584 testOneIfunc(&thisChecker, iatoms, 0.0);
589 TEST_P (ListedForcesTest, Ifunc)
594 //! Function types for testing bonds. Add new terms at the end.
595 std::vector<iListInput> c_InputBonds =
597 { iListInput().setHarmonic(F_BONDS, 0.15, 500.0) },
598 { iListInput(2e-6f, 1e-8).setHarmonic(F_BONDS, 0.15, 500.0, 0.17, 400.0) },
599 { iListInput(1e-4f, 1e-8).setHarmonic(F_G96BONDS, 0.15, 50.0) },
600 { iListInput().setHarmonic(F_G96BONDS, 0.15, 50.0, 0.17, 40.0) },
601 { iListInput().setCubic(0.16, 50.0, 2.0) },
602 { iListInput(2e-6f, 1e-8).setMorse(0.15, 50.0, 2.0, 0.17, 40.0, 1.6) },
603 { iListInput(2e-6f, 1e-8).setMorse(0.15, 30.0, 2.7) },
604 { iListInput().setFene(0.4, 5.0) }
607 //! Constants for Quartic Angles
608 const real cQuarticAngles[5] = { 1.1, 2.3, 4.6, 7.8, 9.2 };
610 //! Function types for testing angles. Add new terms at the end.
611 std::vector<iListInput> c_InputAngles =
613 { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES, 100.0, 50.0) },
614 { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES, 100.15, 50.0, 95.0, 30.0) },
615 { iListInput(8e-3, 1e-8).setHarmonic(F_G96ANGLES, 100.0, 50.0) },
616 { iListInput(8e-3, 1e-8).setHarmonic(F_G96ANGLES, 100.0, 50.0, 95.0, 30.0) },
617 { iListInput().setLinearAngle(50.0, 0.4) },
618 { iListInput().setLinearAngle(50.0, 0.4, 40.0, 0.6) },
619 { iListInput(2e-6, 1e-8).setCrossBondBonds(0.8, 0.7, 45.0) },
620 { iListInput(2e-6, 1e-8).setCrossBondAngles(0.8, 0.7, 0.3, 45.0) },
621 { iListInput(2e-2, 1e-8).setUreyBradley(950.0, 46.0, 0.3, 5.0) },
622 { iListInput(2e-2, 1e-8).setUreyBradley(100.0, 45.0, 0.3, 5.0, 90.0, 47.0, 0.32, 7.0) },
623 { iListInput(2e-3, 1e-8).setQuarticAngles(87.0, cQuarticAngles ) }
626 //! Constants for Ryckaert-Bellemans A
627 const real rbcA[NR_RBDIHS] = { -5.35, 13.6, 8.4, -16.7, 0.3, 12.4 };
629 //! Constants for Ryckaert-Bellemans B
630 const real rbcB[NR_RBDIHS] = { -6.35, 12.6, 8.1, -10.7, 0.9, 15.4 };
632 //! Constants for Ryckaert-Bellemans without FEP
633 const real rbc[NR_RBDIHS] = { -7.35, 13.6, 8.4, -16.7, 1.3, 12.4 };
635 //! Function types for testing dihedrals. Add new terms at the end.
636 std::vector<iListInput> c_InputDihs =
638 { iListInput(1e-4, 1e-8).setProperDihedrals(-100.0, 10.0, 2, -80.0, 20.0) },
639 { iListInput(1e-4, 1e-8).setProperDihedrals(-105.0, 15.0, 2) },
640 { iListInput(2e-4, 1e-8).setHarmonic(F_IDIHS, 100.0, 50.0) },
641 { iListInput(2e-4, 1e-8).setHarmonic(F_IDIHS, 100.15, 50.0, 95.0, 30.0) },
642 { iListInput(3e-4, 1e-8).setRbDihedrals(rbcA, rbcB) },
643 { iListInput(3e-4, 1e-8).setRbDihedrals(rbc) }
646 //! Function types for testing polarization. Add new terms at the end.
647 std::vector<iListInput> c_InputPols =
649 { iListInput(2e-5, 1e-8).setPolarization(0.12) },
650 { iListInput(1.7e-3, 1e-8).setAnharmPolarization(0.0013, 0.02, 1235.6) },
651 { iListInput(1.4e-3, 1e-8).setTholePolarization(0.26, 0.07, 0.09, 1.6) },
652 { iListInput().setWaterPolarization(0.001, 0.0012, 0.0016, 0.095, 0.15, 0.02) },
655 //! Coordinates for testing
656 std::vector<std::vector<gmx::RVec> > c_coordinatesForTests =
658 {{ 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.2 }, { 0.005, 0.0, 0.1 }, { -0.001, 0.1, 0.0 }},
659 {{ 0.5, 0.0, 0.0 }, { 0.5, 0.0, 0.15 }, { 0.5, 0.07, 0.22 }, { 0.5, 0.18, 0.22 }},
660 {{ -0.1143, -0.0282, 0.0 }, { 0.0, 0.0434, 0.0 }, { 0.1185, -0.0138, 0.0 }, { -0.0195, 0.1498, 0.0 }}
663 //! PBC values for testing
664 std::vector<int> c_pbcForTests = { epbcNONE, epbcXY, epbcXYZ };
666 INSTANTIATE_TEST_CASE_P(Bond, ListedForcesTest, ::testing::Combine(::testing::ValuesIn(c_InputBonds), ::testing::ValuesIn(c_coordinatesForTests), ::testing::ValuesIn(c_pbcForTests)));
668 INSTANTIATE_TEST_CASE_P(Angle, ListedForcesTest, ::testing::Combine(::testing::ValuesIn(c_InputAngles), ::testing::ValuesIn(c_coordinatesForTests), ::testing::ValuesIn(c_pbcForTests)));
670 INSTANTIATE_TEST_CASE_P(Dihedral, ListedForcesTest, ::testing::Combine(::testing::ValuesIn(c_InputDihs), ::testing::ValuesIn(c_coordinatesForTests), ::testing::ValuesIn(c_pbcForTests)));
672 INSTANTIATE_TEST_CASE_P(Polarize, ListedForcesTest, ::testing::Combine(::testing::ValuesIn(c_InputPols), ::testing::ValuesIn(c_coordinatesForTests), ::testing::ValuesIn(c_pbcForTests)));