src/gromacs/mdlib/tests/calc_verletbuf.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 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 Tests for the Verlet buffer calculation algorithm.
  37  *
  38  * \author Berk Hess <hess@kth.se>
  39  */
  40 #include "gmxpre.h"
  41
  42 #include "gromacs/mdlib/calc_verletbuf.h"
  43
  44 #include <algorithm>
  45
  46 #include <gtest/gtest.h>
  47
  48 #include "gromacs/math/functions.h"
  49
  50 #include "testutils/testasserts.h"
  51
  52 namespace gmx
  53 {
  54
  55 namespace
  56 {
  57
  58 class VerletBufferConstraintTest : public ::testing::Test
  59 {
  60 };
  61
  62 /* This test covers the displacement correction for constrained atoms.
  63  * This test does not check exact values, but rather checks that the MSD
  64  * estimate for a constrained atom is smaller than that of a free atom
  65  * and checks that the MSD is not smaller and also not much larger
  66  * than the maximum of the exact value for rotational MSD beyond
  67  * the location of the maximum. Furthermore, we check that the MSD estimate
  68  * never decreases, as this is a requirement for the Verlet buffer size
  69  * estimation. Together these criteria provide tight margins on
  70  * the shape and values of the estimate.
  71  *
  72  * Additionally we check the 3D MSD for the COM of the two atoms.
  73  */
  74 TEST_F(VerletBufferConstraintTest, EqualMasses)
  75 {
  76     // The location and value of the MSD maximum for the exact displacement
  77     // is described in the source file. We need to divide the maximum given
  78     // there by 2, since sigma2 is per DOF for the 2 DOF constraint rotation.
  79     const real sigma2RelMaxLocation  = 4.5119;
  80     const real sigma2RelMaxValue     = 2.5695/2;
  81
  82     // Our max of our current estimate is 3% above the exact value.
  83     const real sigma2RelMaxMargin    = 1.04;
  84
  85     // The exact parameter values here don't actually matter.
  86     real mass = 10;
  87     real arm  = 0.1;
  88
  89     atom_nonbonded_kinetic_prop_t prop;
  90     prop.mass     = mass;
  91     prop.type     = -1;
  92     prop.q        = 0;
  93     prop.bConstr  = TRUE;
  94     prop.con_mass = mass;
  95     prop.con_len  = 2*arm;
  96
  97     // We scan a range of rotation distributions by scanning over T.
  98     int  numPointsBeforeMax = 0;
  99     int  numPointsAfterMax  = 0;
 100     real sigma2_2d_prev     = 0;
 101     for (int i = 0; i <= 200; i++)
 102     {
 103         real ktFac = i*0.01;
 104         // The rotational displacement is Gaussian with a sigma^2 of:
 105         real sigma2_rot = ktFac/(2*mass);
 106
 107         // Get the estimate for the Cartesian displacement.
 108         real sigma2_2d, sigma2_3d;
 109         constrained_atom_sigma2(ktFac, &prop, &sigma2_2d, &sigma2_3d);
 110
 111         // Check we are not decreasing sigma2_2d
 112         EXPECT_EQ(std::max(sigma2_2d_prev, sigma2_2d), sigma2_2d);
 113         // Check that sigma2_2d is not larger than sigma2 for free motion.
 114         EXPECT_EQ(std::min(sigma2_rot, sigma2_2d), sigma2_2d);
 115
 116         // Check that we don't underestimate sigma2_rot beyond the real maximum
 117         // and that our overestimate is tight.
 118         real sigma2Rel = sigma2_rot/gmx::square(arm);
 119         if (sigma2Rel >= sigma2RelMaxLocation)
 120         {
 121             EXPECT_EQ(std::max(sigma2_2d, sigma2RelMaxValue*gmx::square(arm)), sigma2_2d);
 122             EXPECT_EQ(std::min(sigma2_2d, sigma2RelMaxMargin*sigma2RelMaxValue*gmx::square(arm)), sigma2_2d);
 123
 124             numPointsAfterMax++;
 125         }
 126         else
 127         {
 128             numPointsBeforeMax++;
 129         }
 130
 131         // Also check sigma2 for the COM of the two atoms
 132         EXPECT_EQ(sigma2_rot, sigma2_3d);
 133     }
 134
 135     GMX_RELEASE_ASSERT(numPointsBeforeMax >= 20 && numPointsAfterMax >= 20, "This test only provides full coverage when we test a sufficient number of points before and after the location of the maximum value for the exact formula.");
 136 }
 137
 138 }      // namespace anonymous
 139
 140 }      // namespace gmx