src/gromacs/listed-forces/tests/bonded.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2016,2017, 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  * Implements test of bonded force routines
  38  *
  39  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  40  * \ingroup module_listed-forces
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "gromacs/listed-forces/bonded.h"
  45
  46 #include <cmath>
  47
  48 #include <memory>
  49
  50 #include <gtest/gtest.h>
  51
  52 #include "gromacs/math/units.h"
  53 #include "gromacs/pbcutil/ishift.h"
  54 #include "gromacs/pbcutil/pbc.h"
  55
  56 #include "testutils/refdata.h"
  57 #include "testutils/testasserts.h"
  58 #include "testutils/testfilemanager.h"
  59
  60 namespace gmx
  61 {
  62 namespace
  63 {
  64
  65 //! Number of atoms used in these tests.
  66 #define NATOMS 4
  67
  68 class BondedTest : public ::testing::Test
  69 {
  70     protected:
  71         rvec   x[NATOMS];
  72         matrix box;
  73         test::TestReferenceData           refData_;
  74         test::TestReferenceChecker        checker_;
  75         BondedTest( ) :
  76             checker_(refData_.rootChecker())
  77         {
  78             test::FloatingPointTolerance tolerance(test::relativeToleranceAsFloatingPoint(1.0, 1e-6));
  79             checker_.setDefaultTolerance(tolerance);
  80             clear_rvecs(NATOMS, x);
  81             x[1][2] = 1;
  82             x[2][1] = x[2][2] = 1;
  83             x[3][0] = x[3][1] = x[3][2] = 1;
  84
  85             clear_mat(box);
  86             box[0][0] = box[1][1] = box[2][2] = 1.5;
  87         }
  88
  89         void testBondAngle(int epbc)
  90         {
  91             rvec  r_ij, r_kj;
  92             real  cosine_angle, angle;
  93             int   t1, t2;
  94             t_pbc pbc;
  95
  96             set_pbc(&pbc, epbc, box);
  97             angle = bond_angle(x[0], x[1], x[2], &pbc,
  98                                r_ij, r_kj, &cosine_angle,
  99                                &t1, &t2);
 100             checker_.checkReal(angle, "angle");
 101             checker_.checkReal(cosine_angle, "cosine_angle");
 102             checker_.checkInteger(t1, "t1");
 103             checker_.checkInteger(t2, "t2");
 104         }
 105
 106         void testDihedralAngle(int epbc)
 107         {
 108             rvec  r_ij, r_kj, r_kl, m, n;
 109             real  cosine_angle, angle;
 110             int   t1, t2, t3;
 111             t_pbc pbc;
 112
 113             set_pbc(&pbc, epbc, box);
 114             angle = dih_angle(x[0], x[1], x[2], x[3], &pbc,
 115                               r_ij, r_kj, r_kl, m, n, &cosine_angle,
 116                               &t1, &t2, &t3);
 117
 118             checker_.checkReal(angle, "angle");
 119             checker_.checkReal(cosine_angle, "cosine_angle");
 120             checker_.checkInteger(t1, "t1");
 121             checker_.checkInteger(t2, "t2");
 122             checker_.checkInteger(t3, "t3");
 123         }
 124
 125         void testIfunc(int                         ftype,
 126                        const std::vector<t_iatom> &iatoms,
 127                        const t_iparams             iparams[],
 128                        int                         epbc)
 129         {
 130             real  lambda = 0;
 131             real  dvdlambda;
 132             rvec4 f[NATOMS];
 133             for (int i = 0; i < NATOMS; i++)
 134             {
 135                 for (int j = 0; j < 4; j++)
 136                 {
 137                     f[i][j] = 0;
 138                 }
 139             }
 140             rvec  fshift[N_IVEC];
 141             clear_rvecs(N_IVEC, fshift);
 142             t_pbc pbc;
 143             set_pbc(&pbc, epbc, box);
 144             int   ddgatindex = 0;
 145             real  energy     = interaction_function[ftype].ifunc(iatoms.size(),
 146                                                                  iatoms.data(),
 147                                                                  iparams,
 148                                                                  x, f, fshift,
 149                                                                  &pbc,
 150                                                                  /* const struct t_graph *g */ nullptr,
 151                                                                  lambda, &dvdlambda,
 152                                                                  /* const struct t_mdatoms *md */ nullptr,
 153                                                                  /* struct t_fcdata *fcd */ nullptr,
 154                                                                  &ddgatindex);
 155             checker_.checkReal(energy, interaction_function[ftype].longname);
 156         }
 157
 158 };
 159
 160 TEST_F (BondedTest, BondAnglePbcNone)
 161 {
 162     testBondAngle(epbcNONE);
 163 }
 164
 165 TEST_F (BondedTest, BondAnglePbcXy)
 166 {
 167     testBondAngle(epbcXY);
 168 }
 169
 170 TEST_F (BondedTest, BondAnglePbcXyz)
 171 {
 172     testBondAngle(epbcXYZ);
 173 }
 174
 175 TEST_F (BondedTest, DihedralAnglePbcNone)
 176 {
 177     testDihedralAngle(epbcNONE);
 178 }
 179
 180 TEST_F (BondedTest, DihedralAnglePbcXy)
 181 {
 182     testDihedralAngle(epbcXY);
 183 }
 184
 185 TEST_F (BondedTest, DihedarlAnglePbcXyz)
 186 {
 187     testDihedralAngle(epbcXYZ);
 188 }
 189
 190 TEST_F (BondedTest, IfuncBondsPbcNo)
 191 {
 192     std::vector<t_iatom> iatoms = { 0, 0, 1, 0, 1, 2, 0, 2, 3 };
 193     t_iparams            iparams;
 194     iparams.harmonic.rA  = iparams.harmonic.rB  = 0.8;
 195     iparams.harmonic.krA = iparams.harmonic.krB = 50;
 196     testIfunc(F_BONDS, iatoms, &iparams, epbcNONE);
 197 }
 198
 199 TEST_F (BondedTest, IfuncBondsPbcXy)
 200 {
 201     std::vector<t_iatom> iatoms = { 0, 0, 1, 0, 1, 2, 0, 2, 3 };
 202     t_iparams            iparams;
 203     iparams.harmonic.rA  = iparams.harmonic.rB  = 0.8;
 204     iparams.harmonic.krA = iparams.harmonic.krB = 50;
 205     testIfunc(F_BONDS, iatoms, &iparams, epbcXY);
 206 }
 207
 208 TEST_F (BondedTest, IfuncBondsPbcXyz)
 209 {
 210     std::vector<t_iatom> iatoms = { 0, 0, 1, 0, 1, 2, 0, 2, 3 };
 211     t_iparams            iparams;
 212     iparams.harmonic.rA  = iparams.harmonic.rB  = 0.8;
 213     iparams.harmonic.krA = iparams.harmonic.krB = 50;
 214     testIfunc(F_BONDS, iatoms, &iparams, epbcXYZ);
 215 }
 216
 217 TEST_F (BondedTest, IfuncAnglesPbcNo)
 218 {
 219     std::vector<t_iatom> iatoms = { 0, 0, 1, 2, 0, 1, 2, 3 };
 220     t_iparams            iparams;
 221     real                 k = 50;
 222     iparams.harmonic.rA  = iparams.harmonic.rB  = 100;
 223     iparams.harmonic.krA = iparams.harmonic.krB = k;
 224     testIfunc(F_ANGLES, iatoms, &iparams, epbcNONE);
 225 }
 226
 227 TEST_F (BondedTest, IfuncAnglesPbcXy)
 228 {
 229     std::vector<t_iatom> iatoms  = { 0, 0, 1, 2, 0, 1, 2, 3 };
 230     t_iparams            iparams;
 231     real                 k = 50;
 232     iparams.harmonic.rA  = iparams.harmonic.rB  = 100;
 233     iparams.harmonic.krA = iparams.harmonic.krB = k;
 234     testIfunc(F_ANGLES, iatoms, &iparams, epbcXY);
 235 }
 236
 237 TEST_F (BondedTest, IfuncAnglesPbcXYZ)
 238 {
 239     std::vector<t_iatom> iatoms = { 0, 0, 1, 2, 0, 1, 2, 3 };
 240     t_iparams            iparams;
 241     real                 k = 50;
 242     iparams.harmonic.rA  = iparams.harmonic.rB  = 100;
 243     iparams.harmonic.krA = iparams.harmonic.krB = k;
 244     testIfunc(F_ANGLES, iatoms, &iparams, epbcXYZ);
 245 }
 246
 247 TEST_F (BondedTest, IfuncProperDihedralsPbcNo)
 248 {
 249     std::vector<t_iatom> iatoms = { 0, 0, 1, 2, 3 };
 250     t_iparams            iparams;
 251     iparams.pdihs.phiA = iparams.pdihs.phiB = -100;
 252     iparams.pdihs.cpA  = iparams.pdihs.cpB  = 10;
 253     iparams.pdihs.mult = 1;
 254     testIfunc(F_PDIHS, iatoms, &iparams, epbcNONE);
 255 }
 256
 257 TEST_F (BondedTest, IfuncProperDihedralsPbcXy)
 258 {
 259     std::vector<t_iatom> iatoms = { 0, 0, 1, 2, 3 };
 260     t_iparams            iparams;
 261     iparams.pdihs.phiA = iparams.pdihs.phiB = -100;
 262     iparams.pdihs.cpA  = iparams.pdihs.cpB  = 10;
 263     iparams.pdihs.mult = 1;
 264     testIfunc(F_PDIHS, iatoms, &iparams, epbcXY);
 265 }
 266
 267 TEST_F (BondedTest, IfuncProperDihedralsPbcXyz)
 268 {
 269     std::vector<t_iatom> iatoms  = { 0, 0, 1, 2, 3 };
 270     t_iparams            iparams;
 271     iparams.pdihs.phiA = iparams.pdihs.phiB = -100;
 272     iparams.pdihs.cpA  = iparams.pdihs.cpB  = 10;
 273     iparams.pdihs.mult = 1;
 274     testIfunc(F_PDIHS, iatoms, &iparams, epbcXYZ);
 275 }
 276
 277 }
 278
 279 }