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], refForcesFloat[i][m], refForcesDouble[i][m],
  87                     // Todo: why does the tolerance need to be so low?
  88                     gmx::test::relativeToleranceAsFloatingPoint(refForcesDouble[i][m], 5e-5));
  89         }
  90     }
  91 }
  92
  93 class ListedExampleData : public ::testing::Test
  94 {
  95 protected:
  96     void SetUp() override
  97     {
  98         // methanol-spc data
  99         HarmonicBondType              bond1{ 376560, 0.136 };
 100         HarmonicBondType              bond2{ 313800, 0.1 };
 101         std::vector<HarmonicBondType> bonds{ bond1, bond2 };
 102         // one bond between atoms 0-1 with bond1 parameters and another between atoms 1-2 with bond2 parameters
 103         std::vector<InteractionIndex<HarmonicBondType>> bondIndices{ { 0, 1, 0 }, { 1, 2, 1 } };
 104
 105         DefaultAngle                                angle(Degrees(108.53), 397.5);
 106         std::vector<DefaultAngle>                   angles{ angle };
 107         std::vector<InteractionIndex<DefaultAngle>> angleIndices{ { 0, 1, 2, 0 } };
 108
 109         pickType<HarmonicBondType>(interactions).indices    = bondIndices;
 110         pickType<HarmonicBondType>(interactions).parameters = bonds;
 111
 112         pickType<DefaultAngle>(interactions).indices    = angleIndices;
 113         pickType<DefaultAngle>(interactions).parameters = angles;
 114
 115         // initial position for the methanol atoms from the spc-water example
 116         x = std::vector<gmx::RVec>{ { 1.97, 1.46, 1.209 }, { 1.978, 1.415, 1.082 }, { 1.905, 1.46, 1.03 } };
 117         forces = std::vector<gmx::RVec>(3, gmx::RVec{ 0, 0, 0 });
 118
 119         refBondForcesFloat =
 120                 std::valarray<gmx::BasicVector<float>>{ { -22.8980637, 128.801575, 363.505951 },
 121                                                         { -43.2698593, -88.0130997, -410.639252 },
 122                                                         { 66.167923, -40.788475, 47.1333084 } };
 123         refAngleForcesFloat =
 124                 std::valarray<gmx::BasicVector<float>>{ { 54.7276611, -40.1688995, 17.6805191 },
 125                                                         { -81.8118973, 86.1988525, 60.1752243 },
 126                                                         { 27.0842342, -46.0299492, -77.8557434 } };
 127
 128         refBondForcesDouble = std::valarray<gmx::BasicVector<double>>{
 129             { -22.89764839974935, 128.79927224858977, 363.50016834602064 },
 130             { -43.24622441913251, -88.025652017772231, -410.61635172385434 },
 131             { 66.14387281888186, -40.773620230817542, 47.116183377833721 }
 132         };
 133         refAngleForcesDouble = std::valarray<gmx::BasicVector<double>>{
 134             { 54.726206806506234, -40.167809526198099, 17.680008528590257 },
 135             { -81.809781666748606, 86.196545126117257, 60.173723525141448 },
 136             { 27.083574860242372, -46.028735599919159, -77.853732053731704 }
 137         };
 138
 139         refBondEnergyFloat  = 0.2113433;
 140         refAngleEnergyFloat = 0.112774156;
 141
 142         refBondEnergyDouble  = 0.2113273434867636;
 143         refAngleEnergyDouble = 0.11276812148357591;
 144
 145         box.reset(new Box(3, 3, 3));
 146         pbc.reset(new PbcHolder(*box));
 147     }
 148
 149     std::vector<gmx::RVec> x;
 150     std::vector<gmx::RVec> forces;
 151
 152     ListedInteractionData interactions;
 153
 154     std::shared_ptr<Box>       box;
 155     std::shared_ptr<PbcHolder> pbc;
 156
 157     // reference values
 158     std::valarray<gmx::BasicVector<float>>  refBondForcesFloat, refAngleForcesFloat;
 159     std::valarray<gmx::BasicVector<double>> refBondForcesDouble, refAngleForcesDouble;
 160
 161     float  refBondEnergyFloat, refAngleEnergyFloat;
 162     double refBondEnergyDouble, refAngleEnergyDouble;
 163 };
 164
 165 TEST_F(ListedExampleData, DISABLED_ComputeHarmonicBondForces)
 166 {
 167     auto indices = pickType<HarmonicBondType>(interactions).indices;
 168     auto bonds   = pickType<HarmonicBondType>(interactions).parameters;
 169     real energy  = computeForces(indices, bonds, x, &forces, *pbc);
 170
 171     EXPECT_FLOAT_DOUBLE_EQ_TOL(energy, refBondEnergyFloat, refBondEnergyDouble,
 172                                gmx::test::relativeToleranceAsFloatingPoint(refBondEnergyDouble, 1e-5));
 173
 174     compareVectors(forces, refBondForcesFloat, refBondForcesDouble);
 175 }
 176
 177 TEST_F(ListedExampleData, ComputeHarmonicAngleForces)
 178 {
 179     auto indices = pickType<DefaultAngle>(interactions).indices;
 180     auto angles  = pickType<DefaultAngle>(interactions).parameters;
 181     real energy  = computeForces(indices, angles, x, &forces, *pbc);
 182
 183     EXPECT_FLOAT_DOUBLE_EQ_TOL(energy, refAngleEnergyFloat, refAngleEnergyDouble,
 184                                gmx::test::relativeToleranceAsFloatingPoint(refAngleEnergyDouble, 1e-5));
 185
 186     compareVectors(forces, refAngleForcesFloat, refAngleForcesDouble);
 187 }
 188
 189 TEST_F(ListedExampleData, DISABLED_CanReduceForces)
 190 {
 191     auto energies    = reduceListedForces(interactions, x, &forces, *pbc);
 192     real totalEnergy = std::accumulate(begin(energies), end(energies), 0.0);
 193
 194     EXPECT_FLOAT_DOUBLE_EQ_TOL(totalEnergy, refBondEnergyFloat + refAngleEnergyFloat,
 195                                refBondEnergyDouble + refAngleEnergyDouble,
 196                                gmx::test::relativeToleranceAsFloatingPoint(refBondEnergyDouble, 1e-5));
 197
 198     compareVectors(forces, refBondForcesFloat + refAngleForcesFloat,
 199                    refBondForcesDouble + refAngleForcesDouble);
 200 }
 201
 202
 203 void compareArray(const ListedForceCalculator::EnergyType& energies,
 204                   const ListedForceCalculator::EnergyType& refEnergies)
 205 {
 206     for (size_t i = 0; i < energies.size(); ++i)
 207     {
 208         EXPECT_REAL_EQ_TOL(energies[i], refEnergies[i],
 209                            gmx::test::relativeToleranceAsFloatingPoint(refEnergies[i], 1e-5));
 210     }
 211 }
 212
 213
 214 //! \brief sets up an interaction tuple for a linear chain with nParticles
 215 class LinearChainDataFixture : public ::testing::Test
 216 {
 217 protected:
 218     void SetUp() override
 219     {
 220         LinearChainData data(430);
 221
 222         x            = data.x;
 223         interactions = data.interactions;
 224         box          = data.box;
 225         refForces    = data.forces;
 226         // pbc.reset(new PbcHolder(*box));
 227
 228         refEnergies = reduceListedForces(interactions, x, &refForces, NoPbc{});
 229     }
 230
 231     void testEnergies(const ListedForceCalculator::EnergyType& energies) const
 232     {
 233         compareArray(energies, refEnergies);
 234     }
 235
 236     void testForces(const std::vector<gmx::RVec>& forces) const
 237     {
 238         compareVectors(forces, refForces, refForces);
 239     }
 240
 241     std::vector<gmx::RVec> x;
 242     ListedInteractionData  interactions;
 243     std::shared_ptr<Box>   box;
 244     // std::shared_ptr<PbcHolder> pbc;
 245
 246 private:
 247     std::vector<gmx::RVec>            refForces;
 248     ListedForceCalculator::EnergyType refEnergies;
 249 };
 250
 251 TEST_F(LinearChainDataFixture, Multithreading)
 252 {
 253     ListedForceCalculator lfCalculator(interactions, x.size(), 4, *box);
 254
 255     std::vector<Vec3>                 forces(x.size(), Vec3{ 0, 0, 0 });
 256     ListedForceCalculator::EnergyType energies;
 257     lfCalculator.compute(x, forces, energies);
 258
 259     testEnergies(energies);
 260     testForces(forces);
 261 }
 262
 263
 264 } // namespace
 265 } // namespace test
 266 } // namespace nblib