api/nblib/listed_forces/tests/calculator.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 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.
   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 /*! \internal \file
  36  * \brief
  37  * This implements basic nblib utility tests
  38  *
  39  * \author Victor Holanda <victor.holanda@cscs.ch>
  40  * \author Joe Jordan <ejjordan@kth.se>
  41  * \author Prashanth Kanduri <kanduri@cscs.ch>
  42  * \author Sebastian Keller <keller@cscs.ch>
  43  */
  44 #include "gmxpre.h"
  45
  46 // includesorter/check-source want this include here. IMO a bug
  47
  48
  49 #include <valarray>
  50
  51 #include <gtest/gtest.h>
  52
  53 #include "nblib/listed_forces/calculator.h"
  54 #include "nblib/listed_forces/dataflow.hpp"
  55 #include "nblib/listed_forces/tests/linear_chain_input.hpp"
  56 #include "nblib/tests/testhelpers.h"
  57
  58 #include "testutils/refdata.h"
  59 #include "testutils/testasserts.h"
  60
  61
  62 namespace nblib
  63 {
  64 namespace test
  65 {
  66 namespace
  67 {
  68
  69 TEST(NBlibTest, ListedForceCalculatorCanConstruct)
  70 {
  71     ListedInteractionData interactions;
  72     Box                   box(1, 1, 1);
  73     EXPECT_NO_THROW(ListedForceCalculator listedForceCalculator(interactions, 2, 1, box));
  74 }
  75
  76 template<class TestSeq, class SeqFloat, class SeqDouble>
  77 void compareVectors(const TestSeq&                    forces,
  78                     [[maybe_unused]] const SeqFloat&  refForcesFloat,
  79                     [[maybe_unused]] const SeqDouble& refForcesDouble)
  80 {
  81     for (size_t i = 0; i < forces.size(); ++i)
  82     {
  83         for (int m = 0; m < dimSize; ++m)
  84         {
  85             EXPECT_FLOAT_DOUBLE_EQ_TOL(
  86                     forces[i][m],
  87                     refForcesFloat[i][m],
  88                     refForcesDouble[i][m],
  89                     // Todo: why does the tolerance need to be so low?
  90                     gmx::test::relativeToleranceAsFloatingPoint(refForcesDouble[i][m], 5e-5));
  91         }
  92     }
  93 }
  94
  95 class ListedExampleData : public ::testing::Test
  96 {
  97 protected:
  98     void SetUp() override
  99     {
 100         // methanol-spc data
 101         HarmonicBondType              bond1{ 376560, 0.136 };
 102         HarmonicBondType              bond2{ 313800, 0.1 };
 103         std::vector<HarmonicBondType> bonds{ bond1, bond2 };
 104         // one bond between atoms 0-1 with bond1 parameters and another between atoms 1-2 with bond2 parameters
 105         std::vector<InteractionIndex<HarmonicBondType>> bondIndices{ { 0, 1, 0 }, { 1, 2, 1 } };
 106
 107         DefaultAngle                                angle(Degrees(108.53), 397.5);
 108         std::vector<DefaultAngle>                   angles{ angle };
 109         std::vector<InteractionIndex<DefaultAngle>> angleIndices{ { 0, 1, 2, 0 } };
 110
 111         pickType<HarmonicBondType>(interactions).indices    = bondIndices;
 112         pickType<HarmonicBondType>(interactions).parameters = bonds;
 113
 114         pickType<DefaultAngle>(interactions).indices    = angleIndices;
 115         pickType<DefaultAngle>(interactions).parameters = angles;
 116
 117         // initial position for the methanol atoms from the spc-water example
 118         x = std::vector<gmx::RVec>{ { 1.97, 1.46, 1.209 }, { 1.978, 1.415, 1.082 }, { 1.905, 1.46, 1.03 } };
 119         forces = std::vector<gmx::RVec>(3, gmx::RVec{ 0, 0, 0 });
 120
 121         refBondForcesFloat =
 122                 std::valarray<gmx::BasicVector<float>>{ { -22.8980637, 128.801575, 363.505951 },
 123                                                         { -43.2698593, -88.0130997, -410.639252 },
 124                                                         { 66.167923, -40.788475, 47.1333084 } };
 125         refAngleForcesFloat =
 126                 std::valarray<gmx::BasicVector<float>>{ { 54.7276611, -40.1688995, 17.6805191 },
 127                                                         { -81.8118973, 86.1988525, 60.1752243 },
 128                                                         { 27.0842342, -46.0299492, -77.8557434 } };
 129
 130         refBondForcesDouble = std::valarray<gmx::BasicVector<double>>{
 131             { -22.89764839974935, 128.79927224858977, 363.50016834602064 },
 132             { -43.24622441913251, -88.025652017772231, -410.61635172385434 },
 133             { 66.14387281888186, -40.773620230817542, 47.116183377833721 }
 134         };
 135         refAngleForcesDouble = std::valarray<gmx::BasicVector<double>>{
 136             { 54.726206806506234, -40.167809526198099, 17.680008528590257 },
 137             { -81.809781666748606, 86.196545126117257, 60.173723525141448 },
 138             { 27.083574860242372, -46.028735599919159, -77.853732053731704 }
 139         };
 140
 141         refBondEnergyFloat  = 0.2113433;
 142         refAngleEnergyFloat = 0.112774156;
 143
 144         refBondEnergyDouble  = 0.2113273434867636;
 145         refAngleEnergyDouble = 0.11276812148357591;
 146
 147         box.reset(new Box(3, 3, 3));
 148         pbc.reset(new PbcHolder(*box));
 149     }
 150
 151     std::vector<gmx::RVec> x;
 152     std::vector<gmx::RVec> forces;
 153
 154     ListedInteractionData interactions;
 155
 156     std::shared_ptr<Box>       box;
 157     std::shared_ptr<PbcHolder> pbc;
 158
 159     // reference values
 160     std::valarray<gmx::BasicVector<float>>  refBondForcesFloat, refAngleForcesFloat;
 161     std::valarray<gmx::BasicVector<double>> refBondForcesDouble, refAngleForcesDouble;
 162
 163     float  refBondEnergyFloat, refAngleEnergyFloat;
 164     double refBondEnergyDouble, refAngleEnergyDouble;
 165 };
 166
 167 TEST_F(ListedExampleData, DISABLED_ComputeHarmonicBondForces)
 168 {
 169     auto indices = pickType<HarmonicBondType>(interactions).indices;
 170     auto bonds   = pickType<HarmonicBondType>(interactions).parameters;
 171     real energy  = computeForces(indices, bonds, x, &forces, *pbc);
 172
 173     EXPECT_FLOAT_DOUBLE_EQ_TOL(energy,
 174                                refBondEnergyFloat,
 175                                refBondEnergyDouble,
 176                                gmx::test::relativeToleranceAsFloatingPoint(refBondEnergyDouble, 1e-5));
 177
 178     compareVectors(forces, refBondForcesFloat, refBondForcesDouble);
 179 }
 180
 181 TEST_F(ListedExampleData, ComputeHarmonicAngleForces)
 182 {
 183     auto indices = pickType<DefaultAngle>(interactions).indices;
 184     auto angles  = pickType<DefaultAngle>(interactions).parameters;
 185     real energy  = computeForces(indices, angles, x, &forces, *pbc);
 186
 187     EXPECT_FLOAT_DOUBLE_EQ_TOL(energy,
 188                                refAngleEnergyFloat,
 189                                refAngleEnergyDouble,
 190                                gmx::test::relativeToleranceAsFloatingPoint(refAngleEnergyDouble, 1e-5));
 191
 192     compareVectors(forces, refAngleForcesFloat, refAngleForcesDouble);
 193 }
 194
 195 TEST_F(ListedExampleData, DISABLED_CanReduceForces)
 196 {
 197     auto energies    = reduceListedForces(interactions, x, &forces, *pbc);
 198     real totalEnergy = std::accumulate(begin(energies), end(energies), 0.0);
 199
 200     EXPECT_FLOAT_DOUBLE_EQ_TOL(totalEnergy,
 201                                refBondEnergyFloat + refAngleEnergyFloat,
 202                                refBondEnergyDouble + refAngleEnergyDouble,
 203                                gmx::test::relativeToleranceAsFloatingPoint(refBondEnergyDouble, 1e-5));
 204
 205     compareVectors(forces, refBondForcesFloat + refAngleForcesFloat, refBondForcesDouble + refAngleForcesDouble);
 206 }
 207
 208
 209 void compareArray(const ListedForceCalculator::EnergyType& energies,
 210                   const ListedForceCalculator::EnergyType& refEnergies)
 211 {
 212     for (size_t i = 0; i < energies.size(); ++i)
 213     {
 214         EXPECT_REAL_EQ_TOL(energies[i],
 215                            refEnergies[i],
 216                            gmx::test::relativeToleranceAsFloatingPoint(refEnergies[i], 1e-5));
 217     }
 218 }
 219
 220
 221 //! \brief sets up an interaction tuple for a linear chain with nParticles
 222 class LinearChainDataFixture : public ::testing::Test
 223 {
 224 protected:
 225     void SetUp() override
 226     {
 227         LinearChainData data(430);
 228
 229         x            = data.x;
 230         interactions = data.interactions;
 231         box          = data.box;
 232         refForces    = data.forces;
 233         // pbc.reset(new PbcHolder(*box));
 234
 235         refEnergies = reduceListedForces(interactions, x, &refForces, NoPbc{});
 236     }
 237
 238     void testEnergies(const ListedForceCalculator::EnergyType& energies) const
 239     {
 240         compareArray(energies, refEnergies);
 241     }
 242
 243     void testForces(const std::vector<gmx::RVec>& forces) const
 244     {
 245         compareVectors(forces, refForces, refForces);
 246     }
 247
 248     std::vector<gmx::RVec> x;
 249     ListedInteractionData  interactions;
 250     std::shared_ptr<Box>   box;
 251     // std::shared_ptr<PbcHolder> pbc;
 252
 253 private:
 254     std::vector<gmx::RVec>            refForces;
 255     ListedForceCalculator::EnergyType refEnergies;
 256 };
 257
 258 TEST_F(LinearChainDataFixture, Multithreading)
 259 {
 260     ListedForceCalculator lfCalculator(interactions, x.size(), 4, *box);
 261
 262     std::vector<Vec3>                 forces(x.size(), Vec3{ 0, 0, 0 });
 263     ListedForceCalculator::EnergyType energies;
 264     lfCalculator.compute(x, forces, energies);
 265
 266     testEnergies(energies);
 267     testForces(forces);
 268 }
 269
 270
 271 } // namespace
 272 } // namespace test
 273 } // namespace nblib