src/gromacs/mdlib/tests/updategroups.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2018,2019,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 for the update groups functionality.
  38  *
  39  * \author berk Hess <hess@kth.se>
  40  * \ingroup module_mdlib
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "gromacs/mdlib/updategroups.h"
  45
  46 #include <gtest/gtest.h>
  47
  48 #include "gromacs/topology/topology.h"
  49
  50 #include "testutils/loggertest.h"
  51 #include "testutils/testasserts.h"
  52
  53 namespace gmx
  54 {
  55
  56 namespace
  57 {
  58
  59 /* TODO: Actually initialize moltype.atoms.atom when this is converted to C++ */
  60
  61 /*! \brief Returns a flexible ethane united-atom molecule */
  62 gmx_moltype_t flexibleEthaneUA()
  63 {
  64     gmx_moltype_t moltype = {};
  65
  66     moltype.atoms.nr              = 2;
  67     moltype.ilist[F_BONDS].iatoms = { 0, 0, 1 };
  68
  69     return moltype;
  70 }
  71
  72 /*! \brief Returns an ethane united-atom molecule */
  73 gmx_moltype_t ethaneUA()
  74 {
  75     gmx_moltype_t moltype = {};
  76
  77     moltype.atoms.nr               = 2;
  78     moltype.ilist[F_CONSTR].iatoms = { 0, 0, 1 };
  79
  80     return moltype;
  81 }
  82
  83 /*! \brief Returns a methane molecule */
  84 gmx_moltype_t methane()
  85 {
  86     gmx_moltype_t moltype = {};
  87
  88     moltype.atoms.nr               = 5;
  89     moltype.ilist[F_CONSTR].iatoms = { 0, 0, 1, 0, 0, 2, 0, 0, 3, 0, 0, 4 };
  90
  91     return moltype;
  92 }
  93
  94 /*! \brief Returns an ethane molecule */
  95 gmx_moltype_t ethane()
  96 {
  97     gmx_moltype_t moltype = {};
  98
  99     moltype.atoms.nr               = 8;
 100     moltype.ilist[F_CONSTR].iatoms = { 0, 0, 1, 0, 0, 2, 0, 0, 3, 0, 4, 5, 0, 4, 6, 0, 4, 7 };
 101     moltype.ilist[F_ANGLES].iatoms = { 1, 1, 0, 2, 1, 1, 0, 3, 1, 2, 0, 3,
 102                                        1, 5, 4, 6, 1, 5, 4, 7, 1, 6, 4, 7 };
 103
 104     return moltype;
 105 }
 106
 107 /*! \brief Returns a butane fully-constrained united-atom molecule */
 108 gmx_moltype_t butaneUA()
 109 {
 110     gmx_moltype_t moltype = {};
 111
 112     moltype.atoms.nr               = 4;
 113     moltype.ilist[F_CONSTR].iatoms = { 0, 0, 1, 0, 1, 2, 0, 2, 3 };
 114
 115     return moltype;
 116 }
 117
 118 /*! \brief Returns a three-site water molecule */
 119 gmx_moltype_t waterThreeSite()
 120 {
 121     gmx_moltype_t moltype = {};
 122
 123     moltype.atoms.nr               = 3;
 124     moltype.ilist[F_SETTLE].iatoms = { 0, 0, 1, 2 };
 125
 126     return moltype;
 127 }
 128
 129 /*! \brief Returns a four-site water molecule with virtual site */
 130 gmx_moltype_t waterFourSite()
 131 {
 132     gmx_moltype_t moltype = {};
 133
 134     moltype.atoms.nr               = 4;
 135     moltype.ilist[F_SETTLE].iatoms = { 0, 1, 2, 3 };
 136     moltype.ilist[F_VSITE3].iatoms = { 1, 0, 1, 2, 3 };
 137
 138     return moltype;
 139 }
 140
 141 /*! \brief Returns a water molecule with flexible angle */
 142 gmx_moltype_t waterFlexAngle()
 143 {
 144     gmx_moltype_t moltype = {};
 145
 146     moltype.atoms.nr               = 3;
 147     moltype.ilist[F_CONSTR].iatoms = {
 148         0, 0, 1, 0, 0, 2,
 149     };
 150     moltype.ilist[F_ANGLES].iatoms = {
 151         1,
 152         1,
 153         0,
 154         2,
 155     };
 156
 157     return moltype;
 158 }
 159
 160 //! Test fixture class
 161 class UpdateGroupsTest : public ::testing::Test
 162 {
 163 public:
 164     //! Global toplogy to use in tests
 165     gmx_mtop_t mtop_;
 166     //! Default temperature for tests
 167     real temperature_ = 298;
 168     //! Logger to use in tests
 169     test::LoggerTestHelper logHelper_;
 170 };
 171
 172 TEST_F(UpdateGroupsTest, WithEthaneUA)
 173 {
 174     mtop_.moltype.emplace_back(ethaneUA());
 175     {
 176         t_iparams iparams;
 177         iparams.constr = { 0.3, 0.3 };
 178         mtop_.ffparams.iparams.push_back(iparams);
 179     }
 180
 181     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 182
 183     ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1);
 184     EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 1);
 185
 186     real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 187     EXPECT_FLOAT_EQ(maxRadius, 0.3 / 2);
 188
 189     logHelper_.expectNoEntries(MDLogger::LogLevel::Info);
 190     UpdateGroups updateGroups = makeUpdateGroups(
 191             logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real);
 192 }
 193
 194 TEST_F(UpdateGroupsTest, WithMethane)
 195 {
 196     mtop_.moltype.emplace_back(methane());
 197     {
 198         t_iparams iparams;
 199         iparams.constr = { 0.1, 0.1 };
 200         mtop_.ffparams.iparams.push_back(iparams);
 201     }
 202
 203     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 204
 205     ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1);
 206     EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 1);
 207
 208     real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 209     EXPECT_FLOAT_EQ(maxRadius, 0.14);
 210
 211     logHelper_.expectNoEntries(MDLogger::LogLevel::Info);
 212     UpdateGroups updateGroups = makeUpdateGroups(
 213             logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real);
 214 }
 215
 216 TEST_F(UpdateGroupsTest, WithEthane)
 217 {
 218     mtop_.moltype.emplace_back(ethane());
 219     {
 220         t_iparams iparams;
 221         iparams.constr = { 0.1, 0.1 };
 222         mtop_.ffparams.iparams.push_back(iparams);
 223         iparams.harmonic = { 107.800, 276.144, 107.800, 276.144 };
 224         mtop_.ffparams.iparams.push_back(iparams);
 225     }
 226
 227     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 228
 229     ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1);
 230     EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 2);
 231
 232     real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 233     EXPECT_FLOAT_EQ(maxRadius, 0.094746813);
 234
 235     logHelper_.expectNoEntries(MDLogger::LogLevel::Info);
 236     UpdateGroups updateGroups = makeUpdateGroups(
 237             logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real);
 238 }
 239
 240 TEST_F(UpdateGroupsTest, CheckRadiusCalculationAtDifferentTemperaturesWithEthane)
 241 {
 242     mtop_.moltype.emplace_back(ethane());
 243     {
 244         t_iparams iparams;
 245         iparams.constr = { 0.1, 0.1 };
 246         mtop_.ffparams.iparams.push_back(iparams);
 247         iparams.harmonic = { 107.800, 276.144, 107.800, 276.144 };
 248         mtop_.ffparams.iparams.push_back(iparams);
 249     }
 250
 251     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 252
 253     ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1);
 254     EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 2);
 255
 256     real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 257     EXPECT_FLOAT_EQ(maxRadius, 0.094746813);
 258
 259     // Observe that the temperature affects the radius only when valid
 260     temperature_ = 0;
 261     maxRadius    = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 262     EXPECT_FLOAT_EQ(maxRadius, 0.10310466);
 263
 264     temperature_ = -1;
 265     maxRadius    = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 266     EXPECT_FLOAT_EQ(maxRadius, 0.125);
 267 }
 268
 269 TEST_F(UpdateGroupsTest, WithButaneUALogsThatUnsuitableForUpdateGroups)
 270 {
 271     mtop_.moltype.emplace_back(butaneUA());
 272     {
 273         t_iparams iparams;
 274         iparams.constr = { 0.3, 0.3 };
 275         mtop_.ffparams.iparams.push_back(iparams);
 276     }
 277
 278     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 279
 280     EXPECT_EQ(updateGroupingsPerMoleculeType.size(), 0);
 281
 282     real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 283     EXPECT_FLOAT_EQ(maxRadius, 0.0);
 284
 285     logHelper_.expectEntryMatchingRegex(
 286             MDLogger::LogLevel::Info,
 287             "At least one moleculetype does not conform to the requirements");
 288     UpdateGroups updateGroups = makeUpdateGroups(
 289             logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real);
 290 }
 291
 292 TEST_F(UpdateGroupsTest, WithWaterThreeSite)
 293 {
 294     mtop_.moltype.emplace_back(waterThreeSite());
 295     {
 296         t_iparams iparams;
 297         iparams.settle = { 0.1, 0.1633 };
 298         mtop_.ffparams.iparams.push_back(iparams);
 299     }
 300
 301     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 302
 303     ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1);
 304     EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 1);
 305
 306     real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 307     EXPECT_FLOAT_EQ(maxRadius, 0.083887339);
 308
 309     logHelper_.expectNoEntries(MDLogger::LogLevel::Info);
 310     UpdateGroups updateGroups = makeUpdateGroups(
 311             logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real);
 312 }
 313
 314 // Tests update group with virtual site
 315 TEST_F(UpdateGroupsTest, WithWaterFourSite)
 316 {
 317     mtop_.moltype.emplace_back(waterFourSite());
 318     {
 319         t_iparams iparams[2];
 320         iparams[0].settle = { 0.1, 0.1633 };
 321         iparams[1].vsite  = { 0.128, 0.128 };
 322         mtop_.ffparams.iparams.push_back(iparams[0]);
 323         mtop_.ffparams.iparams.push_back(iparams[1]);
 324     }
 325
 326     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 327
 328     ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1);
 329     EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 1);
 330
 331     real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 332     EXPECT_FLOAT_EQ(maxRadius, 0.083887339);
 333
 334     logHelper_.expectNoEntries(MDLogger::LogLevel::Info);
 335     UpdateGroups updateGroups = makeUpdateGroups(
 336             logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real);
 337 }
 338
 339 TEST_F(UpdateGroupsTest, WithFourAtomsWithSettle)
 340 {
 341     mtop_.moltype.emplace_back(waterThreeSite());
 342     mtop_.moltype.back().atoms.nr = 4;
 343     {
 344         t_iparams iparams[2];
 345         iparams[0].settle = { 0.1, 0.1633 };
 346         iparams[1].vsite  = { 0.128, 0.128 };
 347         mtop_.ffparams.iparams.push_back(iparams[0]);
 348         mtop_.ffparams.iparams.push_back(iparams[1]);
 349     }
 350
 351     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 352
 353     ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1);
 354     EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 2);
 355
 356     real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 357     EXPECT_FLOAT_EQ(maxRadius, 0.083887339);
 358
 359     logHelper_.expectNoEntries(MDLogger::LogLevel::Info);
 360     UpdateGroups updateGroups = makeUpdateGroups(
 361             logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real);
 362 }
 363
 364 // Tests groups with two constraints and an angle potential
 365 TEST_F(UpdateGroupsTest, WithWaterFlexAngle)
 366 {
 367     mtop_.moltype.emplace_back(waterFlexAngle());
 368     {
 369         t_iparams iparams;
 370         iparams.constr = { 0.1, 0.1 };
 371         mtop_.ffparams.iparams.push_back(iparams);
 372         iparams.harmonic = { 109.47, 383.0, 109.47, 383.0 };
 373         mtop_.ffparams.iparams.push_back(iparams);
 374     }
 375
 376     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 377
 378     ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1);
 379     EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 1);
 380
 381     real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 382     EXPECT_FLOAT_EQ(maxRadius, 0.090824135);
 383
 384     logHelper_.expectNoEntries(MDLogger::LogLevel::Info);
 385     UpdateGroups updateGroups = makeUpdateGroups(
 386             logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real);
 387 }
 388
 389 TEST_F(UpdateGroupsTest, CheckRadiusCalculationAtDifferentTemperaturesWithWaterFlexAngle)
 390 {
 391     mtop_.moltype.emplace_back(waterFlexAngle());
 392     {
 393         t_iparams iparams;
 394         iparams.constr = { 0.1, 0.1 };
 395         mtop_.ffparams.iparams.push_back(iparams);
 396         iparams.harmonic = { 109.47, 383.0, 109.47, 383.0 };
 397         mtop_.ffparams.iparams.push_back(iparams);
 398     }
 399
 400     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 401
 402     ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1);
 403     EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 1);
 404
 405     real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 406     EXPECT_FLOAT_EQ(maxRadius, 0.090824135);
 407
 408     // Observe that the temperature affects the radius only when valid
 409     temperature_ = 0;
 410     maxRadius    = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 411     EXPECT_FLOAT_EQ(maxRadius, 0.1);
 412
 413     temperature_ = -1;
 414     maxRadius    = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 415     EXPECT_FLOAT_EQ(maxRadius, 0.1);
 416 }
 417
 418 TEST_F(UpdateGroupsTest, WithTwoMoltypes)
 419 {
 420     mtop_.moltype.emplace_back(methane());
 421     {
 422         t_iparams iparams;
 423         iparams.constr = { 0.1, 0.1 };
 424         mtop_.ffparams.iparams.push_back(iparams);
 425     }
 426
 427     mtop_.moltype.emplace_back(waterThreeSite());
 428     // Note: iparams not accessed for SETTLE when not computing radius
 429
 430     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 431
 432     ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 2);
 433     EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 1);
 434     EXPECT_EQ(updateGroupingsPerMoleculeType[1].numBlocks(), 1);
 435
 436     real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 437     EXPECT_FLOAT_EQ(maxRadius, 0.14);
 438
 439     logHelper_.expectNoEntries(MDLogger::LogLevel::Info);
 440     UpdateGroups updateGroups = makeUpdateGroups(
 441             logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real);
 442 }
 443
 444 TEST_F(UpdateGroupsTest, LogsWhenSizesAreInvalid)
 445 {
 446     mtop_.moltype.emplace_back(methane());
 447     {
 448         t_iparams iparams;
 449         iparams.constr = { 0.1, 0.1 };
 450         mtop_.ffparams.iparams.push_back(iparams);
 451     }
 452
 453     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 454     logHelper_.expectEntryMatchingRegex(MDLogger::LogLevel::Info,
 455                                         "The combination of rlist and box size prohibits");
 456     UpdateGroups updateGroups = makeUpdateGroups(
 457             logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), 1e9_real, true, true, 1e6_real);
 458 }
 459
 460 TEST_F(UpdateGroupsTest, LogsWhenUpdateGroupsAreNotUseful)
 461 {
 462     mtop_.moltype.emplace_back(flexibleEthaneUA());
 463     {
 464         t_iparams iparams;
 465         iparams.harmonic = { 0.1, 10.0, 0.1, 10.0 };
 466         mtop_.ffparams.iparams.push_back(iparams);
 467     }
 468
 469     auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_);
 470
 471     ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1);
 472     EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 2);
 473
 474     real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_);
 475     EXPECT_FLOAT_EQ(maxRadius, 0);
 476
 477     logHelper_.expectEntryMatchingRegex(MDLogger::LogLevel::Info,
 478                                         "No constraints or virtual sites are in use");
 479     UpdateGroups updateGroups = makeUpdateGroups(
 480             logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, false, 1e6_real);
 481 }
 482
 483
 484 } // namespace
 485
 486 } // namespace gmx