api/nblib/listed_forces/tests/calculator.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020,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 /*! \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         HarmonicAngle                                angle(397.5, Degrees(108.53));
 108         std::vector<HarmonicAngle>                   angles{ angle };
 109         std::vector<InteractionIndex<HarmonicAngle>> angleIndices{ { 0, 1, 2, 0 } };
 110
 111         pickType<HarmonicBondType>(interactions).indices    = bondIndices;
 112         pickType<HarmonicBondType>(interactions).parameters = bonds;
 113
 114         pickType<HarmonicAngle>(interactions).indices    = angleIndices;
 115         pickType<HarmonicAngle>(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         box.reset(new Box(3, 3, 3));
 122         pbc.reset(new PbcHolder(PbcType::Xyz, *box));
 123     }
 124
 125     std::vector<gmx::RVec> x;
 126     std::vector<gmx::RVec> forces;
 127
 128     ListedInteractionData interactions;
 129
 130     std::shared_ptr<Box>       box;
 131     std::shared_ptr<PbcHolder> pbc;
 132 };
 133
 134 TEST_F(ListedExampleData, ComputeHarmonicBondForces)
 135 {
 136     gmx::ArrayRef<const InteractionIndex<HarmonicBondType>> indices =
 137             pickType<HarmonicBondType>(interactions).indices;
 138     gmx::ArrayRef<const HarmonicBondType> bonds = pickType<HarmonicBondType>(interactions).parameters;
 139     computeForces(indices, bonds, x, &forces, *pbc);
 140
 141     RefDataChecker vector3DTest(1e-3);
 142     vector3DTest.testArrays<Vec3>(forces, "Bond forces");
 143 }
 144
 145 TEST_F(ListedExampleData, ComputeHarmonicBondEnergies)
 146 {
 147     gmx::ArrayRef<const InteractionIndex<HarmonicBondType>> indices =
 148             pickType<HarmonicBondType>(interactions).indices;
 149     gmx::ArrayRef<const HarmonicBondType> bonds = pickType<HarmonicBondType>(interactions).parameters;
 150     real                                  energy = computeForces(indices, bonds, x, &forces, *pbc);
 151
 152     RefDataChecker vector3DTest(1e-4);
 153     vector3DTest.testReal(energy, "Bond energy");
 154 }
 155
 156 TEST_F(ListedExampleData, ComputeHarmonicAngleForces)
 157 {
 158     gmx::ArrayRef<const InteractionIndex<HarmonicAngle>> indices =
 159             pickType<HarmonicAngle>(interactions).indices;
 160     gmx::ArrayRef<const HarmonicAngle> angles = pickType<HarmonicAngle>(interactions).parameters;
 161     computeForces(indices, angles, x, &forces, *pbc);
 162
 163     RefDataChecker vector3DTest(1e-4);
 164     vector3DTest.testArrays<Vec3>(forces, "Angle forces");
 165 }
 166
 167 TEST_F(ListedExampleData, CanReduceForces)
 168 {
 169     reduceListedForces(interactions, x, &forces, *pbc);
 170
 171     RefDataChecker vector3DTest(1e-2);
 172     vector3DTest.testArrays<Vec3>(forces, "Reduced forces");
 173 }
 174
 175 TEST_F(ListedExampleData, CanReduceEnergies)
 176 {
 177     auto energies    = reduceListedForces(interactions, x, &forces, *pbc);
 178     real totalEnergy = std::accumulate(begin(energies), end(energies), 0.0);
 179
 180     RefDataChecker vector3DTest(1e-4);
 181     vector3DTest.testReal(totalEnergy, "Reduced energy");
 182 }
 183
 184
 185 void compareArray(const ListedForceCalculator::EnergyType& energies,
 186                   const ListedForceCalculator::EnergyType& refEnergies)
 187 {
 188     for (size_t i = 0; i < energies.size(); ++i)
 189     {
 190         EXPECT_REAL_EQ_TOL(energies[i],
 191                            refEnergies[i],
 192                            gmx::test::relativeToleranceAsFloatingPoint(refEnergies[i], 1e-5));
 193     }
 194 }
 195
 196
 197 //! \brief sets up an interaction tuple for a linear chain with nParticles
 198 class LinearChainDataFixture : public ::testing::Test
 199 {
 200 protected:
 201     void SetUp() override
 202     {
 203         LinearChainData data(430);
 204
 205         x            = data.x;
 206         interactions = data.interactions;
 207         box          = data.box;
 208         refForces    = data.forces;
 209         // pbc.reset(new PbcHolder(*box));
 210
 211         refEnergies = reduceListedForces(interactions, x, &refForces, NoPbc{});
 212     }
 213
 214     void testEnergies(const ListedForceCalculator::EnergyType& energies) const
 215     {
 216         compareArray(energies, refEnergies);
 217     }
 218
 219     void testForces(const std::vector<gmx::RVec>& forces) const
 220     {
 221         compareVectors(forces, refForces, refForces);
 222     }
 223
 224     std::vector<gmx::RVec> x;
 225     ListedInteractionData  interactions;
 226     std::shared_ptr<Box>   box;
 227     // std::shared_ptr<PbcHolder> pbc;
 228
 229 private:
 230     std::vector<gmx::RVec>            refForces;
 231     ListedForceCalculator::EnergyType refEnergies;
 232 };
 233
 234 TEST_F(LinearChainDataFixture, Multithreading)
 235 {
 236     ListedForceCalculator lfCalculator(interactions, x.size(), 4, *box);
 237
 238     std::vector<Vec3>                 forces(x.size(), Vec3{ 0, 0, 0 });
 239     ListedForceCalculator::EnergyType energies;
 240     lfCalculator.compute(x, forces, energies);
 241
 242     testEnergies(energies);
 243     testForces(forces);
 244 }
 245
 246
 247 } // namespace
 248 } // namespace test
 249 } // namespace nblib