From 513813b51406baae3a1b6835daf50017d0f9cb7f Mon Sep 17 00:00:00 2001 From: Mark Abraham Date: Thu, 11 Sep 2014 12:01:14 +0200 Subject: [PATCH] Clean up cshake() and resolve old issues C++ conversion in Id190e36 exposed some "interesting" use of the variable toler caused by an ancient bug that has since been resolved. Investigation showed that variable never meant anything to do with a tolerance! Renamed a bunch of variables for clarity and consistency. Documented the core function. Clarified code comments and manual sections accordingly. Added some unit tests while trying to understand what things were actually going on. Refs #1255 Change-Id: I5b5eee0b8f3f4761ce9cb681dbdb0d4526a6761d --- docs/manual/forcefield.tex | 37 ++- src/gromacs/legacyheaders/constr.h | 6 +- src/gromacs/mdlib/CMakeLists.txt | 3 + src/gromacs/mdlib/clincs.c | 1 + src/gromacs/mdlib/constr.c | 2 +- src/gromacs/mdlib/shakef.c | 300 +++++++++++---------- src/gromacs/mdlib/tests/CMakeLists.txt | 36 +++ src/gromacs/mdlib/tests/shake.cpp | 347 +++++++++++++++++++++++++ 8 files changed, 582 insertions(+), 150 deletions(-) create mode 100644 src/gromacs/mdlib/tests/CMakeLists.txt create mode 100644 src/gromacs/mdlib/tests/shake.cpp diff --git a/docs/manual/forcefield.tex b/docs/manual/forcefield.tex index 70829ad13e..c406e054a6 100644 --- a/docs/manual/forcefield.tex +++ b/docs/manual/forcefield.tex @@ -1810,27 +1810,42 @@ after taking the derivative, we {\em can} insert \ve{p} = m\ve{v}, such that: \subsubsection{Constraints} \label{subsubsec:constraints} -\newcommand{\clam}{C_{\lambda}} The constraints are formally part of the Hamiltonian, and therefore they give a contribution to the free energy. In {\gromacs} this can be calculated using the \normindex{LINCS} or the \normindex{SHAKE} algorithm. -If we have a number of constraint equations $g_k$: +If we have $k = 1 \ldots K$ constraint equations $g_k$ for LINCS, then \beq -g_k = \ve{r}_{k} - d_{k} +g_k = |\ve{r}_{k}| - d_{k} \eeq -where $\ve{r}_k$ is the distance vector between two particles and -$d_k$ is the constraint distance between the two particles, we can write -this using a $\LAM$-dependent distance as +where $\ve{r}_k$ is the displacement vector between two particles and +$d_k$ is the constraint distance between the two particles. We can express +the fact that the constraint distance has a $\LAM$ dependency by \beq -g_k = \ve{r}_{k} - \left(\LL d_{k}^A + \LAM d_k^B\right) +d_k = \LL d_{k}^A + \LAM d_k^B \eeq -the contribution $\clam$ -to the Hamiltonian using Lagrange multipliers $\lambda$: + +Thus the $\LAM$-dependent constraint equation is +\beq +g_k = |\ve{r}_{k}| - \left(\LL d_{k}^A + \LAM d_k^B\right). +\eeq + +The (zero) contribution $G$ to the Hamiltonian from the constraints +(using Lagrange multipliers $\lambda_k$, which are logically distinct +from the free-energy $\LAM$) is \bea -\clam &=& \sum_k \lambda_k g_k \\ -\dvdl{\clam} &=& \sum_k \lambda_k \left(d_k^B-d_k^A\right) +G &=& \sum^K_k \lambda_k g_k \\ +\dvdl{G} &=& \frac{\partial G}{\partial d_k} \dvdl{d_k} \\ + &=& - \sum^K_k \lambda_k \left(d_k^B-d_k^A\right) \eea +For SHAKE, the constraint equations are +\beq +g_k = \ve{r}_{k}^2 - d_{k}^2 +\eeq +with $d_k$ as before, so +\bea +\dvdl{G} &=& -2 \sum^K_k \lambda_k \left(d_k^B-d_k^A\right) +\eea \subsection{Soft-core interactions\index{soft-core interactions}} \begin{figure} diff --git a/src/gromacs/legacyheaders/constr.h b/src/gromacs/legacyheaders/constr.h index 9d01ba8ccc..79ac3a7bf0 100644 --- a/src/gromacs/legacyheaders/constr.h +++ b/src/gromacs/legacyheaders/constr.h @@ -125,9 +125,9 @@ void settle_proj(gmx_settledata_t settled, int econq, * of coordinates working on settle type constraint. */ -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); +void cshake(const atom_id iatom[], int ncon, int *nnit, int maxnit, + const real dist2[], real xp[], const real rij[], const real m2[], real omega, + const real invmass[], const real tt[], real lagr[], int *nerror); /* Regular iterative shake */ void crattle(atom_id iatom[], int ncon, int *nnit, int maxnit, diff --git a/src/gromacs/mdlib/CMakeLists.txt b/src/gromacs/mdlib/CMakeLists.txt index fb9cfad823..d068618f83 100644 --- a/src/gromacs/mdlib/CMakeLists.txt +++ b/src/gromacs/mdlib/CMakeLists.txt @@ -39,3 +39,6 @@ if(GMX_GPU) endif() set(MDLIB_SOURCES ${MDLIB_SOURCES} PARENT_SCOPE) +if (BUILD_TESTING) + add_subdirectory(tests) +endif() diff --git a/src/gromacs/mdlib/clincs.c b/src/gromacs/mdlib/clincs.c index 3423b0811a..a239df22bb 100644 --- a/src/gromacs/mdlib/clincs.c +++ b/src/gromacs/mdlib/clincs.c @@ -1573,6 +1573,7 @@ gmx_bool constrain_lincs(FILE *fplog, gmx_bool bLog, gmx_bool bEner, { real dt_2, dvdl = 0; + /* TODO This should probably use invdt, so that sd integrator scaling works properly */ dt_2 = 1.0/(ir->delta_t*ir->delta_t); for (i = 0; (i < lincsd->nc); i++) { diff --git a/src/gromacs/mdlib/constr.c b/src/gromacs/mdlib/constr.c index fb075d4648..07866e1901 100644 --- a/src/gromacs/mdlib/constr.c +++ b/src/gromacs/mdlib/constr.c @@ -78,7 +78,7 @@ typedef struct gmx_constr { int nblocks; /* The number of SHAKE blocks */ int *sblock; /* The SHAKE blocks */ int sblock_nalloc; /* The allocation size of sblock */ - real *lagr; /* Lagrange multipliers for SHAKE */ + real *lagr; /* -2 times the Lagrange multipliers for SHAKE */ int lagr_nalloc; /* The allocation size of lagr */ int maxwarn; /* The maximum number of warnings */ int warncount_lincs; diff --git a/src/gromacs/mdlib/shakef.c b/src/gromacs/mdlib/shakef.c index 46030d0cbc..6631e59dc5 100644 --- a/src/gromacs/mdlib/shakef.c +++ b/src/gromacs/mdlib/shakef.c @@ -48,9 +48,9 @@ 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; @@ -64,11 +64,11 @@ gmx_shakedata_t shake_init() 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; @@ -91,27 +91,56 @@ static void pv(FILE *log, char *s, rvec x) 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; @@ -121,9 +150,9 @@ void cshake(atom_id iatom[], int ncon, int *nnit, int maxnit, 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; @@ -135,42 +164,48 @@ void cshake(atom_id iatom[], int ncon, int *nnit, int maxnit, 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; } } } @@ -183,36 +218,36 @@ int vec_shakef(FILE *fplog, gmx_shakedata_t shaked, 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) { @@ -220,30 +255,30 @@ int vec_shakef(FILE *fplog, gmx_shakedata_t shaked, 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; } @@ -270,25 +305,28 @@ int vec_shakef(FILE *fplog, gmx_shakedata_t shaked, 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 */ } @@ -298,36 +336,34 @@ int vec_shakef(FILE *fplog, gmx_shakedata_t shaked, { 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; @@ -378,35 +414,34 @@ static void check_cons(FILE *log, int nc, rvec x[], rvec prime[], rvec v[], 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 @@ -423,10 +458,10 @@ gmx_bool bshakef(FILE *log, gmx_shakedata_t shaked, } 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? */ @@ -435,27 +470,20 @@ gmx_bool bshakef(FILE *log, gmx_shakedata_t shaked, 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; } @@ -487,8 +515,9 @@ gmx_bool bshakef(FILE *log, gmx_shakedata_t shaked, } 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 @@ -503,7 +532,8 @@ void crattle(atom_id iatom[], int ncon, int *nnit, int maxnit, 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; @@ -539,19 +569,19 @@ void crattle(atom_id iatom[], int ncon, int *nnit, int maxnit, 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; diff --git a/src/gromacs/mdlib/tests/CMakeLists.txt b/src/gromacs/mdlib/tests/CMakeLists.txt new file mode 100644 index 0000000000..73e80a0bd3 --- /dev/null +++ b/src/gromacs/mdlib/tests/CMakeLists.txt @@ -0,0 +1,36 @@ +# +# 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. + +gmx_add_unit_test(ShakeUnitTests shake-test + shake.cpp) diff --git a/src/gromacs/mdlib/tests/shake.cpp b/src/gromacs/mdlib/tests/shake.cpp new file mode 100644 index 0000000000..4ab751fc30 --- /dev/null +++ b/src/gromacs/mdlib/tests/shake.cpp @@ -0,0 +1,347 @@ +/* + * 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 + +#include + +#include + +#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 +computeDisplacements(const std::vector &iatom, + const std::vector &positions) +{ + assert(0 == iatom.size() % constraintStride); + int numConstraints = iatom.size() / constraintStride; + std::vector 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 +computeHalfOfReducedMasses(const std::vector &iatom, + const std::vector &inverseMasses) +{ + int numConstraints = iatom.size() / constraintStride; + std::vector 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 +computeDistancesSquared(const std::vector &displacements) +{ + assert(0 == displacements.size() % DIM); + int numDistancesSquared = displacements.size() / DIM; + std::vector 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 &iatom, + const std::vector &constrainedDistances, + const std::vector &inverseMasses, + const std::vector &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(numAtoms)); + } + } + std::vector distanceSquaredTolerances; + std::vector lagrangianValues; + std::vector 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 halfOfReducedMasses = computeHalfOfReducedMasses(iatom, inverseMasses); + std::vector initialDisplacements = computeDisplacements(iatom, positions); + + std::vector 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 finalDisplacements = computeDisplacements(iatom, finalPositions); + std::vector 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 inverseMassesDatabase_; + //! Database of atom positions (three reals per atom) + std::vector 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 iatom; + iatom.push_back(-1); // unused + iatom.push_back(0); // i atom index + iatom.push_back(1); // j atom index + + std::vector constrainedDistances; + constrainedDistances.push_back(2.0); + + std::vector inverseMasses(inverseMassesDatabase_.begin(), + inverseMassesDatabase_.begin() + numAtoms); + std::vector 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 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 constrainedDistances; + constrainedDistances.push_back(2.0); + constrainedDistances.push_back(1.0); + + std::vector inverseMasses(inverseMassesDatabase_.begin(), + inverseMassesDatabase_.begin() + numAtoms); + std::vector 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 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 constrainedDistances; + constrainedDistances.push_back(2.0); + constrainedDistances.push_back(1.0); + + std::vector inverseMasses(inverseMassesDatabase_.begin(), + inverseMassesDatabase_.begin() + numAtoms); + std::vector 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 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 constrainedDistances; + constrainedDistances.push_back(2.0); + constrainedDistances.push_back(1.0); + constrainedDistances.push_back(1.0); + + std::vector inverseMasses(inverseMassesDatabase_.begin(), + inverseMassesDatabase_.begin() + numAtoms); + std::vector positions(positionsDatabase_.begin(), + positionsDatabase_.begin() + numAtoms * DIM); + + runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions); +} + +} // namespace -- 2.22.0