*/
-/* 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
-{
- 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.
* 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,
/*
* 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.
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
#
# 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.
# 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)
--- /dev/null
+/*
+ * 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 <hess@kth.se>
+ */
+#include "gmxpre.h"
+
+#include "gromacs/mdlib/calc_verletbuf.h"
+
+#include <algorithm>
+
+#include <gtest/gtest.h>
+
+#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
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff