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