Exact tunneling splittings from symmetrized path integrals

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.³⁴ 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,³⁵ 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.³⁹ 

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,⁵⁰ for which instanton approaches would fail drastically. The torsional barrier for the methyl group in methanol is not quite as low (around 350 cm⁻¹) but still comparable to the harmonic torsional frequency of ∼290 cm⁻¹.⁵¹ 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,⁵⁶ 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,”⁶⁹ 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.⁷⁰ 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,¹⁰⁸ 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.

In addition to discrete symmetries of the kind described in Sec. II A, a molecular Hamiltonian is invariant under rotation.⁷⁶ 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.

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.⁷⁸ 

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, 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.³⁴ 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).⁵⁹ 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, x₁ and x_N, 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.⁷¹ For example, Eq. (10) evaluated for the totally symmetric irrep of the complete nuclear permutation-inversion (CNPI) group⁷³ 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.⁸⁰ 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⁶¹ 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.

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⁽¹⁾ 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.⁸⁶ 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 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.

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

An isolated water molecule (H₂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,⁹² 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⁹⁰ 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 H_N, the standard ring-polymer Hamiltonian. The corresponding prefactor u₀(r) is 1. For $H̃NP1$⁠, we used $H̃N$ from Eq. (41b), and u₁ 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.

The hydronium ion (H₃O⁺) is one of the simplest molecules to exhibit tunneling splitting, similar to the isoelectronic NH₃ 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%,⁴² making hydronium a suitable test system for our method. To describe the interatomic interactions, we have chosen the PES developed by Yu and Bowman,²¹ 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⁻¹. The lowest vibrational frequency of ∼887 cm⁻¹ 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.⁹⁴ 

The molecular symmetry group of hydronium is isomorphic to D_3h,⁹⁵ 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 $σ̂$ 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) quadrature⁹⁰ 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⁻¹ 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⁻¹ 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).

To put this value into context, we compare it against the instanton tunneling splitting, which we calculated following the protocol in Ref. 40.⁹⁶ The calculation converges to 69.8 cm⁻¹ at $β=8000Eh−1$ (T ≈ 40 K) for N = 1024, overestimating the PIMD value by over 30%.⁹⁷ Our other comparison is to the vibrational configuration interaction (VCI) value of 52.48 cm⁻¹ calculated by Yu and Bowman.²¹ 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⁻¹, which both approaches underestimate by 2–3 cm⁻¹. 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 (CH₃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₁ and a doubly degenerate E level by torsional tunneling through a relatively low potential energy barrier.⁹⁸ 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,⁵¹ fitted to CCSD(T)-F12b/aug-cc-pVDZ energies for configurations relevant to the present study, with a torsional barrier of 350 cm⁻¹ and a harmonic torsional frequency of 286 cm⁻¹. Following Sec. II B, symmetry analysis was conducted in the {E, (123), (132)} subgroup, which is isomorphic to the C₃ 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₃ 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 CH₃ group.

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.