typedef struct gmx_shakedata
{
rvec *rij;
- real *M2;
- real *tt;
- real *dist2;
+ real *half_of_reduced_mass;
+ real *distance_squared_tolerance;
+ real *constraint_distance_squared;
int nalloc;
/* SOR stuff */
real delta;
snew(d, 1);
- d->nalloc = 0;
- d->rij = NULL;
- d->M2 = NULL;
- d->tt = NULL;
- d->dist2 = NULL;
+ d->nalloc = 0;
+ d->rij = NULL;
+ d->half_of_reduced_mass = NULL;
+ d->distance_squared_tolerance = NULL;
+ d->constraint_distance_squared = NULL;
/* SOR initialization */
d->delta = 0.1;
fflush(log);
}
-void cshake(atom_id iatom[], int ncon, int *nnit, int maxnit,
- real dist2[], real xp[], real rij[], real m2[], real omega,
- real invmass[], real tt[], real lagr[], int *nerror)
+/*! \brief Inner kernel for SHAKE constraints
+ *
+ * Original implementation from R.C. van Schaik and W.F. van Gunsteren
+ * (ETH Zuerich, June 1992), adapted for GROMACS by David van der
+ * Spoel November 1992.
+ *
+ * The algorithm here is based section five of Ryckaert, Ciccotti and
+ * Berendsen, J Comp Phys, 23, 327, 1977.
+ *
+ * \param[in] iatom Mini-topology of triples of constraint type (unused in this
+ * function) and indices of the two atoms involved
+ * \param[in] ncon Number of constraints
+ * \param[out] nnit Number of iterations performed
+ * \param[in] maxnit Maximum number of iterations permitted
+ * \param[in] constraint_distance_squared The objective value for each constraint
+ * \param[inout] positions The initial (and final) values of the positions of all atoms
+ * \param[in] initial_displacements The initial displacements of each constraint
+ * \param[in] half_of_reduced_mass Half of the reduced mass for each constraint
+ * \param[in] omega SHAKE over-relaxation factor (set non-1.0 by
+ * using shake-sor=yes in the .mdp, but there is no documentation anywhere)
+ * \param[in] invmass Inverse mass of each atom
+ * \param[in] distance_squared_tolerance Multiplicative tolerance on the difference in the
+ * square of the constrained distance (see code)
+ * \param[out] scaled_lagrange_multiplier Scaled Lagrange multiplier for each constraint (-2 * eta from p. 336
+ * of the paper, divided by the constraint distance)
+ * \param[out] nerror Zero upon success, returns one more than the index of the
+ * problematic constraint if the input was malformed
+ *
+ * \todo Make SHAKE use better data structures, in particular for iatom. */
+void cshake(const atom_id iatom[], int ncon, int *nnit, int maxnit,
+ const real constraint_distance_squared[], real positions[],
+ const real initial_displacements[], const real half_of_reduced_mass[], real omega,
+ const real invmass[], const real distance_squared_tolerance[],
+ real scaled_lagrange_multiplier[], int *nerror)
{
- /*
- * r.c. van schaik and w.f. van gunsteren
- * eth zuerich
- * june 1992
- * Adapted for use with Gromacs by David van der Spoel november 92 and later.
- */
/* default should be increased! MRS 8/4/2009 */
- const real mytol = 1e-10;
-
- int ll, i, j, i3, j3, l3;
- int ix, iy, iz, jx, jy, jz;
- real toler, rpij2, rrpr, tx, ty, tz, diff, acor, im, jm;
- real xh, yh, zh, rijx, rijy, rijz;
- real tix, tiy, tiz;
- real tjx, tjy, tjz;
- int nit, error, nconv;
- real iconvf;
+ const real mytol = 1e-10;
+
+ int ll, i, j, i3, j3, l3;
+ int ix, iy, iz, jx, jy, jz;
+ real r_dot_r_prime;
+ real constraint_distance_squared_ll;
+ real r_prime_squared;
+ real scaled_lagrange_multiplier_ll;
+ real r_prime_x, r_prime_y, r_prime_z, diff, im, jm;
+ real xh, yh, zh, rijx, rijy, rijz;
+ real tix, tiy, tiz;
+ real tjx, tjy, tjz;
+ int nit, error, nconv;
+ real iconvf;
error = 0;
nconv = 1;
for (ll = 0; (ll < ncon) && (error == 0); ll++)
{
l3 = 3*ll;
- rijx = rij[l3+XX];
- rijy = rij[l3+YY];
- rijz = rij[l3+ZZ];
+ rijx = initial_displacements[l3+XX];
+ rijy = initial_displacements[l3+YY];
+ rijz = initial_displacements[l3+ZZ];
i = iatom[l3+1];
j = iatom[l3+2];
i3 = 3*i;
jy = j3+YY;
jz = j3+ZZ;
- tx = xp[ix]-xp[jx];
- ty = xp[iy]-xp[jy];
- tz = xp[iz]-xp[jz];
- rpij2 = tx*tx+ty*ty+tz*tz;
- toler = dist2[ll];
- diff = toler-rpij2;
+ /* Compute r prime between atoms i and j, which is the
+ displacement *before* this update stage */
+ r_prime_x = positions[ix]-positions[jx];
+ r_prime_y = positions[iy]-positions[jy];
+ r_prime_z = positions[iz]-positions[jz];
+ r_prime_squared = (r_prime_x * r_prime_x +
+ r_prime_y * r_prime_y +
+ r_prime_z * r_prime_z);
+ constraint_distance_squared_ll = constraint_distance_squared[ll];
+ diff = constraint_distance_squared_ll - r_prime_squared;
/* iconvf is less than 1 when the error is smaller than a bound */
- /* But if tt is too big, then it will result in looping in iconv */
-
- iconvf = fabs(diff)*tt[ll];
+ iconvf = fabs(diff) * distance_squared_tolerance[ll];
- if (iconvf > 1)
+ if (iconvf > 1.0)
{
- nconv = iconvf;
- rrpr = rijx*tx+rijy*ty+rijz*tz;
+ nconv = iconvf;
+ r_dot_r_prime = (rijx * r_prime_x +
+ rijy * r_prime_y +
+ rijz * r_prime_z);
- if (rrpr < toler*mytol)
+ if (r_dot_r_prime < constraint_distance_squared_ll * mytol)
{
error = ll+1;
}
else
{
- acor = omega*diff*m2[ll]/rrpr;
- lagr[ll] += acor;
- xh = rijx*acor;
- yh = rijy*acor;
- zh = rijz*acor;
- im = invmass[i];
- jm = invmass[j];
- xp[ix] += xh*im;
- xp[iy] += yh*im;
- xp[iz] += zh*im;
- xp[jx] -= xh*jm;
- xp[jy] -= yh*jm;
- xp[jz] -= zh*jm;
+ /* The next line solves equation 5.6 (neglecting
+ the term in g^2), for g */
+ scaled_lagrange_multiplier_ll = omega*diff*half_of_reduced_mass[ll]/r_dot_r_prime;
+ scaled_lagrange_multiplier[ll] += scaled_lagrange_multiplier_ll;
+ xh = rijx * scaled_lagrange_multiplier_ll;
+ yh = rijy * scaled_lagrange_multiplier_ll;
+ zh = rijz * scaled_lagrange_multiplier_ll;
+ im = invmass[i];
+ jm = invmass[j];
+ positions[ix] += xh*im;
+ positions[iy] += yh*im;
+ positions[iz] += zh*im;
+ positions[jx] -= xh*jm;
+ positions[jy] -= yh*jm;
+ positions[jz] -= zh*jm;
}
}
}
real invmass[], int ncon,
t_iparams ip[], t_iatom *iatom,
real tol, rvec x[], rvec prime[], real omega,
- gmx_bool bFEP, real lambda, real lagr[],
+ gmx_bool bFEP, real lambda, real scaled_lagrange_multiplier[],
real invdt, rvec *v,
gmx_bool bCalcVir, tensor vir_r_m_dr, int econq,
t_vetavars *vetavar)
{
rvec *rij;
- real *M2, *tt, *dist2;
+ real *half_of_reduced_mass, *distance_squared_tolerance, *constraint_distance_squared;
int maxnit = 1000;
- int nit = 0, ll, i, j, type;
+ int nit = 0, ll, i, j, d, d2, type;
t_iatom *ia;
- real L1, tol2, toler;
+ real L1;
real mm = 0., tmp;
int error = 0;
- real g, vscale, rscale, rvscale;
+ real vscale, rscale, rvscale;
+ real constraint_distance;
if (ncon > shaked->nalloc)
{
shaked->nalloc = over_alloc_dd(ncon);
srenew(shaked->rij, shaked->nalloc);
- srenew(shaked->M2, shaked->nalloc);
- srenew(shaked->tt, shaked->nalloc);
- srenew(shaked->dist2, shaked->nalloc);
+ srenew(shaked->half_of_reduced_mass, shaked->nalloc);
+ srenew(shaked->distance_squared_tolerance, shaked->nalloc);
+ srenew(shaked->constraint_distance_squared, shaked->nalloc);
}
- rij = shaked->rij;
- M2 = shaked->M2;
- tt = shaked->tt;
- dist2 = shaked->dist2;
+ rij = shaked->rij;
+ half_of_reduced_mass = shaked->half_of_reduced_mass;
+ distance_squared_tolerance = shaked->distance_squared_tolerance;
+ constraint_distance_squared = shaked->constraint_distance_squared;
L1 = 1.0-lambda;
- tol2 = 2.0*tol;
ia = iatom;
for (ll = 0; (ll < ncon); ll++, ia += 3)
{
i = ia[1];
j = ia[2];
- mm = 2*(invmass[i]+invmass[j]);
- rij[ll][XX] = x[i][XX]-x[j][XX];
- rij[ll][YY] = x[i][YY]-x[j][YY];
- rij[ll][ZZ] = x[i][ZZ]-x[j][ZZ];
- M2[ll] = 1.0/mm;
+ mm = 2.0*(invmass[i]+invmass[j]);
+ rij[ll][XX] = x[i][XX]-x[j][XX];
+ rij[ll][YY] = x[i][YY]-x[j][YY];
+ rij[ll][ZZ] = x[i][ZZ]-x[j][ZZ];
+ half_of_reduced_mass[ll] = 1.0/mm;
if (bFEP)
{
- toler = sqr(L1*ip[type].constr.dA + lambda*ip[type].constr.dB);
+ constraint_distance = L1*ip[type].constr.dA + lambda*ip[type].constr.dB;
}
else
{
- toler = sqr(ip[type].constr.dA);
+ constraint_distance = ip[type].constr.dA;
}
- dist2[ll] = toler;
- tt[ll] = 1.0/(toler*tol2);
+ constraint_distance_squared[ll] = sqr(constraint_distance);
+ distance_squared_tolerance[ll] = 0.5/(constraint_distance_squared[ll]*tol);
}
switch (econq)
{
case econqCoord:
- cshake(iatom, ncon, &nit, maxnit, dist2, prime[0], rij[0], M2, omega, invmass, tt, lagr, &error);
+ cshake(iatom, ncon, &nit, maxnit, constraint_distance_squared, prime[0], rij[0], half_of_reduced_mass, omega, invmass, distance_squared_tolerance, scaled_lagrange_multiplier, &error);
break;
case econqVeloc:
- crattle(iatom, ncon, &nit, maxnit, dist2, prime[0], rij[0], M2, omega, invmass, tt, lagr, &error, invdt, vetavar);
+ crattle(iatom, ncon, &nit, maxnit, constraint_distance_squared, prime[0], rij[0], half_of_reduced_mass, omega, invmass, distance_squared_tolerance, scaled_lagrange_multiplier, &error, invdt, vetavar);
break;
}
nit = 0;
}
- /* Constraint virial and correct the lagrange multipliers for the length */
+ /* Constraint virial and correct the Lagrange multipliers for the length */
ia = iatom;
for (ll = 0; (ll < ncon); ll++, ia += 3)
{
+ type = ia[0];
+ i = ia[1];
+ j = ia[2];
if ((econq == econqCoord) && v != NULL)
{
/* Correct the velocities */
- mm = lagr[ll]*invmass[ia[1]]*invdt/vetavar->rscale;
- for (i = 0; i < DIM; i++)
+ mm = scaled_lagrange_multiplier[ll]*invmass[i]*invdt/vetavar->rscale;
+ for (d = 0; d < DIM; d++)
{
- v[ia[1]][i] += mm*rij[ll][i];
+ v[ia[1]][d] += mm*rij[ll][d];
}
- mm = lagr[ll]*invmass[ia[2]]*invdt/vetavar->rscale;
- for (i = 0; i < DIM; i++)
+ mm = scaled_lagrange_multiplier[ll]*invmass[j]*invdt/vetavar->rscale;
+ for (d = 0; d < DIM; d++)
{
- v[ia[2]][i] -= mm*rij[ll][i];
+ v[ia[2]][d] -= mm*rij[ll][d];
}
/* 16 flops */
}
{
if (econq == econqCoord)
{
- mm = lagr[ll]/vetavar->rvscale;
+ mm = scaled_lagrange_multiplier[ll]/vetavar->rvscale;
}
if (econq == econqVeloc)
{
- mm = lagr[ll]/(vetavar->vscale*vetavar->vscale_nhc[0]);
+ mm = scaled_lagrange_multiplier[ll]/(vetavar->vscale*vetavar->vscale_nhc[0]);
}
- for (i = 0; i < DIM; i++)
+ for (d = 0; d < DIM; d++)
{
- tmp = mm*rij[ll][i];
- for (j = 0; j < DIM; j++)
+ tmp = mm*rij[ll][d];
+ for (d2 = 0; d2 < DIM; d2++)
{
- vir_r_m_dr[i][j] -= tmp*rij[ll][j];
+ vir_r_m_dr[d][d2] -= tmp*rij[ll][d2];
}
}
/* 21 flops */
}
- /* Correct the lagrange multipliers for the length */
- /* (more details would be useful here . . . )*/
-
- type = ia[0];
+ /* cshake and crattle produce Lagrange multipliers scaled by
+ the reciprocal of the constraint length, so fix that */
if (bFEP)
{
- toler = L1*ip[type].constr.dA + lambda*ip[type].constr.dB;
+ constraint_distance = L1*ip[type].constr.dA + lambda*ip[type].constr.dB;
}
else
{
- toler = ip[type].constr.dA;
- lagr[ll] *= toler;
+ constraint_distance = ip[type].constr.dA;
}
+ scaled_lagrange_multiplier[ll] *= constraint_distance;
}
return nit;
gmx_bool bshakef(FILE *log, gmx_shakedata_t shaked,
real invmass[], int nblocks, int sblock[],
t_idef *idef, t_inputrec *ir, rvec x_s[], rvec prime[],
- t_nrnb *nrnb, real *lagr, real lambda, real *dvdlambda,
+ t_nrnb *nrnb, real *scaled_lagrange_multiplier, real lambda, real *dvdlambda,
real invdt, rvec *v, gmx_bool bCalcVir, tensor vir_r_m_dr,
gmx_bool bDumpOnError, int econq, t_vetavars *vetavar)
{
t_iatom *iatoms;
- real *lam, dt_2, dvdl;
- int i, n0, ncons, blen, type;
+ real dt_2, dvdl;
+ int i, n0, ncon, blen, type, ll;
int tnit = 0, trij = 0;
#ifdef DEBUG
fprintf(log, "nblocks=%d, sblock[0]=%d\n", nblocks, sblock[0]);
#endif
- ncons = idef->il[F_CONSTR].nr/3;
+ ncon = idef->il[F_CONSTR].nr/3;
- for (i = 0; i < ncons; i++)
+ for (ll = 0; ll < ncon; ll++)
{
- lagr[i] = 0;
+ scaled_lagrange_multiplier[ll] = 0;
}
iatoms = &(idef->il[F_CONSTR].iatoms[sblock[0]]);
- lam = lagr;
for (i = 0; (i < nblocks); )
{
blen = (sblock[i+1]-sblock[i]);
blen /= 3;
n0 = vec_shakef(log, shaked, invmass, blen, idef->iparams,
iatoms, ir->shake_tol, x_s, prime, shaked->omega,
- ir->efep != efepNO, lambda, lam, invdt, v, bCalcVir, vir_r_m_dr,
+ ir->efep != efepNO, lambda, scaled_lagrange_multiplier, invdt, v, bCalcVir, vir_r_m_dr,
econq, vetavar);
#ifdef DEBUGSHAKE
}
return FALSE;
}
- tnit += n0*blen;
- trij += blen;
- iatoms += 3*blen; /* Increment pointer! */
- lam += blen;
+ tnit += n0*blen;
+ trij += blen;
+ iatoms += 3*blen; /* Increment pointer! */
+ scaled_lagrange_multiplier += blen;
i++;
}
/* only for position part? */
if (ir->efep != efepNO)
{
real bondA, bondB;
+ /* TODO This should probably use invdt, so that sd integrator scaling works properly */
dt_2 = 1/sqr(ir->delta_t);
dvdl = 0;
- for (i = 0; i < ncons; i++)
+ for (ll = 0; ll < ncon; ll++)
{
- type = idef->il[F_CONSTR].iatoms[3*i];
-
- /* dh/dl contribution from constraint force is dh/dr (constraint force) dot dr/dl */
- /* constraint force is -\sum_i lagr_i* d(constraint)/dr, with constrant = r^2-d^2 */
- /* constraint force is -\sum_i lagr_i* 2 r */
- /* so dh/dl = -\sum_i lagr_i* 2 r * dr/dl */
- /* However, by comparison with lincs and with
- comparison with a full thermodynamics cycle (see
- redmine issue #1255), this is off by a factor of
- two -- the 2r should apparently just be r. Further
- investigation should be done at some point to
- understand why and see if there is something deeper
- we are missing */
+ type = idef->il[F_CONSTR].iatoms[3*ll];
+ /* Per equations in the manual, dv/dl = -2 \sum_ll lagrangian_ll * r_ll * (d_B - d_A) */
+ /* The vector scaled_lagrange_multiplier[ll] contains the value -2 r_ll eta_ll (eta_ll is the
+ estimate of the Langrangian, definition on page 336 of Ryckaert et al 1977),
+ so the pre-factors are already present. */
bondA = idef->iparams[type].constr.dA;
bondB = idef->iparams[type].constr.dB;
- dvdl += lagr[i] * dt_2 * ((1.0-lambda)*bondA + lambda*bondB) * (bondB-bondA);
+ dvdl += scaled_lagrange_multiplier[ll] * dt_2 * (bondB - bondA);
}
*dvdlambda += dvdl;
}
}
void crattle(atom_id iatom[], int ncon, int *nnit, int maxnit,
- real dist2[], real vp[], real rij[], real m2[], real omega,
- real invmass[], real tt[], real lagr[], int *nerror, real invdt, t_vetavars *vetavar)
+ real constraint_distance_squared[], real vp[], real rij[], real m2[], real omega,
+ real invmass[], real distance_squared_tolerance[], real scaled_lagrange_multiplier[],
+ int *nerror, real invdt, t_vetavars *vetavar)
{
/*
* r.c. van schaik and w.f. van gunsteren
int ll, i, j, i3, j3, l3, ii;
int ix, iy, iz, jx, jy, jz;
- real toler, rijd, vpijd, vx, vy, vz, diff, acor, xdotd, fac, im, jm, imdt, jmdt;
+ real constraint_distance_squared_ll;
+ real rijd, vpijd, vx, vy, vz, diff, acor, xdotd, fac, im, jm, imdt, jmdt;
real xh, yh, zh, rijx, rijy, rijz;
real tix, tiy, tiz;
real tjx, tjy, tjz;
vz = vp[iz]-vp[jz];
vpijd = vx*rijx+vy*rijy+vz*rijz;
- toler = dist2[ll];
+ constraint_distance_squared_ll = constraint_distance_squared[ll];
/* this is r(t+dt) \dotproduct \dot{r}(t+dt) */
- xdotd = vpijd*vscale_nhc + veta*toler;
+ xdotd = vpijd*vscale_nhc + veta*constraint_distance_squared_ll;
/* iconv is zero when the error is smaller than a bound */
- iconvf = fabs(xdotd)*(tt[ll]/invdt);
+ iconvf = fabs(xdotd)*(distance_squared_tolerance[ll]/invdt);
if (iconvf > 1)
{
nconv = iconvf;
- fac = omega*2.0*m2[ll]/toler;
+ fac = omega*2.0*m2[ll]/constraint_distance_squared_ll;
acor = -fac*xdotd;
- lagr[ll] += acor;
+ scaled_lagrange_multiplier[ll] += acor;
xh = rijx*acor;
yh = rijy*acor;
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2014, 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.
+ */
+#include "gmxpre.h"
+
+#include <assert.h>
+
+#include <vector>
+
+#include <gtest/gtest.h>
+
+#include "gromacs/legacyheaders/constr.h"
+#include "gromacs/legacyheaders/types/simple.h"
+
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+
+namespace
+{
+
+/*! \brief Stride of the vector of atom_id used to describe each SHAKE
+ * constraint
+ *
+ * Like other such code, SHAKE is hard-wired to use t_ilist.iatoms as
+ * a flat vector of tuples of general data. Here, they are triples
+ * containing the index of the constraint type, and then the indices
+ * of the two atoms involved. So for each constraint, we must stride
+ * this vector by three to get access to its information. */
+const int constraintStride = 3;
+
+/*! \brief Compute the displacements between pairs of constrained
+ * atoms described in the iatom "topology". */
+std::vector<real>
+computeDisplacements(const std::vector<atom_id> &iatom,
+ const std::vector<real> &positions)
+{
+ assert(0 == iatom.size() % constraintStride);
+ int numConstraints = iatom.size() / constraintStride;
+ std::vector<real> displacements;
+
+ for (int ll = 0; ll != numConstraints; ++ll)
+ {
+ int atom_i = iatom[ll*constraintStride + 1];
+ int atom_j = iatom[ll*constraintStride + 2];
+
+ for (int d = 0; d != DIM; d++)
+ {
+ displacements.push_back(positions[atom_i*DIM + d] - positions[atom_j*DIM + d]);
+ }
+ }
+
+ return displacements;
+}
+
+/*! \brief Compute half of the reduced mass of each pair of constrained
+ * atoms in the iatom "topology".
+ *
+ * The reduced mass is m = 1/(1/m_i + 1/m_j)) */
+std::vector<real>
+computeHalfOfReducedMasses(const std::vector<atom_id> &iatom,
+ const std::vector<real> &inverseMasses)
+{
+ int numConstraints = iatom.size() / constraintStride;
+ std::vector<real> halfOfReducedMasses;
+
+ for (int ll = 0; ll != numConstraints; ++ll)
+ {
+ int atom_i = iatom[ll*constraintStride + 1];
+ int atom_j = iatom[ll*constraintStride + 2];
+
+ halfOfReducedMasses.push_back(0.5 / (inverseMasses[atom_i] + inverseMasses[atom_j]));
+ }
+
+ return halfOfReducedMasses;
+}
+
+/*! \brief Compute the distances corresponding to the vector of displacements components */
+std::vector<real>
+computeDistancesSquared(const std::vector<real> &displacements)
+{
+ assert(0 == displacements.size() % DIM);
+ int numDistancesSquared = displacements.size() / DIM;
+ std::vector<real> distanceSquared;
+
+ for (int i = 0; i != numDistancesSquared; ++i)
+ {
+ distanceSquared.push_back(0.0);
+ for (int d = 0; d != DIM; ++d)
+ {
+ real displacement = displacements[i*DIM + d];
+ distanceSquared.back() += displacement * displacement;
+ }
+ }
+
+ return distanceSquared;
+}
+
+/*! \brief Test fixture for testing SHAKE */
+class ShakeTest : public ::testing::Test
+{
+ public:
+ /*! \brief Set up data for test cases to use when constructing
+ their input */
+ void SetUp()
+ {
+ inverseMassesDatabase_.push_back(2.0);
+ inverseMassesDatabase_.push_back(3.0);
+ inverseMassesDatabase_.push_back(4.0);
+ inverseMassesDatabase_.push_back(1.0);
+
+ positionsDatabase_.push_back(2.5);
+ positionsDatabase_.push_back(-3.1);
+ positionsDatabase_.push_back(15.7);
+
+ positionsDatabase_.push_back(0.51);
+ positionsDatabase_.push_back(-3.02);
+ positionsDatabase_.push_back(15.55);
+
+ positionsDatabase_.push_back(-0.5);
+ positionsDatabase_.push_back(-3.0);
+ positionsDatabase_.push_back(15.2);
+
+ positionsDatabase_.push_back(-1.51);
+ positionsDatabase_.push_back(-2.95);
+ positionsDatabase_.push_back(15.05);
+ }
+
+ //! Run the test
+ void runTest(size_t gmx_unused numAtoms,
+ size_t numConstraints,
+ const std::vector<atom_id> &iatom,
+ const std::vector<real> &constrainedDistances,
+ const std::vector<real> &inverseMasses,
+ const std::vector<real> &positions)
+ {
+ // Check the test input is consistent
+ assert(numConstraints * constraintStride == iatom.size());
+ assert(numConstraints == constrainedDistances.size());
+ assert(numAtoms == inverseMasses.size());
+ assert(numAtoms * DIM == positions.size());
+ for (size_t i = 0; i != numConstraints; ++i)
+ {
+ for (size_t j = 1; j < 3; j++)
+ {
+ // Check that the topology refers to atoms that have masses and positions
+ assert(iatom[i*constraintStride + j] >= 0);
+ assert(iatom[i*constraintStride + j] < static_cast<atom_id>(numAtoms));
+ }
+ }
+ std::vector<real> distanceSquaredTolerances;
+ std::vector<real> lagrangianValues;
+ std::vector<real> constrainedDistancesSquared;
+
+ for (size_t i = 0; i != numConstraints; ++i)
+ {
+ constrainedDistancesSquared.push_back(constrainedDistances[i] * constrainedDistances[i]);
+ distanceSquaredTolerances.push_back(1.0 / (constrainedDistancesSquared.back() * ShakeTest::tolerance_));
+ lagrangianValues.push_back(0.0);
+ }
+ std::vector<real> halfOfReducedMasses = computeHalfOfReducedMasses(iatom, inverseMasses);
+ std::vector<real> initialDisplacements = computeDisplacements(iatom, positions);
+
+ std::vector<real> finalPositions = positions;
+ int numIterations = 0;
+ int numErrors = 0;
+
+ cshake(&iatom[0], numConstraints, &numIterations,
+ ShakeTest::maxNumIterations_, &constrainedDistancesSquared[0],
+ &finalPositions[0], &initialDisplacements[0], &halfOfReducedMasses[0],
+ omega_, &inverseMasses[0],
+ &distanceSquaredTolerances[0],
+ &lagrangianValues[0],
+ &numErrors);
+
+ std::vector<real> finalDisplacements = computeDisplacements(iatom, finalPositions);
+ std::vector<real> finalDistancesSquared = computeDistancesSquared(finalDisplacements);
+ assert(numConstraints == finalDistancesSquared.size());
+
+ EXPECT_EQ(0, numErrors);
+ EXPECT_GT(numIterations, 1);
+ EXPECT_LT(numIterations, ShakeTest::maxNumIterations_);
+ // TODO wrap this in a Google Mock matcher if there's
+ // other tests like it some time?
+ for (size_t i = 0; i != numConstraints; ++i)
+ {
+ gmx::test::FloatingPointTolerance constraintTolerance =
+ gmx::test::relativeToleranceAsFloatingPoint(constrainedDistancesSquared[i],
+ ShakeTest::tolerance_);
+ // Assert that the constrained distances are within the required tolerance
+ EXPECT_FLOAT_EQ_TOL(constrainedDistancesSquared[i],
+ finalDistancesSquared[i],
+ constraintTolerance);
+ }
+ }
+
+ //! Tolerance for SHAKE conversion (ie. shake-tol .mdp setting)
+ static const real tolerance_;
+ //! Maximum number of iterations permitted in these tests
+ static const int maxNumIterations_;
+ //! SHAKE over-relaxation (SOR) factor
+ static const real omega_;
+ //! Database of inverse masses of atoms in the topology
+ std::vector<real> inverseMassesDatabase_;
+ //! Database of atom positions (three reals per atom)
+ std::vector<real> positionsDatabase_;
+};
+
+const real ShakeTest::tolerance_ = 1e-5;
+const int ShakeTest::maxNumIterations_ = 30;
+const real ShakeTest::omega_ = 1.0;
+
+TEST_F(ShakeTest, ConstrainsOneBond)
+{
+ int numAtoms = 2;
+ int numConstraints = 1;
+
+ std::vector<atom_id> iatom;
+ iatom.push_back(-1); // unused
+ iatom.push_back(0); // i atom index
+ iatom.push_back(1); // j atom index
+
+ std::vector<real> constrainedDistances;
+ constrainedDistances.push_back(2.0);
+
+ std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
+ inverseMassesDatabase_.begin() + numAtoms);
+ std::vector<real> positions(positionsDatabase_.begin(),
+ positionsDatabase_.begin() + numAtoms * DIM);
+
+ runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
+}
+
+TEST_F(ShakeTest, ConstrainsTwoDisjointBonds)
+{
+ int numAtoms = 4;
+ int numConstraints = 2;
+
+ std::vector<atom_id> iatom;
+ iatom.push_back(-1); // unused
+ iatom.push_back(0); // i atom index
+ iatom.push_back(1); // j atom index
+
+ iatom.push_back(-1); // unused
+ iatom.push_back(2); // i atom index
+ iatom.push_back(3); // j atom index
+
+ std::vector<real> constrainedDistances;
+ constrainedDistances.push_back(2.0);
+ constrainedDistances.push_back(1.0);
+
+ std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
+ inverseMassesDatabase_.begin() + numAtoms);
+ std::vector<real> positions(positionsDatabase_.begin(),
+ positionsDatabase_.begin() + numAtoms * DIM);
+
+ runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
+}
+
+TEST_F(ShakeTest, ConstrainsTwoBondsWithACommonAtom)
+{
+ int numAtoms = 3;
+ int numConstraints = 2;
+
+ std::vector<atom_id> iatom;
+ iatom.push_back(-1); // unused
+ iatom.push_back(0); // i atom index
+ iatom.push_back(1); // j atom index
+
+ iatom.push_back(-1); // unused
+ iatom.push_back(1); // i atom index
+ iatom.push_back(2); // j atom index
+
+ std::vector<real> constrainedDistances;
+ constrainedDistances.push_back(2.0);
+ constrainedDistances.push_back(1.0);
+
+ std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
+ inverseMassesDatabase_.begin() + numAtoms);
+ std::vector<real> positions(positionsDatabase_.begin(),
+ positionsDatabase_.begin() + numAtoms * DIM);
+
+ runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
+}
+
+TEST_F(ShakeTest, ConstrainsThreeBondsWithCommonAtoms)
+{
+ int numAtoms = 4;
+ int numConstraints = 3;
+
+ std::vector<atom_id> iatom;
+ iatom.push_back(-1); // unused
+ iatom.push_back(0); // i atom index
+ iatom.push_back(1); // j atom index
+
+ iatom.push_back(-1); // unused
+ iatom.push_back(1); // i atom index
+ iatom.push_back(2); // j atom index
+
+ iatom.push_back(-1); // unused
+ iatom.push_back(2); // i atom index
+ iatom.push_back(3); // j atom index
+
+ std::vector<real> constrainedDistances;
+ constrainedDistances.push_back(2.0);
+ constrainedDistances.push_back(1.0);
+ constrainedDistances.push_back(1.0);
+
+ std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
+ inverseMassesDatabase_.begin() + numAtoms);
+ std::vector<real> positions(positionsDatabase_.begin(),
+ positionsDatabase_.begin() + numAtoms * DIM);
+
+ runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
+}
+
+} // namespace
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff