Size independent Hessian for normal mode analysis
In order to perform a normal mode analysis GROMACS computes a numerical derivative of the force using
F'(x) = ( F(x+delta)-F(x-delta) ) / ( 2 delta )
The delta up to gromacs 2018 was 10 sqrt(GMX_REAL_EPS) which works for classical models but does not give reproducible results for polarizable models. There is a connection to shell minimization that has a finite precision as well.
Making delta smaller may yield problems for large structures where the delta may be zero in practice.
The general solution for NMA would be to compute a delta that is (just) big enough such that the delta is not zero in any case.
For that it would have to depend on the system size, e.g.
delta = max(box[XX], max(box[YY], box[ZZ])
Added new testing code for normal mode analysis.
Part of #2771
Removed normal modes testing.
Normal mode analysis testing is now part of the gromacs
core code and is not needed here anymore.
Part of #2771
Implemented native normal mode analysis tests.
Moved the normal mode analysis tests from the regression
test set to core gromacs.
In order to do so the tests were converted to use the
Verlet scheme which leads to minor changes in the eigenvalues
produced by the normal modes analysis.
Part of #2771
#1 Updated by Erik Lindahl about 2 years ago
This problem is actually quite a bit harder than it might appear at first sight, and has to do with the fundamental problem of calculating numerical derivatives.
In the previous (and current...) versions of GROMACS, we mostly just fudge it and used a 2-point derivative approximation, initially with the length equal to sqrt(EPS), but then we increased this do be able to handle coordinates with values up to O(10) at least.
I looked into this a bit when designing the new table classes, where we use higher-order approximations. This is described in more detail e.g. at
Briefly, the challenge is to handle the balance between finite precision biting us when the interval is too small, vs. the formula/approximation error biting us when it is too large - but on top of this it's also easy to get problems with compilers "optimizing" away some of the careful logic and cause errors anyway.
I think it might be a good solution to always pick h=sqrt(eps)*x, where x is the coordinate value, but we'll have to do this carefully to make sure it's both consistent between our different minimizers, and test it both on large and small systems close to and far away from minima.
... which in turn means it might not be a small/trivial fix, so if we want that in release-2019, it has to be ready the next week or two :-)
#3 Updated by David van der Spoel about 2 years ago
In order to clarify the connection to shell minimization: there is a precision issue there that the change in coordinates in shell minimization proceeds using Hooke's law F = k dx. This is implemented using the following:
1) compute F on all shells
2) compute dx = F / k
3) update x by dx.
Since k is typically large, e.g. 400000 in the Charmm polarizable force field, a force that is considered small enough for usage in normal mode analysis (1e-3) would give dx = 2.5e-9 nm, but since F and dx are vectors, the dx[m], m in XX,YY,ZZ, will be smaller and precision problems may occur.
#4 Updated by Erik Lindahl about 2 years ago
David van der Spoel wrote:
h = sqrt(eps) * std::max(fabs(x), 0.001)
No, the spacing for the 2-point approximation of a derivative should never be set only based on the function value, and 0.001*sqrt(eps) will be way too small.
1. To handle the rounding (precision) error, you cannot go below x*sqrt(eps). This problem is not merely due to one addition (then it would rather be eps without the square root), but the properties of the quadratic form of the function close to the minimum, which makes it way more sensitive (as I explained earlier in https://redmine.gromacs.org/issues/2641 ).
2. To handle the formula/secant error (this is where the function value is important), you need to derive the proper formulae (depending on the approximation # points) based on the quotient between the function value and second derivative, but that itself gets very complex close to a local minimum.
It is not possible to solve them both by simply picking a value, and "fixing" it by looking at specific numbers for one case will just mean we cause problems for the other case.
I wasn't lying when I said it isn't trivial - if it was, I would already have fixed the corner cases of the table generation code ;-)
#5 Updated by David van der Spoel about 2 years ago
I wasn't accusing you of lying :)
In fact we may be able to do slightly better than we are in NMA at no cost, since we have the force at the minimum as well. In that manner we could determine the F' from three points, but whether it will lead to practical improvement is another matter.
The more I think about it, the more I think that the real problem for normal modes of shells is not the exact value of h, but what I wrote in comment 4. That is, the force is not reproducible enough because of precision issues minimizing shells. The reason why we are seeing this now is that we have not used the code like this. It would be great not to have to code an analytical Hessian ;)