2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2021, 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 Implements test of 1-4 interactions
38 * This test is copied from the bonded interactions test and slightly
39 * modified since 'do_pairs' takes a different set of arguments than
40 * 'calculateSimpleBond'. To keep the test setup uncluttered this test is
41 * therefore not merged into the bonded test but implemented standalone.
43 * The test setup consists of 2 atom pairs that are tested in an fep setting
44 * (vanishing charge and lennard-jones parameters of one atom) and without
45 * fep. Placement of the atoms in the box is such that shift-forces and pbc
46 * paths in do_pairs are covered.
48 * \author Sebastian Kehl <sebastian.kehl@mpcdf.mpg.de>
49 * \ingroup module_listed_forces
53 #include "gromacs/listed_forces/bonded.h"
58 #include <unordered_map>
60 #include <gtest/gtest.h>
62 #include "gromacs/listed_forces/listed_forces.h"
63 #include "gromacs/listed_forces/pairs.h"
64 #include "gromacs/math/paddedvector.h"
65 #include "gromacs/math/units.h"
66 #include "gromacs/math/vec.h"
67 #include "gromacs/math/vectypes.h"
68 #include "gromacs/mdtypes/mdatom.h"
69 #include "gromacs/mdtypes/simulation_workload.h"
70 #include "gromacs/mdtypes/enerdata.h"
71 #include "gromacs/mdtypes/forcerec.h"
72 #include "gromacs/mdtypes/inputrec.h"
73 #include "gromacs/mdtypes/interaction_const.h"
74 #include "gromacs/mdtypes/nblist.h"
75 #include "gromacs/tables/forcetable.h"
76 #include "gromacs/pbcutil/ishift.h"
77 #include "gromacs/pbcutil/pbc.h"
78 #include "gromacs/topology/idef.h"
79 #include "gromacs/utility/enumerationhelpers.h"
80 #include "gromacs/utility/strconvert.h"
81 #include "gromacs/utility/stringstream.h"
82 #include "gromacs/utility/textwriter.h"
83 #include "gromacs/utility/fatalerror.h"
85 #include "testutils/refdata.h"
86 #include "testutils/testasserts.h"
95 //! Number of atoms used in these tests.
96 constexpr int c_numAtoms = 3;
98 /*! \brief Output from pairs kernels
101 struct OutputQuantities
103 OutputQuantities(int energyGroup) :
104 energy(energyGroup), dvdLambda(static_cast<int>(FreeEnergyPerturbationCouplingType::Count), 0.0)
108 //! Energy of this interaction
109 gmx_grppairener_t energy;
110 //! Derivative with respect to lambda
111 std::vector<real> dvdLambda;
113 rvec fShift[gmx::detail::c_numIvecs] = { { 0 } };
115 alignas(GMX_REAL_MAX_SIMD_WIDTH * sizeof(real)) rvec4 f[c_numAtoms] = { { 0 } };
118 /*! \brief Utility to check the output from pairs tests
120 * \param[in] checker Reference checker
121 * \param[in] output The output from the test to check
122 * \param[in] bondedKernelFlavor Flavor for determining what output to check
123 * \param[in] functionType type of the interaction
125 void checkOutput(TestReferenceChecker* checker,
126 const OutputQuantities& output,
127 const BondedKernelFlavor bondedKernelFlavor,
128 const int functionType)
130 if (computeEnergy(bondedKernelFlavor))
132 switch (functionType)
136 checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::Coulomb14][0],
138 checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::LJ14][0],
142 checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR][0],
144 checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::LJSR][0],
147 default: gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", functionType);
149 checker->checkReal(output.dvdLambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)],
151 checker->checkReal(output.dvdLambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)],
154 checker->checkSequence(std::begin(output.f), std::end(output.f), "Forces");
157 /* \brief Utility class to setup forcerec and interaction parameters
159 * Data is only initialized as necessary for the 1-4 interactions to work!
166 fepVals_.sc_alpha = 0.3;
167 fepVals_.sc_power = 1;
168 fepVals_.sc_r_power = 6.0;
169 fepVals_.sc_sigma = 0.3;
170 fepVals_.sc_sigma_min = 0.3;
171 fepVals_.bScCoul = true;
172 fepVals_.scGapsysScaleLinpointLJ = 0.85;
173 fepVals_.scGapsysScaleLinpointQ = 0.3;
174 fepVals_.scGapsysSigmaLJ = 0.3;
175 fepVals_.softcoreFunction = SoftcoreType::Beutler;
177 fr_.ic = std::make_unique<interaction_const_t>();
179 fr_.ic->softCoreParameters = std::make_unique<interaction_const_t::SoftCoreParameters>(fepVals_);
182 real tableRange = 2.9;
183 fr_.pairsTable = make_tables(nullptr, fr_.ic.get(), nullptr, tableRange, GMX_MAKETABLES_14ONLY);
186 fr_.use_simd_kernels = useSimd_;
187 fr_.bMolPBC = haveMolPBC_;
190 t_forcerec* get() { return &fr_; }
193 FreeEnergyPerturbationType haveFep_ = FreeEnergyPerturbationType::No;
194 bool useSimd_ = false;
195 bool haveMolPBC_ = false;
201 ForcerecHelper frHelper;
204 /*! \brief Input structure for listed forces tests
213 //! Tolerance for float evaluation
214 float floatToler = 1e-6;
215 //! Tolerance for double evaluation
216 double doubleToler = 1e-8;
217 //! Interaction parameters
218 t_iparams iparams = { { 0 } };
223 /*! \brief Constructor with tolerance
225 * \param[in] ftol Single precision tolerance
226 * \param[in] dtol Double precision tolerance
228 ListInput(float ftol, double dtol)
234 /*! \brief Set parameters for lj14 interaction
236 * Fep is used if either c6A != c6B or c12A != c12B.
238 * \param[in] c6A lj-c6 of state A
239 * \param[in] c12A lj-c12 of state A
240 * \param[in] c6B lj-c6 of state B
241 * \param[in] c12B lj-c12 of state B
243 ListInput setLj14Interaction(real c6A, real c12A, real c6B, real c12B)
246 fep = (c6A != c6B || c12A != c12B);
247 iparams.lj14.c6A = c6A;
248 iparams.lj14.c12A = c12A;
249 iparams.lj14.c6B = c6B;
250 iparams.lj14.c12B = c12B;
255 /*! \brief Set parameters for ljc14 interaction
257 * \param[in] qi charge i
258 * \param[in] qj charge j
259 * \param[in] c6 lj-c6
260 * \param[in] c12 lj-c12
261 * \param[in] fudgeQ fudge factor
263 ListInput setLjc14Interaction(real qi, real qj, real c6, real c12, real fudgeQ)
266 iparams.ljc14.qi = qi;
267 iparams.ljc14.qj = qj;
268 iparams.ljc14.c6 = c6;
269 iparams.ljc14.c12 = c12;
270 iparams.ljc14.fqq = fudgeQ;
275 /*! \brief Set parameters for ljcnb interaction
277 * \param[in] qi charge i
278 * \param[in] qj charge j
279 * \param[in] c6 lj-c6
280 * \param[in] c12 lj-c12
282 ListInput setLjcnbInteraction(real qi, real qj, real c6, real c12)
284 fType = F_LJC_PAIRS_NB;
285 iparams.ljcnb.qi = qi;
286 iparams.ljcnb.qj = qj;
287 iparams.ljcnb.c6 = c6;
288 iparams.ljcnb.c12 = c12;
294 class ListedForcesPairsTest :
295 public ::testing::TestWithParam<std::tuple<ListInput, PaddedVector<RVec>, PbcType>>
300 PaddedVector<RVec> x_;
303 TestReferenceData refData_;
304 TestReferenceChecker checker_;
306 ListedForcesPairsTest() : checker_(refData_.rootChecker())
308 input_ = std::get<0>(GetParam());
309 x_ = std::get<1>(GetParam());
310 pbcType_ = std::get<2>(GetParam());
312 box_[0][0] = box_[1][1] = box_[2][2] = 1.0;
313 set_pbc(&pbc_, pbcType_, box_);
315 FloatingPointTolerance tolerance = relativeToleranceAsPrecisionDependentFloatingPoint(
316 1.0, input_.floatToler, input_.doubleToler);
318 checker_.setDefaultTolerance(tolerance);
321 void testOneIfunc(TestReferenceChecker* checker, const real lambda, const SoftcoreType softcoreType)
323 SCOPED_TRACE(std::string("Testing PBC type: ") + c_pbcTypeNames[pbcType_]);
325 // 'definition of pairs' is a concatenation of #npairs (here 2)
326 // 'nAtomsPerPair+1'-tuples (fType a_0 a_i ... a_nAtomsPerPair)
327 std::vector<t_iatom> iatoms = { 0, 1, 2, 0, 0, 2 };
329 std::vector<int> ddgatindex = { 0, 1, 2 };
330 std::vector<real> chargeA = { 1.0, -0.5, -0.5 };
331 std::vector<real> chargeB = { 0.0, 0.0, 0.0 };
332 std::vector<unsigned short> egrp = { 0, 0, 0 };
333 t_mdatoms mdatoms = { 0 };
335 mdatoms.chargeA = chargeA.data();
336 mdatoms.chargeB = chargeB.data();
337 mdatoms.cENER = egrp.data();
338 // nPerturbed is not decisive for fep to be used; it is overruled by
339 // other conditions in do_pairs_general; just here to not segfault
341 mdatoms.nPerturbed = 0;
343 t_forcerec* fr = frHelper.get();
344 fr->efep = input_.fep ? FreeEnergyPerturbationType::Yes : FreeEnergyPerturbationType::No;
345 fr->ic->softCoreParameters->softcoreType = softcoreType;
346 if (pbcType_ != PbcType::No)
351 std::vector<BondedKernelFlavor> flavors = { BondedKernelFlavor::ForcesAndVirialAndEnergy };
353 if (!input_.fep || lambda == 0)
355 fr->use_simd_kernels = true;
356 flavors.push_back(BondedKernelFlavor::ForcesSimdWhenAvailable);
359 for (const auto flavor : flavors)
361 SCOPED_TRACE("Testing bonded kernel flavor: " + c_bondedKernelFlavorStrings[flavor]);
363 StepWorkload stepWork;
364 stepWork.computeVirial = computeVirial(flavor);
365 stepWork.computeEnergy = computeEnergy(flavor);
367 bool havePerturbedInteractions = input_.fep;
368 if (flavor == BondedKernelFlavor::ForcesSimdWhenAvailable)
370 havePerturbedInteractions = false;
373 int numEnergyTerms = static_cast<int>(NonBondedEnergyTerms::Count);
374 int numFepCouplingTerms = static_cast<int>(FreeEnergyPerturbationCouplingType::Count);
375 OutputQuantities output(numEnergyTerms);
376 std::vector<real> lambdas(numFepCouplingTerms, lambda);
378 do_pairs(input_.fType,
382 as_rvec_array(x_.data()),
387 output.dvdLambda.data(),
388 gmx::arrayRefFromArray(mdatoms.chargeA, mdatoms.nr),
389 gmx::arrayRefFromArray(mdatoms.chargeB, mdatoms.nr),
390 gmx::arrayRefFromArray(mdatoms.bPerturbed, mdatoms.nr),
391 gmx::arrayRefFromArray(mdatoms.cENER, mdatoms.nr),
394 havePerturbedInteractions,
399 checkOutput(checker, output, flavor, input_.fType);
400 auto shiftForcesChecker = checker->checkCompound("Shift-Forces", "Shift-forces");
402 if (computeVirial(flavor))
404 shiftForcesChecker.checkVector(output.fShift[gmx::c_centralShiftIndex], "Central");
408 // Permit omitting to compare shift forces with
409 // reference data when that is useless.
410 shiftForcesChecker.disableUnusedEntriesCheck();
417 TestReferenceChecker thisChecker =
418 checker_.checkCompound("FunctionType", interaction_function[input_.fType].name)
419 .checkCompound("FEP", (input_.fep ? "Yes" : "No"));
423 const int numLambdas = 3;
424 for (int i = 0; i < numLambdas; ++i)
426 const real lambda = i / (numLambdas - 1.0);
427 for (SoftcoreType c : EnumerationWrapper<SoftcoreType>{})
429 auto lambdaChecker = thisChecker.checkCompound("Lambda", toString(lambda));
430 auto softcoreChecker = lambdaChecker.checkCompound("Sofcore", enumValueToString(c));
431 testOneIfunc(&softcoreChecker, lambda, c);
437 testOneIfunc(&thisChecker, 0.0, SoftcoreType::Beutler);
442 TEST_P(ListedForcesPairsTest, Ifunc)
447 //! Function types for testing 1-4 interaction. Add new terms at the end.
448 std::vector<ListInput> c_14Interaction = {
449 { ListInput(1e-5, 1e-7).setLj14Interaction(0.001458, 1.0062882e-6, 0.0, 0.0) },
450 { ListInput(1e-5, 1e-7).setLj14Interaction(0.001458, 1.0062882e-6, 0.001458, 1.0062882e-6) },
451 { ListInput(1e-5, 1e-7).setLjc14Interaction(1.0, -1.0, 0.001458, 1.0062882e-6, 0.5) },
452 { ListInput(1e-5, 1e-7).setLjcnbInteraction(1.0, -1.0, 0.001458, 1.0062882e-6) }
455 //! PBC values for testing
456 std::vector<PbcType> c_pbcForTests = { PbcType::No, PbcType::XY, PbcType::Xyz };
458 /*! \brief Coordinates for testing 1-4 interaction
460 * Define coordinates for 3 atoms here, which will be used in 2 interactions.
462 std::vector<PaddedVector<RVec>> c_coordinatesFor14Interaction = {
463 { { 0.0, 0.0, 0.0 }, { 1.0, 1.0, 1.0 }, { 1.1, 1.2, 1.3 } }
466 INSTANTIATE_TEST_SUITE_P(14Interaction,
467 ListedForcesPairsTest,
468 ::testing::Combine(::testing::ValuesIn(c_14Interaction),
469 ::testing::ValuesIn(c_coordinatesFor14Interaction),
470 ::testing::ValuesIn(c_pbcForTests)));