GROMACS

The parade of mdrun fixes for 4.6.2 continues :-)

But everybody loves a parade!

Is this really a "normal" priority bug? Should it not be at least high priority and perhaps announced so people know about it and can redo their simulations if necessary?

Just to clarify about whether the approximation for md-vv (doubling the position restraint correction term) is correct,
here's the test with changing the dipole length from 0.5 to 0.2 nm with md, md-vv, and sd, 10 ns simulatiom. Note there is one sd simulation (lambda state = 1) that is running really slowly for some pinning issue, affecting the overall uncertainty-- I will fix it when it finishes. I don't think it affects the conclusions.

                   md:                  md-vv:                sd:
   0 -- 1         1.760  +-  0.018     1.752  +-  0.018      1.680  +-  0.078
   1 -- 2         1.831  +-  0.018     1.805  +-  0.018      1.764  +-  0.067
   2 -- 3         2.307  +-  0.019     2.267  +-  0.019      2.315  +-  0.024
   3 -- 4         2.887  +-  0.019     2.831  +-  0.019      2.887  +-  0.022
   4 -- 5         3.412  +-  0.019     3.386  +-  0.019      3.415  +-  0.021
   5 -- 6         4.225  +-  0.018     4.235  +-  0.018      4.226  +-  0.021
   6 -- 7         5.246  +-  0.017     5.239  +-  0.018      5.268  +-  0.020
   7 -- 8         6.361  +-  0.017     6.341  +-  0.017      6.335  +-  0.020
   8 -- 9         7.340  +-  0.016     7.288  +-  0.016      7.287  +-  0.019
   9 -- 10        8.085  +-  0.016     8.032  +-  0.019      8.044  +-  0.018
---------- ---------------- -----------------  ------------------
 TOTAL:          43.455  +-  0.069    43.178  +-  0.070     43.221  +-  0.154

(can't figure out how to fix spacing)
So, they appear to be almost statistically indistinguishable. The approximation for md-vv (if it is an approximation) gives results that appear to be well within noise for sd, and close to noise to md-vv (2.6 sigma). This should be good enough for essentially all situations.

Tabs are evil. Fixed :-)

While converting files to C++, I unearthed the cause of the bug that led to this fix.

Michael inserted some braces wrongly in http://redmine.gromacs.org/projects/gromacs/repository/revisions/a6ee084bcf3ac87d9dd493702b1c737d86506bdc/diff/src/mdlib/shakef.c that broke some FEP functionality at the end of vec_shake, which led to him later re-implementing that in bshakef (see http://redmine.gromacs.org/projects/gromacs/repository/revisions/714b43e235bf3ca67ebbe73bf5d267e1bf352b7b/diff/src/mdlib/shakef.c)

Berk's original implementation re-used the toler variable to post-multiply the Lagrangian multipliers by the constraint distance (weighted by FE-lambda in that case), but those multipliers are unused in the non-FE case. I don't think there is any use for that other than computing dH/dl, but we should clear that up rather than returning a Lagrangian that is scaled by the constraint distance. Michael's subsequent changes have converged on an implementation equivalent to Berk's original one, but with clarity in what value lagr[] actually holds.

Michael also observed that there is a factor of two that seems to be missing in the SHAKE case. The original Ryckaert SHAKE paper derives the constraints based on (r_i-r_j)^2-d^2=0, but the LINCS paper uses |r_i-r_j|-d=0. GROMACS documentation is inconsistent, with neither manual section 3.6.1 and 3.6.2 directly agreeing with section 4.5. The implementation seems to use the LINCS definition in the SHAKE case for computing dH/dl, which offhand seems OK so long as it is consistently implemented (and documented). So, is it? And how should we fix the docs?

Gerrit received a related patchset '1' for Issue #1255.
Uploader: Mark Abraham (mark.j.abraham@gmail.com)
Change-Id: I5b5eee0b8f3f4761ce9cb681dbdb0d4526a6761d
Gerrit URL: https://gerrit.gromacs.org/4057

Mark Abraham wrote:

While converting files to C++, I unearthed the cause of the bug that led to this fix.

Michael inserted some braces wrongly in http://redmine.gromacs.org/projects/gromacs/repository/revisions/a6ee084bcf3ac87d9dd493702b1c737d86506bdc/diff/src/mdlib/shakef.c that broke some FEP functionality at the end of vec_shake, which led to him later re-implementing that in bshakef (see http://redmine.gromacs.org/projects/gromacs/repository/revisions/714b43e235bf3ca67ebbe73bf5d267e1bf352b7b/diff/src/mdlib/shakef.c)

Berk's original implementation re-used the toler variable to post-multiply the Lagrangian multipliers by the constraint distance (weighted by FE-lambda in that case), but those multipliers are unused in the non-FE case. I don't think there is any use for that other than computing dH/dl, but we should clear that up rather than returning a Lagrangian that is scaled by the constraint distance. Michael's subsequent changes have converged on an implementation equivalent to Berk's original one, but with clarity in what value lagr[] actually holds.

Clarified - the cshake/crattle routines return the Lagrangian multipliers divided by their respective constraint distances, so the post-multiplication was simply finishing the job of computing the Lagrangian multipliers.

Michael also observed that there is a factor of two that seems to be missing in the SHAKE case. The original Ryckaert SHAKE paper derives the constraints based on (r_i-r_j)^2-d^2=0, but the LINCS paper uses |r_i-r_j|-d=0. GROMACS documentation is inconsistent, with neither manual section 3.6.1 and 3.6.2 directly agreeing with section 4.5. The implementation seems to use the LINCS definition in the SHAKE case for computing dH/dl, which offhand seems OK so long as it is consistently implemented (and documented). So, is it? And how should we fix the docs?

Fixed docs, and made derivation in 4.5 make more sense, which clarifies why the code was correct all along. Note the distinction in the case of SHAKE between the quadratic constraint equation, and the linear lambda-interpolation.