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  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,
 116                               &t1, &t2, &t3);
 117
 118             checker_.checkReal(angle, "angle");
 119             checker_.checkInteger(t1, "t1");
 120             checker_.checkInteger(t2, "t2");
 121             checker_.checkInteger(t3, "t3");
 122         }
 123
 124         void testIfunc(int                         ftype,
 125                        const std::vector<t_iatom> &iatoms,
 126                        const t_iparams             iparams[],
 127                        int                         epbc)
 128         {
 129             real  lambda    = 0;
 130             real  dvdlambda = 0;
 131             rvec4 f[NATOMS];
 132             for (int i = 0; i < NATOMS; i++)
 133             {
 134                 for (int j = 0; j < 4; j++)
 135                 {
 136                     f[i][j] = 0;
 137                 }
 138             }
 139             rvec  fshift[N_IVEC];
 140             clear_rvecs(N_IVEC, fshift);
 141             t_pbc pbc;
 142             set_pbc(&pbc, epbc, box);
 143             int   ddgatindex = 0;
 144             real  energy     = interaction_function[ftype].ifunc(iatoms.size(),
 145                                                                  iatoms.data(),
 146                                                                  iparams,
 147                                                                  x, f, fshift,
 148                                                                  &pbc,
 149                                                                  /* const struct t_graph *g */ nullptr,
 150                                                                  lambda, &dvdlambda,
 151                                                                  /* const struct t_mdatoms *md */ nullptr,
 152                                                                  /* struct t_fcdata *fcd */ nullptr,
 153                                                                  &ddgatindex);
 154             checker_.checkReal(energy, interaction_function[ftype].longname);
 155         }
 156
 157 };
 158
 159 TEST_F (BondedTest, BondAnglePbcNone)
 160 {
 161     testBondAngle(epbcNONE);
 162 }
 163
 164 TEST_F (BondedTest, BondAnglePbcXy)
 165 {
 166     testBondAngle(epbcXY);
 167 }
 168
 169 TEST_F (BondedTest, BondAnglePbcXyz)
 170 {
 171     testBondAngle(epbcXYZ);
 172 }
 173
 174 TEST_F (BondedTest, DihedralAnglePbcNone)
 175 {
 176     testDihedralAngle(epbcNONE);
 177 }
 178
 179 TEST_F (BondedTest, DihedralAnglePbcXy)
 180 {
 181     testDihedralAngle(epbcXY);
 182 }
 183
 184 TEST_F (BondedTest, DihedralAnglePbcXyz)
 185 {
 186     testDihedralAngle(epbcXYZ);
 187 }
 188
 189 TEST_F (BondedTest, IfuncBondsPbcNo)
 190 {
 191     std::vector<t_iatom> iatoms = { 0, 0, 1, 0, 1, 2, 0, 2, 3 };
 192     t_iparams            iparams;
 193     iparams.harmonic.rA  = iparams.harmonic.rB  = 0.8;
 194     iparams.harmonic.krA = iparams.harmonic.krB = 50;
 195     testIfunc(F_BONDS, iatoms, &iparams, epbcNONE);
 196 }
 197
 198 TEST_F (BondedTest, IfuncBondsPbcXy)
 199 {
 200     std::vector<t_iatom> iatoms = { 0, 0, 1, 0, 1, 2, 0, 2, 3 };
 201     t_iparams            iparams;
 202     iparams.harmonic.rA  = iparams.harmonic.rB  = 0.8;
 203     iparams.harmonic.krA = iparams.harmonic.krB = 50;
 204     testIfunc(F_BONDS, iatoms, &iparams, epbcXY);
 205 }
 206
 207 TEST_F (BondedTest, IfuncBondsPbcXyz)
 208 {
 209     std::vector<t_iatom> iatoms = { 0, 0, 1, 0, 1, 2, 0, 2, 3 };
 210     t_iparams            iparams;
 211     iparams.harmonic.rA  = iparams.harmonic.rB  = 0.8;
 212     iparams.harmonic.krA = iparams.harmonic.krB = 50;
 213     testIfunc(F_BONDS, iatoms, &iparams, epbcXYZ);
 214 }
 215
 216 TEST_F (BondedTest, IfuncAnglesPbcNo)
 217 {
 218     std::vector<t_iatom> iatoms = { 0, 0, 1, 2, 0, 1, 2, 3 };
 219     t_iparams            iparams;
 220     real                 k = 50;
 221     iparams.harmonic.rA  = iparams.harmonic.rB  = 100;
 222     iparams.harmonic.krA = iparams.harmonic.krB = k;
 223     testIfunc(F_ANGLES, iatoms, &iparams, epbcNONE);
 224 }
 225
 226 TEST_F (BondedTest, IfuncAnglesPbcXy)
 227 {
 228     std::vector<t_iatom> iatoms  = { 0, 0, 1, 2, 0, 1, 2, 3 };
 229     t_iparams            iparams;
 230     real                 k = 50;
 231     iparams.harmonic.rA  = iparams.harmonic.rB  = 100;
 232     iparams.harmonic.krA = iparams.harmonic.krB = k;
 233     testIfunc(F_ANGLES, iatoms, &iparams, epbcXY);
 234 }
 235
 236 TEST_F (BondedTest, IfuncAnglesPbcXYZ)
 237 {
 238     std::vector<t_iatom> iatoms = { 0, 0, 1, 2, 0, 1, 2, 3 };
 239     t_iparams            iparams;
 240     real                 k = 50;
 241     iparams.harmonic.rA  = iparams.harmonic.rB  = 100;
 242     iparams.harmonic.krA = iparams.harmonic.krB = k;
 243     testIfunc(F_ANGLES, iatoms, &iparams, epbcXYZ);
 244 }
 245
 246 TEST_F (BondedTest, IfuncProperDihedralsPbcNo)
 247 {
 248     std::vector<t_iatom> iatoms = { 0, 0, 1, 2, 3 };
 249     t_iparams            iparams;
 250     iparams.pdihs.phiA = iparams.pdihs.phiB = -100;
 251     iparams.pdihs.cpA  = iparams.pdihs.cpB  = 10;
 252     iparams.pdihs.mult = 1;
 253     testIfunc(F_PDIHS, iatoms, &iparams, epbcNONE);
 254 }
 255
 256 TEST_F (BondedTest, IfuncProperDihedralsPbcXy)
 257 {
 258     std::vector<t_iatom> iatoms = { 0, 0, 1, 2, 3 };
 259     t_iparams            iparams;
 260     iparams.pdihs.phiA = iparams.pdihs.phiB = -100;
 261     iparams.pdihs.cpA  = iparams.pdihs.cpB  = 10;
 262     iparams.pdihs.mult = 1;
 263     testIfunc(F_PDIHS, iatoms, &iparams, epbcXY);
 264 }
 265
 266 TEST_F (BondedTest, IfuncProperDihedralsPbcXyz)
 267 {
 268     std::vector<t_iatom> iatoms  = { 0, 0, 1, 2, 3 };
 269     t_iparams            iparams;
 270     iparams.pdihs.phiA = iparams.pdihs.phiB = -100;
 271     iparams.pdihs.cpA  = iparams.pdihs.cpB  = 10;
 272     iparams.pdihs.mult = 1;
 273     testIfunc(F_PDIHS, iatoms, &iparams, epbcXYZ);
 274 }
 275
 276 }
 277
 278 }