of 10
Supplementary material for: Equilibrium-nonequilibrium ring-polymer molecular
dynamics for nonlinear spectroscopy
Tomislav Beguˇsi ́c,
1,
a)
Xuecheng Tao,
1
Geoffrey A. Blake,
1, 2
and Thomas F. Miller III
1
1)
Division of Chemistry and Chemical Engineering, California Institute of
Technology, Pasadena, California 91125, USA
2)
Division of Geological and Planetary Sciences, California Institute of Technology,
Pasadena, California 91125, USA
(Dated: 13 March 2022)
a)
Electronic mail: tbegusic@caltech.edu
1
I. COMPUTATIONAL DETAILS
A. One-dimensional model
0
0.05
0.1
0.15
0.2
0.25
0.3
-
10
-
5
0
5
10
0
20
40
60
80
100
q
V
(
q
)
FIG. 1. Potentials
V
(
q
) = 0
.
5
q
2
+
aq
3
+
a
2
q
4
of varying degree of anharmonicity, controlled by the
parameter
a
, plotted on the position grid used in quantum calculations.
Here, the response function
R
(
t
2
,t
1
) was computed for 120 steps along both time coordi-
nates with a time step of 0.25, with the exact result obtained by solving the time-independent
Schr ̈odinger equation on a position grid (see Fig. 1) ranging from -10 to 10 with a step of
0.01 and transforming the operators
ˆ
A
,
ˆ
B
, and
ˆ
C
to the eigenstate basis. Depending on the
temperature, between 10 and 60 states were used in the evaluation of the response function.
The same approach was used for computing the exact DKT correlation function. Since
the relation between the DKT correlation function
K
(
t
2
,t
1
) and
R
(
t
2
,t
1
) assumes a Fourier
transform involving both positive and negative
t
1
and
t
2
, DKT correlation functions were
computed for positive and negative times, i.e., from -30 to 30. In addition, the Fourier trans-
form assumes infinite integral boundaries, which can be replaced by final times only if the
integrand is zero beyond some time. This is usually the case because correlation functions
are typically damped by the environment. In our simple system, no damping is applied to the
correlation functions, which is why the relation between the computed
K
(
t
2
,t
1
) and
R
(
t
2
,t
1
)
holds only approximately. To avoid artifacts due to boundary effects, the Fourier transform
of the DKT correlation function,
K
(
ω
2
1
), was multiplied by exp[
0
.
02(
ω
2
1
+
ω
2
2
)]. For the
same reason, only
t
1
,t
2
<
27
.
5 are used to offer a fair comparison between DKT-based and
2
other methods. Figure 2 shows that the numerical error introduced in this way is negligible.
Exact
Exact DKT
Symmetrized DKT
RPMD DKT
β
=
1
0
5
10
15
20
25
0.0
0.5
1.0
1.5
0
5
10
15
20
25
-
1.0
-
0.5
0.0
0.5
1.0
β
=
8
0
5
10
15
20
25
0.0
0.5
1.0
1.5
Time
t
2
0
5
10
15
20
25
-
1.0
-
0.5
0.0
0.5
1.0
Time
t
1
R
(
t
2
,
t
1
=
10
)
R
(
t
2
=
10,
t
1
)
FIG. 2. Same as Fig. 2 of the main text, but comparing two exact quantum simulations of the
response function, using either the direct evaluation or the DKT correlation function, and two
approximate calculations based on the symmetrized DKT correlation function and its RPMD
approximation.
To obtain classical and RPMD results,
K
= 2
17
10
5
initial conditions (
q
(
k
)
,p
(
k
)
) were
sampled from a thermal trajectory computed using the Andersen thermostat.
1
Then, four
trajectories were run for each initial condition: two nonequilibrium trajectories starting from
(
q
(
k
)
,p
(
k
)
±
(
ε
2
/
2)
N
[
∂B
N
(
q
(
k
)
)
/∂q
(
k
)
]) with
ε
= 0
.
01 (see Fig. 3), a forward trajectory starting
from (
q
(
k
)
,p
(
k
)
), and a backward trajectory that is computed by propagating (
q
(
k
)
,
p
(
k
)
)
forward in time. Trajectories were propagated using a second-order symplectic algorithm
2
with a time step of 0.05 for 600 steps, while the observables were evaluated every 5 steps of the
dynamics. Classical results were obtained by setting the number of beads
N
= 1, while
N
=
64 was used for the converged RPMD results. The nonequilibrium and backward trajectories
were used for evaluating the equilibrium-nonequilibrium response function according to
R
(
t
2
,t
1
) =
∂t
1
(
β
2
K
X
k
=1
[
C
N
(
q
(
k
)
+
,t
2
)
C
N
(
q
(
k
)
,t
2
)]
A
N
(
q
(
k
)
t
1
)
)
,
(1)
where the time derivative was computed numerically. The forward and backward equilibrium
trajectories were used for evaluating the RPMD DKT response function for positive and
3
0.005
0.01
0.02
0.04
0.08
0.16
0.32
0.64
1.28
2.56
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
ε
2
parameter
Error
FIG. 3. Effect of the
ε
2
parameter on the equilibrium-nonequilibrium RPMD response function.
Error is measured by Eq. (38) of the main text with
R
exact
(
t
2
,t
1
) replaced by the RPMD response
function computed with
ε
2
= 0
.
01. Figure shows that the results are not affected by the choice of
ε
2
for a wide range of values.
negative
t
1
and
t
2
times.
B. Two-dimensional model
The two-time response function of the two-dimensional model was evaluated for 200 steps
(time step 0.25) in
t
1
and
t
2
. Exact results were obtained in the same way as for the one-
dimensional system, using a two-dimensional grid between
20 and 20 (step 0
.
01) along
each coordinate and the first 50 eigenstates. For the RPMD (
N
= 64) and classical (
N
= 1)
simulations, we averaged the response function over 2
17
samples and the trajectories were
run for 1000 steps (time step 0.05, observables evaluated every 5 steps). External field
parameter was
ε
2
= 0
.
01.
Spectra for the two-dimensional model were evaluated as
S
C
(
ω
2
1
) =
Z
0
dt
2
Z
0
dt
1
R
(
t
2
,t
1
)
e
(
t
1
+
t
2
)
cos(
ω
1
t
1
) cos(
ω
2
t
2
)
,
(2)
where
τ
= 7
.
5. Before taking the discrete cosine transform, the damped response function
was padded with zeros up to
t
1
,t
2
= 1300. The double cosine transform of Eq. (2) is related
to the Fourier transform
S
F
(
ω
2
1
) =
Z
0
dt
2
Z
0
dt
1
R
(
t
2
,t
1
)
e
(
t
1
+
t
2
)
e
1
t
1
2
t
2
(3)
4
by
S
(
ω
2
1
) =
1
2
Re[
S
F
(
ω
2
1
) +
S
F
(
ω
2
1
)]
.
(4)
An alternative is to use the sine transform
S
S
(
ω
2
1
) =
Z
0
dt
2
Z
0
dt
1
R
(
t
2
,t
1
)
e
(
t
1
+
t
2
)
sin(
ω
1
t
1
) sin(
ω
2
t
2
)
(5)
=
1
2
Re[
S
F
(
ω
2
1
)
S
F
(
ω
2
1
)]
,
(6)
which was proposed in Ref. 3 as a way to extract the purely absorptive component of the
two-dimensional spectrum. We show these different transforms in Fig. 4. In this case, the
main off-diagonal peak at (
ω
1
2
)
(Ω
1
= 0
.
5
,
2
= 2
.
0) does not appear in the sine
transform, but is well resolved in the cosine transform. The peaks appearing in the Fourier
transform spectrum (Fig. 4, middle) can be assigned to the Feynman diagrams of Ref. 4.
FIG. 4. Cosine [Eq. (2)], Fourier [Eq. (3), real part], and sine [Eq. (5)] transforms of the damped
two-time response function of the two-mode system, presented in the main text, evaluated with
exact quantum dynamics at
β
= 8.
II. HARMONIC LIMIT
In this Section, we show that the RPMD method proposed in the main text is exact in
the limit
N
→∞
for the harmonic oscillator,
V
(
q
) =
1
2
2
q
2
,
(7)
and when two of the operators
ˆ
A
,
ˆ
B
, and
ˆ
C
are linear in position.
5
A. Quantum response functions
To derive the exact quantum results, we first rewrite
R
(
t
2
,t
1
) =
1
2
Tr
n
ˆ
C
(
t
2
+
t
1
)[
ˆ
B
(
t
1
)
,
[
ˆ
A,
ˆ
ρ
]]
o
(8)
=
1
2
Tr
n
[[
ˆ
C
(
t
2
+
t
1
)
,
ˆ
B
(
t
1
)]
,
ˆ
A
] ˆ
ρ
o
(9)
=
1
2
D
[[
ˆ
C
(
t
2
+
t
1
)
,
ˆ
B
(
t
1
)]
,
ˆ
A
]
E
(10)
and note that
q
(
t
)
,f
q
)] =
1
[ ˆ
p,f
q
)] =
i
f
q
) sin
ωt,
(11)
where
f
(
q
) =
∂f/∂q
.
Then, for
ˆ
A
=
ˆ
B
= ˆ
q
and
ˆ
C
=
C
q
), we have
R
(
t
2
,t
1
) =
1
2
[[
C
q
)
,
ˆ
q
(
t
2
)]
,
ˆ
q
(
t
1
t
2
)]
(12)
=
1
2

i

sin
ωt
2
[
C
q
)
,
ˆ
q
(
t
1
t
2
)]
(13)
=
i
m
ω
sin
ωt
2

i

sin
ω
(
t
1
+
t
2
)
C
′′
q
)
(14)
=
C
′′
q
)
(
)
2
sin
ωt
2
sin
ω
(
t
1
+
t
2
)
,
(15)
where we applied Eq. (11) twice. For
ˆ
A
=
ˆ
C
= ˆ
q
and
ˆ
B
=
B
q
):
R
(
t
2
,t
1
) =
1
2
[[ˆ
q
(
t
2
)
,B
q
)]
,
ˆ
q
(
t
1
)]
(16)
=
1
2

i

sin
ωt
2
[
B
q
)
,
ˆ
q
(
t
1
)]
(17)
=
B
′′
q
)
(
)
2
sin
ωt
1
sin
ωt
2
.
(18)
For
ˆ
B
=
ˆ
C
= ˆ
q
and
ˆ
A
=
A
q
):
R
(
t
2
,t
1
) =
1
2
D
[[ˆ
q
(
t
2
)
,
ˆ
q
]
,
ˆ
A
(
t
1
)]
E
(19)
=
1
2

i

sin
ωt
2
,
ˆ
A
(
t
1
)

(20)
= 0
,
(21)
because scalars commute with operators.
6
B. RPMD response functions
Assuming
ε
2
is small, we can derive the fully equilibrium RPMD from the equilibrium-
nonequilibrium RPMD:
R
(
t
2
,t
1
) =
β
ε
2
[
C
N
(
q
+
,t
2
)
C
N
(
q
,t
2
)]
̇
A
N
(
q
t
1
)
,
(22)
=
*

∂C
N
(
q
t
2
)
∂q
t
2

T
·
M
qp,t
2
·
∂B
N
(
q
)
∂q
̇
A
N
(
q
t
1
)
+
,
(23)
where we used
C
N
(
q
±
,t
2
)
C
N
(
q
t
2
)
±
ε
2
2

∂C
N
(
q
t
2
)
∂p

T
·
N
∂B
N
(
q
)
∂q
(24)
C
N
(
q
t
2
)
±
ε
2
2
N

∂C
N
(
q
t
2
)
∂q
t
2

T
·
∂q
t
2
∂p
·
∂B
N
(
q
)
∂q
(25)
and defined
M
xy,t
=
∂x
t
/∂y
. Equation (23) provides yet another way to employ RPMD in
the simulation of two-time response function
R
(
t
2
,t
1
). Similar to the hybrid equilibrium-
nonequilibrium approach, the equilibrium method also reduces to its classical MD analogue.
5–7
However, this method is less practical due to the need for computing and converging the
stability matrix element
M
qp,t
2
, which has been deemed difficult even in the classical MD
setting. Finally, we note that the equilibrium-nonequilibrium and fully equilibrium methods
give equal results when
ε
2
is small.
Before continuing with specific cases, we rewrite Eq. (23)
R
(
t
2
,t
1
) =
*"
̃
C
N
( ̃
q
t
2
)
̃
q
t
2
#
T
·
M
̃
q
̃
p,t
2
·
̃
B
N
( ̃
q
)
̃
q
̇
̃
A
N
( ̃
q
t
1
)
+
,
(26)
in terms of the ring-polymer normal mode coordinates ̃
q
=
O
·
q
, where
O
is an orthogonal
N
×
N
matrix that diagonalizes the ring-polymer force-constant and
̃
f
( ̃
q
) =
f
(
q
). Specifically,
the first normal mode ̃
q
1
=
Nq
c
is a scalar multiple of the ring-polymer centroid
q
c
=
1
N
P
N
i
=1
q
i
and evolves according to
̃
q
1
,t
= ̃
q
1
cos
ωt
+
1
̃
p
1
sin
ωt.
(27)
7
Then, for
A
(
q
) =
B
(
q
) =
q
and arbitrary
C
(
q
), we have:
R
(
t
2
,t
1
) =
*
̃
C
N
( ̃
q
t
2
)
̃
q
1
,t
2
M
̃
q
1
̃
p
1
,t
2
1
N
̇
̃
A
N
( ̃
q
t
1
)
+
(28)
=
sin
ωt
2
*
̃
C
N
( ̃
q
t
2
)
̃
q
1
,t
2
̇
̃
A
N
( ̃
q
t
1
)
+
(29)
=
sin
ωt
2
*
̃
C
N
( ̃
q
)
̃
q
1
̇
̃
A
N
( ̃
q
t
1
t
2
)
+
(30)
=
β
m
sin
ωt
2
sin
ω
(
t
1
+
t
2
)
*
̃
C
N
( ̃
q
)
̃
q
1
̃
q
1
+
(31)
=
N
(
)
2
*
2
̃
C
N
( ̃
q
)
̃
q
2
1
+
sin
ωt
2
sin
ω
(
t
1
+
t
2
)
(32)
=
1
(
)
2
*
1
N
N
X
i
=1
2
C
(
q
i
)
∂q
2
i
+
sin
ωt
2
sin
ω
(
t
1
+
t
2
)
,
(33)
=
C
′′
RP
(
)
2
sin
ωt
2
sin
ω
(
t
1
+
t
2
)
,
(34)
where we used
̃
A
N
( ̃
q
) =
̃
B
N
( ̃
q
) =
q
c
,
∂q
c
/∂
̃
q
1
= 1
/
N
, and
̃
q
i,t
/∂
̃
p
1
=
δ
i
1
sin
ωt/
(
) in
Eqs. (28) and (29),
̇
̃
A
N
( ̃
q
t
) =
(
ω/
N
) ̃
q
1
sin
ωt
+ (1
/m
N
) ̃
p
1
cos
ωt
to derive Eq. (31),
in which we immediately dropped the second term that vanishes, because
p
1
= 0, and
̃
q
1
exp(
β
N
H
N
) =
(
N/βmω
2
)(
∂/∂
̃
q
1
) exp(
β
N
H
N
) together with the integration by parts
to obtain Eq. (32). Finally, Eq. (33) follows from Eq. (32) because
2
̃
C
N
( ̃
q
)
̃
q
2
1
=

∂q
̃
q
1

T
·
2
C
N
(
q
)
∂q
2
·
∂q
̃
q
1
=
1
N
e
T
·
2
C
N
(
q
)
∂q
2
·
e
=
1
N
2
N
X
i
=1
2
C
(
q
i
)
∂q
i
,
(35)
where we used
∂q/∂
̃
q
1
=
e/
N
and defined the
N
-dimensional vector
e
= (1
...
1)
T
. In
the limit of
N
→ ∞
, Eq. (33) converges to the quantum result (15). Similarly, for
A
(
q
) =
C
(
q
) =
q
and arbitrary
B
(
q
):
R
(
t
2
,t
1
) =
*
1
N
1
sin
ωt
2
̃
B
N
( ̃
q
)
̃
q
1
ω
N
̃
q
1
sin
ωt
1
+
(36)
=
β
m
*
̃
B
N
( ̃
q
)
̃
q
1
̃
q
1
+
sin
ωt
1
sin
ωt
2
(37)
=
B
′′
RP
(
)
2
sin
ωt
1
sin
ωt
2
,
(38)
which reduces to the corresponding quantum expression (18) when converged with respect
8
to the number of beads
N
. Finally, for
B
(
q
) =
C
(
q
) =
q
,
R
(
t
2
,t
1
) =

1
N
1
sin
ωt
2
1
N
̇
̃
A
N
( ̃
q
t
1
)

(39)
=
1
sin
ωt
2
D
̇
̃
A
N
( ̃
q
t
1
)
E
= 0
,
(40)
because equilibrium expectation values are time-independent.
Interestingly, in certain cases, even the classical method (
N
= 1) is exact. This is trivially
true for the last studied case [
B
(
q
) =
C
(
q
) =
q
], but also for certain non-trivial cases, e.g.,
for
C
(
q
)
q
2
in the first studied example [
A
(
q
) =
B
(
q
) =
q
] and
B
(
q
)
q
2
in the second
studied case [
A
(
q
) =
C
(
q
) =
q
] because the corresponding final expressions are independent
of temperature. This can be seen in Fig. 3a of the main text for the anharmonicity parameter
a
= 0,
A
(
q
) =
B
(
q
) =
q
, and
C
(
q
) =
q
2
/
2, where both classical and RPMD results are exact.
III. TWO-MODE SYSTEM IN THE ABSENCE OF COUPLING OR
ANHARMONICITY
In the main text we showed spectra for the system described by the Hamiltonian
ˆ
H
2D
[see Eq. (37) of the main text] and operators
ˆ
A
= ˆ
q
1
and
ˆ
B
=
ˆ
C
= ˆ
q
2
. The corresponding
two-time response function
R
(
t
2
,t
1
) =
1
2
Tr
{
ˆ
q
2
(
t
2
+
t
1
)[ˆ
q
2
(
t
1
)
,
q
1
,
ˆ
ρ
]]
}
(41)
can be shown to vanish in the absence of coupling (
λ
= 0) or anharmonicity (
a
= 0). In the
limit
λ
= 0, ˆ
ρ
= ˆ
ρ
1
ˆ
ρ
2
and, therefore,
R
(
t
2
,t
1
) =
1
2
Tr
q
2
{
ˆ
q
2
(
t
2
+
t
1
)[ˆ
q
2
(
t
1
)
,
ˆ
ρ
2
]
}
Tr
q
1
{
q
1
,
ˆ
ρ
1
]
}
= 0
,
(42)
because the trace of a commutator is zero. If
a
= 0, the system is harmonic, which means
it can be decomposed into normal modes. Since all operators are linear in
q
1
and
q
2
, they
will also be linear in the normal-mode coordinates and, as shown in the previous Section,
the response function is zero for linear operators in a harmonic potential.
REFERENCES
1
H. C. Andersen, J. Chem. Phys.
72
, 2384 (1980).
9
2
M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopoulos, J. Chem. Phys.
133
,
124104 (2010).
3
P. Hamm and J. Savolainen, J. Chem. Phys.
136
, 094516 (2012).
4
L. Vietze, E. H. G. Backus, M. Bonn, and M. Grechko, J. Chem. Phys.
154
, 174201 (2021).
5
S. Saito and I. Ohmine, Phys. Rev. Lett.
88
, 207401 (2002).
6
S. Saito and I. Ohmine, J. Chem. Phys.
119
, 9073 (2003).
7
T. Hasegawa and Y. Tanimura, J. Chem. Phys.
125
, 074512 (2006).
10