From 7861aae58fea4aaa1dfbfe4ef26e94e9b3e55ff9 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Tue, 21 Mar 2017 16:30:10 +0100 Subject: [PATCH] Improve Verlet buffer constraint estimate The displacement estimate for a constrained atom (typically H) rotating around the COM with a partner atom was not thoroughly derived and the code documentation was not fully correct. Now there is an exact MSD integral. Approximating that gives the same v/(1+av+bv^2) variance scaling formula with a=1/6, but b changed from 5/180 to 8/180, slightly lowering the buffer size. Note that we (still) use a Gaussian with matched variance, which results in a much larger buffer than necessary, since the tail of the distribution sets the buffer size and the Gaussian has a long tail whereas the actual distribution has no tail. Also added a test for this estimate. Change-Id: Ie6fb817dcd55177e5992facc7b68616f318572a3 --- src/gromacs/mdlib/calc_verletbuf.cpp | 107 +++++++++------------- src/gromacs/mdlib/calc_verletbuf.h | 30 ++++++- src/gromacs/mdlib/tests/CMakeLists.txt | 3 +- src/gromacs/mdlib/tests/calc_verletbuf.cpp | 140 +++++++++++++++++++++++++++++ 4 files changed, 213 insertions(+), 67 deletions(-) create mode 100644 src/gromacs/mdlib/tests/calc_verletbuf.cpp diff --git a/src/gromacs/mdlib/calc_verletbuf.cpp b/src/gromacs/mdlib/calc_verletbuf.cpp index 19f8ee73ed..fcc6b690e5 100644 --- a/src/gromacs/mdlib/calc_verletbuf.cpp +++ b/src/gromacs/mdlib/calc_verletbuf.cpp @@ -95,20 +95,6 @@ */ typedef struct { - real mass; /* mass */ - int type; /* type (used for LJ parameters) */ - real q; /* charge */ - gmx_bool bConstr; /* constrained, if TRUE, use #DOF=2 iso 3 */ - real con_mass; /* mass of heaviest atom connected by constraints */ - real con_len; /* constraint length to the heaviest atom */ -} atom_nonbonded_kinetic_prop_t; - -/* Struct for unique atom type for calculating the energy drift. - * The atom displacement depends on mass and constraints. - * The energy jump for given distance depend on LJ type and q. - */ -typedef struct -{ atom_nonbonded_kinetic_prop_t prop; /* non-bonded and kinetic atom prop. */ int n; /* #atoms of this type in the system */ } verletbuf_atomtype_t; @@ -479,69 +465,60 @@ static void get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop, * into account. If an atom has multiple constraints, this will result in * an overestimate of the displacement, which gives a larger drift and buffer. */ -static void constrained_atom_sigma2(real kT_fac, - const atom_nonbonded_kinetic_prop_t *prop, - real *sigma2_2d, - real *sigma2_3d) +void constrained_atom_sigma2(real kT_fac, + const atom_nonbonded_kinetic_prop_t *prop, + real *sigma2_2d, + real *sigma2_3d) { - real sigma2_rot; - real com_dist; - real sigma2_rel; - real scale; - /* Here we decompose the motion of a constrained atom into two * components: rotation around the COM and translation of the COM. */ - /* Determine the variance for the displacement of the rotational mode */ - sigma2_rot = kT_fac/(prop->mass*(prop->mass + prop->con_mass)/prop->con_mass); + /* Determine the variance of the arc length for the two rotational DOFs */ + real massFraction = prop->con_mass/(prop->mass + prop->con_mass); + real sigma2_rot = kT_fac*massFraction/prop->mass; /* The distance from the atom to the COM, i.e. the rotational arm */ - com_dist = prop->con_len*prop->con_mass/(prop->mass + prop->con_mass); + real comDistance = prop->con_len*massFraction; /* The variance relative to the arm */ - sigma2_rel = sigma2_rot/(com_dist*com_dist); - /* At 6 the scaling formula has slope 0, - * so we keep sigma2_2d constant after that. + real sigma2_rel = sigma2_rot/gmx::square(comDistance); + + /* For sigma2_rel << 1 we don't notice the rotational effect and + * we have a normal, Gaussian displacement distribution. + * For larger sigma2_rel the displacement is much less, in fact it can + * not exceed 2*comDistance. We can calculate MSD/arm^2 as: + * integral_x=0-inf distance2(x) x/sigma2_rel exp(-x^2/(2 sigma2_rel)) dx + * where x is angular displacement and distance2(x) is the distance^2 + * between points at angle 0 and x: + * distance2(x) = (sin(x) - sin(0))^2 + (cos(x) - cos(0))^2 + * The limiting value of this MSD is 2, which is also the value for + * a uniform rotation distribution that would be reached at long time. + * The maximum is 2.5695 at sigma2_rel = 4.5119. + * We approximate this integral with a rational polynomial with + * coefficients from a Taylor expansion. This approximation is an + * overestimate for all values of sigma2_rel. Its maximum value + * of 2.6491 is reached at sigma2_rel = sqrt(45/2) = 4.7434. + * We keep the approximation constant after that. + * We use this approximate MSD as the variance for a Gaussian distribution. + * + * NOTE: For any sensible buffer tolerance this will result in a (large) + * overestimate of the buffer size, since the Gaussian has a long tail, + * whereas the actual distribution can not reach values larger than 2. */ - if (sigma2_rel < 6) - { - /* A constrained atom rotates around the atom it is constrained to. - * This results in a smaller linear displacement than for a free atom. - * For a perfectly circular displacement, this lowers the displacement - * by: 1/arcsin(arc_length) - * and arcsin(x) = 1 + x^2/6 + ... - * For sigma2_rel<<1 the displacement distribution is erfc - * (exact formula is provided below). For larger sigma, it is clear - * that the displacement can't be larger than 2*com_dist. - * It turns out that the distribution becomes nearly uniform. - * For intermediate sigma2_rel, scaling down sigma with the third - * order expansion of arcsin with argument sigma_rel turns out - * to give a very good approximation of the distribution and variance. - * Even for larger values, the variance is only slightly overestimated. - * Note that the most relevant displacements are in the long tail. - * This rotation approximation always overestimates the tail (which - * runs to infinity, whereas it should be <= 2*com_dist). - * Thus we always overestimate the drift and the buffer size. - */ - scale = 1/(1 + sigma2_rel/6); - *sigma2_2d = sigma2_rot*scale*scale; - } - else - { - /* sigma_2d is set to the maximum given by the scaling above. - * For large sigma2 the real displacement distribution is close - * to uniform over -2*con_len to 2*com_dist. - * Our erfc with sigma_2d=sqrt(1.5)*com_dist (which means the sigma - * of the erfc output distribution is con_dist) overestimates - * the variance and additionally has a long tail. This means - * we have a (safe) overestimation of the drift. - */ - *sigma2_2d = 1.5*com_dist*com_dist; - } + /* Coeffients obtained from a Taylor expansion */ + const real a = 1.0/3.0; + const real b = 2.0/45.0; + + /* Our approximation is constant after sigma2_rel = 1/sqrt(b) */ + sigma2_rel = std::min(sigma2_rel, 1/std::sqrt(b)); + + /* Compute the approximate sigma^2 for 2D motion due to the rotation */ + *sigma2_2d = gmx::square(comDistance)* + sigma2_rel/(1 + a*sigma2_rel + b*gmx::square(sigma2_rel)); /* The constrained atom also moves (in 3D) with the COM of both atoms */ - *sigma2_3d = kT_fac/(prop->mass + prop->con_mass); + *sigma2_3d = kT_fac/(prop->mass + prop->con_mass); } static void get_atom_sigma2(real kT_fac, diff --git a/src/gromacs/mdlib/calc_verletbuf.h b/src/gromacs/mdlib/calc_verletbuf.h index 5ceccb8404..097e900d78 100644 --- a/src/gromacs/mdlib/calc_verletbuf.h +++ b/src/gromacs/mdlib/calc_verletbuf.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by + * Copyright (c) 2012,2013,2014,2015,2017, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -91,6 +91,34 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, int *n_nonlin_vsite, real *rlist); +/* Struct for unique atom type for calculating the energy drift. + * The atom displacement depends on mass and constraints. + * The energy jump for given distance depend on LJ type and q. + */ +struct atom_nonbonded_kinetic_prop_t +{ + real mass; /* mass */ + int type; /* type (used for LJ parameters) */ + real q; /* charge */ + gmx_bool bConstr; /* constrained, if TRUE, use #DOF=2 iso 3 */ + real con_mass; /* mass of heaviest atom connected by constraints */ + real con_len; /* constraint length to the heaviest atom */ +}; + +/* This function computes two components of the estimate of the variance + * in the displacement of one atom in a system of two constrained atoms. + * Returns in sigma2_2d the variance due to rotation of the constrained + * atom around the atom to which it constrained. + * Returns in sigma2_3d the variance due to displacement of the COM + * of the whole system of the two constrained atoms. + * + * Only exposed here for testing purposes. + */ +void constrained_atom_sigma2(real kT_fac, + const atom_nonbonded_kinetic_prop_t *prop, + real *sigma2_2d, + real *sigma2_3d); + #ifdef __cplusplus } #endif diff --git a/src/gromacs/mdlib/tests/CMakeLists.txt b/src/gromacs/mdlib/tests/CMakeLists.txt index 0ad14a2204..fbcf73f88f 100644 --- a/src/gromacs/mdlib/tests/CMakeLists.txt +++ b/src/gromacs/mdlib/tests/CMakeLists.txt @@ -1,7 +1,7 @@ # # This file is part of the GROMACS molecular simulation package. # -# Copyright (c) 2014,2016, by the GROMACS development team, led by +# Copyright (c) 2014,2016,2017, by the GROMACS development team, led by # Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, # and including many others, as listed in the AUTHORS file in the # top-level source directory and at http://www.gromacs.org. @@ -33,6 +33,7 @@ # the research papers on the package. Check out http://www.gromacs.org. gmx_add_unit_test(MdlibUnitTest mdlib-test + calc_verletbuf.cpp settle.cpp shake.cpp simulationsignal.cpp) diff --git a/src/gromacs/mdlib/tests/calc_verletbuf.cpp b/src/gromacs/mdlib/tests/calc_verletbuf.cpp new file mode 100644 index 0000000000..2e1531efe1 --- /dev/null +++ b/src/gromacs/mdlib/tests/calc_verletbuf.cpp @@ -0,0 +1,140 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2017, by the GROMACS development team, led by + * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, + * and including many others, as listed in the AUTHORS file in the + * top-level source directory and at http://www.gromacs.org. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ +/*! \internal \file + * \brief Tests for the Verlet buffer calculation algorithm. + * + * \author Berk Hess + */ +#include "gmxpre.h" + +#include "gromacs/mdlib/calc_verletbuf.h" + +#include + +#include + +#include "gromacs/math/functions.h" + +#include "testutils/testasserts.h" + +namespace gmx +{ + +namespace +{ + +class VerletBufferConstraintTest : public ::testing::Test +{ +}; + +/* This test covers the displacement correction for constrained atoms. + * This test does not check exact values, but rather checks that the MSD + * estimate for a constrained atom is smaller than that of a free atom + * and checks that the MSD is not smaller and also not much larger + * than the maximum of the exact value for rotational MSD beyond + * the location of the maximum. Furthermore, we check that the MSD estimate + * never decreases, as this is a requirement for the Verlet buffer size + * estimation. Together these criteria provide tight margins on + * the shape and values of the estimate. + * + * Additionally we check the 3D MSD for the COM of the two atoms. + */ +TEST_F(VerletBufferConstraintTest, EqualMasses) +{ + // The location and value of the MSD maximum for the exact displacement + // is described in the source file. We need to divide the maximum given + // there by 2, since sigma2 is per DOF for the 2 DOF constraint rotation. + const real sigma2RelMaxLocation = 4.5119; + const real sigma2RelMaxValue = 2.5695/2; + + // Our max of our current estimate is 3% above the exact value. + const real sigma2RelMaxMargin = 1.04; + + // The exact parameter values here don't actually matter. + real mass = 10; + real arm = 0.1; + + atom_nonbonded_kinetic_prop_t prop; + prop.mass = mass; + prop.type = -1; + prop.q = 0; + prop.bConstr = TRUE; + prop.con_mass = mass; + prop.con_len = 2*arm; + + // We scan a range of rotation distributions by scanning over T. + int numPointsBeforeMax = 0; + int numPointsAfterMax = 0; + real sigma2_2d_prev = 0; + for (int i = 0; i <= 200; i++) + { + real ktFac = i*0.01; + // The rotational displacement is Gaussian with a sigma^2 of: + real sigma2_rot = ktFac/(2*mass); + + // Get the estimate for the Cartesian displacement. + real sigma2_2d, sigma2_3d; + constrained_atom_sigma2(ktFac, &prop, &sigma2_2d, &sigma2_3d); + + // Check we are not decreasing sigma2_2d + EXPECT_EQ(std::max(sigma2_2d_prev, sigma2_2d), sigma2_2d); + // Check that sigma2_2d is not larger than sigma2 for free motion. + EXPECT_EQ(std::min(sigma2_rot, sigma2_2d), sigma2_2d); + + // Check that we don't underestimate sigma2_rot beyond the real maximum + // and that our overestimate is tight. + real sigma2Rel = sigma2_rot/gmx::square(arm); + if (sigma2Rel >= sigma2RelMaxLocation) + { + EXPECT_EQ(std::max(sigma2_2d, sigma2RelMaxValue*gmx::square(arm)), sigma2_2d); + EXPECT_EQ(std::min(sigma2_2d, sigma2RelMaxMargin*sigma2RelMaxValue*gmx::square(arm)), sigma2_2d); + + numPointsAfterMax++; + } + else + { + numPointsBeforeMax++; + } + + // Also check sigma2 for the COM of the two atoms + EXPECT_EQ(sigma2_rot, sigma2_3d); + } + + 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."); +} + +} // namespace anonymous + +} // namespace gmx -- 2.11.4.GIT