src/gromacs/mdlib/tests/updategroups.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2018,2019, 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/testasserts.h"
  51
  52 namespace gmx
  53 {
  54
  55 namespace
  56 {
  57
  58 /* TODO: Actually initialize moltype.atoms.atom when this is converted to C++ */
  59
  60 /*! \brief Returns an ethane united-atom molecule */
  61 gmx_moltype_t ethaneUA()
  62 {
  63     gmx_moltype_t moltype = {};
  64
  65     moltype.atoms.nr               = 2;
  66     moltype.ilist[F_CONSTR].iatoms = { 0, 0, 1 };
  67
  68     return moltype;
  69 }
  70
  71 /*! \brief Returns a methane molecule */
  72 gmx_moltype_t methane()
  73 {
  74     gmx_moltype_t moltype = {};
  75
  76     moltype.atoms.nr               = 5;
  77     moltype.ilist[F_CONSTR].iatoms = { 0, 0, 1, 0, 0, 2, 0, 0, 3, 0, 0, 4 };
  78
  79     return moltype;
  80 }
  81
  82 /*! \brief Returns an ethane molecule */
  83 gmx_moltype_t ethane()
  84 {
  85     gmx_moltype_t moltype = {};
  86
  87     moltype.atoms.nr               = 8;
  88     moltype.ilist[F_CONSTR].iatoms = { 0, 0, 1, 0, 0, 2, 0, 0, 3, 0, 4, 5, 0, 4, 6, 0, 4, 7 };
  89     moltype.ilist[F_ANGLES].iatoms = { 1, 1, 0, 2, 1, 1, 0, 3, 1, 2, 0, 3,
  90                                        1, 5, 4, 6, 1, 5, 4, 7, 1, 6, 4, 7 };
  91
  92     return moltype;
  93 }
  94
  95 /*! \brief Returns a butane united-atom molecule */
  96 gmx_moltype_t butaneUA()
  97 {
  98     gmx_moltype_t moltype = {};
  99
 100     moltype.atoms.nr               = 4;
 101     moltype.ilist[F_CONSTR].iatoms = { 0, 0, 1, 0, 1, 2, 0, 2, 3 };
 102
 103     return moltype;
 104 }
 105
 106 /*! \brief Returns a three-site water molecule */
 107 gmx_moltype_t waterThreeSite()
 108 {
 109     gmx_moltype_t moltype = {};
 110
 111     moltype.atoms.nr               = 3;
 112     moltype.ilist[F_SETTLE].iatoms = { 0, 0, 1, 2 };
 113
 114     return moltype;
 115 }
 116
 117 /*! \brief Returns a four-site water molecule with virtual site */
 118 gmx_moltype_t waterFourSite()
 119 {
 120     gmx_moltype_t moltype = {};
 121
 122     moltype.atoms.nr               = 4;
 123     moltype.ilist[F_SETTLE].iatoms = { 0, 1, 2, 3 };
 124     moltype.ilist[F_VSITE3].iatoms = { 1, 0, 1, 2, 3 };
 125
 126     return moltype;
 127 }
 128
 129 /*! \brief Returns a water molecule with flexible angle */
 130 gmx_moltype_t waterFlexAngle()
 131 {
 132     gmx_moltype_t moltype = {};
 133
 134     moltype.atoms.nr               = 3;
 135     moltype.ilist[F_CONSTR].iatoms = {
 136         0, 0, 1, 0, 0, 2,
 137     };
 138     moltype.ilist[F_ANGLES].iatoms = {
 139         1,
 140         1,
 141         0,
 142         2,
 143     };
 144
 145     return moltype;
 146 }
 147
 148 TEST(UpdateGroups, ethaneUA)
 149 {
 150     gmx_mtop_t mtop;
 151
 152     mtop.moltype.emplace_back(ethaneUA());
 153     t_iparams iparams;
 154     iparams.constr = { 0.3, 0.3 };
 155     mtop.ffparams.iparams.push_back(iparams);
 156
 157     auto updateGroups = gmx::makeUpdateGroups(mtop);
 158
 159     ASSERT_EQ(updateGroups.size(), 1);
 160     EXPECT_EQ(updateGroups[0].numBlocks(), 1);
 161
 162     real temperature = 298;
 163     real maxRadius   = computeMaxUpdateGroupRadius(mtop, updateGroups, temperature);
 164     EXPECT_FLOAT_EQ(maxRadius, 0.3 / 2);
 165 }
 166
 167 TEST(UpdateGroups, methane)
 168 {
 169     gmx_mtop_t mtop;
 170
 171     mtop.moltype.emplace_back(methane());
 172     t_iparams iparams;
 173     iparams.constr = { 0.1, 0.1 };
 174     mtop.ffparams.iparams.push_back(iparams);
 175
 176     auto updateGroups = gmx::makeUpdateGroups(mtop);
 177
 178     ASSERT_EQ(updateGroups.size(), 1);
 179     EXPECT_EQ(updateGroups[0].numBlocks(), 1);
 180
 181     real temperature = 298;
 182     real maxRadius   = computeMaxUpdateGroupRadius(mtop, updateGroups, temperature);
 183     EXPECT_FLOAT_EQ(maxRadius, 0.14);
 184 }
 185 TEST(UpdateGroups, ethane)
 186 {
 187     gmx_mtop_t mtop;
 188
 189     mtop.moltype.emplace_back(ethane());
 190     t_iparams iparams;
 191     iparams.constr = { 0.1, 0.1 };
 192     mtop.ffparams.iparams.push_back(iparams);
 193     iparams.harmonic = { 107.800, 276.144, 107.800, 276.144 };
 194     mtop.ffparams.iparams.push_back(iparams);
 195
 196     auto updateGroups = gmx::makeUpdateGroups(mtop);
 197
 198     ASSERT_EQ(updateGroups.size(), 1);
 199     EXPECT_EQ(updateGroups[0].numBlocks(), 2);
 200
 201     real temperature = 298;
 202     real maxRadius   = computeMaxUpdateGroupRadius(mtop, updateGroups, temperature);
 203     EXPECT_FLOAT_EQ(maxRadius, 0.094746813);
 204
 205     temperature = 0;
 206     maxRadius   = computeMaxUpdateGroupRadius(mtop, updateGroups, temperature);
 207     EXPECT_FLOAT_EQ(maxRadius, 0.10310466);
 208
 209     temperature = -1;
 210     maxRadius   = computeMaxUpdateGroupRadius(mtop, updateGroups, temperature);
 211     EXPECT_FLOAT_EQ(maxRadius, 0.125);
 212 }
 213
 214 TEST(UpdateGroups, butaneUA)
 215 {
 216     gmx_mtop_t mtop;
 217
 218     mtop.moltype.emplace_back(butaneUA());
 219     t_iparams iparams;
 220     iparams.constr = { 0.3, 0.3 };
 221     mtop.ffparams.iparams.push_back(iparams);
 222
 223     auto updateGroups = gmx::makeUpdateGroups(mtop);
 224
 225     EXPECT_EQ(updateGroups.size(), 0);
 226 }
 227
 228 TEST(UpdateGroups, waterThreeSite)
 229 {
 230     gmx_mtop_t mtop;
 231
 232     mtop.moltype.emplace_back(waterThreeSite());
 233     t_iparams iparams;
 234     iparams.settle = { 0.1, 0.1633 };
 235     mtop.ffparams.iparams.push_back(iparams);
 236
 237     auto updateGroups = gmx::makeUpdateGroups(mtop);
 238
 239     ASSERT_EQ(updateGroups.size(), 1);
 240     EXPECT_EQ(updateGroups[0].numBlocks(), 1);
 241
 242     real temperature = 298;
 243     real maxRadius   = computeMaxUpdateGroupRadius(mtop, updateGroups, temperature);
 244     EXPECT_FLOAT_EQ(maxRadius, 0.083887339);
 245 }
 246
 247 // Tests update group with virtual site
 248 TEST(UpdateGroups, waterFourSite)
 249 {
 250     gmx_mtop_t mtop;
 251
 252     mtop.moltype.emplace_back(waterFourSite());
 253     t_iparams iparams[2];
 254     iparams[0].settle = { 0.1, 0.1633 };
 255     iparams[1].vsite  = { 0.128, 0.128 };
 256     mtop.ffparams.iparams.push_back(iparams[0]);
 257     mtop.ffparams.iparams.push_back(iparams[1]);
 258
 259     auto updateGroups = gmx::makeUpdateGroups(mtop);
 260
 261     ASSERT_EQ(updateGroups.size(), 1);
 262     EXPECT_EQ(updateGroups[0].numBlocks(), 1);
 263 }
 264
 265 TEST(UpdateGroups, fourAtomsWithSettle)
 266 {
 267     gmx_mtop_t mtop;
 268
 269     mtop.moltype.emplace_back(waterThreeSite());
 270     mtop.moltype.back().atoms.nr = 4;
 271
 272     auto updateGroups = gmx::makeUpdateGroups(mtop);
 273
 274     ASSERT_EQ(updateGroups.size(), 1);
 275     EXPECT_EQ(updateGroups[0].numBlocks(), 2);
 276 }
 277
 278 // Tests groups with two constraints and an angle potential
 279 TEST(UpdateGroups, waterFlexAngle)
 280 {
 281     gmx_mtop_t mtop;
 282
 283     mtop.moltype.emplace_back(waterFlexAngle());
 284     t_iparams iparams;
 285     iparams.constr = { 0.1, 0.1 };
 286     mtop.ffparams.iparams.push_back(iparams);
 287     iparams.harmonic = { 109.47, 383.0, 109.47, 383.0 };
 288     mtop.ffparams.iparams.push_back(iparams);
 289
 290     auto updateGroups = gmx::makeUpdateGroups(mtop);
 291
 292     ASSERT_EQ(updateGroups.size(), 1);
 293     EXPECT_EQ(updateGroups[0].numBlocks(), 1);
 294
 295     real temperature = 298;
 296     real maxRadius   = computeMaxUpdateGroupRadius(mtop, updateGroups, temperature);
 297     EXPECT_FLOAT_EQ(maxRadius, 0.090824135);
 298
 299     temperature = 0;
 300     maxRadius   = computeMaxUpdateGroupRadius(mtop, updateGroups, temperature);
 301     EXPECT_FLOAT_EQ(maxRadius, 0.1);
 302
 303     temperature = -1;
 304     maxRadius   = computeMaxUpdateGroupRadius(mtop, updateGroups, temperature);
 305     EXPECT_FLOAT_EQ(maxRadius, 0.1);
 306 }
 307
 308 TEST(UpdateGroups, twoMoltypes)
 309 {
 310     gmx_mtop_t mtop;
 311
 312     mtop.moltype.emplace_back(methane());
 313     t_iparams iparams;
 314     iparams.constr = { 0.1, 0.1 };
 315     mtop.ffparams.iparams.push_back(iparams);
 316
 317     mtop.moltype.emplace_back(waterThreeSite());
 318     // Note: iparams not accessed for SETTLE when not computing radius
 319
 320     auto updateGroups = gmx::makeUpdateGroups(mtop);
 321
 322     ASSERT_EQ(updateGroups.size(), 2);
 323     EXPECT_EQ(updateGroups[0].numBlocks(), 1);
 324     EXPECT_EQ(updateGroups[1].numBlocks(), 1);
 325 }
 326
 327 } // namespace
 328
 329 } // namespace gmx