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