src/gromacs/pbcutil/tests/com.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  * Tests COM handling code.
  38  *
  39  * \author Paul Bauer <paul.bauer.q@gmail.com>
  40  * \ingroup module_pbcutil
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "gromacs/pbcutil/com.h"
  45
  46 #include <vector>
  47
  48 #include <gtest/gtest.h>
  49
  50 #include "gromacs/math/vectypes.h"
  51 #include "gromacs/pbcutil/pbc.h"
  52 #include "gromacs/topology/mtop_util.h"
  53 #include "gromacs/topology/topology.h"
  54
  55 #include "testutils/refdata.h"
  56 #include "testutils/testasserts.h"
  57
  58 namespace gmx
  59 {
  60
  61 namespace test
  62 {
  63
  64 namespace
  65 {
  66
  67 /*! \brief Populates a molectype for generate a graph
  68  *
  69  * This moleculetype has the first and last atom not connected to the other atoms
  70  * so we test optimizations in the shift code that generates the actual graph only
  71  * for the atom range from the first to the last connected atoms.
  72  *
  73  * Defines three residues to test residue and molecule COM handling.
  74  *
  75  * \todo Return moltype once t_atoms is a proper class.
  76  */
  77 void populateMoleculeType(gmx_moltype_t* moltype)
  78 {
  79     constexpr int atomNumber        = 5;
  80     constexpr int residueNumber     = 3;
  81     moltype->atoms.nr               = atomNumber;
  82     moltype->ilist[F_CONSTR].iatoms = { 0, 1, 2 };
  83     moltype->ilist[F_ANGLES].iatoms = { 1, 2, 1, 3 };
  84
  85     snew(moltype->atoms.atom, atomNumber);
  86     snew(moltype->atoms.resinfo, residueNumber);
  87     moltype->atoms.atom[0].resind = 0;
  88     moltype->atoms.atom[1].resind = 1;
  89     moltype->atoms.atom[2].resind = 1;
  90     moltype->atoms.atom[3].resind = 1;
  91     moltype->atoms.atom[4].resind = 2;
  92
  93     moltype->atoms.atom[0].m = 42;
  94     moltype->atoms.atom[1].m = 13;
  95     moltype->atoms.atom[2].m = 7;
  96     moltype->atoms.atom[3].m = 23;
  97     moltype->atoms.atom[4].m = 2;
  98 }
  99
 100 //! Set up initial coordinates.
 101 std::vector<RVec> initialCoordinates()
 102 {
 103     std::vector<RVec> coordinates;
 104     coordinates.emplace_back(-1, 0, 3);
 105     coordinates.emplace_back(1.5, 1.5, 4.5);
 106     coordinates.emplace_back(1.6, 1.5, 4.5);
 107     coordinates.emplace_back(1.6, 1.7, 4.5);
 108     coordinates.emplace_back(1, 4, 2);
 109
 110     coordinates.emplace_back(1, 0, -3);
 111     coordinates.emplace_back(-1.5, -1.5, -4.5);
 112     coordinates.emplace_back(-1.6, -1.5, -4.5);
 113     coordinates.emplace_back(-1.6, -1.7, -4.5);
 114     coordinates.emplace_back(-1, -4, -2);
 115
 116     return coordinates;
 117 }
 118
 119 using COMInPlaceTestParams = std::tuple<UnitCellType, CenteringType, PbcType>;
 120 using test::FloatingPointTolerance;
 121
 122 /*! \brief
 123  * Test fixture for checking correct molecule COM treatment.
 124  */
 125 class COMInPlaceTest : public ::testing::Test, public ::testing::WithParamInterface<COMInPlaceTestParams>
 126 {
 127 public:
 128     COMInPlaceTest();
 129
 130     //! Run the test with the given input for molecule COM.
 131     void runTestMolecule(const matrix box);
 132     //! Run the test with the given input for residue COM.
 133     void runTestResidue(const matrix box);
 134
 135 private:
 136     //! Coordinates to use.
 137     std::vector<RVec> testCoordinates_;
 138     //! Dummy topology to use.
 139     gmx_mtop_t testTopology_;
 140     //! Reference data holder.
 141     TestReferenceData data_;
 142     //! Checker for data.
 143     TestReferenceChecker checker_;
 144 };
 145
 146 COMInPlaceTest::COMInPlaceTest() :
 147     testCoordinates_(initialCoordinates()), checker_(data_.rootChecker())
 148 {
 149     auto& moltype = testTopology_.moltype.emplace_back();
 150     populateMoleculeType(&moltype);
 151     auto& molblock       = testTopology_.molblock.emplace_back();
 152     molblock.nmol        = 2;
 153     molblock.type        = 0;
 154     testTopology_.natoms = moltype.atoms.nr * molblock.nmol;
 155     testTopology_.finalize();
 156     FloatingPointTolerance tolerance(
 157             FloatingPointTolerance(1.0e-6, 1.0e-6, 1.0e-8, 1.0e-12, 10000, 100, false));
 158     checker_.setDefaultTolerance(tolerance);
 159 }
 160
 161 void COMInPlaceTest::runTestMolecule(const matrix box)
 162 {
 163     auto params   = GetParam();
 164     auto unitcell = std::get<0>(params);
 165     auto center   = std::get<1>(params);
 166     auto pbcType  = std::get<2>(params);
 167     placeCoordinatesWithCOMInBox(
 168             pbcType, unitcell, center, box, testCoordinates_, testTopology_, COMShiftType::Molecule);
 169     std::string testString = "Molecule " + std::string(unitCellTypeNames(unitcell))
 170                              + std::string(centerTypeNames(center)) + c_pbcTypeNames[pbcType]
 171                              + c_pbcTypeNames[guessPbcType(box)];
 172     checker_.checkSequence(std::begin(testCoordinates_), std::end(testCoordinates_), testString.c_str());
 173 }
 174
 175 void COMInPlaceTest::runTestResidue(const matrix box)
 176 {
 177     auto params   = GetParam();
 178     auto unitcell = std::get<0>(params);
 179     auto center   = std::get<1>(params);
 180     auto pbcType  = std::get<2>(params);
 181     placeCoordinatesWithCOMInBox(
 182             pbcType, unitcell, center, box, testCoordinates_, testTopology_, COMShiftType::Residue);
 183     std::string testString = "Residue " + std::string(unitCellTypeNames(unitcell))
 184                              + std::string(centerTypeNames(center)) + c_pbcTypeNames[pbcType]
 185                              + c_pbcTypeNames[guessPbcType(box)];
 186     checker_.checkSequence(std::begin(testCoordinates_), std::end(testCoordinates_), testString.c_str());
 187 }
 188
 189 TEST(ShiftTest, CoordinateShiftWorks)
 190 {
 191     auto       coords = initialCoordinates();
 192     const RVec shift(1.5, -2.2, 0);
 193
 194     shiftAtoms(shift, coords);
 195     EXPECT_FLOAT_EQ(coords[4][0], 2.5);
 196     EXPECT_FLOAT_EQ(coords[3][1], -0.5);
 197     EXPECT_FLOAT_EQ(coords[9][2], -2);
 198 }
 199
 200
 201 TEST_P(COMInPlaceTest, MatrixDefault)
 202 {
 203     const matrix box = { { 3, 0, 0 }, { 0, 3, 0 }, { 0, 0, 3 } };
 204
 205     runTestMolecule(box);
 206     runTestResidue(box);
 207 }
 208
 209 INSTANTIATE_TEST_CASE_P(CorrectCoordinates,
 210                         COMInPlaceTest,
 211                         testing::Combine(::testing::Values(UnitCellType::Compact,
 212                                                            UnitCellType::Rectangular,
 213                                                            UnitCellType::Triclinic),
 214                                          ::testing::Values(CenteringType::Rectangular,
 215                                                            CenteringType::Triclinic,
 216                                                            CenteringType::Zero),
 217                                          ::testing::Values(PbcType::No, PbcType::Xyz, PbcType::XY)));
 218
 219 // TODO add PbcType::Screw once it is fully supported.
 220
 221 } // namespace
 222
 223 } // namespace test
 224
 225 } // namespace gmx