## Task #2771

### Size independent Hessian for normal mode analysis

**Description**

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])

### Associated revisions

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

Change-Id: I82df371feeabc40a0951bf21afa09a50c08e045f

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

Change-Id: Ice2e4423d8f48f1656aaa296e9207fda3bc4d233

### History

#### #1 Updated by Erik Lindahl almost 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

https://en.wikipedia.org/wiki/Numerical_differentiation

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 :-)

#### #2 Updated by David van der Spoel almost 2 years ago

Thanks for looking into this.

I guess that would have to combined with a minimum coordinate value of e.g. 1 pm as well, something like:

h = sqrt(eps) * std::max(fabs(x), 0.001)

I could add a couple of tests to begin with and then implement this.

#### #3 Updated by David van der Spoel almost 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 almost 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.

In summary:

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 almost 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 ;)

#### #6 Updated by Gerrit Code Review Bot almost 2 years ago

Gerrit received a related patchset '3' for Issue #2771.

Uploader: David van der Spoel (spoel@xray.bmc.uu.se)

Change-Id: regressiontests~master~I046a5069f24ef65077690838ff6f5ec9e2a88573

Gerrit URL: https://gerrit.gromacs.org/8743

#### #7 Updated by Gerrit Code Review Bot almost 2 years ago

Gerrit received a related patchset '1' for Issue #2771.

Uploader: David van der Spoel (spoel@xray.bmc.uu.se)

Change-Id: gromacs~master~Ifb44841b36be5076fe8d77d5982876f881fdc69b

Gerrit URL: https://gerrit.gromacs.org/8747

#### #8 Updated by David van der Spoel almost 2 years ago

Once the first testing patch is in, all the normal mode analysis tests will be moved from the regression tests to the core code and removed from the regression tests.

#### #9 Updated by Paul Bauer almost 2 years ago

**Tracker**changed from*Bug*to*Task***Target version**changed from*2019*to*2020***Affected version**deleted ()*2018*

I changed this to 2020, because I don't think this will happen in time for 2019 and would also mean some substantial changes that we shouldn't risk.

#### #10 Updated by Gerrit Code Review Bot almost 2 years ago

Gerrit received a related patchset '1' for Issue #2771.

Uploader: David van der Spoel (spoel@xray.bmc.uu.se)

Change-Id: gromacs~master~Ice2e4423d8f48f1656aaa296e9207fda3bc4d233

Gerrit URL: https://gerrit.gromacs.org/8795

#### #11 Updated by Gerrit Code Review Bot almost 2 years ago

Gerrit received a related patchset '1' for Issue #2771.

Uploader: David van der Spoel (spoel@xray.bmc.uu.se)

Change-Id: regressiontests~master~I82df371feeabc40a0951bf21afa09a50c08e045f

Gerrit URL: https://gerrit.gromacs.org/8796

#### #12 Updated by Paul Bauer 9 months ago

**Target version**changed from*2020*to*future*

Added new testing code for normal mode analysis.

Part of #2771

Change-Id: Ifb44841b36be5076fe8d77d5982876f881fdc69b