Reference for Suzuki-Yoshida weights

[thomasloux]

While making the batched version of Nosé Hoover NVT integrator, I discovered the higher order integration using Suzuki Yoshida weights. I was a bit surprised to see float values hard coded, instead of the numerical expression. Despite at many references given in torchsim, jax-md and openMM and gromacs, I did not find any papers reporting the exact same values in these softwares:

SUZUKI_YOSHIDA_WEIGHTS = {
    1: [1],
    3: [0.828981543588751, -0.657963087177502, 0.828981543588751],
    5: [0.2967324292201065, 0.2967324292201065, -0.186929716880426,
        0.2967324292201065, 0.2967324292201065],
    7: [0.784513610477560, 0.235573213359357, -1.17767998417887,
        1.31518632068391, -1.17767998417887, 0.235573213359357,
        0.784513610477560]
}

Only in Gromacs they actually provide the expression:

/* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration  */
/* for n=1, w0 = 1 */
/* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
/* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */

https://github.com/gromacs/gromacs/blob/8bed29775ddcd00d2ac3ae5602f79b74ab079af4/src/gromacs/mdlib/coupling.cpp#L89

What is surprising is that ASE and related papers rather provide these values:

Coefficients for the fourth-order Suzuki-Yoshida integration scheme
Ref: H. Yoshida, Phys. Lett. A 150, 5-7, 262-268 (1990).
https://doi.org/10.1016/0375-9601(90)90092-3

FOURTH_ORDER_COEFFS = [
    1 / (2 - 2 ** (1 / 3)),
    -(2 ** (1 / 3)) / (2 - 2 ** (1 / 3)),
    1 / (2 - 2 ** (1 / 3)),
]

So everything is similar in the expression up to a negative sign in the exponent, which changes everything. In the end, I don't think they are any difference in the implementation (jax-md and torch sim are similar as it was an inspiration, ase is similar although I may be wrong). Is there any reason for this difference? In the paper it's ASE version that is given.

Actually if one closely looks at Suzuki-Yoshida method, it's supposed to provide around terms 3**n (actually less as some terms can be factorised) to get a 2n order integrator. Here it seems to be only 2n-1 terms.

Do you have the derivation of these weights?

References:

• Canonical dynamics: Equilibrium phase-space distributions, Hoover 85
• Glenn J. Martyna , Mark E. Tuckerman , Douglas J. Tobias & Michael L. Klein
  (1996) Explicit reversible integrators for extended systems dynamics, Molecular Physics, 87:5,
  1117-1157, DOI: 10.1080/00268979600100761
• Construction of higher order symplectic integrators, Haruo Yoshida, 90

[CompRhys]

I believe these were copied from jax-md which we're copied from open-mm.

Neither of those references fully assign back to literature. I will send off a message to the person on the OpenMM team who originally implemented their nose-hoover and see if they can help trace this back to an earlier ground truth

[CompRhys]

Andreas Kraemer responded to me and the memory of the original source is lost to time, he says that they're likely sourced from the Tuckerman book although I do not have access to it to check. There are fairly extensive tests of the OpenMM implementation and they checked that they do work with the equations in the Yoshida 1990 paper.

[abhijeetgangan]

@CompRhys is right. I won't say it's lost in time. The derivation and weights can be found in - Statistical Mechanics: Theory and Molecular Simulation Section 4.12 by Tuckerman and they match with what all the packages have (TorchSim, OpenMM etc).

[thomasloux]

Ok thanks for the reference. Looking at it at Section 4.11 (Integrating the Nose-Hoover chain equations), equation (4.11.12) indeed provide the decomposition SUZUKI_YOSHIDA_WEIGHTS[7] in Torch Sim. The book provide a decomposition in 4th order with 3 or 5 terms but they are not the ones used in openMM, jax-md, torchsim or gromacs:
The book uses
for n=3, w0 = w2 = 1 /(2−2^(1/3)), w1 = 1 − 2w0
for n=5, w0 = w1 = w3 = w4 = 1/(4−4^(1/3)), w1 = 1-4w0

While, these constants are used in these MD codes:
for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2w0
for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4w0

By the way, I actually did not notice but as you said the 6th order derivation was provided in Yoshida 1990 paper and is the one used in the MD codes. I was confused as it's written in scientific notation and only write w_1, w_2, w_3 (See Solution A of Table 1).

Also, while I understand that SUZUKI_YOSHIDA_WEIGHTS[7] is a 6th order integrator for L_NH, as far as I understand, SUZUKI_YOSHIDA_WEIGHTS 3 and 5 should both be 4th order integrator. Are they any other difference then? It seems like 4th order with 3 decomposition should always be preferred compared to the 5 terms for efficiency reasons. For instance, in gromacs, they actually force the use of the 5 terms. In ASE, they force the use of 3 terms (with actually the weights provided by Yoshida).

[thomasloux]

@abhijeetgangan let me know what you think about that when you have time!

[abhijeetgangan]

I will post a longer answer after I check again with literature if I am missing something but tldr these do not correspond to Yoshida's weights as per the paper in the sense the 3 stage and 5 stage schemes do not satisfy all the moment conditions. The 3 stage and 5 stage equations (4.3) and (4.11) respectively have unique real solutions and the weights here are not one of them. 7 stage is fine as those are the same weights. It should be safe to switch all weights to the original ones from Yoshida's paper.

RE default scheme: We could use the 3 stage scheme because doing for loops in python might be slow. 1 stage should possibly be avoided. Would need to look into literature or test it ourself the optimal for speed vs accuracy.