src/gromacs/gmxlib/nonbonded/tests/nb_free_energy.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35
  36 /*! \internal \file
  37  * \brief Implements test of nonbonded fep kernel
  38  *
  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.
  42  *
  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.
  46  *
  47  * \author Sebastian Kehl <sebastian.kehl@mpcdf.mpg.de>
  48  * \ingroup module_gmxlib_nonbonded
  49  */
  50 #include "gmxpre.h"
  51
  52 #include "gromacs/gmxlib/nonbonded/nb_free_energy.h"
  53 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
  54
  55 #include <cmath>
  56
  57 #include <gtest/gtest.h>
  58
  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"
  85
  86 #include "testutils/refdata.h"
  87 #include "testutils/testasserts.h"
  88
  89 namespace gmx
  90 {
  91 namespace test
  92 {
  93 namespace
  94 {
  95
  96 //! Number of atoms used in these tests.
  97 constexpr int c_numAtoms     = 4;
  98 constexpr int c_numAtomTypes = 3;
  99
 100 /*! \brief Output from nonbonded fep kernel
 101  *
 102  */
 103 struct OutputQuantities
 104 {
 105     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 })
 110     {
 111     }
 112
 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;
 121 };
 122
 123 /*! \brief Utility to check the output from nonbonded test
 124  *
 125  * \param[in] checker Reference checker
 126  * \param[in] output  The output from the test to check
 127  */
 128 void checkOutput(TestReferenceChecker* checker, const OutputQuantities& output)
 129 {
 130     checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::LJSR][0], "EVdw ");
 131     checker->checkReal(output.energy.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR][0],
 132                        "ECoul ");
 133     checker->checkReal(output.dvdLambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)],
 134                        "dVdlCoul ");
 135     checker->checkReal(output.dvdLambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)],
 136                        "dVdlVdw ");
 137
 138     checker->checkSequence(std::begin(output.f), std::end(output.f), "Forces");
 139
 140     auto shiftForcesChecker = checker->checkCompound("Shift-Forces", "Shift-forces");
 141     shiftForcesChecker.checkVector(output.fShift[0], "Central");
 142 }
 143
 144 class InteractionConstHelper
 145 {
 146 public:
 147     InteractionConstHelper() {}
 148
 149     //! init data to construct interaction_const
 150     void initInteractionConst(CoulombInteractionType coulType, VanDerWaalsType vdwType, InteractionModifiers vdwMod)
 151     {
 152         coulType_ = coulType;
 153         vdwType_  = vdwType;
 154         vdwMod_   = vdwMod;
 155
 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>();
 166
 167         init_interaction_const_tables(nullptr, &tmp, 1.0, 0.0);
 168         coulombTables_ = *tmp.coulombEwaldTables;
 169         vdwTables_     = *tmp.vdwEwaldTables;
 170     }
 171
 172     /*! \brief Setup interaction_const_t
 173      *
 174      * \param[in]  fepVals t_lambda struct of fep values
 175      * \param[out] ic      interaction_const_t pointer with data
 176      */
 177     void getInteractionConst(const t_lambda& fepVals, interaction_const_t* ic)
 178     {
 179         ic->softCoreParameters = std::unique_ptr<interaction_const_t::SoftCoreParameters>(
 180                 new interaction_const_t::SoftCoreParameters(fepVals));
 181
 182         ic->coulombEwaldTables  = std::unique_ptr<EwaldCorrectionTables>(new EwaldCorrectionTables);
 183         *ic->coulombEwaldTables = coulombTables_;
 184
 185         ic->vdwEwaldTables  = std::unique_ptr<EwaldCorrectionTables>(new EwaldCorrectionTables);
 186         *ic->vdwEwaldTables = vdwTables_;
 187
 188         // set coulomb and vdw types
 189         ic->eeltype      = coulType_;
 190         ic->vdwtype      = vdwType_;
 191         ic->vdw_modifier = vdwMod_;
 192
 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;
 203     }
 204
 205 private:
 206     //! correction tables
 207     EwaldCorrectionTables coulombTables_;
 208     EwaldCorrectionTables vdwTables_;
 209
 210     //! coulomb and vdw type specifiers
 211     CoulombInteractionType coulType_;
 212     real                   coulEwaldCoeff_;
 213     VanDerWaalsType        vdwType_;
 214     InteractionModifiers   vdwMod_;
 215     real                   vdwEwaldCoeff_;
 216 };
 217
 218
 219 /* \brief Utility class to setup forcerec
 220  *
 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.
 224  *
 225  * Data is only initialized as necessary for the nonbonded kernel to work!
 226  */
 227 class ForcerecHelper
 228 {
 229 public:
 230     ForcerecHelper()
 231     {
 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_.scScaleLinpointLJGapsys = 0.3;
 239         fepVals_.scScaleLinpointQGapsys  = 0.3;
 240         fepVals_.scSigmaLJGapsys         = 0.3;
 241         fepVals_.softcoreFunction        = SoftcoreType::Beutler;
 242     }
 243
 244     //! initialize data structure to construct forcerec
 245     void initForcerec(const gmx_ffparams_t&  idef,
 246                       CoulombInteractionType coulType,
 247                       VanDerWaalsType        vdwType,
 248                       InteractionModifiers   vdwMod)
 249     {
 250         icHelper_.initInteractionConst(coulType, vdwType, vdwMod);
 251         nbfp_        = makeNonBondedParameterLists(idef.atnr, idef.iparams, false);
 252         ljPmeC6Grid_ = makeLJPmeC6GridCorrectionParameters(idef.atnr, idef.iparams, LongRangeVdW::Geom);
 253     }
 254
 255     void setSoftcoreAlpha(const real scAlpha) { fepVals_.sc_alpha = scAlpha; }
 256     void setSoftcoreCoulomb(const bool scCoulomb) { fepVals_.bScCoul = scCoulomb; }
 257
 258     //! get forcerec data as wanted by the nonbonded kernel
 259     void getForcerec(t_forcerec* fr)
 260     {
 261         fr->ic = std::make_unique<interaction_const_t>();
 262
 263         // set data in ic
 264         icHelper_.getInteractionConst(fepVals_, fr->ic.get());
 265
 266         // set data in fr
 267         fr->ljpme_c6grid = ljPmeC6Grid_;
 268         fr->nbfp         = nbfp_;
 269         fr->shift_vec    = { { 0.0, 0.0, 0.0 } };
 270         fr->rlist        = rlist_;
 271         fr->ntype        = c_numAtomTypes;
 272
 273         // simd
 274         fr->use_simd_kernels = GMX_USE_SIMD_KERNELS;
 275     }
 276
 277 private:
 278     InteractionConstHelper icHelper_;
 279     std::vector<real>      ljPmeC6Grid_;
 280     std::vector<real>      nbfp_;
 281     t_lambda               fepVals_;
 282     real                   rlist_ = 1.0;
 283 };
 284
 285 /*! \brief Utility structure to hold atoms data
 286  *
 287  * A system having 4 interactions with the following perturbation pattern:
 288  *  - no perturbation
 289  *  - vdw- and coulomb-perturbation
 290  *  - coulomb-perturbation only
 291  *  - vdw-perturbation only
 292  *
 293  * This is realized by defining 3 different atom types that control
 294  * the vdw-perturbation. The coulomb-perturbation is controlled by directly
 295  * setting the charge of the atoms at the lambda states A/B.
 296  */
 297 struct AtomData
 298 {
 299     AtomData()
 300     {
 301         idef.atnr = c_numAtomTypes;
 302         idef.iparams.resize(2 * c_numAtomTypes * c_numAtomTypes);
 303
 304         // set interaction parameters for different combinations of types
 305         idef.iparams[0].lj = { 0.001458, 1.0062882e-6 }; // 0-0
 306         idef.iparams[1].lj = { 0.0, 0.0 };               // 0-1
 307         idef.iparams[2].lj = { 0.001458, 1.0062882e-6 }; // 0-2
 308         idef.iparams[3].lj = { 0.0, 0.0 };               // 1-0
 309         idef.iparams[4].lj = { 0.0, 0.0 };               // 1-1
 310         idef.iparams[5].lj = { 0.0, 0.0 };               // 1-2
 311         idef.iparams[6].lj = { 0.001458, 1.0062882e-6 }; // 2-0
 312         idef.iparams[7].lj = { 0.0, 0.0 };               // 2-1
 313         idef.iparams[8].lj = { 0.001458, 1.0062882e-6 }; // 2-2
 314
 315         GMX_ASSERT(chargeA.size() == c_numAtoms, "This test wants 4 atoms");
 316         GMX_ASSERT(chargeB.size() == c_numAtoms, "This test wants 4 atoms");
 317         GMX_ASSERT(typeA.size() == c_numAtoms, "This test wants 4 atoms");
 318         GMX_ASSERT(typeB.size() == c_numAtoms, "This test wants 4 atoms");
 319     }
 320
 321     // forcefield parameters
 322     gmx_ffparams_t idef;
 323
 324     // atom data
 325     std::vector<real> chargeA = { 1.0, -1.0, -1.0, 1.0 };
 326     std::vector<real> chargeB = { 1.0, 0.0, 0.0, 1.0 };
 327     std::vector<int>  typeA   = { 0, 0, 0, 0 };
 328     std::vector<int>  typeB   = { 0, 1, 2, 1 };
 329     // perturbation pattern: {no-pert, vdw- and coul-pert, coul-pert, vdw-pert}
 330
 331     // neighbourhood information
 332     std::vector<int> iAtoms  = { 0 };
 333     std::vector<int> jAtoms  = { 0, 1, 2, 3 };
 334     std::vector<int> jIndex  = { 0, 4 };
 335     std::vector<int> shift   = { 0 };
 336     std::vector<int> gid     = { 0 };
 337     std::vector<int> exclFep = { 0, 1, 1, 1 };
 338
 339     // construct t_nblist
 340     t_nblist getNbList()
 341     {
 342         t_nblist nbl;
 343         nbl.nri      = 1;
 344         nbl.nrj      = 4;
 345         nbl.iinr     = iAtoms;
 346         nbl.jindex   = jIndex;
 347         nbl.jjnr     = jAtoms;
 348         nbl.shift    = shift;
 349         nbl.gid      = gid;
 350         nbl.excl_fep = exclFep;
 351         return nbl;
 352     }
 353 };
 354
 355 /*! \brief Input structure for nonbonded fep kernel
 356  */
 357 struct ListInput
 358 {
 359 public:
 360     //! Function type
 361     int fType = F_LJ;
 362     //! Tolerance for float evaluation
 363     float floatToler = 1e-6;
 364     //! Tolerance for double evaluation
 365     double doubleToler = 1e-8;
 366     //! atom parameters
 367     AtomData atoms;
 368     //! forcerec helper
 369     ForcerecHelper frHelper;
 370
 371     //! Constructor
 372     ListInput() {}
 373
 374     /*! \brief Constructor with tolerance
 375      *
 376      * \param[in] ftol Single precision tolerance
 377      * \param[in] dtol Double precision tolerance
 378      */
 379     ListInput(float ftol, double dtol)
 380     {
 381         floatToler  = ftol;
 382         doubleToler = dtol;
 383     }
 384
 385     /*! \brief Set parameters for nonbonded interaction
 386      *
 387      * \param[in] coulType coulomb type
 388      * \param[in] vdwType  vdw type
 389      * \param[in] vdwMod   vdw potential modifier
 390      */
 391     ListInput setInteraction(CoulombInteractionType coulType, VanDerWaalsType vdwType, InteractionModifiers vdwMod)
 392     {
 393         frHelper.initForcerec(atoms.idef, coulType, vdwType, vdwMod);
 394         return *this;
 395     }
 396 };
 397
 398 class NonbondedFepTest :
 399     public ::testing::TestWithParam<std::tuple<ListInput, PaddedVector<RVec>, real, real, bool>>
 400 {
 401 protected:
 402     PaddedVector<RVec>   x_;
 403     ListInput            input_;
 404     real                 lambda_;
 405     real                 softcoreAlpha_;
 406     bool                 softcoreCoulomb_;
 407     TestReferenceData    refData_;
 408     TestReferenceChecker checker_;
 409
 410     NonbondedFepTest() : checker_(refData_.rootChecker())
 411     {
 412         input_           = std::get<0>(GetParam());
 413         x_               = std::get<1>(GetParam());
 414         lambda_          = std::get<2>(GetParam());
 415         softcoreAlpha_   = std::get<3>(GetParam());
 416         softcoreCoulomb_ = std::get<4>(GetParam());
 417
 418         // Note that the reference data for Ewald type interactions has been generated
 419         // with accurate analytical approximations for the long-range corrections.
 420         // When the free-energy kernel switches from tabulated to analytical corrections,
 421         // the double precision tolerance can be tightend to 1e-11.
 422         test::FloatingPointTolerance tolerance(
 423                 input_.floatToler, input_.doubleToler, 1.0e-6, 1.0e-11, 10000, 100, false);
 424         checker_.setDefaultTolerance(tolerance);
 425     }
 426
 427     void testKernel()
 428     {
 429         input_.frHelper.setSoftcoreAlpha(softcoreAlpha_);
 430         input_.frHelper.setSoftcoreCoulomb(softcoreCoulomb_);
 431
 432         // get forcerec and interaction_const
 433         t_forcerec fr;
 434         input_.frHelper.getForcerec(&fr);
 435
 436         // t_nblist
 437         t_nblist nbl = input_.atoms.getNbList();
 438
 439         // output buffers
 440         OutputQuantities output;
 441
 442         // lambda vector
 443         int numFepCouplingTerms = static_cast<int>(FreeEnergyPerturbationCouplingType::Count);
 444         std::vector<real> lambdas(numFepCouplingTerms, lambda_);
 445
 446         // fep kernel data
 447         int doNBFlags = 0;
 448         doNBFlags |= GMX_NONBONDED_DO_FORCE;
 449         doNBFlags |= GMX_NONBONDED_DO_SHIFTFORCE;
 450         doNBFlags |= GMX_NONBONDED_DO_POTENTIAL;
 451
 452         // force buffers
 453         bool                      unusedBool = true; // this bool has no effect in the kernel
 454         gmx::ForceWithShiftForces forces(output.f.arrayRefWithPadding(), unusedBool, output.fShift);
 455
 456         // dummy counter
 457         t_nrnb nrnb;
 458
 459         // run fep kernel
 460         gmx_nb_free_energy_kernel(nbl,
 461                                   x_.arrayRefWithPadding(),
 462                                   fr.use_simd_kernels,
 463                                   fr.ntype,
 464                                   fr.rlist,
 465                                   *fr.ic,
 466                                   fr.shift_vec,
 467                                   fr.nbfp,
 468                                   fr.ljpme_c6grid,
 469                                   input_.atoms.chargeA,
 470                                   input_.atoms.chargeB,
 471                                   input_.atoms.typeA,
 472                                   input_.atoms.typeB,
 473                                   doNBFlags,
 474                                   lambdas,
 475                                   &nrnb,
 476                                   output.f.arrayRefWithPadding(),
 477                                   as_rvec_array(output.fShift.data()),
 478                                   output.energy.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR],
 479                                   output.energy.energyGroupPairTerms[NonBondedEnergyTerms::LJSR],
 480                                   output.dvdLambda);
 481
 482         checkOutput(&checker_, output);
 483     }
 484 };
 485
 486 TEST_P(NonbondedFepTest, testKernel)
 487 {
 488     testKernel();
 489 }
 490
 491 //! configurations to test
 492 std::vector<ListInput> c_interaction = {
 493     { ListInput(1e-6, 1e-8).setInteraction(CoulombInteractionType::Cut, VanDerWaalsType::Cut, InteractionModifiers::None) },
 494     { ListInput(1e-6, 1e-8).setInteraction(CoulombInteractionType::Cut, VanDerWaalsType::Cut, InteractionModifiers::PotSwitch) },
 495     { ListInput(1e-6, 1e-8).setInteraction(CoulombInteractionType::Pme, VanDerWaalsType::Pme, InteractionModifiers::None) }
 496 };
 497
 498 //! test parameters
 499 std::vector<real> c_fepLambdas      = { 0.0, 0.5, 1.0 };
 500 std::vector<real> c_softcoreAlphas  = { 0.0, 0.3 };
 501 std::vector<bool> c_softcoreCoulomb = { true, false };
 502
 503 //! Coordinates for testing
 504 std::vector<PaddedVector<RVec>> c_coordinates = {
 505     { { 1.0, 1.0, 1.0 }, { 1.1, 1.15, 1.2 }, { 0.9, 0.85, 0.8 }, { 1.1, 1.15, 0.8 } }
 506 };
 507
 508 INSTANTIATE_TEST_SUITE_P(NBInteraction,
 509                          NonbondedFepTest,
 510                          ::testing::Combine(::testing::ValuesIn(c_interaction),
 511                                             ::testing::ValuesIn(c_coordinates),
 512                                             ::testing::ValuesIn(c_fepLambdas),
 513                                             ::testing::ValuesIn(c_softcoreAlphas),
 514                                             ::testing::ValuesIn(c_softcoreCoulomb)));
 515
 516 } // namespace
 517
 518 } // namespace test
 519
 520 } // namespace gmx