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-2w0 /
/ for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4w0 */
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