src/gromacs/mdlib/tests/settle.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014,2015,2016, 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 #include "gmxpre.h"
  36
  37 #include <tuple>
  38 #include <vector>
  39
  40 #include <gtest/gtest.h>
  41
  42 #include "gromacs/math/vec.h"
  43 #include "gromacs/math/vectypes.h"
  44 #include "gromacs/mdlib/constr.h"
  45 #include "gromacs/mdtypes/mdatom.h"
  46 #include "gromacs/pbcutil/pbc.h"
  47 #include "gromacs/topology/idef.h"
  48 #include "gromacs/topology/topology.h"
  49 #include "gromacs/utility/scoped_cptr.h"
  50 #include "gromacs/utility/smalloc.h"
  51 #include "gromacs/utility/stringutil.h"
  52
  53 #include "testutils/testasserts.h"
  54
  55 namespace gmx
  56 {
  57
  58 namespace test
  59 {
  60
  61 //! Database of 51 water atom input positions (DIM reals per atom, taken from spc216.gro) for use as test inputs.
  62 const real g_positions[] = {
  63     .130, -.041, -.291,
  64     .120, -.056, -.192,
  65     .044, -.005, -.327,
  66     -.854, -.406, .477,
  67     -.900, -.334, .425,
  68     -.858, -.386, .575,
  69     .351, -.061, .853,
  70     .401, -.147, .859,
  71     .416, .016, .850,
  72     -.067, -.796, .873,
  73     -.129, -.811, .797,
  74     -.119, -.785, .958,
  75     -.635, -.312, -.356,
  76     -.629, -.389, -.292,
  77     -.687, -.338, -.436,
  78     .321, -.919, .242,
  79     .403, -.880, .200,
  80     .294, -1.001, .193,
  81     -.404, .735, .728,
  82     -.409, .670, .803,
  83     -.324, .794, .741,
  84     .461, -.596, -.135,
  85     .411, -.595, -.221,
  86     .398, -.614, -.059,
  87     -.751, -.086, .237,
  88     -.811, -.148, .287,
  89     -.720, -.130, .152,
  90     .202, .285, -.364,
  91     .122, .345, -.377,
  92     .192, .236, -.278,
  93     -.230, -.485, .081,
  94     -.262, -.391, .071,
  95     -.306, -.548, .069,
  96     .464, -.119, .323,
  97     .497, -.080, .409,
  98     .540, -.126, .258,
  99     -.462, .107, .426,
 100     -.486, .070, .336,
 101     -.363, .123, .430,
 102     .249, -.077, -.621,
 103     .306, -.142, -.571,
 104     .233, -.110, -.714,
 105     -.922, -.164, .904,
 106     -.842, -.221, .925,
 107     -.971, -.204, .827,
 108     .382, .700, .480,
 109     .427, .610, .477,
 110     .288, .689, .513,
 111     .781, .264, -.113,
 112     .848, .203, -.070,
 113     .708, .283, -.048
 114 };
 115
 116 //! Simple cubic simulation box to use in tests
 117 matrix g_box = {{1.86206, 0, 0}, {0, 1.86206, 0}, {0, 0, 1.86206}};
 118
 119 //! Convenience typedef
 120 typedef std::tuple<int, bool, bool, bool> SettleTestParameters;
 121
 122 /*! \brief Test fixture for testing SETTLE position updates
 123  *
 124  * \todo This also tests that if the calling code requires velocities
 125  * and virial updates, that those outputs do change, but does not test
 126  * that those changes are correct.
 127  *
 128  * \todo Only no-PBC and cubic-PBC are tested here, but the correct
 129  * function of the SIMD version of set_pbx_auic in all cases should be
 130  * tested elsewhere.
 131  */
 132 class SettleTest : public ::testing::TestWithParam<SettleTestParameters>
 133 {
 134     public:
 135         //! Updated water atom positions to constrain (DIM reals per atom)
 136         std::vector<real> updatedPositions_;
 137         //! Water atom velocities to constrain (DIM reals per atom)
 138         std::vector<real> velocities_;
 139         //! PBC option to test
 140         t_pbc             pbcNone_;
 141         //! PBC option to test
 142         t_pbc             pbcXYZ_;
 143
 144         //! Constructor
 145         SettleTest() :
 146             updatedPositions_(std::begin(g_positions), std::end(g_positions)),
 147             velocities_(updatedPositions_.size(), 0)
 148         {
 149             set_pbc(&pbcNone_, epbcNONE, g_box);
 150             set_pbc(&pbcXYZ_, epbcXYZ, g_box);
 151
 152             // Perturb the atom positions, to appear like an
 153             // "update," and where there is definitely constraining
 154             // work to do.
 155             for (size_t i = 0; i != updatedPositions_.size(); ++i)
 156             {
 157                 if (i % 4 == 0)
 158                 {
 159                     updatedPositions_[i] += 0.01;
 160                 }
 161                 else if (i % 4 == 1)
 162                 {
 163                     updatedPositions_[i] -= 0.01;
 164                 }
 165                 else if (i % 4 == 2)
 166                 {
 167                     updatedPositions_[i] += 0.02;
 168                 }
 169                 else if (i % 4 == 3)
 170                 {
 171                     updatedPositions_[i] -= 0.02;
 172                 }
 173             }
 174         }
 175 };
 176
 177 TEST_P(SettleTest, SatisfiesConstraints)
 178 {
 179     int  numSettles;
 180     bool usePbc, useVelocities, calcVirial;
 181     // Make some symbolic names for the parameter combination under
 182     // test.
 183     std::tie(numSettles, usePbc, useVelocities, calcVirial) = GetParam();
 184
 185     // Make a string that describes which parameter combination is
 186     // being tested, to help make failing tests comprehensible.
 187     std::string testDescription = formatString("while testing %d SETTLEs, %sPBC, %svelocities and %scalculating the virial",
 188                                                numSettles,
 189                                                usePbc ? "with " : "without ",
 190                                                useVelocities ? "with " : "without ",
 191                                                calcVirial ? "" : "not ");
 192
 193     const int settleType     = 0;
 194     const int atomsPerSettle = NRAL(F_SETTLE);
 195     ASSERT_LE(numSettles, updatedPositions_.size() / (atomsPerSettle * DIM)) << "cannot test that many SETTLEs " << testDescription;
 196
 197     // Set up the topology. We still have to make some raw pointers,
 198     // but they are put into scope guards for automatic cleanup.
 199     gmx_mtop_t             *mtop;
 200     snew(mtop, 1);
 201     scoped_cptr<gmx_mtop_t> mtopGuard(mtop);
 202     mtop->mols.nr  = 1;
 203     mtop->nmoltype = 1;
 204     snew(mtop->moltype, mtop->nmoltype);
 205     scoped_cptr<gmx_moltype_t> moltypeGuard(mtop->moltype);
 206     mtop->nmolblock = 1;
 207     snew(mtop->molblock, mtop->nmolblock);
 208     scoped_cptr<gmx_molblock_t> molblockGuard(mtop->molblock);
 209     mtop->molblock[0].type = 0;
 210     std::vector<int>            iatoms;
 211     for (int i = 0; i < numSettles; ++i)
 212     {
 213         iatoms.push_back(settleType);
 214         iatoms.push_back(i*atomsPerSettle+0);
 215         iatoms.push_back(i*atomsPerSettle+1);
 216         iatoms.push_back(i*atomsPerSettle+2);
 217     }
 218     mtop->moltype[0].ilist[F_SETTLE].iatoms = iatoms.data();
 219     mtop->moltype[0].ilist[F_SETTLE].nr     = iatoms.size();
 220
 221     // Set up the SETTLE parameters.
 222     mtop->ffparams.ntypes = 1;
 223     snew(mtop->ffparams.iparams, mtop->ffparams.ntypes);
 224     scoped_cptr<t_iparams> iparamsGuard(mtop->ffparams.iparams);
 225     const real             dOH = 0.09572;
 226     const real             dHH = 0.15139;
 227     mtop->ffparams.iparams[settleType].settle.doh = dOH;
 228     mtop->ffparams.iparams[settleType].settle.dhh = dHH;
 229
 230     // Set up the masses.
 231     t_mdatoms         mdatoms;
 232     std::vector<real> mass, massReciprocal;
 233     const real        oxygenMass = 15.9994, hydrogenMass = 1.008;
 234     for (int i = 0; i < numSettles; ++i)
 235     {
 236         mass.push_back(oxygenMass);
 237         mass.push_back(hydrogenMass);
 238         mass.push_back(hydrogenMass);
 239         massReciprocal.push_back(1./oxygenMass);
 240         massReciprocal.push_back(1./hydrogenMass);
 241         massReciprocal.push_back(1./hydrogenMass);
 242     }
 243     mdatoms.massT   = mass.data();
 244     mdatoms.invmass = massReciprocal.data();
 245     mdatoms.homenr  = numSettles * atomsPerSettle;
 246
 247     // Finally make the settle data structures
 248     gmx_settledata_t settled = settle_init(mtop);
 249     settle_set_constraints(settled, &mtop->moltype[0].ilist[F_SETTLE], &mdatoms);
 250
 251     // Run the test
 252     bool       errorOccured;
 253     tensor     virial             = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
 254     int        numThreads         = 1, threadIndex = 0;
 255     const real reciprocalTimeStep = 1.0/0.002;
 256     csettle(settled, numThreads, threadIndex,
 257             usePbc ? &pbcXYZ_ : &pbcNone_,
 258             std::begin(g_positions), updatedPositions_.data(), reciprocalTimeStep,
 259             useVelocities ? velocities_.data() : nullptr,
 260             calcVirial, virial, &errorOccured);
 261     EXPECT_FALSE(errorOccured) << testDescription;
 262
 263     // The necessary tolerances for the test to pass were determined
 264     // empirically. This isn't nice, but the required behaviour that
 265     // SETTLE produces constrained coordinates consistent with
 266     // sensible sampling needs to be tested at a much higher level.
 267     FloatingPointTolerance tolerance =
 268         relativeToleranceAsPrecisionDependentUlp(dOH*dOH, 80, 380);
 269
 270     // Verify the updated coordinates match the requirements
 271     int positionIndex = 0, velocityIndex = 0;
 272     for (int i = 0; i < numSettles; ++i)
 273     {
 274         rvec positionO  {
 275             updatedPositions_[positionIndex++], updatedPositions_[positionIndex++], updatedPositions_[positionIndex++]
 276         };
 277         rvec positionH1 {
 278             updatedPositions_[positionIndex++], updatedPositions_[positionIndex++], updatedPositions_[positionIndex++]
 279         };
 280         rvec positionH2 {
 281             updatedPositions_[positionIndex++], updatedPositions_[positionIndex++], updatedPositions_[positionIndex++]
 282         };
 283
 284         EXPECT_REAL_EQ_TOL(dOH*dOH, distance2(positionO, positionH1), tolerance) << formatString("for water %d ", i) << testDescription;
 285         EXPECT_REAL_EQ_TOL(dOH*dOH, distance2(positionO, positionH2), tolerance) << formatString("for water %d ", i) << testDescription;
 286         EXPECT_REAL_EQ_TOL(dHH*dHH, distance2(positionH1, positionH2), tolerance) << formatString("for water %d ", i) << testDescription;
 287
 288         // This merely tests whether the velocities were
 289         // updated from the starting values of zero (or not),
 290         // but not whether the update was correct.
 291         for (int j = 0; j < atomsPerSettle * DIM; ++j, ++velocityIndex)
 292         {
 293             EXPECT_TRUE(useVelocities == (0. != velocities_[velocityIndex])) << formatString("for water %d velocity coordinate %d ", i, j) << testDescription;
 294         }
 295     }
 296
 297     // This merely tests whether the viral was updated from
 298     // the starting values of zero (or not), but not whether
 299     // any update was correct.
 300     for (int d = 0; d < DIM; ++d)
 301     {
 302         for (int dd = 0; dd < DIM; ++dd)
 303         {
 304             EXPECT_TRUE(calcVirial == (0. != virial[d][dd])) << formatString("for virial component[%d][%d] ", d, dd) << testDescription;
 305         }
 306     }
 307 }
 308
 309 using ::testing::Bool;
 310 // Scan the full Cartesian product of numbers of SETTLE interactions
 311 // (4 and 17 are chosen to test cases that do and do not match
 312 // hardware SIMD widths), and whether or not we use PBC, velocities or
 313 // calculate the virial contribution. It would be nicer to generate
 314 // these combinations with ::testing::Combine, but gcc 4.6 can't cope
 315 // with the template meta-programming required to generate the tuples.
 316 INSTANTIATE_TEST_CASE_P(WithParameters, SettleTest,
 317                             ::testing::Values(SettleTestParameters(1,  true,  true,  true),
 318                                               SettleTestParameters(4,  true,  true,  true),
 319                                               SettleTestParameters(17, true,  true,  true),
 320                                               SettleTestParameters(1,  false, true,  true),
 321                                               SettleTestParameters(4,  false, true,  true),
 322                                               SettleTestParameters(17, false, true,  true),
 323                                               SettleTestParameters(1,  true,  false, true),
 324                                               SettleTestParameters(4,  true,  false, true),
 325                                               SettleTestParameters(17, true,  false, true),
 326                                               SettleTestParameters(1,  false, false, true),
 327                                               SettleTestParameters(4,  false, false, true),
 328                                               SettleTestParameters(17, false, false, true),
 329                                               SettleTestParameters(1,  true,  true,  false),
 330                                               SettleTestParameters(4,  true,  true,  false),
 331                                               SettleTestParameters(17, true,  true,  false),
 332                                               SettleTestParameters(1,  false, true,  false),
 333                                               SettleTestParameters(4,  false, true,  false),
 334                                               SettleTestParameters(17, false, true,  false),
 335                                               SettleTestParameters(1,  true,  false, false),
 336                                               SettleTestParameters(4,  true,  false, false),
 337                                               SettleTestParameters(17, true,  false, false),
 338                                               SettleTestParameters(1,  false, false, false),
 339                                               SettleTestParameters(4,  false, false, false),
 340                                               SettleTestParameters(17, false, false, false)));
 341
 342
 343 } // namespace
 344 } // namespace