of 9
View
Online
Export
Citation
CrossMark
RAPID COMMUNICATION
|
APRIL 07 2022
Equilibrium–nonequilibrium ring-polymer molecular
dynamics for nonlinear spectroscopy
Tomislav Begušić
;
Xuecheng T
ao
;
Geoffrey A. Blake
;
Thomas F
. Miller
, III
J. Chem. Phys.
156, 131
102 (2022)
https://doi.org/10.1063/5.0087156
CHORUS
Articles Y
ou May Be Interested In
Two-dimensional Raman spectroscopy of Lennard-Jones liquids via ring-polymer molecular dynamics
J. Chem. Phys.
(July 2020)
Drafting, kissing, and tumbling of a pair of particles settling in non-Newtonian fluids
Physics of Fluids
(February 2022)
The effects of channel width on particle sedimentation in fluids using a coupled lattice Boltzmann-discrete
element model
Physics of Fluids
(May 2023)
06 October 2023 18:05:37
The Journal
of Chemical Physics
COMMUNICATION
scitation.org/journal/jcp
Equilibrium–nonequilibrium ring-polymer
molecular dynamics for nonlinear spectroscopy
Cite as: J. Chem. Phys.
156
, 131102 (2022); doi: 10.1063/5.0087156
Submitted: 2 February 2022
Accepted: 15 March 2022
Published Online: 7 April 2022
Tomislav Begu
ˇ
si
́
c,
1, a)
Xuecheng Tao,
1
Geoffrey A. Blake,
1
,
2
and Thomas F. Miller III
1
AFFILIATIONS
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
a)
Author to whom correspondence should be addressed:
tbegusic@caltech.edu
ABSTRACT
Two-dimensional Raman and hybrid terahertz-Raman spectroscopic techniques provide invaluable insight into molecular structures and
dynamics of condensed-phase systems. However, corroborating experimental results with theory is difficult due to the high computational
cost of incorporating quantum-mechanical effects in the simulations. Here, we present the equilibrium–nonequilibrium ring-polymer molec-
ular dynamics (RPMD), a practical computational method that can account for nuclear quantum effects on the two-time response function
of nonlinear optical spectroscopy. Unlike a recently developed approach based on the double Kubo transformed (DKT) correlation func-
tion, our method is exact in the classical limit, where it reduces to the established equilibrium-nonequilibrium classical molecular dynamics
method. Using benchmark model calculations, we demonstrate the advantages of the equilibrium–nonequilibrium RPMD over classical
and DKT-based approaches. Importantly, its derivation, which is based on the nonequilibrium RPMD, obviates the need for identifying an
appropriate Kubo transformed correlation function and paves the way for applying real-time path-integral techniques to multidimensional
spectroscopy.
Published under an exclusive license by AIP Publishing.
https://doi.org/10.1063/5.0087156
I. INTRODUCTION
Two-dimensional vibrational spectroscopy is a versatile tech-
nique to study microscopic interactions at their natural, fem-
tosecond time scales.
1
Recently, a series of hybrid spectroscopic
experiments
2–5
involving mid-infrared, far-infrared (or terahertz),
and visible (Raman) pulses, have been developed to study electri-
cal and mechanical anharmonicities,
6–9
structural heterogeneities
of liquids,
10–12
and the couplings between intermolecular and
intramolecular vibrational modes.
13,14
However, the interpretation
of such spectra is still an open question and requires adequate
simulation methods.
12,15–17
Several computational methods have been proposed for
simulating two-dimensional off-resonant Raman and hybrid
terhertz–Raman spectra. Model-based approaches aim at construct-
ing simplified few-dimensional model systems that can be solved
fully quantum mechanically in the presence of a harmonic oscil-
lator bath.
18,19
These approaches can disentangle contributions to
spectra and relate them to different physical effects encoded in the
model parameters. The parameters can be obtained by fitting to the
existing experiment
8
or molecular dynamics (MD) simulations.
20
Alternatively, MD can be used to simulate spectra directly,
6,21–31
as a way to cross-validate the simplified model and further test
the molecular mechanics force field or the
ab initio
quantum
chemistry method used for describing the forces between atoms,
dipole moments, and polarizabilities.
32–34
However, one improve-
ment of MD is highly desired—MD simulations are based on classi-
cal atomic nuclei and neglect nuclear quantum effects, which can
significantly modify the spectra. Indeed, recent two-dimensional
Raman–THz–THz spectra of H
2
O and D
2
O revealed the effect of
isotopic substitution beyond what would be expected from classical
mechanics, pointing at nuclear quantum effects of the light hydrogen
atoms.
12
Quantum simulations of the condensed phase have been
enabled by various semiclassical methods based on the classi-
cal MD framework, including linearized semiclassical initial value
representation (LSC-IVR),
35,36
path-integral Liouville dynamics,
37
centroid molecular dynamics (CMD),
38
and ring-polymer molec-
ular dynamics (RPMD).
39
These methods have proven useful for
the computation of one-time correlation functions
40
related to
J. Chem. Phys.
156
, 131102 (2022); doi: 10.1063/5.0087156
156
, 131102-1
Published under an exclusive license by AIP Publishing
06 October 2023 18:05:37
The Journal
of Chemical Physics
COMMUNICATION
scitation.org/journal/jcp
various dynamical properties, including reaction rates,
35,41–45
diffu-
sion coefficients,
46,47
and one-dimensional vibrational spectra.
48–53
Some of them have also been applied to two-dimensional
electronic
54–58
and infrared spectroscopies,
59–62
in which the quan-
tum subsystem, consisting of, e.g., electronic or high-frequency
vibrational degrees of freedom, can be well defined.
63–65
In contrast,
their application to two-dimensional Raman and hybrid terahertz-
Raman spectroscopic techniques, in which all vibrational degrees of
freedom should be treated on an equal footing, has been limited.
Recently, the group of Batista
66–68
proposed a set of methods that
approximates the symmetric contribution to the so-called double
Kubo transformed (DKT) correlation function, which is a two-time
extension of the original one-time Kubo transformed correlation
function.
40
However, the relation between the symmetrized DKT
correlation function and the spectroscopic response function is only
approximate.
Here, we present an alternative RPMD approach to two-
dimensional spectroscopy, which directly considers the two-time
response function and does not rely on the additional approx-
imation related to the DKT correlation function. To derive the
proposed method, we start from the recently developed nonequi-
librium RPMD
69–73
and employ classical response theory. We
then explore its validity and limitations both theoretically and
numerically.
II. THEORY
In two-dimensional Raman
74,75
and hybrid terahertz-Raman
spectroscopies,
3–5
the signal is measured as a function of two delay
times between three ultrashort light pulses, whose specific sequence
determines the type of spectroscopy. To keep the discussion general,
we consider the time-dependent Hamiltonian
ˆ
H
tot
=
ˆ
H
ˆ
AF
1
(
t
)
ˆ
BF
2
(
t
)
,
(1)
comprised of the system’s field-free Hamiltonian
ˆ
H
and two inter-
actions with the first two ultrashort pulses, controlled, respectively,
by coordinate-dependent operators
ˆ
A
and
ˆ
B
and by time-dependent
functions
F
1,2
(
t
)
representing the pulse shapes. For terahertz pulses,
the interaction operators are the system’s dipole moments and
functions
F
(
t
)
are the electric fields of the light, whereas for the
visible/near-infrared (Raman) pulses or sum-frequency terahertz
excitation,
9
the operators are the system’s polarizabilities and the
time-dependent functions correspond to the squares of the electric
fields.
F
1
is centered at
t
=
t
1
,
F
2
is centered at
t
=
0, and the sig-
nal is measured at
t
=
t
2
, where
t
1
and
t
2
represent the time delays
between the three light pulses. The recorded signal is proportional
to the expectation value
7
ˆ
C
(
t
2
)⟩
=
Tr
[
ˆ
C
ˆ
ρ
(
t
2
)]
(2)
of another coordinate-dependent operator
ˆ
C
, which is again
either the polarizability or dipole moment operator of the
system. In Eq. (2), the system’s equilibrium density operator
ˆ
ρ
=
exp
(
β
ˆ
H
)/
Tr
[
exp
(
β
ˆ
H
)]
evolves under the time-dependent
Hamiltonian (1), i.e.,
ˆ
ρ
(
t
)
=
ˆ
U
tot
(
t
,
−∞
)
ˆ
ρ
ˆ
U
tot
(
t
,
−∞
)
,
(3)
ˆ
U
tot
(
t
f
,
t
i
)
=
T
exp
[
i
̵
h
t
f
t
i
d
τ
ˆ
H
tot
(
τ
)]
.
(4)
Here,
β
is the inverse temperature, and
T
is the time ordering
operator. In practice, the experiments are designed to extract the
components of the signal that are proportional to
F
1
and
F
2
, which
can be accounted for by evaluating the signal as
8
S
(
t
2
)
=
ˆ
C
(
t
2
)⟩
++
ˆ
C
(
t
2
)⟩
+−
ˆ
C
(
t
2
)⟩
−+
+
ˆ
C
(
t
2
)⟩
−−
,
(5)
where
jk
represents the expectation value (2) of the system
evolved under the Hamiltonian (1) with
F
1
(
t
)
(
j
/
2
)
F
1
(
t
)
and
F
2
(
t
)
(
k
/
2
)
F
2
(
t
)
. For sufficiently weak external fields, we can
further invoke the second-order time-dependent perturbation the-
ory and recover the well-known result
74
S
(
t
2
)
=
0
d
τ
2
0
d
τ
1
R
(
τ
2
,
τ
1
)
F
2
(
t
2
τ
2
)
F
1
(
t
2
τ
2
τ
1
)
, (6)
where
R
(
τ
2
,
τ
1
)
=
1
̵
h
2
Tr
{
ˆ
C
(
τ
2
+
τ
1
)[
ˆ
B
(
τ
1
)
,
[
ˆ
A
,
ˆ
ρ
]]
}
(7)
is the corresponding response function, which depends only on the
system’s properties. Because it can be easily convolved with differ-
ent experimental pulses to recover the observed experimental signals
(6), the response function is often the direct target of many theoret-
ical studies. Here, we also intend to focus on the response function
but take a slightly different route in developing the computational
approach. Specifically, we note that for delta pulses,
F
1
(
t
)
=
ε
1
δ
(
t
+
t
1
)
and
F
2
(
t
)
=
ε
2
δ
(
t
)
, the signal measured after the second time
delay
t
2
is proportional to the response function,
S
(
t
2
)
=
ε
1
ε
2
R
(
t
2
,
t
1
)
,
(8)
where
ε
1,2
control the amplitude of the external fields. In the follow-
ing, we derive an approximation to
S
(
t
2
)
under weak delta pulses
and relate it to
R
(
t
2
,
t
1
)
using Eq. (8).
To develop a tractable theory for condensed-phase simulations
of spectra, we replace the exact quantum-mechanical expression (2)
by its nonequilibrium RPMD
69
approximation,
ˆ
C
(
t
)⟩
jk
C
(
t
)⟩
RP
jk
=
dq
dp C
N
(
q
jk
,
t
)
ρ
N
(
q
,
p
)
,
(9)
where
q
and
p
are the positions and momenta of the extended system
comprised of
N
replicas of the original
D
-dimensional system,
ρ
N
(
q
,
p
)
=
e
β
N
H
N
(
q
,
p
)
dq
dp e
β
N
H
N
(
q
,
p
)
(10)
is the corresponding phase-space distribution,
J. Chem. Phys.
156
, 131102 (2022); doi: 10.1063/5.0087156
156
, 131102-2
Published under an exclusive license by AIP Publishing
06 October 2023 18:05:37