We develop a new simulation technique based on path-integral molecular dynamics for calculating ground-state tunneling splitting patterns from ratios of symmetrized partition functions. In particular, molecular systems are rigorously projected onto their J = 0 rotational state by an “Eckart spring” that connects two adjacent beads in a ring polymer. Using this procedure, the tunneling splitting can be obtained from thermodynamic integration at just one (sufficiently low) temperature. Converged results are formally identical to the values that would have been obtained by solving the full rovibrational Schrödinger equation on a given Born–Oppenheimer potential energy surface. The new approach is showcased with simulations of hydronium and methanol, which are in good agreement with wavefunction-based calculations and experimental measurements. The method will be of particular use for the study of low-barrier methyl rotations and other floppy modes, where instanton theory is not valid.

Quantum tunneling leads to a splitting of molecular energy levels that is very sensitive to the height and shape of the underlying potential energy barrier. This makes the tunneling splitting an effective experimental observable for probing the non-equilibrium parts of a potential energy surface (PES).1–4 It can help uncover structural and energetic details in molecular solids, where splittings are typically measured by inelastic neutron scattering5,6 or nuclear magnetic resonance,7,8 and in the gas phase, where it is extracted from high-resolution microwave9–12 and infrared2,13–17 spectra.

These experimental measurements are particularly effective when coupled with high-accuracy calculations, which facilitate the interpretation of spectroscopic data and can be used to test the accuracy of potential energy surfaces. Some of the theoretical frameworks used to this end involve solving the Schrödinger equation,18–24 resulting in a comprehensive set of the system’s energy levels. They must, however, rely on carefully tailored representations of the underlying wavefunctions, PESs, and operators in order to combat the exponential scaling of computational costs with system size (the curse of dimensionality).25–28 One approach to overcome this difficulty is diffusion Monte Carlo (DMC)—a grid-free method that samples the ground-state wavefunction stochastically.29–33 It readily yields zero-point energies but is challenging to use for tunneling splitting calculations, which are especially difficult to converge if the splitting is small. Furthermore, such calculations may contain an uncontrolled error since sampling excited states requires specifying a nodal surface, which is not known a priori.

Another class of readily scalable methods is based on the path-integral formulation of quantum mechanics.34 In such approaches, quantum-mechanical properties are recast as an integral over the configuration space of an extended classical system. In some cases, this integral can be evaluated by steepest descent,35 leading to a family of “instanton” theories36–44 that constitute rigorous asymptotic approximations to the exact quantum result. Instantons correctly capture multidimensional tunneling and zero-point energy effects but are relatively computationally inexpensive since they only require local information along the optimal tunneling pathway. This makes them an attractive alternative to numerically exact techniques when studying complex molecular systems.12,45–49 Despite these advantages, instantons only yield approximate tunneling splitting values, whose accuracy deteriorates in low-barrier anharmonic systems.39 

In this paper, we will focus on molecular systems with low barriers and anharmonic vibrations, such as methyl rotations, for which instanton theory is not appropriate. Historically, studies of large-amplitude torsional modes have helped make the connection between the electronic structure and conformational/reaction energy barriers. Spectroscopic features associated with methyl torsion remain a sensitive probe of intermolecular interactions in solids and of rotation–vibration coupling in isolated molecules. The low barriers can, however, pose a challenge. Methyl groups connected to planar frameworks (e.g., aryl, nitryl, and boryl) exhibit barriers of only a few wavenumbers,50 for which instanton approaches would fail drastically. The torsional barrier for the methyl group in methanol is not quite as low (around 350 cm−1) but still comparable to the harmonic torsional frequency of ∼290 cm−1.51 Such “floppiness” poses a problem not only for instantons but also for wavefunction-based computations, since large-amplitude motions must be represented with care to keep computational expense to a minimum.27,52,53

Developments in sparse-grid methods have pushed the limit on the size of molecules for which full-dimensional wavefunction calculations are feasible (e.g., malonaldehyde in Ref. 54), but further extensions to larger systems, especially ones with multiple large-amplitude modes, remain a challenge. Such systems can be successfully tackled by standard instanton theory.12,55 However, instanton theory relies on a semiclassical approximation, and although its accuracy can, in some cases, be improved with perturbative corrections,56 these are only applicable up to a point and will fail in extreme cases. Nevertheless, strongly anharmonic systems with low barriers can be successfully simulated if one abandons the asymptotic approximation, choosing instead to evaluate the (discretized) path integral nonperturbatively with numerical sampling methods. This brings us, then, to path-integral Monte Carlo (PIMC) and the main focus of our paper, path-integral molecular dynamics (PIMD).57–61 Early work employed path-integral sampling to extract splitting values from low-temperature thermal density matrix elements of model condensed-phase systems.60–63 The approach was adapted by Mátyus and co-workers to study molecular rovibrational tunneling,64,65 and continued methodological developments have enabled its application to water clusters.66–68 

To further expand the scope and accuracy of PIMD tunneling splitting calculations, we will address two complications associated with this approach. First, the original PIMD procedure requires fitting simulation results to a function not just of the tunneling splitting but also of a second “tunneling time” parameter. At the very least, then, calculations must be conducted at two different temperatures in order to obtain a value for the tunneling splitting. Second, Mátyus et al. accomplish rigorous projection onto states of definite angular momentum by performing numerical quadrature over rotations. Hence, in principle, a complete calculation at one temperature itself comprises multiple “sub-calculations,”69 one for every point on the rotational quadrature grid. At the time of writing, we are not aware of any applications to molecular systems that employ the full procedure as outlined above.70 Instead, previous PIMD calculations chose to reduce the computational cost by using only one orientation as defined by the instanton pathway.64,66–68 This neglects rovibrational coupling,108 which, although it may often be reasonably accurate, formally constitutes an uncontrolled approximation.

Here, we derive an alternative means of calculating tunneling splittings that uses a modified PIMD formulation based on symmetrized partition functions rather than density matrix elements. The dynamics involve a new type of spring that accounts for permutational and rotational symmetry in an apparently simple modification that eliminates both the “tunneling time” parameter and the numerical quadrature over orientations. This results in a “one-shot” procedure that is formally exact and only calls for a single (thermodynamic integration) calculation to yield the tunneling splitting. We lay down the theoretical foundations for the approach in Sec. II, which describes how to extract tunneling splittings from symmetrized partition functions. In Sec. III, we give a brief summary of standard PIMD, introduce the necessary modifications for dealing with symmetry projections and rotational degrees of freedom, and summarize how these are used to yield ratios of partition functions. Section IV discusses the numerical tests of the new rotational projection procedure on water and the application of our method to splitting calculations in hydronium and methanol. Section V concludes the article.

Let us consider the tunneling splitting between the totally symmetric ground state and a set of energy levels that all have different symmetries and that are all below the first excited totally symmetric state.65 For reasons explained in Sec. II C, we refer to systems with this property as “rotationless.” One such system is the symmetric one-dimensional double well, with quantum energy levels ɛn (see Fig. 1). In the low-temperature limit, only the lowest-energy states contribute significantly to the partition function
Z=n=0n|eβĤ|n=eβε0+eβε1+.
(1)
In order to extract the tunneling splitting, defined as Δ = ɛ1ɛ0, we introduce the inversion operator σ̂, which acts on coordinates as σ̂x=x and maps the two degenerate wells onto each other. The ground state is symmetric under inversion, while the first excited state is antisymmetric, meaning σ̂|0=|0 and σ̂|1=|1. We can use this property to obtain the low-temperature limit of the “symmetrized” partition function61,71
Zσ=n=0n|eβĤσ̂|n=eβε0eβε1+.
(2)
The tunneling splitting is then readily obtained from
ZZσZ+Zσ=1ZσZ1+ZσZ=eβΔ,
(3)
in the limit as β → ∞.
FIG. 1.

A double-well potential (blue) with a mirror plane (segment drawn in gray) centered on the origin. The dashed lines denote the two potential energy surfaces in a hypothetical system with no tunneling between the minima at ±x0. The corresponding unsplit energy levels are shown as thick horizontal lines on the left and right sides of the plot. Thick horizontal lines in the center depict energy levels of the double well for which the tunneling splitting is observed. Small arrows indicate the ground-state splitting Δ.

FIG. 1.

A double-well potential (blue) with a mirror plane (segment drawn in gray) centered on the origin. The dashed lines denote the two potential energy surfaces in a hypothetical system with no tunneling between the minima at ±x0. The corresponding unsplit energy levels are shown as thick horizontal lines on the left and right sides of the plot. Thick horizontal lines in the center depict energy levels of the double well for which the tunneling splitting is observed. Small arrows indicate the ground-state splitting Δ.

Close modal
The argument can be extended to multidimensional systems with other symmetries in a similar manner to Ref. 65. Consider a Hamiltonian whose symmetries form the group G of order g. Let P̂ be some symmetry operation in this group, for which we define the symmetrized partition function
ZP=TreβĤP̂.
(4)
As in Eq. (2), we can expand the trace in the energy eigenbasis
ZP=nlnl|eβĤP̂|nl,
(5)
where n ranges over all the energy levels, and l ranges over the degenerate states comprising level n. These states form an orthonormal basis for the irreducible representation (irrep) Γ(n) and transform according to72 
P̂|nl=j|njDjl(n)(P̂).
(6)
Here D(n)(P̂) is a representation matrix, whose trace
χP(n)=jDjj(n)(P̂)
(7)
is known as the character of P̂ in Γ(n). By virtue of Eqs. (6) and (7), and the orthonormality of |nl⟩, we can rewrite Eq. (5) as
ZP=n,l,jnl|eβĤ|njDjl(n)(P̂)=n,l,jeβεnδljDjl(n)(P̂)=neβεnχP(n).
(8)
Note that if two or more operations belong to the same conjugacy class α of order gα, they have identical characters.72,73 It is, therefore, sufficient to calculate ZP = Zα for a single member of a conjugacy class. To extract individual energy levels from Eq. (8), we can, in general, use the character orthogonality relation
αgαχα(m)*χα(n)=gδΓ(m)Γ(n),
(9)
where the summation is over all the conjugacy classes of G, χα(m)* denotes the complex conjugate of χα(m), and δΓ(m)Γ(n) is 1 if energy levels m and n correspond to the same irrep, and 0 otherwise. From Eq. (9) it follows that
αgαgχα(m)*Zα=neβεnδΓ(m)Γ(n)=eβεm,
(10)
for β → ∞, with ɛm denoting the lowest energy level whose states transform as Γ(m). The tunneling splitting pattern (Δm,n = ɛmɛn) may, therefore, be deduced by evaluating expressions of the form
eβΔm,n=αgαχα(m)*Zααgαχα(n)*Zα,
(11)
using tabulated values of gα and χα(m). The expression includes Zα for every conjugacy class, but typically one needs only to calculate a subset of those, as illustrated in the following example.
Consider the C3v point group, which describes the symmetry of a (2π/3)-periodic potential, sketched in Fig. 2, and whose characters are given in Table I. Cs is a subgroup of C3v describing the feasible symmetries74 of the hypothetical non-tunneling version of the system confined to a single potential energy well (our “reference system”). Following Bunker and Jensen,73 we determine the representation in C3v that is induced by the irrep of a reference state in Cs. For the ground state, we must consider the totally symmetric irrep (A′ in the Cs point group). Its induced representation in C3v is A1E, which gives the symmetries of the tunneling-split states, namely, a totally symmetric A1 ground state and a doubly degenerate E excited state. There is, therefore, just one tunneling splitting to consider, Δ1,0 ≡ Δ, given by Eq. (11) as
eβΔ=2Z2ZC3Z+2ZC3+3Zσv.
(12)
Taking the low-temperature limit of Eq. (8), we find that
Z=eβεA1+2eβεE,
(13a)
ZC3=eβεA1eβεE,
(13b)
Zσv=eβεA1,
(13c)
where ɛΓ denotes the lowest energy level that transforms according to the irreducible representation Γ. Equation (13) does not contain any terms in εA2, since, based on our symmetry analysis, this energy derives from an excited vibrational state of the reference system. This is assumed to be well separated from the tunneling-split ground state, such that the Boltzmann factor eβεA2 becomes negligible compared to the other terms in Eq. (13) as β → ∞.
FIG. 2.

(a) Contour plot of a model potential with C3v symmetry. The solid black line denotes the symmetry plane of the reference non-tunneling system confined to the top well and described by the Cs point group. The full C3v group includes two further symmetry planes (dotted black lines) and two three-fold rotations (curved arrows). (b) The tunneling splitting pattern for the A′ ground state and the A″ vibrationally excited state of the reference system (the energy diagram does not include any excitations in the radial mode). Thick black lines denote the energy levels before (Cs) and after (C3v) tunneling splitting.

FIG. 2.

(a) Contour plot of a model potential with C3v symmetry. The solid black line denotes the symmetry plane of the reference non-tunneling system confined to the top well and described by the Cs point group. The full C3v group includes two further symmetry planes (dotted black lines) and two three-fold rotations (curved arrows). (b) The tunneling splitting pattern for the A′ ground state and the A″ vibrationally excited state of the reference system (the energy diagram does not include any excitations in the radial mode). Thick black lines denote the energy levels before (Cs) and after (C3v) tunneling splitting.

Close modal
TABLE I.

Character table for the C3v point group. Column headings refer to conjugacy classes, with leading integers giving the corresponding order gα, understood to be 1 if not stated explicitly. Conventional irreducible representation labels are used as row headings.

C3vE2C33σv
A1 
A2 −1 
E −1 
C3vE2C33σv
A1 
A2 −1 
E −1 
Given this, there are several possible ways to derive a simplified expression for Δ. The first possibility is to re-evaluate Eq. (13) in the smallest subgroup that still assigns different symmetry labels to the different energy levels of the tunneling-split ground state. Here, the appropriate subgroup is C3. The second possibility is to express one of the symmetrized partition functions in terms of the others, e.g., Zσv=13[Z+2ZC3], followed by substitution into Eq. (12). Yet another option is to pick as many expressions as there are non-negligible Boltzmann factors, e.g., Eqs. (13a) and (13b), and to solve the simultaneous linear equations for these factors (and, hence, for the tunneling splittings) directly. The three approaches all give the same result,75 which for our example reads
eβΔ=ZZC3Z+2ZC3=1ZC3Z1+2ZC3Z
(14)
and is subtly different from the simple double-well case described by Eq. (3). The general procedure for obtaining such expressions can be summarized as follows:
  1. Determine the symmetries of the tunneling-split states.

  2. Obtain the low-temperature limit of Eq. (8) for a subset of all conjugacy classes, keeping only the Boltzmann factors that correspond to the states in step 1.

  3. Once the set contains as many linearly independent expressions as the number of levels in step 1, solve for the Boltzmann factors.

  4. Express the tunneling splitting constant(s) in terms of the symmetrized partition functions by taking ratios of the expressions in step 3.

For this introductory example, we have considered a point group such as one might encounter in discussions of rigid molecules. However, the results of this section also apply to non-rigid molecules, for which the P̂ operators represent the permutations of indistinguishable atoms and inversion through the molecular center of mass (CoM). Feasible combinations of these operations belong to the molecular symmetry group (MS),74 which replaces the point group in subsequent sections. For completeness, we note that it is with the help of the MS group that nuclear spin-statistical weights can be assigned to tunneling-split energy levels.73 

In addition to discrete symmetries of the kind described in Sec. II A, a molecular Hamiltonian is invariant under rotation.76 Consequently, its eigenfunctions can be labeled by two further quantum numbers: the total angular momentum quantum number, J, and the projection quantum number, M. The latter gives the projection of the angular momentum onto the space-fixed z-axis and ranges in integer steps from −J to J. As we consider systems in the absence of external electromagnetic fields, states that only differ by the value of M are degenerate.

The energy separation between states of different J can, in some cases, be similar to or even smaller than a rovibrational tunneling splitting. This means that simply taking the low-temperature limit of Eq. (11), as outlined in Secs. II A and II B, will not be sufficient to isolate contributions from the tunneling-split ground state. We must instead tackle this problem by explicitly selecting for a particular J, i.e., projecting onto states with a specific total angular momentum.

From here on, we assume that the molecule in question is closed-shell, so there are no interactions between rotational angular momentum and electron spin. We also neglect any hyperfine structure due to coupling to nuclear spin. Under these assumptions, we can denote the eigenfunctions of the molecular Hamiltonian as |n, J, M⟩, where n is shorthand for any other relevant quantum numbers or symmetry labels. These eigenfunctions transform as
R̂(Ω)|n,J,M=L=JJ|n,J,LDLM(J)(Ω),
(15)
where R̂(Ω) is a rotation specified by Ω (which could be a set of Euler angles, or a unit quaternion, or any other valid representation), and DLM(J)(Ω) is a Wigner D-matrix.77 The latter satisfy an orthogonality relation
dΩDLM(J)(Ω)*DLM(J)(Ω)=δJJδLLδMM8π22J+1,
(16)
which can be used to construct a projection onto a state with a specified J. By analogy with Sec. II A, we define a “rotated” partition function
Z(Ω)=TreβĤR̂(Ω)
(17a)
=n,J,Mn,J,M|eβĤR̂(Ω)|n,J,M
(17b)
=n,JM,LDLM(J)(Ω)n,J,M|eβĤ|n,J,L,
(17c)
and use Eq. (16) to get
DMM(J)(Ω)*Z(Ω)dΩ=8π22J+1nn,J,M|eβĤ|n,J,M
(18a)
=8π2(2J+1)2ZJ=J.
(18b)
On going from Eq. (18a) to Eq. (18b), we recall the (2J′ + 1)-fold degeneracy of the states that only differ by the projection quantum number and recognize that the resulting sum is the contribution to the total partition function from states with angular momentum J′. Given that D00(0)(Ω)=1, Eq. (18) can be rearranged to read
ZJ=0=18π2Tr[eβĤR̂(Ω)]dΩ,
(19)
i.e., states with J = 0 can be isolated by taking an unweighted average over all “rotated” partition functions.
Having treated rotational symmetries, we can extract ground-state tunneling splitting constants using the same approach as in Sec. II B. In a straightforward generalization, we define
ZP(Ω)=TreβĤP̂R̂(Ω),
(20a)
ZPJ=0=18π2ZP(Ω)dΩ.
(20b)

Given a ZαJ=0 for every conjugacy class, α, of the MS group (or a subgroup thereof), Eq. (11) directly yields the desired splitting constants.78 

For Eq. (11) to be of practical use, we need a means of calculating partition functions without having to solve the Schrödinger equation. The path-integral formulation of quantum mechanics is particularly suited for this purpose and is easily extended to obtain the symmetrized partition function, ZP. First, we apply the symmetric Trotter product formula,79 
eβĤ=limNeβNV̂/2eβNT̂eβNV̂/2N,
(21)
where βN = β/N, V̂ is the potential energy operator, and T̂ is the kinetic energy operator. Denoting the bracketed product of operators as ÔN, in the limit as N → ∞,
ZP=dx1x1|(ÔN)NP̂|x1=dxxN|ÔNP̂|x1×i=1N1xi|ÔN|xi+1.
(22)
To arrive at the right-hand side of Eq. (22), we inserted N − 1 position resolutions of the identity. The bold cursive x refers to the set of N vectors {xi|i = 1…N}, where xi=(x1(i),x2(i),,xf(i)), and f is the number of configurational degrees of freedom (DoF). For the integration, we use the shorthand notation
dxi=1Ndxi.
(23)
Upon evaluating the matrix elements in Eq. (22), we can express the partition function as
ZP=(2π)NfdpdxeβNHNP(p,x),
(24a)
HNP=HN(open)+k=1fmkωN22xk(N)P̂xk(1)2,
(24b)
HN(open)=i=1Nk=1f(pk(i))22mk+V(xi)+i=1N1k=1fmkωN22xk(i+1)xk(i)2,
(24c)
where mk is the mass of the k-th DoF, and p refers to the set of momenta {pi|i = 1…N}. The Hamiltonian in this expression describes an extended classical system comprised of N replicas of the original set of particles that reside on the potential energy surface V(x). Replicas with adjacent indices (i and i ± 1) are connected by harmonic springs of frequency ωN = 1/βN, so that the term HN(open) corresponds to a linear polymer or an open necklace in which replicas are the “beads.”

The second term in Eq. (24b) introduces a spring between the ends of the necklace (beads 1 and N). When the symmetry operator P̂ is taken to be the identity Ê, where Êx=x, the additional spring connects the extended classical system into a “ring polymer.” In this case, ZE refers to the standard quantum partition function, Z, and Eq. (24a) represents the discretized version of a Feynman path integral expression for the trace of the Boltzmann operator.34 Crucially, Eq. (24a) is in the form of a classical partition function (for the ring polymer at a temperature T × N). As such, it can be evaluated using well-established techniques from classical statistical mechanics. Because of the connection to the path-integral formalism, simulations that sample Eq. (24a) using classical dynamics generated by HNP (in combination with a suitable thermostat to ensure canonical sampling) are known as path-integral molecular dynamics (PIMD).59 Many PIMD techniques are easily compatible with other choices of P̂, for which the cyclic symmetry of the ring polymer is broken. As shown in Secs. IV B and IV C, the ends of the necklace, x1 and xN, tend to explore degenerate minima related by the symmetry operation P̂, and the intermediate beads sample the potential energy barrier connecting the two low-energy regions.

We note in passing that the symmetrized partition functions used here also arise in the context of bosonic and fermionic quantum statistics.71 For example, Eq. (10) evaluated for the totally symmetric irrep of the complete nuclear permutation-inversion (CNPI) group73 is the contribution to the bosonic partition function from the states that have a totally symmetric nuclear spin part; evaluated for the antisymmetric irrep, it gives the contribution from the states with an antisymmetric nuclear spin part. The contributions to the fermionic partition function can be obtained in an equivalent way.80 Other irreps could be paired with nuclear spin states of different symmetries or ascribed to hypothetical particles with more exotic quantum statistics. Due to this formal connection, we could in principle make use of the PIMC61 and PIMD methods81,82 developed for simulating the quantum thermodynamics of bosonic and fermionic systems. For now, however, it seems more pragmatic to use the PIMD formulation presented in this section, combined with thermodynamic integration (TI) as described in Sec. III C.

The distinguishing feature of the J = 0 partition function in Eq. (19) is the additional integral over orientations, Ω. Repeating the manipulations from Sec. III A results in the following “polymer” partition function:
ZJ=0=(2π)Nf18π2dpdreβNHN(open)(p,r)×dΩexpβNωN22a=1natmara(N)R̂(Ω)ra(1)2,
(25)
where nat is the number of atoms in the molecule, f = 3nat, and ra(i) is a three-dimensional coordinate vector xa(i),ya(i),za(i) describing the a-th atom of mass ma in the i-th system replica (bead). The integral over Ω can be expressed as
I(r)=18π2dΩexpNfΩ;r(N),r(1)2β2,
(26)
where
fΩ;r(N),r(1)=a=1natmara(N)R̂(Ω)ra(1)2.
(27)
For the purposes of evaluating Eq. (26), we can approximate fΩ;r(N),r(1) by a Taylor series expansion in Ω about its minimum, truncated at second order. In doing so, we replace the (rather non-trivial) integral over Ω with a three-dimensional Gaussian integral, which can be evaluated analytically. It can be shown that the relative error introduced by this approximation vanishes in the many-bead (N → ∞) limit, which is the same limit as we generally aim for in PIMD simulations. This way of approximating integrals is known as Laplace’s method or, more generally, the method of steepest descent.35 It is a standard technique of asymptotic analysis and forms the foundation of instanton theory.40 However, unlike in instanton theory, the approximation here can be made exact by increasing N. The error scales as O(N−1), in principle decaying more slowly than the O(N−2) error from ring-polymer discretization.83 Despite the formally slower scaling, we observe that the Eckart spring does not negatively impact the rate at which tunneling splitting calculations converge with respect to the number of beads.

For the steepest-descent approach to be practical, we need an efficient way of locating the minimum of f(Ω) (from here on, the parametric dependence on r(N) and r(1) is dropped for brevity). The problem amounts to rotating replica 1 about its CoM such that the mass-weighted square distance to replica N is minimized. This concept is closely tied to the Eckart frame, which is central to the construction of the rovibrational Hamiltonian84,85 and has previously featured in another PIMD method developed by one of the authors.86 Note, however, that the Eckart frame is not used to define a set of curvilinear coordinates in the present work; all our expressions are in a Cartesian frame.

The algorithm for finding the Eckart rotation is comprehensively laid out in Ref. 87. Here, we summarize the key details of the procedure. First, the rotation is represented by a unit quaternion q = (q0, q1, q2, q3) with a real part q0. This is related to the rotation matrix via
R(Ω)=12q222q322q1q2+2q0q32q1q32q0q22q1q22q0q312q122q322q2q3+2q0q12q1q3+2q0q22q2q32q0q112q122q22,
(28)
such that
R̂(Ω)ra(1)=R(Ω)ra(1)rcom(1)+rcom(1),
(29)
with
rcom(1)=a=1natmara(1)/a=1natma
(30)
denoting the CoM of replica 1. For the rest of this section, it is assumed that we are working in the CoM frame of replica 1, such that rcom(1)=0. With the definitions in Eqs. (28)(30), it can be shown that the function in Eq. (27) is a quadratic form in q
f(Ω)=qCq,
(31)
with C denoting a 4 × 4 real symmetric matrix that is a simple function of ma, ra(1), and ra(N) (a = 1, …, nat). The explicit form of C is given in Eq. (24) of Ref. 87 where, in our case,
x+a=xa(N)+xa(1),xa=xa(N)xa(1),
(32)
and similarly for y and z coordinates. The quadratic form can be diagonalized by finding the eigenvalues and eigenvectors of C. The smallest eigenvalue λ0 corresponds to the minimum of f(Ω), and the associated normalized eigenvector q(0) is the quaternion that accomplishes the rotation onto the minimum. In our case, the ring-polymer potential ensures that r(N) and r(1) have a strong structural similarity and can be rotated onto each other such that they are nearly coincident. The similarity of the structures means that the optimal rotation is unique, unless the molecules fluctuate about a (nearly) linear geometry.88 The construction and diagonalization of C are numerically inexpensive relative to the calculation of the PES, so the outlined procedure can be readily incorporated into a PIMD simulation.
To proceed with the steepest-descent approximation, we denote the optimal rotation with R(Ω̃), and the corresponding rotated structure with r̃a(1)=R(Ω̃)ra(1). Using R(Ω)=R(δΩ)R(Ω̃), the function f(Ω) is rewritten as
f(Ω)=a=1natmara(N)r̃a(1)2+2a=1natmara(N)[IR(δΩ)]r̃a(1),
(33)
where I is the 3 × 3 identity matrix. Next we express the rotation R(δΩ) about the minimum in terms of Tait–Bryan angles ϕ,θ,χ,89 
R(δΩ)=cθcχcθsχsθsϕsθcχ+cϕsχsϕsθsχ+cϕcχsϕcθcϕsθcχ+sϕsχcϕsθsχ+sϕcχcϕcθ,
(34)
where cϕ = cos ϕ, sϕ = sin ϕ, etc. The corresponding integral over orientations is
dΩ=ππdϕπ/2π/2cosθdθππdχ=ππdϕ11dzππdχ,
(35)
with z = cos θ. Now we expand IR(δΩ) in a Taylor series about ϕ,z,χ=0. Up to second order,
IR(δΩ)0χzχ0ϕzϕ0+12z2+χ2002zϕϕ2+χ202ϕχ2zχz2+ϕ2.
(36)
The first (linear) term of Eq. (36) contributes
2a=1natmara(N)0χzχ0ϕzϕ0r̃a(1)
(37)
to f(Ω). This evaluates to zero because r̃a(1) satisfies the rotational Eckart conditions.85 Hence, only the second- (and higher-) order terms make a non-zero contribution around the stationary point
f(Ω)a=1natmara(N)r̃a(1)2+δΩΘδΩ,
(38)
where δΩ stands for ϕ,z,χ, and Θ = Θ[r(N), r(1)] is a 3 × 3 matrix with elements
Θij[r(N),r(1)]=a=1natma(ra(N)r̃a(1))δijra,i(N)r̃a,j(1),
(39)
which is symmetric because of the rotational Eckart conditions. This finally allows us to approximate the integrand in Eq. (26) by a Gaussian, so that
dΩexpβNωN22a=1natmara(N)R̂(Ω)ra(1)22πβNωN23/21detΘ1/2expβNωN22a=1natmara(N)r̃a(1)2.
(40)
Under this steepest-descent approximation (which becomes exact as N → ∞),
ZJ=0=(2π)Nfdpdru(r)eβNH̃N(p,r),
(41a)
H̃N=HN(open)+a=1natmaωN22ra(N)r̃a(1)2,
(41b)
u(r)=18π22πβNωN23/21detΘ[r(N),r(1)]1/2,
(41c)
giving the final form of the rotationally projected partition function.

The projection introduces into the integrand a prefactor u, proportional to det Θ−1/2. The matrix Θ is approximately equal to the molecule’s inertia tensor, and its determinant will fluctuate with a relatively small amplitude over the course of a simulation. The thermal expectation value of u can, therefore, be calculated efficiently using standard PIMD techniques.

The other change is the replacement of the harmonic spring between replica 1 and replica N with an “Eckart spring,” as defined by the sum in Eq. (41b). This modified spring does not specifically favor r1 and rN that are close to each other. Instead, low energies are achieved for replica pairs whose CoMs and shapes are near coincident, regardless of relative orientation. The two end beads will tend to explore the same potential energy well but at a range of relative orientations. This stems from the rotational averaging that led us to ZJ=0 in the first place. The forces due to the Eckart spring, required to simulate such trajectories, can be readily obtained from the C-matrix and the normalized eigenvector q(0) associated with its smallest eigenvalue:
rb(i)a=1natmara(N)r̃a(1)2=μ,ν=03qμ(0)Cμνrb(i)qν(0).
(42)
The same line of reasoning can be applied to the symmetrized partition function in Eq. (20b). Going through the derivation, we find that Eq. (41) still applies, but with a Hamiltonian
H̃NP=HN(open)+a=1natmaωN22P̂1ra(N)r̃a(1)2,
(43)
where P̂1 is the inverse of P̂. The prefactor is similarly modified, now depending on Θ[P̂1r(N),r(1)] defined as in Eq. (39). The inclusion of an Eckart spring with a permutation–inversion operation similarly allows the end beads to assume a range of relative orientations, except now they tend to explore different potential energy wells, with minima related by the operation P̂.
The general expression for the tunneling splitting, Eq. (11), and its particular realizations for the double [Eq. (3)] and triple well [Eq. (14)] are (or can be) written in terms of ratios of symmetrized partition functions. These functions are of the form
ZPJ=0=(2π)Nfdpdru(r)eβNH̃NP(p,r).
(44)
The ratio of two partition functions that are symmetrized by the operators P̂1 and P̂0, respectively, can be expressed as ZP1J=0/ZP0J=0=AeβNΔF, where
eβNΔF=dpdreβNH̃NP1p,rdpdreβNH̃NP0p,r
(45)
and A=u1(r)P1/u0(r)P0. Here, ui(r) refers to Eq. (41c) evaluated with Θ[P̂i1r(N),r(1)], and
u(r)P=dpdru(r)eβNH̃NP(p,r)dpdreβNH̃NP(p,r)
(46)
denotes the thermal average of u for a system described by the Hamiltonian H̃NP. Such averages are directly available from standard PIMD simulations. The term in Eq. (45) presents more of a challenge; it is a ratio of partition functions, which we get by calculating the associated free energy difference, ΔF, using thermodynamic integration (TI). To this end, we define
Zξ=(2π)NfdpdreβNH̃ξ,
(47a)
H̃ξ=(1ξ)H̃NP0+ξH̃NP1,
(47b)
and use the relation
ΔF=1βN01dlnZξdξdξ=01H̃NP1H̃NP0ξdξ.
(48)
In practice, the integral over ξ is approximated using a numerical quadrature scheme, which prescribes a set of grid points at which to calculate the integrand. As in previous work,64,66–68 we find that Clenshaw–Curtis (CC) and Gauss–Legendre (GL) quadrature90 both offer a substantial improvement over simpler schemes, e.g., Simpson’s rule. At every grid point, the integrand is in the form of a thermal average for a system described by the mixed Hamiltonian H̃ξ. Most of the terms in the average cancel exactly; the only remaining contribution is the difference of potential energies in H̃NP1 and H̃NP0 due to the springs connecting replica 1 and replica N.

Having obtained all the necessary averages, we calculate their weighted sum with the weights prescribed by the quadrature scheme. This yields a numerical approximation to the ξ-integral in Eq. (48) and, hence, a value for ΔF. PIMD simulations at ξ = 0 and 1 produce the prefactor A, which combines with ΔF to give us the target ratio ZP1J=0/ZP0J=0=AeβNΔF. This can be substituted into Eq. (11), to give the desired tunneling splitting constant.

The entire process consists of running PIMD simulations (no wavefunction calculations are involved). To propagate the dynamics and ensure proper sampling, we use a multiple-time-step (MTS) integrator combined with an Andersen thermostat.59,91 The “slow” force, used to update the momenta at time intervals of Δt, derives from the “external” potential iV(ri). The “fast” force, used to update momenta at intervals of δt, derives from the spring potential, which includes the Eckart spring and a permutation of atomic indices. The procedure used to obtain confidence intervals for the reported tunneling splittings is described in the  Appendix.

An isolated water molecule (H2O) does not exhibit a tunneling splitting. However, we can still use this system as a test of our procedure for projecting onto states of zero angular momentum (J = 0), as described in Sec. III B. A water molecule is small enough that its quantum energy levels can be calculated numerically using only modest computational resources. To this end, we used the DVR3D program suite developed by Tennyson and co-workers,92 from which we obtained the ratio ZJ=0/Z at two different temperatures. As it was not our goal to reproduce experimental data, we used the simple qTIP4P/F potential to model interatomic interactions, with the same parameters as in Ref. 93 and with atomic masses set to 1.007 825 Da and 15.994 915 Da for H and O, respectively.

PIMD simulations were all performed with a conservative integration time step of 0.1 fs and an MTS factor of Δt/δt = 8. A CC quadrature grid90 with nξ = 33 points was used for thermodynamic integration. This scheme is particularly convenient for running convergence tests since a 17-point CC grid can be constructed by taking every other point in a 33-point grid. In Eq. (45), we replaced H̃NP0 with HN, the standard ring-polymer Hamiltonian. The corresponding prefactor u0(r) is 1. For H̃NP1, we used H̃N from Eq. (41b), and u1 was defined as in Eq. (41c). Every set of TI calculations consisted of 32 independent trajectories initialized at the equilibrium geometry and propagated in parallel. Each one was thermalized for tth = 5 ps under HN, subject to an Andersen thermostat with an average resampling time of 25 fs for every particle. The same propagator and thermostat were used for production simulations, which ran sequentially over Hξ [Eq. (47b)] in order of increasing ξ, starting at ξ = 0. At every new grid point, we ran a trajectory of trun = 25 ps, starting from the configuration at which the simulation for the preceding grid point had completed. The energy differences in Eq. (48) were recorded at intervals of 2 fs, and the values of the prefactor in Eq. (41c) were recorded every 5 fs. Data from the first 500 fs of every production trajectory were discarded.

Our simulation results are summarized in Table II, from which we can see that nξ = 17 CC quadrature points and N = 48 replicas are sufficient to converge the PIMD results, with an estimated statistical error <1%. At both simulation temperatures, the results are within one standard deviation from the exact (DVR3D) values, verifying our procedure for projecting onto J = 0 states.

TABLE II.

ZJ=0/Z × 102 for an isolated molecule of qTIP4P/F water. The number of replicas used in PIMD calculations is indicated in curly brackets, and the TI grid sizes are given by nξ. Errors in the last digit (in parentheses) are reported at 2σ. The exact results were obtained by running DVR3D92 and summing the resulting Boltzmann factors.

T (K)PIMD {N = 48}PIMD {N = 64}
nξ = 17nξ = 33nξ = 17nξ = 33Exact
300 1.21(2) 1.21(1) 1.20(1) 1.21(1) 1.204 
200 2.20(2) 2.20(1) 2.20(3) 2.20(2) 2.200 
T (K)PIMD {N = 48}PIMD {N = 64}
nξ = 17nξ = 33nξ = 17nξ = 33Exact
300 1.21(2) 1.21(1) 1.20(1) 1.21(1) 1.204 
200 2.20(2) 2.20(1) 2.20(3) 2.20(2) 2.200 

The hydronium ion (H3O+) is one of the simplest molecules to exhibit tunneling splitting, similar to the isoelectronic NH3 but with a lower barrier and, correspondingly, larger splitting. The inversion barrier is sufficiently low and anharmonic for instanton theory to be in error by over 20%,42 making hydronium a suitable test system for our method. To describe the interatomic interactions, we have chosen the PES developed by Yu and Bowman,21 which was fitted to a set of electronic energies obtained using coupled-cluster theory with singles, doubles, and perturbative triples [CCSD(T)-F12b/aug-cc-pVQZ], and has a barrier of 674.3 cm−1. The lowest vibrational frequency of ∼887 cm−1 corresponds to the tunneling umbrella mode, from which it can be deduced that only one pair of tunneling-split J = 0 states reside below the barrier.94 

The molecular symmetry group of hydronium is isomorphic to D3h,95 and its characters can be directly substituted into Eq. (11) to produce an expression for the tunneling splitting. However, as noted in the discussion after Eq. (11), the analysis can also be performed using the smallest subgroup that distinguishes between the tunneling-split energy levels. In the present case, we use the {E, (12)} subgroup (isomorphic to Cs), whose (12) operation permutes hydrogens 1 and 2 (see Fig. 3 for the atomic indices). The general expression for the tunneling splitting then reduces to the double-well formula in Eq. (3), where σ̂ is replaced with P̂12, the (12) permutation. TI was performed using H̃N [Eq. (41b)] and H̃N(12) [Eq. (43) with P̂1=P̂12] for the H̃NP0 and H̃NP1 Hamiltonians, respectively. The corresponding prefactor functions were u0=detΘ[r(N),r(1)]1/2 and u1=detΘ[P̂12r(N),r(1)]1/2. The same PIMD simulation parameters as in Sec. IV A were used for this section, except that the average resampling time in the Andersen thermostat was increased to 100 fs and Gauss–Legendre (GL) quadrature90 was used instead of CC, as it was found to converge more rapidly with grid size.

FIG. 3.

(a) A low-energy hydronium (H3O+) “polymer” subject to an Eckart spring with a (12) permutation [final term in Eq. (43)]. The hydrogen atoms of beads 1 and N are labeled with the corresponding indices, and the intermediate beads, spread out for clarity, are represented with coiled blue lines; in reality, the oxygen atoms are nearly coincident. A gray rectangle shows the mirror plane that approximately maps the end beads onto each other. This operation is equivalent to the permutation of atoms and rotation shown on the right. (b) A schematic illustrating the calculation of the Eckart spring potential. In the starting configuration (top), the atoms of the Nth replica are permuted (middle), and the 1st replica is rotated about its center of mass (bottom) to match the permuted replica as closely as possible. A good match implies low energy in the Eckart spring.

FIG. 3.

(a) A low-energy hydronium (H3O+) “polymer” subject to an Eckart spring with a (12) permutation [final term in Eq. (43)]. The hydrogen atoms of beads 1 and N are labeled with the corresponding indices, and the intermediate beads, spread out for clarity, are represented with coiled blue lines; in reality, the oxygen atoms are nearly coincident. A gray rectangle shows the mirror plane that approximately maps the end beads onto each other. This operation is equivalent to the permutation of atoms and rotation shown on the right. (b) A schematic illustrating the calculation of the Eckart spring potential. In the starting configuration (top), the atoms of the Nth replica are permuted (middle), and the 1st replica is rotated about its center of mass (bottom) to match the permuted replica as closely as possible. A good match implies low energy in the Eckart spring.

Close modal

For Eq. (3) to give an accurate estimate of the tunneling splitting, the dynamics have to be simulated at a sufficiently low temperature. The lowest vibrational frequency of 887 cm−1 corresponds to a temperature of 1280 K. At this temperature, we can expect some excited vibrational states to have a significant thermal population. Hence, we reduced this temperature by a factor of approximately six, taking 200 K as our starting point. We then ran PIMD simulations for a series of progressively lower temperatures, taking care to converge the estimated tunneling splittings with respect to the number of beads, N, and the number of GL grid points, nξ, at each temperature (see Table III). Convergence is already satisfactory at 150 K for N = 96 and nξ = 15, and the final value of Δ = 53.2(4) cm−1 is taken from the simulation at T = 100 K, N = 128, and nξ = 20 (since Δ was judged converged at this grid size and the simulation with nξ = 30 involved fewer, shorter trajectories and, therefore, resulted in a larger statistical error).

TABLE III.

Ground-state tunneling splittings calculated with PIMD for the hydronium potential from Ref. 21. T specifies the temperature, N gives the number of beads, and nξ is the number of Gauss–Legendre quadrature points. Errors in the last digit (in parentheses) are reported at 2σ. In the bottom part of the table, our most reliable result (PIMD; T = 100 K, N = 128, and nξ = 20 — see main text) is compared against instanton (Inst.), vibrational configuration interaction (VCI), and experimentally measured (Expt.) tunneling splittings.

T (K)Δ (cm−1)
Nnξ = 15nξ = 20nξ = 30
200 48 55.2(7) 55.7(6) 55.2(5) 
64 54.5(8) 54.6(8) 54.2(6) 
150a 64 53.7(3) 53.8(3) 53.4(5) 
96 53.1(3) 52.9(3) 52.8(6) 
100a 96 53.3(4) 53.7(4) 53.3(6) 
128 53.2(4) 53.2(4) 53.2(7) 
Method Inst.b VCI21  PIMDb Expt.14  
Δ (cm−169.8 52.48 53.2(4) 55.35 
T (K)Δ (cm−1)
Nnξ = 15nξ = 20nξ = 30
200 48 55.2(7) 55.7(6) 55.2(5) 
64 54.5(8) 54.6(8) 54.2(6) 
150a 64 53.7(3) 53.8(3) 53.4(5) 
96 53.1(3) 52.9(3) 52.8(6) 
100a 96 53.3(4) 53.7(4) 53.3(6) 
128 53.2(4) 53.2(4) 53.2(7) 
Method Inst.b VCI21  PIMDb Expt.14  
Δ (cm−169.8 52.48 53.2(4) 55.35 
a

For nξ = 15 and 20, the number of independent trajectories was increased to 64, and their duration was doubled to 50 ps per quadrature point.

b

This work.

To put this value into context, we compare it against the instanton tunneling splitting, which we calculated following the protocol in Ref. 40.96 The calculation converges to 69.8 cm−1 at β=8000Eh1 (T ≈ 40 K) for N = 1024, overestimating the PIMD value by over 30%.97 Our other comparison is to the vibrational configuration interaction (VCI) value of 52.48 cm−1 calculated by Yu and Bowman.21 There is a small but statistically significant difference between that and our PIMD result, which may indicate that VCI is not sufficiently converged with basis set size or with respect to the number of modes used to represent the vibrational angular momentum terms. Even so, this discrepancy is small compared to the deviation from the experimentally measured splitting of 55.35 cm−1, which both approaches underestimate by 2–3 cm−1. It is likely that the deviation is caused by imperfections in the PES, and we provide evidence in the supplementary material that it arises from a combination of fitting and basis-set truncation errors.

Methanol (CH3OH) is a multi-well system with a MS group that is isomorphic to C3v (see Sec. II B and Fig. 2). Its ground state is split into an A1 and a doubly degenerate E level by torsional tunneling through a relatively low potential energy barrier.98 With six atoms, methanol is large enough to present a challenge to wavefunction-based approaches24,27 and sufficiently anharmonic to introduce a sizable error into the instanton approximation.

To run our PIMD simulations of methanol, we turned to the potential by Qu and Bowman,51 fitted to CCSD(T)-F12b/aug-cc-pVDZ energies for configurations relevant to the present study, with a torsional barrier of 350 cm−1 and a harmonic torsional frequency of 286 cm−1. Following Sec. II B, symmetry analysis was conducted in the {E, (123), (132)} subgroup, which is isomorphic to the C3 group. The tunneling splitting can, therefore, be calculated using Eq. (14), where it is understood that all constituent partition functions are projected onto J = 0, and the C3 operation is replaced by the (123) cyclic permutation. Our indexing convention is shown in Fig. 4, along with the dominant polymer configuration, corresponding to the minimum of H̃N(132) [Eq. (43) with P̂1=P̂123]. The tunneling pathway traced out by the polymer beads highlights the large extent to which the OH hydrogen participates in the torsional motion of the CH3 group.

FIG. 4.

(a) A minimum-energy methanol (CH3OH) “polymer” subject to an Eckart spring with a (123) permutation. The hydrogen atoms of the first replica are labeled with the corresponding indices; taken together, the beads trace out the torsional tunneling pathway, and the direction of motion is indicated with blue arrows. (b) Newman projections of the first and final replicas (beads). Here, the replicas are mapped onto each other exactly by a (123) permutation and rotation about the center of mass and, therefore, the Eckart spring is fully contracted [see also Fig. 3(b)].

FIG. 4.

(a) A minimum-energy methanol (CH3OH) “polymer” subject to an Eckart spring with a (123) permutation. The hydrogen atoms of the first replica are labeled with the corresponding indices; taken together, the beads trace out the torsional tunneling pathway, and the direction of motion is indicated with blue arrows. (b) Newman projections of the first and final replicas (beads). Here, the replicas are mapped onto each other exactly by a (123) permutation and rotation about the center of mass and, therefore, the Eckart spring is fully contracted [see also Fig. 3(b)].

Close modal

Thermodynamic integration for this system was performed between H̃N and H̃N(132), with prefactor functions u0=detΘ[r(N),r(1)]1/2 and u1=detΘ[P̂123r(N),r(1)]1/2, using a GL quadrature grid. We used a simulation setup analogous to Sec. IV B, with a propagation time step of 0.2 fs, an MTS factor of 8, and a mean resampling time of 100 fs. Further simulation details and calculated tunneling splitting values are given in Table IV.

TABLE IV.

Ground-state tunneling splittings calculated with PIMD for the methanol potential from Ref. 51. Column headings and error specifications follow the convention used in Table III. In the bottom part of the table, the final PIMD result is compared against instanton (Inst.), wavefunction (WF), and experimental (Expt.) tunneling splittings.

T (K)Δ (cm−1)
Nnξ = 5nξ = 10nξ = 15
75a 64 8.4(2) 8.39(18) ⋯ 
96 8.2(3) 8.28(14) 8.33(11) 
50b 96 9.18(15) 9.15(13) 9.21(8) 
128 9.08(12) 9.05(11) 9.03(8) 
35b 128 9.2(2) 9.10(13) 9.09(11) 
196 9.00(19) 9.1(3) 9.14(13) 
Method Inst.c WFd PIMDc Expt.16  
Δ (cm−111.3 9.15 9.14(13) 9.12 
T (K)Δ (cm−1)
Nnξ = 5nξ = 10nξ = 15
75a 64 8.4(2) 8.39(18) ⋯ 
96 8.2(3) 8.28(14) 8.33(11) 
50b 96 9.18(15) 9.15(13) 9.21(8) 
128 9.08(12) 9.05(11) 9.03(8) 
35b 128 9.2(2) 9.10(13) 9.09(11) 
196 9.00(19) 9.1(3) 9.14(13) 
Method Inst.c WFd PIMDc Expt.16  
Δ (cm−111.3 9.15 9.14(13) 9.12 
a

Thermalized for 5 ps, propagated for 25 ps at every grid point, 32 independent trajectories.

b

Thermalized for 10 ps, propagated for 50 ps at every grid point, 64 independent trajectories.

c

This work.

d

Value by Lauvergnat and Nauts27 obtained from sparse-grid wavefunction calculations on the PES from Ref. 19.

Given the lowest vibrational frequency of 286 cm−1 (411 K), we started our simulations at 75 K and progressively lowered the temperature until reaching convergence. Our most reliable result was calculated at a temperature of 35 K with N = 196, giving a final value of Δ = 9.14(13) cm−1. The corresponding instanton tunneling splitting calculation40,96 converged to 11.3 cm−1 at β=20000Eh1 (T ≈ 16 K) and N = 1024.

We are not aware of any wavefunction-based tunneling splitting calculations for the PES employed in this work. However, an earlier PES by Bowman and co-workers,19 fitted to CCSD(T)/aug-cc-pVTZ energies at many of the same molecular configurations, is believed to achieve similar accuracy. This PES was used by Lauvergnat and Nauts in Ref. 27 to calculate the splitting by solving the Schrödinger equation on a sparse grid. Our PIMD estimate is in excellent agreement with their result of 9.15 cm−1, and both are very close to the experimentally measured splitting16 of 9.12 cm−1. Compared to the results in Sec. IV B, this level of agreement is surprising given the relatively small basis set used to calculate the underlying ab initio energies. It seems that the PES benefits from a fortuitous cancellation between fitting and basis-set truncation errors, as shown in the supplementary material. Even so, the key point is that the PIMD method has achieved its goal of probing the accuracy of the PES.

We have presented an exact path-integral method for calculating ground-state rovibrational tunneling splittings from ratios of symmetrized partition functions. By construction, splitting values are the only unknown parameters in our approach, and a rigorous projection onto the J = 0 state is achieved by means of an Eckart spring (avoiding the need to integrate over orientations explicitly). The first feature makes it possible to extract tunneling splittings from just a single thermodynamic integration, offering a simplification over the procedure used by Mátyus and co-workers.64,65 The second feature means that rotational projection is performed in an exact and computationally efficient manner, unlike in previous PIMD studies,64,66–68,108 such that we can fully account for rovibrational coupling effects.

A perturbatively corrected instanton (RPI+PC) calculation56 would likely approach the accuracy of our PIMD method for the molecular systems in this study. However, whereas the perturbative correction will eventually fail, the PIMD calculation should remain robust even in extremely low-barrier systems, for which large-amplitude motions correspond to above-barrier states.99,100 In addition, unlike the current formulation of RPI+PC, our method can be readily extended to J > 0 states, so that in the future we hope to use Eckart-spring PIMD to study the dependence of tunneling splitting on the rotational state.14 We should also be able to calculate the energy levels of fluxional systems, where rotational and vibrational motion cannot be separated even to a first approximation, as in the case of the CH5+ molecular ion.101 

As expected, we have observed that tunneling splittings are very sensitive to the underlying PES. Using potentials fitted to high-accuracy ab initio energies will, therefore, be crucial if we are to interpret and compare tunneling splitting patterns with experimental data in upcoming studies. While for some molecules, including the ones discussed here, such potentials have already been developed, the availability of an accurate PES can, unfortunately, limit the choice of system. At least in part, this problem can be alleviated with the help of machine learning (ML).102 In contrast to wavefunction-based approaches, which essentially require a global PES, PIMD only samples configurations around the instanton(s). One can, therefore, fit to high-accuracy ab initio data around the instanton path,103,104 expanding the scope of our method.

Instantons may also help us optimize the PIMD procedure itself. A major goal is to minimize the sampling time needed to achieve a specified level of statistical accuracy. This may be accomplished with the help of more sophisticated thermostats66,105 or higher-order propagator splittings, beyond the Trotter product formula.79,106 Alternatively, one could use the system-specific information available from an instanton, which, in a sense, describes the optimal tunneling pathway. This may help us devise the optimal quadrature scheme for thermodynamic integration or otherwise improve the efficiency of our free-energy calculations. It should also be noted that, unlike PIMD methods utilizing ring polymers,59,61 we have not fully exploited the cyclic symmetry of the trace. Doing so should lead to better statistical behavior compared to approaches based on density-matrix elements.64 The development of this idea is left for future work.

See the supplementary material for thermodynamic integration curves, intermediate values used in calculating the results in Tables II–IV, and estimates of systematic errors due to inaccuracies in the PES.

The authors thank Joel Bowman, Qi Yu, and Chen Qu for providing the potential energy surfaces for hydronium and methanol. They also acknowledge financial support from the Swiss National Science Foundation through SNSF Project No. 207772.

The authors have no conflicts to disclose.

George Trenins: Formal analysis (equal); Investigation (lead); Methodology (lead); Software (lead); Supervision (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Lars Meuser: Methodology (supporting); Validation (equal). Hannah Bertschi: Methodology (supporting); Software (supporting). Odysseas Vavourakis: Methodology (supporting); Validation (equal); Visualization (supporting). Reto Flütsch: Formal analysis (equal); Investigation (supporting); Software (supporting). Jeremy O. Richardson: Conceptualization (lead); Supervision (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.

The standard error, σu, of the mean of u(r)P̂ [Eq. (46)], is estimated from the variance of the values obtained from independent trajectories. For ΔF [Eq. (48)], the error σΔF is estimated using the bootstrap method.90 These estimates are combined to get the standard error in the ratio Φ=ZP1J=0/ZP0J=0 according to
σΦ=Φ×σu0/u0(r)P02+σu1/u1(r)P12+βNσΔF21/2.
(A1)
Then, for a double-well tunneling splitting derived from Eq. (3) (as in Sec. IV B), the standard error is
σΔ=2σΦβ(1Φ2).
(A2)
For a triple-well splitting derived from Eq. (14) (as in Sec. IV C), the error is
σΔ=3σΦβ(1Φ)(1+2Φ).
(A3)
1.
V. A.
Benderskii
,
D. E.
Makarov
, and
C. A.
Wight
, in
Chemical Dynamics at Low Temperatures
,
Adv. Chem. Phys.
(
Wiley
,
New York
,
1994
), Vol.
88
.
2.
F. N.
Keutsch
and
R. J.
Saykally
, “
Water clusters: Untangling the mysteries of the liquid, one molecule at a time
,”
Proc. Natl. Acad. Sci. U. S. A.
98
,
10533
10540
(
2001
).
3.
A.
Stone
,
The Theory of Intermolecular Forces
, 2nd ed. (
Oxford
,
2013
).
4.
S.
Khazaei
and
D.
Sebastiani
, “
Methyl rotor quantum states and the effect of chemical environment in organic crystals: γ-Picoline and toluene
,”
J. Chem. Phys.
145
,
234506
(
2016
).
5.
M.
Prager
and
A.
Heidemann
, “
Rotational tunneling and neutron spectroscopy: A compilation
,”
Chem. Rev.
97
,
2933
2966
(
1997
).
6.
J.
Colmenero
,
A. J.
Moreno
, and
A.
Alegría
, “
Neutron scattering investigations on methyl group dynamics in polymers
,”
Prog. Polym. Sci.
30
,
1147
1184
(
2005
).
7.
A. J.
Horsewill
, “
Quantum tunnelling aspects of methyl group rotation studied by NMR
,”
Prog. Nucl. Magn. Reson. Spectrosc.
35
,
359
389
(
1999
).
8.
M.
Šimėnas
,
D.
Klose
,
M.
Ptak
,
K.
Aidas
,
M.
Maczka
,
J.
Banys
,
A.
Pöppl
, and
G.
Jeschke
, “
Magnetic excitation and readout of methyl group tunnel coherence
,”
Sci. Adv.
6
,
eaba1517
(
2020
).
9.
G. T.
Fraser
,
R. D.
Suenram
, and
L. H.
Coudert
, “
Microwave electric-resonance optothermal spectroscopy of (H2O)2
,”
J. Chem. Phys.
90
,
6077
(
1989
).
10.
L. H.
Xu
and
J. T.
Hougen
, “
Global fit of rotational transitions in the ground torsional state of methanol
,”
J. Mol. Spectrosc.
169
,
396
409
(
1995
).
11.
T.
Baba
,
T.
Tanaka
,
I.
Morino
,
K. M. T.
Yamada
, and
K.
Tanaka
, “
Detection of the tunneling-rotation transitions of malonaldehyde in the submillimeter-wave region
,”
J. Chem. Phys.
110
,
4131
4133
(
1999
).
12.
J. O.
Richardson
,
C.
Pérez
,
S.
Lobsiger
,
A. A.
Reid
,
B.
Temelso
,
G. C.
Shields
,
Z.
Kisiel
,
D. J.
Wales
,
B. H.
Pate
, and
S. C.
Althorpe
, “
Concerted hydrogen-bond breaking by quantum tunneling in the water hexamer prism
,”
Science
351
,
1310
1313
(
2016
).
13.
D.-J.
Liu
and
T.
Oka
, “
Experimental determination of the ground-state inversion splitting in H3O+
,”
Phys. Rev. Lett.
54
,
1787
1789
(
1985
).
14.
P.
Verhoeve
,
M.
Versluis
,
J. J.
Ter Meulen
,
W. L.
Meerts
, and
A.
Dymanus
, “
Far infrared laser sideband spectroscopy of H3O+: The pure inversion spectrum around 55 cm−1
,”
Chem. Phys. Lett.
161
,
195
201
(
1989
).
15.
N.
Pugliano
and
R. J.
Saykally
, “
Measurement of quantum tunneling between chiral isomers of the cyclic water trimer
,”
Science
257
,
1937
(
1992
).
16.
L.-H.
Xu
,
X.
Wang
,
T. J.
Cronin
,
D. S.
Perry
,
G. T.
Fraser
, and
A. S.
Pine
, “
Sub-Doppler infrared spectra and torsion–rotation energy manifold of methanol in the CH-stretch fundamental region
,”
J. Mol. Spectrosc.
185
,
158
172
(
1997
).
17.
Ö.
Birer
and
M.
Havenith
, “
High-resolution infrared spectroscopy of the formic acid dimer
,”
Annu. Rev. Phys. Chem.
60
,
263
275
(
2009
).
18.
J. M.
Bowman
, “
The self-consistent-field approach to polyatomic vibrations
,”
Acc. Chem. Res.
19
,
202
208
(
1986
).
19.
J. M.
Bowman
,
X.
Huang
,
N. C.
Handy
, and
S.
Carter
, “
Vibrational levels of methanol calculated by the reaction path version of MULTIMODE, using an ab initio, full-dimensional potential
,”
J. Phys. Chem. A
111
,
7317
7321
(
2007
).
20.
A. G.
Császár
,
C.
Fábri
,
T.
Szidarovszky
,
E.
Mátyus
,
T.
Furtenbacher
, and
G.
Czakó
, “
The fourth age of quantum chemistry: Molecules in motion
,”
Phys. Chem. Chem. Phys.
14
,
1085
1106
(
2012
).
21.
Q.
Yu
and
J. M.
Bowman
, “
Ab initio potential for H3O+ → H+ + H2O: A step to a many-body representation of the hydrated proton?
,”
J. Chem. Theory Comput.
12
,
5284
5292
(
2016
).
22.
M. D.
Coutinho-Neto
,
A.
Viel
, and
U.
Manthe
, “
The ground state tunneling splitting of malonaldehyde: Accurate full dimensional quantum dynamics calculations
,”
J. Chem. Phys.
121
,
9207
9210
(
2004
).
23.
S. N.
Yurchenko
,
J.
Tennyson
,
S.
Miller
,
V. V.
Melnikov
,
J.
O’Donoghue
, and
L.
Moore
, “
ExoMol line lists—XL. Rovibrational molecular line list for the hydronium ion (H3O+)
,”
Mon. Not. R. Astron. Soc.
497
,
2340
2351
(
2020
).
24.
T.
Carrington
, Jr.
, “
Perspective: Computing (ro-) vibrational spectra of molecules with more than four atoms
,”
J. Chem. Phys.
146
,
120902
(
2017
).
25.
P. S.
Thomas
,
T.
Carrington
,
J.
Agarwal
, and
H. F.
Schaefer
, “
Using an iterative eigensolver and intertwined rank reduction to compute vibrational spectra of molecules with more than a dozen atoms: Uracil and naphthalene
,”
J. Chem. Phys.
149
,
064108
(
2018
).
26.
H.
Wang
and
M.
Thoss
, “
Multilayer formulation of the multiconfiguration time-dependent Hartree theory
,”
J. Chem. Phys.
119
,
1289
(
2003
).
27.
D.
Lauvergnat
and
A.
Nauts
, “
Quantum dynamics with sparse grids: A combination of Smolyak scheme and cubature. Application to methanol in full dimensionality
,”
Spectrochim. Acta, Part A
119
,
18
25
(
2014
).
28.
S.
Carter
,
J. M.
Bowman
, and
N. C.
Handy
, “
Extensions and tests of ‘multimode’: A code to obtain accurate vibration/rotation energies of many-mode molecules
,”
Theor. Chem. Acc.
100
,
191
198
(
1998
).
29.
D. F.
Coker
and
R. O.
Watts
, “
Structure and vibrational spectroscopy of the water dimer using quantum simulation
,”
J. Phys. Chem.
91
,
2513
(
1987
).
30.
M. A.
Suhm
and
R. O.
Watts
, “
Quantum Monte Carlo studies of vibrational states in molecules and clusters
,”
Phys. Rep.
204
,
293
329
(
1991
).
31.
J. K.
Gregory
and
D. C.
Clary
, “
Calculations of the tunneling splittings in water dimer and trimer using diffusion Monte Carlo
,”
J. Chem. Phys.
102
,
7817
(
1995
).
32.
M.
Quack
and
M. A.
Suhm
, “
Accurate quantum Monte Carlo calculations of the tunneling splitting in (HF)2 on a six-dimensional potential hypersurface
,”
Chem. Phys. Lett.
234
,
71
76
(
1995
).
33.
Y.
Wang
,
B. J.
Braams
,
J. M.
Bowman
,
S.
Carter
, and
D. P.
Tew
, “
Full-dimensional quantum calculations of ground-state tunneling splitting of malonaldehyde using an accurate ab initio potential energy surface
,”
J. Chem. Phys.
128
,
224314
(
2008
).
34.
R. P.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
(
McGraw-Hill
,
New York
,
1965
).
35.
C. M.
Bender
and
S. A.
Orszag
,
Advanced Mathematical Methods for Scientists and Engineers
(
McGraw-Hill
,
New York
,
1978
).
36.
J. O.
Richardson
, “
Perspective: Ring-polymer instanton theory
,”
J. Chem. Phys.
148
,
200901
(
2018
).
37.
S.
Coleman
, “
The uses of instantons
,” in
Aspects of Symmetry
, edited by
S.
Coleman
(
Cambridge University Press
,
1985
), pp.
265
350
, Chap. 7, Proc. Int. School of Subnuclear Physics (Erice, 1977).
38.
G. V.
Mil’nikov
,
K.
Yagi
,
T.
Taketsugu
,
H.
Nakamura
, and
K.
Hirao
, “
Simple and accurate method to evaluate tunneling splitting in polyatomic molecules
,”
J. Chem. Phys.
120
,
5036
5045
(
2004
).
39.
J. O.
Richardson
and
S. C.
Althorpe
, “
Ring-polymer instanton method for calculating tunneling splittings
,”
J. Chem. Phys.
134
,
054109
(
2011
).
40.
J. O.
Richardson
, “
Ring-polymer instanton theory
,”
Int. Rev. Phys. Chem.
37
,
171
216
(
2018
).
41.
C. L.
Vaillant
and
M. T.
Cvitaš
, “
Rotation-tunneling spectrum of the water dimer from instanton theory
,”
Phys. Chem. Chem. Phys.
20
,
26809
26813
(
2018
).
42.
M. T.
Cvitaš
, “
Quadratic string method for locating instantons in tunneling splitting calculations
,”
J. Chem. Theory Comput.
14
,
1487
1500
(
2018
).
43.
M.
Eraković
and
M. T.
Cvitaš
, “
Tunneling splittings of vibrationally excited states using general instanton paths
,”
J. Chem. Phys.
153
,
134106
(
2020
).
44.
E.
Jahr
,
G.
Laude
, and
J. O.
Richardson
, “
Instanton theory of tunnelling in molecules with asymmetric isotopic substitutions
,”
J. Chem. Phys.
153
,
094101
(
2020
); arXiv:2007.07540 [physics.chem-ph].
45.
J. O.
Richardson
,
S. C.
Althorpe
, and
D. J.
Wales
, “
Instanton calculations of tunneling splittings for water dimer and trimer
,”
J. Chem. Phys.
135
,
124109
(
2011
).
46.
M. T.
Cvitaš
and
J. O.
Richardson
, “
Quantum dynamics in water clusters
,” in
Molecular Spectroscopy and Quantum Dynamics
, edited by
M.
Quack
and
R.
Marquardt
(
Elsevier
,
2021
).
47.
J. O.
Richardson
, “
Full- and reduced-dimensionality instanton calculations of the tunnelling splitting in the formic acid dimer
,”
Phys. Chem. Chem. Phys.
19
,
966
970
(
2017
); arXiv:1611.04816 [physics.chem-ph].
48.
N.
Sahu
,
J. O.
Richardson
, and
R.
Berger
, “
Instanton calculations of tunneling splittings in chiral molecules
,”
J. Comput. Chem.
42
,
210
221
(
2021
).
49.
A.
Nandi
,
G.
Laude
,
S. S.
Khire
,
N. D.
Gurav
,
C.
Qu
,
R.
Conte
,
Q.
Yu
,
S.
Li
,
P. L.
Houston
,
S. R.
Gadre
,
J. O.
Richardson
,
F. A.
Evangelista
, and
J. M.
Bowman
, “
Ring-polymer instanton tunneling splittings of tropolone and isotopomers using a Δ-machine learned CCSD(T) potential. Theory and experiment shake hands
,”
J. Am. Chem. Soc.
145
,
9655
9664
(
2023
).
50.
D.
Lister
,
J.
Macdonald
, and
N.
Owen
,
Internal Rotation and Inversion: An Introduction to Large Amplitude Motions in Molecules
(
Academic Press
,
1978
).
51.
C.
Qu
and
J. M.
Bowman
, “
Full-dimensional, ab initio potential energy surface for CH3OH → CH3+OH
,”
Mol. Phys.
111
,
1964
1971
(
2013
).
52.
S.
Carter
and
N. C.
Handy
, “
The vibrations of H2O2, studied by ‘multimode,’ with a large amplitude motion
,”
J. Chem. Phys.
113
,
987
993
(
2000
).
53.
X.-G.
Wang
and
T.
Carrington
, Jr.
, “
Vibrational energy levels of CH5+
,”
J. Chem. Phys.
129
,
234102
(
2008
).
54.
D.
Lauvergnat
and
A.
Nauts
, “
Smolyak scheme for solving Schrödinger equation: Application to malonaldehyde in full dimensionality
,” chemRxiv:10.26434/chemrxiv-2023-l7npc (
2023
).
55.
Y.
Litman
,
J. O.
Richardson
,
T.
Kumagai
, and
M.
Rossi
, “
Elucidating the quantum dynamics of intramolecular double hydrogen transfer in porphycene
,”
J. Am. Chem. Soc.
141
,
2526
2534
(
2019
); arXiv:1810.05681 [physics.chem-ph].
56.
J. E.
Lawrence
,
J.
Dušek
, and
J. O.
Richardson
, “
Perturbatively corrected ring-polymer instanton theory for accurate tunneling splittings
,”
J. Chem. Phys.
(to be published); arXiv:2304.10963
[physics.chem-ph]
(
2023
).
57.
D.
Chandler
and
P. G.
Wolynes
, “
Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids
,”
J. Chem. Phys.
74
,
4078
4095
(
1981
).
58.
M.
Parrinello
and
A.
Rahman
, “
Study of an F center in molten KCl
,”
J. Chem. Phys.
80
,
860
867
(
1984
).
59.
M. E.
Tuckerman
,
Statistical Mechanics: Theory and Molecular Simulation
(
Oxford University Press
,
Oxford
,
2010
).
60.
D. M.
Ceperley
and
G.
Jacucci
, “
Calculation of exchange frequencies in bcc 3He with the path-integral Monte Carlo method
,”
Phys. Rev. Lett.
58
,
1648
1651
(
1987
).
61.
D. M.
Ceperley
, “
Path integrals in the theory of condensed helium
,”
Rev. Mod. Phys.
67
,
279
(
1995
).
62.
C.
Alexandrou
and
J. W.
Negele
, “
Stochastic calculation of tunneling in systems with many degrees of freedom
,”
Phys. Rev. C
37
,
1513
(
1988
).
63.
M.
Marchi
and
D.
Chandler
, “
Path-integral calculation of the tunnel splitting in aqueous ferrous–ferric electron transfer
,”
J. Chem. Phys.
95
,
889
(
1991
).
64.
E.
Mátyus
,
D. J.
Wales
, and
S. C.
Althorpe
, “
Quantum tunneling splittings from path-integral molecular dynamics
,”
J. Chem. Phys.
144
,
114108
(
2016
).
65.
E.
Mátyus
and
S. C.
Althorpe
, “
Calculating splittings between energy levels of different symmetry using path-integral methods
,”
J. Chem. Phys.
144
,
114109
(
2016
).
66.
C. L.
Vaillant
,
D. J.
Wales
, and
S. C.
Althorpe
, “
Tunneling splittings from path-integral molecular dynamics using a Langevin thermostat
,”
J. Chem. Phys.
148
,
234102
(
2018
); arXiv:1803.04433 [physics.chem-ph].
67.
C. L.
Vaillant
,
D. J.
Wales
, and
S. C.
Althorpe
, “
Tunneling splittings in water clusters from path integral molecular dynamics
,”
J. Phys. Chem. Lett.
10
,
7300
7304
(
2019
).
68.
Y.-C.
Zhu
,
S.
Yang
,
J.-X.
Zeng
,
W.
Fang
,
L.
Jiang
,
D. H.
Zhang
, and
X.-Z.
Li
, “
Torsional tunneling splitting in a water trimer
,”
J. Am. Chem. Soc.
144
,
21356
21362
(
2022
).
69.

By a single “sub-calculation” we mean a complete thermodynamic integration as described in Sec. III C.

70.

A study by Vaillant and Cvitaš41 did perform rotational projection by quadrature, but based on instanton sub-calculations rather than PIMD.

71.
R. P.
Feynman
,
Statistical Mechanics: A Set of Lectures
(
W. A. Benjamin
,
Reading, Mass
,
1972
).
72.
K. F.
Riley
,
M. P.
Hobson
, and
S. J.
Bence
,
Mathematical Methods for Physics and Engineering
, 3rd ed. (
Cambridge University Press
,
2006
).
73.
P. R.
Bunker
and
P.
Jensen
,
Molecular Symmetry and Spectroscopy
, 2nd ed. (
NRC Research Press
,
Ottawa
,
2006
).
74.
H. C.
Longuet-Higgins
, “
The symmetry groups of non-rigid molecules
,”
Mol. Phys.
6
,
445
460
(
1963
).
75.

Alternative expressions in term of Z and Zσv or Zσv and ZC3 can be derived, but ours is the better choice for numerical calculations.

76.

It is also invariant under translation, but the translational contribution is exactly factorizable and can, therefore, be easily removed.73 

77.
R.
Zare
,
Angular Momentum: Understanding Spatial Aspects in Chemistry and Physics
(
Wiley
,
1988
).
78.

As an aside, we note that transitions between rovibronic states with J = 0 are forbidden.73 However, the corresponding tunneling splitting can be extracted from experimental data by fitting to a spectroscopic Hamiltonian, as in Ref. 16.

79.
M.
Suzuki
, “
Decomposition formulas of exponential operators and Lie exponentials with some applications to quantum mechanics and statistical physics
,”
J. Math. Phys.
26
,
601
612
(
1985
).
80.

Here, we are not required to subtract terms of a similar magnitude, and so do not suffer from the fermionic sign problem.

81.
J.
Runeson
,
M.
Nava
, and
M.
Parrinello
, “
Quantum symmetry from enhanced sampling methods
,”
Phys. Rev. Lett.
121
,
140602
(
2018
).
82.
B.
Hirshberg
,
V.
Rizzi
, and
M.
Parrinello
, “
Path integral molecular dynamics for bosons
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
21445
21449
(
2019
).
83.
K. E.
Schmidt
and
M. A.
Lee
, “
High-accuracy Trotter-formula method for path integrals
,”
Phys. Rev. E
51
,
5495
5498
(
1995
).
84.
C.
Eckart
, “
Some studies concerning rotating axes and polyatomic molecules
,”
Phys. Rev.
47
,
552
558
(
1935
).
85.
E.
Bright Wilson
, Jr.
,
J. C.
Decius
, and
P. C.
Cross
,
Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra
(
McGraw-Hill
,
1955
).
86.
G.
Trenins
,
M. J.
Willatt
, and
S. C.
Althorpe
, “
Path-integral dynamics of water using curvilinear centroids
,”
J. Chem. Phys.
151
,
054109
(
2019
).
87.
S. V.
Krasnoshchekov
,
E. V.
Isayeva
, and
N. F.
Stepanov
, “
Determination of the Eckart molecule-fixed frame by use of the apparatus of quaternion algebra
,”
J. Chem. Phys.
140
,
154104
(
2014
).
88.
F.
Jørgensen
, “
Orientation of the Eckart frame in a polyatomic molecule by symmetric orthonormalization
,”
Int. J. Quantum Chem.
14
,
55
63
(
1978
).
89.
H.
Goldstein
,
C.
Poole
, and
J.
Safko
,
Classical Mechanics
, 3rd ed. (
Addison Wesley
,
San Francisco
,
2002
).
90.
W. H.
Press
,
S. A.
Teukolsky
,
W. T.
Vetterling
, and
B. P.
Flannery
,
Numerical Recipes: The Art of Scientific Computing
, 3rd ed. (
Cambridge University Press
,
Cambridge
,
2007
).
91.

Some of the previous tunneling splitting PIMD calculations used a normal-mode propagator64 and the more efficient path-integral Langevin equation thermostat.66 The Eckart spring prevents us from directly employing these in our calculations, but suitable analogous algorithms can probably be developed.

92.
J.
Tennyson
,
M. A.
Kostin
,
P.
Barletta
,
G. J.
Harris
,
O. L.
Polyansky
,
J.
Ramanlal
, and
N. F.
Zobov
, “
DVR3D: A program suite for the calculation of rotation–vibration spectra of triatomic molecules
,”
Comput. Phys. Commun.
163
,
85
116
(
2004
).
93.
S.
Habershon
,
T. E.
Markland
, and
D. E.
Manolopoulos
, “
Competing quantum effects in the dynamics of a flexible water model
,”
J. Chem. Phys.
131
,
024501
(
2009
).
94.

In fact, the lower of the two states is forbidden by nuclear spin symmetry.73 Of course, we can still calculate the corresponding tunneling splitting and compare to the result in Ref. 14, obtained by fitting a rovibrational Hamiltonian to experimental data.

95.
M.
Schnell
, “
Group theory for high-resolution spectroscopy of nonrigid molecules
,” in
Handbook of High-resolution Spectroscopy
, edited by
M.
Quack
and
F.
Merkt
(
John Wiley & Sons, Ltd.
,
2011
), pp.
607
632
.
96.

The calculation is sensitive to the accuracy of the potential energy gradient and Hessian. To achieve satisfactory numerical stability, fourth-order central differences90 were used to calculate first and second derivatives.

97.

Note that the instanton calculation requires a lower temperature and larger number of beads to reach convergence. This is necessary to satisfy the assumptions involved in integrating over the permutational mode and neglecting the tunneling time parameter.39,40,64,107

98.

Both states are allowed by nuclear spin symmetry and have equal spin statistical weights.73 

99.
J.
Sarka
and
A. G.
Császár
, “
Interpretation of the vibrational energy level structure of the astructural molecular ion H5+ and all of its deuterated isotopomers
,”
J. Chem. Phys.
144
,
154309
(
2016
).
100.
X.
Wang
,
S.
Carter
, and
J. M.
Bowman
, “
Pruning the Hamiltonian matrix in MULTIMODE: Test for C2H4 and application to CH3NO2 using a new ab initio potential energy surface
,”
J. Phys. Chem. A
119
,
11632
11640
(
2015
).
101.
R.
Wodraszka
and
U.
Manthe
, “CH5+
: Symmetry and the entangled rovibrational quantum states of a fluxional molecule
,”
J. Phys. Chem. Lett.
6
,
4229
4232
(
2015
).
102.
J. M.
Bowman
,
C.
Qu
,
R.
Conte
,
A.
Nandi
,
P. L.
Houston
, and
Q.
Yu
, “
Δ-machine learned potential energy surfaces and force fields
,”
J. Chem. Theory Comput.
19
,
1
17
(
2023
).
103.
G.
Laude
,
D.
Calderini
,
D. P.
Tew
, and
J. O.
Richardson
, “
Ab initio instanton rate theory made efficient using Gaussian process regression
,”
Faraday Discuss.
212
,
237
258
(
2018
); arXiv:1805.02589 [physics.chem-ph].
104.
S.
Käser
,
J. O.
Richardson
, and
M.
Meuwly
, “
Transfer learning for affordable and high-quality tunneling splittings from instanton calculations
,”
J. Chem. Theory Comput.
18
,
6840
6850
(
2022
); arXiv:2208.01315 [physics.chem-ph].
105.
M.
Ceriotti
,
D. E.
Manolopoulos
, and
M.
Parrinello
, “
Accelerating the convergence of path integral dynamics with a generalized Langevin equation
,”
J. Chem. Phys.
134
,
084104
(
2011
).
106.
V.
Kapil
,
J.
Behler
, and
M.
Ceriotti
, “
High order path integrals made easy
,”
J. Chem. Phys.
145
,
234103
(
2016
).
107.
J. O.
Richardson
and
S. C.
Althorpe
, “
Ring-polymer molecular dynamics rate-theory in the deep-tunneling regime: Connection with semiclassical instanton theory
,”
J. Chem. Phys.
131
,
214106
(
2009
).
108.
Y.-C.
Zhu
,
S.
Yang
,
J.-X.
Zeng
,
W.
Fang
,
L.
Jiang
,
D. H.
Zhang
, and
X.-Z.
Li
Accurate calculation of tunneling splittings in water clusters using path-integral based methods
,”
J. Chem. Phys.
158
,
220901
(
2023
).

Supplementary Material