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.

## I. INTRODUCTION

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 scattering^{5,6} or nuclear magnetic resonance,^{7,8} and in the gas phase, where it is extracted from high-resolution microwave^{9–12} and infrared^{2,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” theories^{36–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.

## II. THEORY

### A. Rotationless double-well systems

^{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

*ɛ*

_{1}−

*ɛ*

_{0}, we introduce the inversion operator $\sigma \u0302$, which acts on coordinates as $\sigma \u0302x=\u2212x$ and maps the two degenerate wells onto each other. The ground state is symmetric under inversion, while the first excited state is antisymmetric, meaning $\sigma \u0302|0\u3009=|0\u3009$ and $\sigma \u0302|1\u3009=\u2212|1\u3009$. We can use this property to obtain the low-temperature limit of the “symmetrized” partition function

^{61,71}

*β*→ ∞.

### B. Generalization to finite symmetries

*g*. Let $P\u0302$ be some symmetry operation in this group, for which we define the symmetrized partition function

*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 to

^{72}

^{(n)}. By virtue of Eqs. (6) and (7), and the orthonormality of |

*nl*⟩, we can rewrite Eq. (5) as

*α*of order

*g*

_{α}, they have identical characters.

^{72,73}It is, therefore, sufficient to calculate

*Z*

_{P}=

*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

*m*and

*n*correspond to the same irrep, and 0 otherwise. From Eq. (9) it follows that

*β*→ ∞, 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

*g*

_{α}and $\chi \alpha (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.

*C*

_{3v}point group, which describes the symmetry of a (2

*π*/3)-periodic potential, sketched in Fig. 2, and whose characters are given in Table I.

*C*

_{s}is a subgroup of

*C*

_{3v}describing the feasible symmetries

^{74}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

*C*

_{3v}that is induced by the irrep of a reference state in

*C*

_{s}. For the ground state, we must consider the totally symmetric irrep (

*A*′ in the

*C*

_{s}point group). Its induced representation in

*C*

_{3v}is

*A*

_{1}⊕

*E*, which gives the symmetries of the tunneling-split states, namely, a totally symmetric

*A*

_{1}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

*ɛ*

_{Γ}denotes the lowest energy level that transforms according to the irreducible representation Γ. Equation (13) does not contain any terms in $\epsilon 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\u2212\beta \epsilon A2$ becomes negligible compared to the other terms in Eq. (13) as

*β*→ ∞.

C_{3v}
. | E
. | 2C_{3}
. | 3σ_{v}
. |
---|---|---|---|

A_{1} | 1 | 1 | 1 |

A_{2} | 1 | 1 | −1 |

E | 2 | −1 | 0 |

C_{3v}
. | E
. | 2C_{3}
. | 3σ_{v}
. |
---|---|---|---|

A_{1} | 1 | 1 | 1 |

A_{2} | 1 | 1 | −1 |

E | 2 | −1 | 0 |

*C*

_{3}. The second possibility is to express one of the symmetrized partition functions in terms of the others, e.g., $Z\sigma 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

Determine the symmetries of the tunneling-split states.

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.

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

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

^{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}

### C. Rotating systems

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.

*n*,

*J*,

*M*⟩, where

*n*is shorthand for any other relevant quantum numbers or symmetry labels. These eigenfunctions transform as

**Ω**(which could be a set of Euler angles, or a unit quaternion, or any other valid representation), and $DLM(J)(\Omega )$ is a Wigner

*D*-matrix.

^{77}The latter satisfy an orthogonality relation

*J*. By analogy with Sec. II A, we define a “rotated” partition function

*J*′ + 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)(\Omega )=1$, Eq. (18) can be rearranged to read

*J*= 0 can be isolated by taking an unweighted average over all “rotated” partition functions.

## III. METHODS

### A. Path-integral molecular dynamics

*Z*

_{P}. First, we apply the symmetric Trotter product formula,

^{79}

*β*

_{N}=

*β*/

*N*, $V\u0302$ is the potential energy operator, and $T\u0302$ is the kinetic energy operator. Denoting the bracketed product of operators as $O\u0302N$, in the limit as

*N*→ ∞,

*N*− 1 position resolutions of the identity. The bold cursive

**refers to the set of**

*x**N*vectors {

**x**

_{i}|

*i*= 1…

*N*}, where $xi=(x1(i),x2(i),\u2026,xf(i))$, and

*f*is the number of configurational degrees of freedom (DoF). For the integration, we use the shorthand notation

*m*

_{k}is the mass of the

*k*-th DoF, and

**refers to the set of momenta {**

*p***p**

_{i}|

*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\u0302$ is taken to be the identity $E\u0302$, where $E\u0302x=x$, the additional spring connects the extended classical system into a “ring polymer.” In this case, *Z*_{E} 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\u0302$, 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, **x**_{1} and **x**_{N}, tend to explore degenerate minima related by the symmetry operation $P\u0302$, 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) group^{73} 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 PIMC^{61} and PIMD methods^{81,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.

### B. Eckart spring

*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:

*n*

_{at}is the number of atoms in the molecule,

*f*= 3

*n*

_{at}, and $ra(i)$ is a three-dimensional coordinate vector $xa(i),ya(i),za(i)$ describing the

*a*-th atom of mass

*m*

_{a}in the

*i*-th system replica (bead). The integral over

**Ω**can be expressed as

**Ω**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 Hamiltonian^{84,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.

**q**= (

*q*

_{0},

*q*

_{1},

*q*

_{2},

*q*

_{3}) with a real part

*q*

_{0}. This is related to the rotation matrix via

**q**

**C**denoting a 4 × 4 real symmetric matrix that is a simple function of

*m*

_{a}, $ra(1)$, and $ra(N)$ (

*a*= 1, …,

*n*

_{at}). The explicit form of

**C**is given in Eq. (24) of Ref. 87 where, in our case,

*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.

*f*(

**Ω**) is rewritten as

**I**is the 3 × 3 identity matrix. Next we express the rotation

**R**(

*δ*

**Ω**) about the minimum in terms of Tait–Bryan angles $\varphi ,\theta ,\chi $,

^{89}

*c*

_{ϕ}= cos

*ϕ*,

*s*

_{ϕ}= sin

*ϕ*, etc. The corresponding integral over orientations is

*z*= cos

*θ*. Now we expand

**I**−

**R**(

*δ*

**Ω**) in a Taylor series about $\varphi ,z,\chi =0$. Up to second order,

*f*(

**Ω**). This evaluates to zero because $r\u0303a(1)$ satisfies the rotational Eckart conditions.

^{85}Hence, only the second- (and higher-) order terms make a non-zero contribution around the stationary point

*δ*

**Ω**stands for $\varphi ,z,\chi $, and

**Θ**=

**Θ**[

**r**

^{(N)},

**r**

^{(1)}] is a 3 × 3 matrix with elements

*N*→ ∞),

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.

*N*with an “Eckart spring,” as defined by the sum in Eq. (41b). This modified spring does not specifically favor

**r**

_{1}and

**r**

_{N}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

*Z*

^{J=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:

*different*potential energy wells, with minima related by the operation $P\u0302$.

### C. Thermodynamic integration

*u*

_{i}(

**) refers to Eq. (41c) evaluated with $\Theta [P\u0302i\u22121r(N),r(1)]$, and**

*r**u*for a system described by the Hamiltonian $H\u0303NP$. 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

*ξ*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) quadrature

^{90}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\u0303\xi $. Most of the terms in the average cancel exactly; the only remaining contribution is the difference of potential energies in $H\u0303NP1$ and $H\u0303NP0$ 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\u2212\beta N\Delta 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 *∑*_{i}*V*(**r**_{i}). 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.

## IV. RESULTS

### A. Water

An isolated water molecule (H_{2}O) 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 *Z*^{J=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 grid^{90} 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\u0303NP0$ with *H*_{N}, the standard ring-polymer Hamiltonian. The corresponding prefactor *u*_{0}(** r**) is 1. For $H\u0303NP1$, we used $H\u0303N$ from Eq. (41b), and

*u*

_{1}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

*t*

_{th}= 5 ps under

*H*

_{N}, 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

*t*

_{run}= 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.

T (K)
. | PIMD {N = 48}
. | PIMD {N = 64}
. | . | ||
---|---|---|---|---|---|

n_{ξ} = 17
. | n_{ξ} = 33
. | n_{ξ} = 17
. | n_{ξ} = 33
. | Exact . | |

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_{ξ} = 17
. | n_{ξ} = 33
. | n_{ξ} = 17
. | n_{ξ} = 33
. | Exact . | |

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 |

### B. Hydronium

The hydronium ion (H_{3}O^{+}) is one of the simplest molecules to exhibit tunneling splitting, similar to the isoelectronic NH_{3} 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 *D*_{3h},^{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 *C*_{s}), 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 $\sigma \u0302$ is replaced with $P\u030212$, the (12) permutation. TI was performed using $H\u0303N$ [Eq. (41b)] and $H\u0303N(12)$ [Eq. (43) with $P\u0302\u22121=P\u030212$] for the $H\u0303NP0$ and $H\u0303NP1$ Hamiltonians, respectively. The corresponding prefactor functions were $u0=det\Theta [r(N),r(1)]\u22121/2$ and $u1=det\Theta [P\u030212r(N),r(1)]\u22121/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) quadrature^{90} was used instead of CC, as it was found to converge more rapidly with grid size.

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).

T (K)
. | . | Δ (cm^{−1})
. | ||
---|---|---|---|---|

N
. | n_{ξ} = 15
. | n_{ξ} = 20
. | n_{ξ} = 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 | VCI^{21} | PIMDb | Expt.^{14} |

Δ (cm^{−1}) | 69.8 | 52.48 | 53.2(4) | 55.35 |

T (K)
. | . | Δ (cm^{−1})
. | ||
---|---|---|---|---|

N
. | n_{ξ} = 15
. | n_{ξ} = 20
. | n_{ξ} = 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 | VCI^{21} | PIMDb | Expt.^{14} |

Δ (cm^{−1}) | 69.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 $\beta =8000Eh\u22121$ (*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.

### C. Methanol

Methanol (CH_{3}OH) is a multi-well system with a MS group that is isomorphic to *C*_{3v} (see Sec. II B and Fig. 2). Its ground state is split into an *A*_{1} 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 approaches^{24,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 *C*_{3} 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 *C*_{3} 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\u0303N(132)$ [Eq. (43) with $P\u0302\u22121=P\u0302123$]. 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 CH_{3} group.

Thermodynamic integration for this system was performed between $H\u0303N$ and $H\u0303N(132)$, with prefactor functions $u0=det\Theta [r(N),r(1)]\u22121/2$ and $u1=det\Theta [P\u0302123r(N),r(1)]\u22121/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.

T (K)
. | . | Δ (cm^{−1})
. | ||
---|---|---|---|---|

N
. | n_{ξ} = 5
. | n_{ξ} = 10
. | n_{ξ} = 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^{−1}) | 11.3 | 9.15 | 9.14(13) | 9.12 |

T (K)
. | . | Δ (cm^{−1})
. | ||
---|---|---|---|---|

N
. | n_{ξ} = 5
. | n_{ξ} = 10
. | n_{ξ} = 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^{−1}) | 11.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.

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 calculation^{40,96} converged to 11.3 cm^{−1} at $\beta =20000Eh\u22121$ (*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 splitting^{16} 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.

## V. CONCLUSIONS

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) calculation^{56} 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 thermostats^{66,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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX: ERROR PROPAGATION

*σ*

_{u}, of the mean of $\u27e8u(r)\u27e9P\u0302$ [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 $\Phi =ZP1J=0/ZP0J=0$ according to

## REFERENCES

*Chemical Dynamics at Low Temperatures*

*Adv. Chem. Phys.*

*The Theory of Intermolecular Forces*

*γ*-Picoline and toluene

_{2}O)

_{2}

_{3}O

^{+}

_{3}O

^{+}: The pure inversion spectrum around 55 cm

^{−1}

_{3}O

^{+}→ H

^{+}+ H

_{2}O: A step to a many-body representation of the hydrated proton?

_{3}O

^{+})

_{2}on a six-dimensional potential hypersurface

*ab initio*potential energy surface

*Quantum Mechanics and Path Integrals*

*Advanced Mathematical Methods for Scientists and Engineers*

*Aspects of Symmetry*

*Molecular Spectroscopy and Quantum Dynamics*

*Internal Rotation and Inversion: An Introduction to Large Amplitude Motions in Molecules*

_{3}OH → CH

_{3}+OH

_{2}O

_{2}, studied by ‘multimode,’ with a large amplitude motion

_{5}

^{+}

*Statistical Mechanics: Theory and Molecular Simulation*

^{3}He with the path-integral Monte Carlo method

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

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

*Statistical Mechanics: A Set of Lectures*

*Mathematical Methods for Physics and Engineering*

*Molecular Symmetry and Spectroscopy*

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

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

*Angular Momentum: Understanding Spatial Aspects in Chemistry and Physics*

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

*Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra*

*Classical Mechanics*

*Numerical Recipes: The Art of Scientific Computing*

Some of the previous tunneling splitting PIMD calculations used a normal-mode propagator^{64} 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.

*Handbook of High-resolution Spectroscopy*

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

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}

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

_{2}H

_{4}and application to CH

_{3}NO

_{2}using a new ab initio potential energy surface

*Ab initio*instanton rate theory made efficient using Gaussian process regression