1549380323-Statistical Mechanics Theory and Molecular Simulation

Integrating Nos ́e–Hoover chains 199

algorithm could improve the accuracy of RESPA without adding significantly to the
computational overhead. Fortunately, high order methods suitable for our purposes
exist and are straightforward to apply. One scheme in particular, due to Suzuki (1991a,
1991 b) and Yoshida (1990), has proved particularly useful for the Nos ́e–Hoover chain
system.
The Suzuki–Yoshida scheme works as follows: LetS(λ) be aprimitive factorization
of the operator exp[λ(A 1 +A 2 )]. For example, a primitive factorization could be the
simple Trotter schemeS(λ) = exp(λA 2 /2) exp(λA 1 ) exp(λA 2 /2). Next, introduce a
set ofnsyweightswαsuch that
∑nsy

α=1

wα= 1. (4.11.9)

These weights are chosen in such a way that error terms up to a certain order 2sare
eliminated in a general factorization of exp[λ(A 1 +A 2 )], yielding a high order scheme.
In the original Suzuki scheme, it was shown that

nsy= 5s−^1 , (4.11.10)

so that a fourth-order scheme would require 5 weights, a sixth-order scheme would
require 25 weights, etc., with all weights having a simple analytical form. For example,
for 2s= 4, the five weights are

w 1 =w 2 =w 4 =w 5 =

1

4 − 41 /^3

w 3 = 1−(w 1 +w 2 +w 4 +w 5 ).

Since the number of weights grows exponentially with the order, an alternative set of
weights, introduced by Yoshida, proves beneficial. In the Yoshida scheme, a numerical
procedure for obtaining the weights is introduced, leading to a muchsmaller number
of weights. For example, only three weights are needed for a fourth-order scheme, and
these are given by

w 1 =w 3 =

1

2 − 21 /^3

w 2 = 1−w 1 −w 3. (4.11.11)

For a sixth-order scheme, seven weights are needed, and these are obtained numerically
with the result

w 1 =w 7 = 0. 784513610477560

w 2 =w 6 = 0. 235573213359357

w 3 =w 5 =− 1. 17767998417887

w 4 = 1−w 1 −w 2 −w 3 −w 5 −w 6 −w 7. (4.11.12)

Once a set of weights is chosen, the factorization of the operatoris then expressed as

eλ(A^1 +A^2 )≈

n∏sy

α=1

S(wαλ). (4.11.13)