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.
37 * \brief Implements test of nonbonded fep kernel
39 * Implements the test logic from the bonded interactions also for the
40 * nonbonded fep kernel. This requires setting up some more input
41 * structures that in the bonded case.
43 * The test setup consists of an atom pair that is evaluated in an fep setting
44 * (vanishing charge and lennard-jones parameters of atom #2) with and without
45 * softcore Potentials.
47 * \author Sebastian Kehl <sebastian.kehl@mpcdf.mpg.de>
48 * \ingroup module_gmxlib_nonbonded
52 #include "gromacs/gmxlib/nonbonded/nb_free_energy.h"
53 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
57 #include <gtest/gtest.h>
59 #include "gromacs/math/paddedvector.h"
60 #include "gromacs/math/units.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/math/vectypes.h"
63 #include "gromacs/math/arrayrefwithpadding.h"
64 #include "gromacs/mdtypes/mdatom.h"
65 #include "gromacs/mdtypes/enerdata.h"
66 #include "gromacs/mdtypes/forcerec.h"
67 #include "gromacs/mdtypes/inputrec.h"
68 #include "gromacs/mdtypes/interaction_const.h"
69 #include "gromacs/mdtypes/nblist.h"
70 #include "gromacs/mdtypes/forceoutput.h"
71 #include "gromacs/mdlib/forcerec.h"
72 #include "gromacs/tables/forcetable.h"
73 #include "gromacs/pbcutil/ishift.h"
74 #include "gromacs/pbcutil/pbc.h"
75 #include "gromacs/topology/idef.h"
76 #include "gromacs/topology/forcefieldparameters.h"
77 #include "gromacs/utility/enumerationhelpers.h"
78 #include "gromacs/utility/strconvert.h"
79 #include "gromacs/utility/stringstream.h"
80 #include "gromacs/utility/textwriter.h"
81 #include "gromacs/utility/arrayref.h"
82 #include "gromacs/utility/smalloc.h"
83 #include "gromacs/ewald/ewald_utils.h"
84 #include "gromacs/gmxlib/nrnb.h"
86 #include "testutils/refdata.h"
87 #include "testutils/testasserts.h"
96 //! Number of atoms used in these tests.
97 constexpr int c_numAtoms = 4;
98 constexpr int c_numAtomTypes = 3;
100 /*! \brief Output from nonbonded fep kernel
103 struct OutputQuantities
106 energy(static_cast<int>(NonBondedEnergyTerms::Count)),
107 dvdLambda(static_cast<int>(FreeEnergyPerturbationCouplingType::Count), 0.0),
108 fShift(1, { 0.0, 0.0, 0.0 }),
109 f(c_numAtoms, { 0.0, 0.0, 0.0 })
113 //! Energies of this interaction (size EgNR)
114 gmx_grppairener_t energy;
115 //! Derivative with respect to lambda (size efptNR)
116 std::vector<real> dvdLambda;
117 //! Shift force vectors (size N_IVEC but in this test only 1)
118 std::vector<RVec> fShift;
119 //! Forces (size c_numAtoms)
120 PaddedVector<RVec> f;
123 /*! \brief Utility to check the output from nonbonded test
125 * \param[in] checker Reference checker
126 * \param[in] output The output from the test to check
128 void checkOutput(TestReferenceChecker* checker, const OutputQuantities& output)
130 checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::LJSR][0], "EVdw ");
131 checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR][0],
133 checker->checkReal(output.dvdLambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)],
135 checker->checkReal(output.dvdLambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)],
138 checker->checkSequence(std::begin(output.f), std::end(output.f), "Forces");
140 auto shiftForcesChecker = checker->checkCompound("Shift-Forces", "Shift-forces");
141 shiftForcesChecker.checkVector(output.fShift[0], "Central");
144 class InteractionConstHelper
147 InteractionConstHelper() {}
149 //! init data to construct interaction_const
150 void initInteractionConst(CoulombInteractionType coulType, VanDerWaalsType vdwType, InteractionModifiers vdwMod)
152 coulType_ = coulType;
156 // initialize correction tables
157 interaction_const_t tmp;
158 tmp.ewaldcoeff_q = calc_ewaldcoeff_q(1.0, 1.0e-5);
159 coulEwaldCoeff_ = tmp.ewaldcoeff_q;
160 tmp.ewaldcoeff_lj = calc_ewaldcoeff_lj(1.0, 1.0e-5);
161 vdwEwaldCoeff_ = tmp.ewaldcoeff_lj;
162 tmp.eeltype = coulType;
163 tmp.vdwtype = vdwType;
164 tmp.coulombEwaldTables = std::make_unique<EwaldCorrectionTables>();
165 tmp.vdwEwaldTables = std::make_unique<EwaldCorrectionTables>();
167 init_interaction_const_tables(nullptr, &tmp, 1.0, 0.0);
168 coulombTables_ = *tmp.coulombEwaldTables;
169 vdwTables_ = *tmp.vdwEwaldTables;
172 /*! \brief Setup interaction_const_t
174 * \param[in] fepVals t_lambda struct of fep values
175 * \param[out] ic interaction_const_t pointer with data
177 void getInteractionConst(const t_lambda& fepVals, interaction_const_t* ic)
179 ic->softCoreParameters = std::unique_ptr<interaction_const_t::SoftCoreParameters>(
180 new interaction_const_t::SoftCoreParameters(fepVals));
182 ic->coulombEwaldTables = std::unique_ptr<EwaldCorrectionTables>(new EwaldCorrectionTables);
183 *ic->coulombEwaldTables = coulombTables_;
185 ic->vdwEwaldTables = std::unique_ptr<EwaldCorrectionTables>(new EwaldCorrectionTables);
186 *ic->vdwEwaldTables = vdwTables_;
188 // set coulomb and vdw types
189 ic->eeltype = coulType_;
190 ic->vdwtype = vdwType_;
191 ic->vdw_modifier = vdwMod_;
193 // some non default parameters used in this testcase
194 ic->epsfac = gmx::c_one4PiEps0 * 0.25;
195 ic->reactionFieldCoefficient = 0.0; // former k_rf
196 ic->reactionFieldShift = 1.0; // former c_rf
197 ic->ewaldcoeff_q = coulEwaldCoeff_;
198 ic->ewaldcoeff_lj = vdwEwaldCoeff_;
199 ic->sh_ewald = 1.0e-5;
200 ic->sh_lj_ewald = -1.0;
201 ic->dispersion_shift.cpot = -1.0;
202 ic->repulsion_shift.cpot = -1.0;
206 //! correction tables
207 EwaldCorrectionTables coulombTables_;
208 EwaldCorrectionTables vdwTables_;
210 //! coulomb and vdw type specifiers
211 CoulombInteractionType coulType_;
212 real coulEwaldCoeff_;
213 VanDerWaalsType vdwType_;
214 InteractionModifiers vdwMod_;
219 /* \brief Utility class to setup forcerec
221 * This helper takes care of handling the neccessary data that are kept
222 * at various places and in various forms in the forcerec hierarchy such
223 * that this class can safely be used.
225 * Data is only initialized as necessary for the nonbonded kernel to work!
232 fepVals_.sc_alpha = 0.3;
233 fepVals_.sc_power = 1;
234 fepVals_.sc_r_power = 6.0;
235 fepVals_.sc_sigma = 0.3;
236 fepVals_.sc_sigma_min = 0.3;
237 fepVals_.bScCoul = true;
238 fepVals_.scGapsysScaleLinpointLJ = 0.85;
239 fepVals_.scGapsysScaleLinpointQ = 0.3;
240 fepVals_.scGapsysSigmaLJ = 0.3;
241 fepVals_.softcoreFunction = SoftcoreType::Beutler;
244 //! initialize data structure to construct forcerec
245 void initForcerec(const gmx_ffparams_t& idef,
246 CoulombInteractionType coulType,
247 VanDerWaalsType vdwType,
248 InteractionModifiers vdwMod)
250 icHelper_.initInteractionConst(coulType, vdwType, vdwMod);
251 nbfp_ = makeNonBondedParameterLists(idef.atnr, idef.iparams, false);
252 ljPmeC6Grid_ = makeLJPmeC6GridCorrectionParameters(idef.atnr, idef.iparams, LongRangeVdW::Geom);
255 void setSoftcoreAlpha(const real scBeutlerAlphaOrGapsysLinpointScaling)
257 fepVals_.sc_alpha = scBeutlerAlphaOrGapsysLinpointScaling;
258 fepVals_.scGapsysScaleLinpointLJ = scBeutlerAlphaOrGapsysLinpointScaling;
259 fepVals_.scGapsysScaleLinpointQ = scBeutlerAlphaOrGapsysLinpointScaling;
261 void setSoftcoreCoulomb(const bool scCoulomb) { fepVals_.bScCoul = scCoulomb; }
262 void setSoftcoreType(const SoftcoreType softcoreType)
264 fepVals_.softcoreFunction = softcoreType;
268 //! get forcerec data as wanted by the nonbonded kernel
269 void getForcerec(t_forcerec* fr)
271 fr->ic = std::make_unique<interaction_const_t>();
274 icHelper_.getInteractionConst(fepVals_, fr->ic.get());
277 fr->ljpme_c6grid = ljPmeC6Grid_;
279 fr->shift_vec = { { 0.0, 0.0, 0.0 } };
281 fr->ntype = c_numAtomTypes;
284 fr->use_simd_kernels = GMX_USE_SIMD_KERNELS;
288 InteractionConstHelper icHelper_;
289 std::vector<real> ljPmeC6Grid_;
290 std::vector<real> nbfp_;
295 /*! \brief Utility structure to hold atoms data
297 * A system having 4 interactions with the following perturbation pattern:
299 * - vdw- and coulomb-perturbation
300 * - coulomb-perturbation only
301 * - vdw-perturbation only
303 * This is realized by defining 3 different atom types that control
304 * the vdw-perturbation. The coulomb-perturbation is controlled by directly
305 * setting the charge of the atoms at the lambda states A/B.
311 idef.atnr = c_numAtomTypes;
312 idef.iparams.resize(2 * c_numAtomTypes * c_numAtomTypes);
314 // set interaction parameters for different combinations of types
315 idef.iparams[0].lj = { 0.001458, 1.0062882e-6 }; // 0-0
316 idef.iparams[1].lj = { 0.0, 0.0 }; // 0-1
317 idef.iparams[2].lj = { 0.001458, 1.0062882e-6 }; // 0-2
318 idef.iparams[3].lj = { 0.0, 0.0 }; // 1-0
319 idef.iparams[4].lj = { 0.0, 0.0 }; // 1-1
320 idef.iparams[5].lj = { 0.0, 0.0 }; // 1-2
321 idef.iparams[6].lj = { 0.001458, 1.0062882e-6 }; // 2-0
322 idef.iparams[7].lj = { 0.0, 0.0 }; // 2-1
323 idef.iparams[8].lj = { 0.001458, 1.0062882e-6 }; // 2-2
325 GMX_ASSERT(chargeA.size() == c_numAtoms, "This test wants 4 atoms");
326 GMX_ASSERT(chargeB.size() == c_numAtoms, "This test wants 4 atoms");
327 GMX_ASSERT(typeA.size() == c_numAtoms, "This test wants 4 atoms");
328 GMX_ASSERT(typeB.size() == c_numAtoms, "This test wants 4 atoms");
331 // forcefield parameters
335 std::vector<real> chargeA = { 1.0, -1.0, -1.0, 1.0 };
336 std::vector<real> chargeB = { 1.0, 0.0, 0.0, 1.0 };
337 std::vector<int> typeA = { 0, 0, 0, 0 };
338 std::vector<int> typeB = { 0, 1, 2, 1 };
339 // perturbation pattern: {no-pert, vdw- and coul-pert, coul-pert, vdw-pert}
341 // neighbourhood information
342 std::vector<int> iAtoms = { 0 };
343 std::vector<int> jAtoms = { 0, 1, 2, 3 };
344 std::vector<int> jIndex = { 0, 4 };
345 std::vector<int> shift = { 0 };
346 std::vector<int> gid = { 0 };
347 std::vector<int> exclFep = { 0, 1, 1, 1 };
349 // construct t_nblist
360 nbl.excl_fep = exclFep;
365 /*! \brief Input structure for nonbonded fep kernel
372 //! Tolerance for float evaluation
373 float floatToler = 1e-6;
374 //! Tolerance for double evaluation
375 double doubleToler = 1e-8;
379 ForcerecHelper frHelper;
384 /*! \brief Constructor with tolerance
386 * \param[in] ftol Single precision tolerance
387 * \param[in] dtol Double precision tolerance
389 ListInput(float ftol, double dtol)
395 /*! \brief Set parameters for nonbonded interaction
397 * \param[in] coulType coulomb type
398 * \param[in] vdwType vdw type
399 * \param[in] vdwMod vdw potential modifier
401 ListInput setInteraction(CoulombInteractionType coulType, VanDerWaalsType vdwType, InteractionModifiers vdwMod)
403 frHelper.initForcerec(atoms.idef, coulType, vdwType, vdwMod);
408 class NonbondedFepTest :
409 public ::testing::TestWithParam<std::tuple<SoftcoreType, ListInput, PaddedVector<RVec>, real, real, bool>>
412 PaddedVector<RVec> x_;
416 bool softcoreCoulomb_;
417 SoftcoreType softcoreType_;
418 TestReferenceData refData_;
419 TestReferenceChecker checker_;
421 NonbondedFepTest() : checker_(refData_.rootChecker())
423 softcoreType_ = std::get<0>(GetParam());
424 input_ = std::get<1>(GetParam());
425 x_ = std::get<2>(GetParam());
426 lambda_ = std::get<3>(GetParam());
427 softcoreAlpha_ = std::get<4>(GetParam());
428 softcoreCoulomb_ = std::get<5>(GetParam());
430 // Note that the reference data for Ewald type interactions has been generated
431 // with accurate analytical approximations for the long-range corrections.
432 // When the free-energy kernel switches from tabulated to analytical corrections,
433 // the double precision tolerance can be tightend to 1e-11.
434 test::FloatingPointTolerance tolerance(
435 input_.floatToler, input_.doubleToler, 1.0e-6, 1.0e-11, 10000, 100, false);
436 checker_.setDefaultTolerance(tolerance);
441 input_.frHelper.setSoftcoreAlpha(softcoreAlpha_);
442 input_.frHelper.setSoftcoreCoulomb(softcoreCoulomb_);
443 input_.frHelper.setSoftcoreType(softcoreType_);
445 // get forcerec and interaction_const
447 input_.frHelper.getForcerec(&fr);
450 t_nblist nbl = input_.atoms.getNbList();
453 OutputQuantities output;
456 int numFepCouplingTerms = static_cast<int>(FreeEnergyPerturbationCouplingType::Count);
457 std::vector<real> lambdas(numFepCouplingTerms, lambda_);
461 doNBFlags |= GMX_NONBONDED_DO_FORCE;
462 doNBFlags |= GMX_NONBONDED_DO_SHIFTFORCE;
463 doNBFlags |= GMX_NONBONDED_DO_POTENTIAL;
466 bool unusedBool = true; // this bool has no effect in the kernel
467 gmx::ForceWithShiftForces forces(output.f.arrayRefWithPadding(), unusedBool, output.fShift);
473 gmx_nb_free_energy_kernel(nbl,
474 x_.arrayRefWithPadding(),
482 input_.atoms.chargeA,
483 input_.atoms.chargeB,
489 output.f.arrayRefWithPadding(),
490 as_rvec_array(output.fShift.data()),
491 output.energy.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR],
492 output.energy.energyGroupPairTerms[NonBondedEnergyTerms::LJSR],
495 checkOutput(&checker_, output);
499 TEST_P(NonbondedFepTest, testKernel)
504 //! configurations to test
505 std::vector<ListInput> c_interaction = {
506 { ListInput(1e-6, 1e-8).setInteraction(CoulombInteractionType::Cut, VanDerWaalsType::Cut, InteractionModifiers::None) },
507 { ListInput(1e-6, 1e-8).setInteraction(CoulombInteractionType::Cut, VanDerWaalsType::Cut, InteractionModifiers::PotSwitch) },
508 { ListInput(1e-6, 1e-8).setInteraction(CoulombInteractionType::Pme, VanDerWaalsType::Pme, InteractionModifiers::None) }
512 std::vector<real> c_fepLambdas = { 0.0, 0.5, 1.0 };
513 std::vector<real> c_softcoreBeutlerAlphaOrGapsysLinpointScaling = { 0.0, 0.3 };
514 std::vector<bool> c_softcoreCoulomb = { true, false };
515 std::vector<SoftcoreType> c_softcoreType = { SoftcoreType::Beutler, SoftcoreType::Gapsys };
517 //! Coordinates for testing
518 std::vector<PaddedVector<RVec>> c_coordinates = {
519 { { 1.0, 1.0, 1.0 }, { 1.1, 1.15, 1.2 }, { 0.9, 0.85, 0.8 }, { 1.1, 1.15, 0.8 } }
522 INSTANTIATE_TEST_SUITE_P(NBInteraction,
524 ::testing::Combine(::testing::ValuesIn(c_softcoreType),
525 ::testing::ValuesIn(c_interaction),
526 ::testing::ValuesIn(c_coordinates),
527 ::testing::ValuesIn(c_fepLambdas),
528 ::testing::ValuesIn(c_softcoreBeutlerAlphaOrGapsysLinpointScaling),
529 ::testing::ValuesIn(c_softcoreCoulomb)));