Tunneling splittings in the energetically low-lying structural isomers of the water hexamer: the prism, the cage and the book

Abstract

Tunneling splitting (TS) patterns of the energetically low-lying structural isomers of the water hexamer are calculated using the modified WKB (Wentzel–Kramers–Brillouin) method in full dimensionality. TSs in the water hexamer prism are determined for a number of excited low-frequency vibrational modes. Internal rotation of a double-donor water monomer is identified as the mechanism that potentially plays a role in the appearance of the TS pattern in vibrationally excited states in addition to the mechanisms that shape the TS pattern in the ground state. The ground-state TSs of the water hexamer cage were found to form a doublet of doublets. The finer splitting is two orders of magnitude smaller due to a stark difference in the barrier heights for bifurcations of the water monomers at the two opposite vertices of the cage. We also give the first estimates of the ground-state TSs in the water hexamer book structure. The TS pattern is again a doublet of doublets, with the wider doublet of similar size to that in the cage and the narrower doublets an order of magnitude larger than that in the cage. The case study of the cage and the book represents the first realistic application of instanton theories to systems with symmetrically inequivalent wells.

---

1 Introduction

Understanding water from first principles requires a detailed knowledge of its molecular interactions.^1,2 These interactions can be analysed on a fundamental level by studying the properties of size-selected water clusters.^3–5 Rearrangements of water clusters between the symmetry-related, equivalent minima produce measurable energy splittings of the otherwise degenerate vibrational states. These spectroscopic signatures of tunneling at zero temperature are highly sensitive to the molecular interactions along the rearrangement pathways. They can vary over many orders of magnitude depending on the cluster structure, the mechanism of the rearrangement and the excitations of the vibrational modes of the cluster.

Tunneling splittings (TSs) of vibrational states in water dimer have been measured in experiments^6–9 and found to vary over four orders of magnitude depending on the rearrangement mechanism.¹⁰ More recent experiments^11,12 found that the vibrational excitation of a librational mode magnifies the ground-state (GS) splitting pattern up to 40 times. Water trimer, tetramer and pentamer form hydrogen-bonded ring structures in their minimum-energy conformations. Each water monomer acts as a hydrogen-bond donor and acceptor with a ‘free’ hydrogen on each monomer pointing either above or below the ring plane. The GS splitting pattern in water trimer has been measured experimentally^13–15 and rationalized in terms of six different rearrangement mechanisms.^16–18 The mechanisms include torsional motion (or flip of ‘free’ hydrogen), whereby a water monomer rotates about its in-ring hydrogen bond, while the ‘free’ hydrogen moves from above to below the ring plane. Torsional states are further split by bifurcations, in which the in-ring hydrogen bond is broken and the monomer rotates about the in-plane axis passing through its oxygen such that the ‘free’ hydrogen rotates into the ring plane to replace the in-bond hydrogen and reforms the hydrogen bond. Bifurcations were found to be accompanied by one or more flips of ‘free’ hydrogens on other monomers (in the remaining five mechanisms). Similar mechanisms are found to be responsible for the TS pattern formation in tetramer^19–21 and pentamer^22–26 (using five dominant rearrangement pathways). Vibrational excitation of a librational mode in water trimer and pentamer was also found to have a pronounced effect on the sizes of the splittings, resulting in 400-fold^12,27 and 4000-fold²⁸ increases, respectively.

Larger water clusters form three-dimensional (3D) structures. The number of these structures is vast²⁹ with many of them lying close in energy. Strong directionality of hydrogen bonds results in many competing effects affecting their energetics, which places high demands on the accuracy of the potential energy surface (PES) for a reliable prediction of stable structures. In a water hexamer, a number of low-energy isomeric structures have been located computationally.^30,31 A planar ring structure, similar in nature to those found in the smaller clusters, has been detected in experiments.³² However, 3D structures, the so-called cage, prism and book isomers, which were simultaneously detected in experiments,³³ have lower energies. Whether the hexamer cage or the prism isomer lies lowest in energy has been a subject of debate.^34,35 It has now been established through experiment^33,36 and theoretical analysis³¹ that the cage lies the lowest. It is only after the inclusion of the nuclear zero-point vibrational motion that the cage structure becomes more favourable than the prism.³¹ Interestingly, the energetic order reverses upon deuteration and the prism isomer becomes more stable.³⁵ With growing temperature, due to entropic effects, the book structure becomes more likely than either the cage or the prism.³⁴

The TS pattern was measured for the hexamer cage and the hexamer prism. It was found that the TS pattern in the hexamer cage in the excited state at 83 cm⁻¹ forms a triplet.^36,37 This was rationalized^37,38 in terms of two rearrangement mechanisms and the assumption of an accidental degeneracy of a pair of states. The TS pattern in the GS of the hexamer prism was found to form a doublet of triplets.³⁹ Computational work using instanton theory identified two rearrangement mechanisms responsible for the formation of the spectral pattern including one which involves a simultaneous breaking of two hydrogen bonds.³⁹ In the vibrationally excited librational modes, large enhancements of TSs, of three orders of magnitude, were found for both the hexamer cage⁴⁰ and the prism,⁴¹ with the accidental degeneracy in the cage removed. Attempts to detect TSs in the water heptamer,⁴² octamer⁴³ and nonamer⁴⁴ did not yield evidence of tunneling. Instanton calculations confirmed that the splittings in the octamer are below the detection limit.⁴³ TSs in the water decamers have been observed for two different isomers⁴⁴ and the widths were found to qualitatively match the calculations using instanton theory.⁴⁵

A water dimer is the only water cluster that has been solved using exact quantum mechanics^46,47 in full dimensionality. Torsional rearrangements in water trimer^48,49 and OH flipping vibrations in the water hexamer cage⁵⁰ have been studied using variational methods in 3D and 2D reduced-dimensional models. Recently, vibrational states of the water trimer, including the tunneling splittings due to bifurcations, were determined in nine dimensions,⁵¹ with all angular degrees of freedom of rigid water monomers treated explicitly, in the calculations that reached good agreement with experiment. The GS splitting patterns can also be determined using diffusion Monte Carlo (DMC) methods in reduced dimensionality.^52–54 A full-dimensional treatment of the fine splittings due to bifurcations has recently been accomplished using a path-integral molecular dynamics (PIMD) approach^55,56 for the water dimer, trimer and hexamer. Other recent studies of the GS splitting patterns in water clusters, from dimers to decamers, have employed semiclassical approaches, instanton theory^18,39,43,45 or the modified WKB^26,57,58 methods. Both PIMD and the semiclassical approaches determine tunneling matrix elements in the so-called tunneling matrix approach, developed in the earlier studies of water clusters.^14,16,17,38 Tunneling matrix elements are calculated in separate calculations in the limited regions of configuration space.

Modified WKB (M-WKB)^59,60 and instanton methods⁶¹ for calculating TSs are closely related. Both use the minimum action path connecting two minima on the PES and quadratic expansions of the potential around it in order to evaluate the energy shifts in full dimensionality. We modified the original formulation^59,60 of the M-WKB method to treat asymmetric paths (with an asymmetry in the potential energy profile along the path),^62,63 which regularly arise in systems with multiple minima, such as the water clusters,⁶⁴ and showed that it is equivalent to the instanton theory in the vibrational GS.⁶² The M-WKB method has recently been tested on vinyl radical⁶⁵ and it was shown that it gives the TSs in quantitative agreement, within a factor of 2–3, with the exact quantum results. In the vibrationally excited states, it reproduces the enhancements of TSs of three orders in magnitude in all cases where they occur. This method was also applied to calculate TSs in the low-lying vibrational states of water trimer,⁵⁸ including the fine structure due to bifurcations. We extended the theory⁶⁶ to calculate TSs in asymmetric systems, having minima that differ in energy or the shape of the potential well. This allowed us to estimate the TS patterns in the partially deuterated water trimers⁵⁷ for the states that lie in the symmetry-related minima but possess an energy asymmetry due to a difference in their zero-point energies (ZPEs).

The objective of the present study is to apply the M-WKB method to study the TS patterns of the low-lying isomers of water hexamers. The GS TS pattern of the hexamer prism has already been studied using instanton theory and the relevant rearrangement mechanisms determined that rationalize the experiment.³⁹ A PIMD⁵⁵ treatment of the hexamer prism showed that the instanton results lie within a factor of two relative to those obtained in a formally exact approach. M-WKB allows us to calculate the TS patterns in the low-lying vibrationally excited states at a similar level of accuracy and to study the effect of modal excitations on the sizes of the splittings below. We also apply the M-WKB method to study the GS TS pattern in the hexamer cage and book isomers. These isomers possess nearly degenerate minima connected by the low-barrier OH flips,⁶⁷ which we treat below using the recent modification of the M-WKB method⁶⁶ to provide numerical estimates of the sizes of the splittings. Our calculations below present the first application of the M-WKB method to a realistic system with inequivalent wells having both the energy and the shape asymmetry, in full dimensionality. A more accurate PIMD approach⁶⁸ cannot be applied to any of the calculations below as it is, at present, limited to the GS and the degenerate potential wells. Water hexamer is a 48-dimensional system and the variational methods could only be applied to it in reduced-dimensional models. The choice of the appropriate coordinates to describe the dynamics accurately then is not straight-forward, especially so in the vibrationally excited states, which involve motion of the whole cluster.

All the above methods for treating the rearrangement dynamics in water clusters rely on an accurate PES. They are thus made possible by the recent developments in constructing accurate water potentials. The commonly used water potentials are the so-called CC-pol,⁶⁹ WHBB,^3,70,71 MB-pol^72–74 and, more recently, a new version of MB-pol⁷⁵ and q-AQUA,⁷⁶ a fully ab initio potential. They all employ a many-body expansion of the monomer potentials, up to 4-body terms, that are fits to high-accuracy ab initio electronic energies. We use MB-pol^72–74 potential below.

The paper is organized as follows. Section 2 summarizes the M-WKB theory for calculating tunneling splittings and the computational methods that we use to evaluate them. In Section 3, we summarize previous work on the GS TS in water hexamer prism and proceed to study the effect of vibrational excitation on TSs. This is followed by the study of TSs due to bifurcations in the GS of the hexamer cage, where we give the first theoretical estimates of their size. We then identify the relevant rearrangements for the appearance of the TS pattern in the hexamer book and estimate the sizes of the splittings in its GS. Section 4 concludes the presented results.

2 Methodology

2.1 Theoretical

In molecular systems with multiple minima, the sets of local vibrational states of each minimum interact by tunneling via feasible rearrangements. The tunneling interaction causes energy shifts in their vibrational spectra. The local vibrational states can be used as a basis set to construct the nuclear Hamiltonian. If tunneling interaction is small compared to the energy differences of the single-well localized states, it is sufficient to include one such state per minimum in the basis set. Hamiltonian is then represented by a square matrix of the dimension equal to the number of minima that are accessible via tunneling and is referred to as the tunneling matrix (TM). Localized state energies of each well lie on its diagonal and pairwise tunneling interactions are off-diagonal elements connecting the corresponding pair of minima that interact via tunneling in the overlap region inside the barrier.

Tunneling interaction between the localized states of two wells, i and j, can be evaluated using the Herring formula,^77,78

where ϕ^(i/j) are the localized wavefunctions in wells i and j, Σ is the dividing plane, placed inside the barrier which separates the two wells, with the unit normal n, oriented in the i → j direction. Eqn (1) assumes the use of Cartesian mass-scaled coordinates and atomic units (ħ = 1). We note that the Herring formula is valid even in the case where the wells i and j are inequivalent and the localized states are not degenerate.⁶⁶

Localized wavefunctions ϕ^(i/j) are obtained using the M-WKB approach of ref. 62 and 63. The wavefunctions are constructed along the characteristics of the M-WKB equations, which follow the classical imaginary time trajectories on the inverted potential that start out at a minimum with zero total energy and run towards the dividing plane Σ. The classical path is a continuous sequence of molecular geometries x(S) in the f-dimensional mass-scaled coordinate space, parametrized by the arc-length distance S from the minimum. The f-dimensional momentum is

where τ is the imaginary time, t is the tangent to the path at S, V(S) is the potential at x(S) and V_min = V(0) is the potential at minimum. The action integral from the minimum at S = 0 to the dividing surface at S = S_Σ is thenwhere p(S) = dS/dτ is the magnitude of the momentum. The two paths from minima i and j are connected at Σ in such a way that the sum of their associated actions in eqn (3), A⁽ⁱ⁾ + A^(j), is minimal. This path produces the largest amplitude of the wavefunction at Σ in the M-WKB approach or, in path-integral formalism, gives the dominant contribution to the double-well partition function, the instanton. Instanton theory^61,79 and M-WKB^59,62 lead to the equivalent expression for the GS TS. In the case of inequivalent wells, the MAP has a tangent discontinuity at S_Σ due to the different V_min, in eqn (2), in wells i and j.

The M-WKB wavefunction for the GS (ν = 0) and the first excited state (ν = 1) of the normal mode U_e with frequency ω_e at a general geometry x then reads,⁶³

where Δx = x − x(S) is a displacement from the reference point on the path x(S) and we drop the indices (i/j) that label minima. The wavefunction in eqn (4) at minimum has the form of the harmonic oscillator wavefunction. The matching to harmonic oscillator determines the first two factors on the right-hand side (r.h.s.) in eqn (4) as normalization constants for the GS and the first excited state, the normal mode vector U(0) as U_e and the f × f matrix A(0) = A₀ as the square root of hessian H₀ at minimum. The evolution of F(S), U(S) and A(S) along the path is determined by the following equations,⁶³

Function F can also be determined⁶³ using F(S) = Up/ω_e. A(S) in eqn (5) describes the Gaussian shape of the wavefunction in eqn (4). Its evolution is governed by H(S), Hessian of the potential at x(S). At each reference point on the path x(S), the wavefunction in eqn (4) represents the M-WKB solution of the Schrödinger equation in the full-dimensional space. The underlying potential is a quadratic expansion of the potential at x(S) on the path. In the third factor on the r.h.s. of eqn (4) (for ν = 1), F(S) describes the shift of the nodal plane away from the minimum action path (MAP), while U(S) traces the direction of the nodal plane along the MAP. The exponential factors that follow describe the change of the amplitude due to the effect of the change of ZPEs of normal modes along the MAP and the Gaussian shape of the wavefunction in all space. The wavefunction amplitude decays as exp(−A) along S.

TSs are evaluated by substituting the wavefunctions ϕ⁽ⁱ⁾ and ϕ^(j), for the wells i and j, of the form given in eqn (4), into the Herring formula, eqn (1). Dividing surface Σ is set to be a plane normal to n. n is fixed along the (p⁽ⁱ⁾ − p^(j)) vector evaluated at S_Σ. For inequivalent wells, there is a MAP tangent discontinuity at S_Σ and momenta p⁽ⁱ⁾ and p^(j) have different magnitudes and directions at S = S_Σ. n is parallel to the weighted mean of the MAP tangents at the two sides of Σ. We note that for the equivalent wells p⁽ⁱ⁾ = −p^(j) and n = t(S_Σ). Directional derivatives, (n∇), in eqn (1) bring down a constant factor −|p⁽ⁱ⁾(S_Σ) − p^(j)(S_Σ)|/2 as the leading order contribution in front of the product ϕ⁽ⁱ⁾ϕ^(j) in the integrand. The product of Gaussians is then integrated over the (f − 1)-dimensional Σ analytically, one dimension at a time, in the diagonal representation of . The symbol ⊥ means that the direction n has explicitly been projected out from Ā. The final formula for the TM element h_ij is given by

where all quantities outside integrals are evaluated at S = S_Σ, p̄ = (p⁽ⁱ⁾ + p^(j))/2, F̃⁽ⁱ⁾ = F⁽ⁱ⁾ − U^(i)TĀ⁻¹p̄, and det′ denotes the product of all non-zero eigenvalues. Zero eigenvalues are associated with the overall translations and rotations, while Ā has an additional one associated with n. Ā⁻¹ in eqn (7) is therefore a pseudoinverse of Ā. In the case of equivalent minima, p̄ = 0. In eqn (7), the number of vibrational quanta in states (i/j) are ν^(i/j) = 0 or 1.

2.2 Computational

Numerical evaluation of the TM element h_ij starts with the locating of the MAP that connects minima i and j. The molecular geometries at the two minima are aligned by minimizing the Euclidean distance between them in mass-scaled Cartesian coordinates using a quaternion-based algorithm.⁸⁰ The initial path is discretized using N system replicae with the geometries given by a linear interpolation at equally-spaced distances S between two minimum geometries. If the linear path passes through physically inadmissible geometries, as it does for water monomer rotations below, we create an intermediate geometry explicitly, e.g., by rotating one water monomer around its C₂ axis by π/2, and perform a linear interpolation at equidistant S values using a ‘three-point’ path. The MAP is then determined iteratively using string method^81,82 that requires the evaluation of potential gradients at all discretization points at each iteration. End points of the string are fixed at minima and their orientation, as well as all intermediate geometries are adjusted to minimize the action integral. The convergence criterion is set to 10⁻⁶ a.u. for the largest magnitude of the perpendicular-to-path Nf-dimensional action gradient at a discretization point along the string. MAP is determined by progressively increasing the number of discretization points at subsequent optimizations. Up to N = 1000 is used in this work. For the water hexamer cage and book, where MAP connects minima at different energies, we first set both V_min's in eqn (3) (for the wells i and j) to the lower energy of the two minima to determine a MAP. We then use this MAP as the initial path in the optimization in which we minimize the sum of actions in eqn (3) evaluated on the two sides relative to the connection point S_Σ which, initially, is chosen at the discretization point with the highest potential V. V_min's are thereby set to different values in the two action integrals. The definition of the path tangent at the connection point, used for the action gradient projections in the string method,^81,82 is adjusted to the difference in the momenta on the two sides of the connection point.

After locating the MAP, cubic spline interpolants for molecular geometries x(S), potential V(S) and Hessians H(S) are constructed, element-by-element, from their values at the discretization points. For MAPs connecting inequivalent minima, two interpolants are constructed on the two sides of S_Σ. These are used to solve for A^(i/j)(S) by solving eqn (5). Cubic spline interpolants of A⁽ⁱ⁾(S) and A^(j)(S) are then constructed and used to solve eqn (6) for each excited mode U_e of interest. Due to p(0) = 0, there is a singularity at S = 0 in eqn (5) and (6) and the solution at S = ε is obtained using a polynomial expansion⁶³ in H, A and U. We varied ε in the interval from 0.1m_e^1/2a₀ to 10m_e^1/2a₀ in our calculations below. TM elements are next calculated using eqn (7). Dependence of TM elements on N and ε is used to assess their convergence.

Finally, TM elements for all relevant rearrangements are inserted into a TM, which is then diagonalized to give tunneling splittings. This is performed separately for each vibrationally excited mode of interest. The symmetry of TM eigenvectors determines the symmetry of vibrational states, which allows one to determine the nuclear spin degeneracies and thus the allowed transitions and their intensity pattern. In a system with multiple symmetry-related minima, diagonal state energies are left out of the TM and the eigenvalues of TM give the relative energy shifts due to tunneling. In systems with inequivalent minima, such as in the hexamer cage and book, we use harmonic energies of the localized states on the diagonal of TM. These are considerably less accurate than the M-WKB TM elements. We thus concentrate only on the splittings due to bifurcations in the lowest torsional branch of the GS of the hexamer cage and book below.

3 Results

3.1 Prism

3.1.1 Tunneling pathways in D_2d. The lowest-energy isomer of the water hexamer is the prism shown in its equilibrium geometry in Fig. 1. The labeling of atoms is introduced to distinguish between different versions, the permutational isomers that differ only in the labeling of atoms. Rearrangements that interconvert versions and produce tunneling splittings of vibrational states are termed feasible. From a large number of possible rearrangement processes, those that break covalent bonds and those that break many of its hydrogen bonds simultaneously result in negligible splittings and can be excluded from consideration. Ref. 39, in a joint experimental and theoretical study, identified two processes that are responsible for the formation of the GS splitting pattern, a doublet of triplets, in the water hexamer prism. By means of isotopic substitutions of ¹⁶O by ¹⁸O, it was found that monomers A and D, labeled by their oxygen atoms, must be involved in the relevant rearrangements. Isotopic substitutions of A and D, B and F, or C and E, all preserve the splitting pattern. The substitution of the heavier isotopologues in the place of A and D resulted in a reduced splitting. Under these constraints, instanton theory was used to estimate the sizes of the splittings for the viable rearrangements to find those that are responsible.

Fig. 1 Minimum-energy geometry of the water hexamer prism labeled in its reference version.

The mechanism responsible for the formation of the doublet splitting in the hexamer prism is the antigeared double flip³⁹ of monomers A and D. In a simultaneous action, the hydrogen atoms 1 and 7 on monomers A and D move from below the basis of the prism, formed by oxygens ABC and DEF, to above, by rotating around the bonds A-2 and D-8, respectively, in opposite directions, as shown in Fig. 2(a). The mechanism is labeled AD by the symbols of oxygen atoms of the monomers involved in the motion. The motion of each monomer is similar to the single monomer flip in water trimer, but, contrary to it, the process breaks and reforms a hydrogen bond. It connects the reference version in Fig. 1 to the version labeled by and related to it by the symmetry operation (A D)(B F)(C E)(1 7)(2 8)(3 11)(4 12)(5 9)(6 10). The operation is its own inverse and the two minima form a symmetric double-well system. The potential profile along the MAP is shown in Fig. 3. The MAP length is 294.3m_e^1/2a₀, which is four times shorter than the minimum-energy path (MEP),³⁹ due to a large corner-cutting effect. The barrier height along the MAP is 1602 cm⁻¹, while it is 1041 cm⁻¹ at the transition state. The barrier from the vibrational ground state, corrected for the harmonic ZPE of orthogonal modes, is 636 cm⁻¹.

Fig. 2 Minimum action paths of the four rearrangements in the water hexamer prism that are responsible for the formation of its tunneling splitting pattern: (a) AD, (b) ÃD/D̃A, (c) rot C, and (d) rot E (see text for notation).

Fig. 3 Potential energy curves along the minimum action paths for different rearrangement mechanisms in the water hexamer prism. See the text for notation.

The mechanism that is responsible for the triplet splitting of each doublet branch is a geared bifurcation of monomer A and a flip of monomer D,³⁹ shown in Fig. 2(b). The mechanism is labeled ÃD, whereby the atom capped by a tilde refers to the monomer undergoing bifurcation, a similar dynamics to the bifurcation mechanism in water trimer.¹⁸ Monomer A rotates about the axis that is perpendicular to the monomer plane and passing through the oxygen A, while monomer D executes similar motion to that in the antigeared double flip AD. The process breaks two hydrogen bonds and produces a smaller splitting than the antigeared double flip AD. The symmetry operation that links the two minima connected by this rearrangement is (A D)(B F)(C E)(1 7 2 8)(3 11)(4 12)(5 9)(6 10). Four consecutive applications of the symmetry operation lead back to the reference version. The four state energies form a triplet in each doublet branch, with the mid level doubly degenerate. The reverse motion is associated with D̃A and the symmetry operation (A D)(B F)(C E)(1 8 2 7)(3 11)(4 12)(5 9)(6 10). The MAP length is 320.7m_e^1/2a₀, which is, according to ref. 39, two times shorter than the MEP for this rearrangement. The potential along the MAP, shown in Fig. 3, is asymmetric (with a barrier top at 52% of its length). The MAP barrier height is 1550 cm⁻¹, lower than that for the double flip AD, with a transition state at 1339 cm⁻¹. The corrected barrier height is, on the other hand, higher, at 952 cm⁻¹. We also note that the MAP length is 9% longer than that of the double flip AD.

Combined together, the two symmetry operations (for AD and ÃD) generate a group of order 8 that is isomorphic to the D_2d point group.³⁹ The remaining symmetry elements (apart from the identity and D̃A) are the rotation of monomer A, denoted ‘rot A’ and associated with the element (1 2), the rotation of monomer D or ‘rot D’, associated with (7 8), the double bifurcation ÃD̃, associated with (A D)(B F)(C E)(1 8)(2 7)(3 11)(4 12)(5 9)(6 10), as well as the element (1 2)(7 8). We determined the MAPs for these rearrangements and found that the last one decomposes into a two-step process ÃD + D̃A (or vice versa), achieved by applying the symmetry element for ÃD twice. Potential curves along MAPs for the remaining rearrangements are shown in Fig. 3. They result in higher barriers and longer paths and, therefore, larger action in eqn (3). Their contribution is neglected, as in ref. 39.

3.1.2 Ground-state tunneling splittings in D_2d. TM in the D_2d group is an 8 × 8 matrix which has three non-zero TM elements in each of its rows/columns (associated with AD, ÃD and D̃A), with two of them equal in size (see Fig. 4B in ref. 39). The resulting TS pattern is shown in Fig. 4. The left-most set of energy levels displays the doublet splitting due to the antigeared double flip AD of width 2h_a, where h_a is the associated TM element (adopting notation from ref. 39). In the mid set, the geared ÃD/D̃A mechanism is introduced with its TM element h_g. Each doublet branch resolves into a triplet with relative energies 2h_g, 0 and −2h_g. TM eigenvectors allow one to determine the symmetries of vibrational states in D_2d, given in Fig. 4, and the nuclear-spin statistical weights of the sextet states, given in ref. 39 and, for the fully deuterated hexamer-d₁₂, in ref. 41.

Fig. 4 Tunneling splitting pattern of the vibrational ground state of the water hexamer prism on MB-pol PES.^72–74 The vertical sets of levels, from the left to right, are the energy splittings due to the mechanisms: AD, AD + ÃD/D̃A, and AD + ÃD/D̃A + rot E + rot C, respectively. See the text for notation.

TM elements h_a and h_g in the GS have been determined using instanton theory and experiment in ref. 39 and, later, using formally exact PIMD in ref. 55. They are compared on the MB-pol PES in Table 1 of ref. 55. Our M-WKB calculation of the TM element for the antigeared AD flip gives h_a = −0.74 MHz, which is in agreement with −0.75 MHz of ref. 39 and it overestimates the PIMD result by a factor of 1.85. The PIMD result lies within 5% of the experimental value.³⁹ This likely means that the error in our calculation originates mainly from the anharmonicity of the potential in directions orthogonal to the MAP, and not from the approximations introduced by using the TM model. For the geared ÃD/D̃A mechanism, we obtain h_g = −0.10 MHz, in good agreement with both, −0.11 MHz of ref. 39 and −0.12 MHz obtained using PIMD.⁵⁵ The experimental value is −0.073 MHz, or ≈30% smaller. The difference between the PIMD and experiment could potentially be attributed to the accuracy of the PES. Instanton calculations³⁹ of TM elements using the potential of ref. 3 differ from those obtained using MB-pol by 36% for h_g and 17% for h_a, which suggest that the potential is possibly more accurate in the region probed by the mechanism AD.

3.1.3 Tunneling pathways and symmetry analysis in G₃₂. Vibrational excitations can have a drastic effect on the sizes of the splittings. Ref. 41 measured a more than 1000-fold increase in the TM element for the geared ÃD/D̃A mechanism in the fully deuterated water hexamer-d₁₂ prism for the excited librational mode at 510 cm⁻¹. GS showed no measurable splittings, which indicates that |h_g| < 1 MHz, while |h_g| = 1720 MHz, in the excited state. In systems with multiple tunneling pathways, the sizes of the TM elements are mode specific; the excitation of a particular vibrational mode can increase the size of one TM element and decrease the size of another, or even change its sign.⁵⁸ For that reason, we calculated the TM elements in the low-lying vibrational modes of the hexamer prism for all single-step mechanisms discussed above, as well as the rotations of water monomers B, C, E and F, associated with the symmetry elements (3 4), (5 6), (9 10) and (11 12), respectively. The MAP barriers for mechanisms ‘rot B’ and ‘rot F’ (not shown in Fig. 3) are higher, at 3333 cm⁻¹ and 3901 cm⁻¹ with the MAP lengths of 316.2m_e^1/2a₀ and 322.5m_e^1/2a₀, respectively. They can be neglected along with the single-step processes of ‘rot A’, ‘rot D’ and ÃD̃. ‘rot E’ and ‘rot C’ mechanisms have higher barriers along their MAPs than mechanisms AD and ÃD/D̃A, 1782 cm⁻¹ and 1910 cm⁻¹, but they pass near transition states at 1135 cm⁻¹ and 1310 cm⁻¹, which are comparable. The corrected barrier heights from the GS, including harmonic contributions of orthogonal modes at transition state, are 755 cm⁻¹ and 898 cm⁻¹, lower than for the geared ÃD/D̃A motion, but higher than for the antigeared double flip AD. The MAP potentials are shown in Fig. 3. We note that the tunneling in ‘rot E’ and ‘rot C’ proceeds over wider barriers than AD and ÃD/D̃A, with MAP lengths of 341.9m_e^1/2a₀ and 357.7m_e^1/2a₀, and that the paths are asymmetric, with the barrier top at 50.8% and 55.6% of the full pathway, respectively. The relevance of these processes to the splitting pattern is discussed next.

Inclusion of ‘rot E’ and ‘rot C’ mechanisms through TM elements h_E and h_C, respectively, generates a group G₃₂ of order 32. TM is then a 32 × 32 matrix with five non-zero elements in each row/column; one associated with h_a, h_E and h_C and two with h_g. Any minimum can be accessed in a maximum of five single-step processes. The group elements divide into 14 classes, defined in Table 9. We devised its character table in Table 10 and therewith defined the names of irreducible representations (irreps) of the group. Symmetries of tunneling states, Γ_tun, are found by reducing the representation of the full set of localized single-well states,⁶⁴ which has character 32 for identity and zero for all other symmetry elements. One obtains Γ_tun = A_1a ⊕ A_1b ⊕ A_2a ⊕ A_2b ⊕ B_1a ⊕ B_1b ⊕ B_2a ⊕ B_2b ⊕ 2E_1a ⊕ 2E_1b ⊕ 2E₂ ⊕ 2E₃ ⊕ 2E₄ ⊕ 2E₅. All one-dimensional representations appear only once in our set of tunneling states, which allows one to construct a symmetry-adapted linear combination of localized states for each of them uniquely. Expectation value of TM in symbolic form in a state of a particular symmetry gives us an analytical expression for its energy. Two-dimensional representations appear twice in the set of tunneling states. We generate four different symmetry-adapted vectors of the particular symmetry and use them to construct a 4 × 4 representation of the TM in this subspace in symbolic form. Diagonalization of this matrix leads to a biquadratic equation and gives us pairs of energies of two doubly-degenerate states. The analytic expressions for energy levels in terms of TM elements are given in Table 1 in the order in which they appear in the GS of water hexamer prism, i.e., assuming that |h_a| ≫ |h_g| ≫ |h_E| ≈ |h_C| and that all TM elements are negative. Setting h_E = h_C = 0 reproduces the set of D_2d energies. Each branch of the D_2d sextet now splits into a triplet of equal width 2|h_C + h_E|. The triplet is equally spaced only if h_C = h_E in this TM model.

Table 1 Analytic expressions for the tunneling energy levels of the water hexamer prism. Energy levels are labeled by the irreducible representations of the G₃₂ group

E(A_1a) = h_a + 2h_g + h_E + h_C

E(A_1b) = h_a + 2h_g − (h_E + h_C)

E(E_1a) = h_a + h_E + h_C

E(E_1b) = h_a − (h_E + h_C)

E(B_1a) = h_a − 2h_g + h_E + h_C

E(B_1b) = h_a − 2h_g − (h_E + h_C)

E(A_2a) = −(h_a − 2h_g) + h_E + h_C

E(A_2b) = −(h_a − 2h_g) − (h_E + h_C)

E(E_1a) = −h_a + h_E + h_C

E(E_1b) = −h_a − (h_E + h_C)

E(B_2a) = −(h_a + 2h_g) + h_E + h_C

E(B_2b) = −(h_a + 2h_g) − (h_E + h_C)

---

Total internal wavefunction is the product of rovibrational state and the nuclear spin state, Γ_rovib ⊗ Γ_nuc ⊃ Γ_int. The wavefunction must be antisymmetric under an exchange of hydrogen nuclei and symmetric under an exchange of oxygen nuclei, and its symmetry is thus Γ_int = B_1b. Inversion is not a symmetry element of the group and there are thus two sets of degenerate prism states. Nuclear spin states span Γ_nuc = 666A_1a ⊕ 78A_1b ⊕ 66A_2a ⊕ 6A_2b ⊕ 78B_1a ⊕ 48B_1b ⊕ 630B_2a ⊕ 66B_2b ⊕ 432E_1a ⊕ 48E_1b ⊕ 144E₂ ⊕ 144E₃ ⊕ 432E₄ ⊕ 48E₅. In the rotational ground state, the rotations do not affect the symmetry, so the vibrational tunneling states have the following statistical weights: Γ_vib = 10A_1a ⊕ 78A_1b ⊕ 66A_2a ⊕ 630A_2b ⊕ 78B_1a ⊕ 666B_1b ⊕ 6B_2a ⊕ 66B_2b ⊕ 48E_1a ⊕ 432E_1b ⊕ 144E₂ ⊕ 144E₃ ⊕ 48E₄ ⊕ 432E₅. For the fully deuterated hexamer-d₁₂ prism, internal wavefunction is of Γ_int = A_1a symmetry and the statistical weights are then Γ_vib = 52 650A_1a ⊕ 13 203A_1b ⊕ 13 041A_2a ⊕ 3240A_2b ⊕ 13 203B_1a ⊕ 3321B_1b ⊕ 52 326B_2a ⊕ 13 041B_2b ⊕ 52 488E_1a ⊕ 13 122E_1b ⊕ 26 244E₂ ⊕ 26 244E₃ ⊕ 52 488E₄ ⊕ 13 122E₅. None of the states have a statistical weight zero. Following the analysis of ref. 39 in the D_2d group, the representation of the dipole moment operator has character 3 under all group operations that do not switch the oxygens A and D and character −1 for the operations that do and thus require a rotation R^π₀ to return the oxygen framework into its initial configuration. Therefore, Γ_dip = A_1a ⊕ 2B_2a. The largest dipole component lies along the principal axis that points approximately perpendicular to the two triangular bases of the prism and is of B_2a symmetry.

TM elements for the four rearrangement mechanisms discussed above and computed using the M-WKB method are listed in Table 2 for the GS and the 13 excited lowest-frequency normal modes of the prism. TM elements h_E and h_C include a multiplicative factor of two to take into account the clockwise and counter-clockwise monomer rotation, in an analogous way to the acceptor tunneling motion in water dimer.¹⁸ The convergence properties of the method degrade for the states with higher excitation energies. We use a heuristic definition of the uncertainty in TM element as the maximum deviation in its size when ε is varied in the interval 0.1–10m_e^1/2a₀ and N between 100–1000. The uncertainty in the states GS-5 is less than 2%. The uncertainty in the TM elements for the excited modes 6–13 goes up to 40%, with several exceptions where the error is larger. The excited mode 10 for ÃD/D̃A, 9, 10 and 12 for ‘rot E’, and 12 and 13 for ‘rot C’ come with relative errors of 50–100%. The sizes of the TM elements for ‘rot E’ in the excited modes 11 and 13 and for ‘rot C’ in the excited mode 8 are smaller than 1.4 × 10⁻⁶ cm⁻¹, 4.0 × 10⁻⁵ cm⁻¹ and 1.0 × 10⁻⁶ cm⁻¹, respectively, with its sign uncertain. They are thus set to zero in Table 2.

Table 2 Tunneling matrix (TM) elements of the water hexamer prism (H₂O)₆ on MB-pol PES^72–74 for rearrangement paths AD (h_a), ÃD/D̃A (h_g), rot E (h_E) and rot C (h_C) in the vibrational ground state and the 13 lowest-frequency excited vibrational modes in cm⁻¹

Mode	ω e	h a/10^−5	h g/10^−6	h E/10^−6	h C/10^−7
GS	0	−2.47	−3.44	−2.92	−2.09
1	62.1	−2.29	−3.53	−3.44	−2.71
2	72.0	−2.95	−4.48	−3.29	−2.33
3	74.4	−0.759	0.0639	−2.96	−2.38
4	110	−0.814	1.56	−2.89	−6.20
5	120	−1.20	−0.351	4.17	−2.09
6	150	−1.5	−2.6	−13.6	−0.60
7	172	−2.0	−1.8	−17.2	7.7
8	179	−5.6	−3.0	−10	0
9	204	−2.7	−10	−4.0	−14
10	218	−19	33	12	−2.8
11	230	−1.9	10	0	−26
12	240	−0.13	2.5	−20	6.0
13	271	−2.7	7.0	0	−22

---

TM elements in the excited states are determined by the projections of U(S_Σ) on the eigenvectors of A(S) with eigenvalues ω_i(S). These projections at Σ are obtained by integration of eqn (6) and exhibit, adiabatically, an exponential growth ∝ exp[(ω_e − ω_i)τ] for ω_i < ω_e. Tiny projections of U_e at S = ε on lower frequency eigenvectors of A can lead to significant contributions at S = S_Σ and this creates numerical noise that leads to uncertainty in the results⁵⁸ and finally to the failure of the method for the higher excited modes.^58,65,83,84

3.1.4 Ground-state tunneling splittings in G₃₂. The TS pattern in the GS of the water hexamer prism is obtained by inserting the TM elements in Table 2 into expressions in Table 1. It is displayed in the right-most set of lines in Fig. 4. The calculated energies are given as shifts relative to the GS energy of the single-well localized GS. Each sextet branch is split into a triplet with a width of 6.3 × 10⁻⁶ cm⁻¹ (or 0.2 MHz). This is comparable to the ÃD/D̃A splitting. The largest B_2a component of the dipole induces the transitions between A_1a ↔ B_2a, A_1b ↔ B_2b and E₄ ↔ E₄ which overlap the pair of lines for the A₁ ↔ B₂ transition in the D_2d group analysis of the spectrum in ref. 39. The analogous analysis applies to the transitions B₁ ↔ A₂ and E ↔ E in D_2d. Ref. 39 does not report a change in the splitting pattern when monomers C and E are substituted by the heavier ¹⁸O isotopes. This means that h_E is overestimated by a factor of more than two, because the effect of this mechanism was not captured in the experiment. h_C is smaller by an order of magnitude than h_E and is not manifested in the appearance of the TS pattern.

A way to reveal whether the ‘rot E’ mechanism is relevant for the GS TS pattern of the hexamer prism would be to perform measurements on the hexamer prism with the deuterated monomers A and D. TM elements for that isotopologue, obtained using M-WKB, are reported in the ESI† in Table S1. Our calculations estimate that the GS splitting due to ‘rot E’ is then reduced by 9% to 5.32 × 10⁻⁶ cm⁻¹ (or 0.16 MHz), which is smaller than the observed ÃD/D̃A splitting with line separation of 0.29 MHz. Nevertheless, it could possibly be observed if the calculations do not overestimate it. In the excited modes 1–5, the ‘rot E’ splitting is up to 16% enhanced (in the excited mode 2) and does not change appreciably. The TM elements for other mechanisms are 1–3 orders of magnitude smaller.

3.1.5 Excited-state tunneling splittings in G₃₂. In Fig. 5, we plot the widths of the AD doublet, the ÃD/D̃A triplets and the ‘rot E + C’ triplets, that are equal to 2|h_a|, 4|h_g| and 2|h_E + h_C|, respectively, for the GS and the excited modes 1–13. As evident in the figure, TS patterns in the excited modes 1, 2, 9, 10 and 13 will have a similar appearance to the GS TS pattern (noting also that |h_E| ≫ |h_C|). TM element h_a has the same sign in all the studied excited states. h_g changes sign in the excited modes 10 and 13 and thus the order of the states in each triplet in the D_2d group in Fig. 4 is reversed. The same holds true for the triplet levels in the G₃₂ group for the excited mode 10, where h_E has the opposite sign to the GS and the order of levels is reversed. The size of h_a varies from that in the excited mode 12, where it is reduced by a factor of ≈20 compared to the GS (uncertainty of 30%), to the excited mode 10, where it is increased by a factor of 7.7 (uncertainty of 5%). The limits of variation in h_g relative to the GS are between the reduced size, by a factor of ≈50 in the excited mode 3, to the increased size, by a factor of ≈10 in the excited mode 10.

Fig. 5 Widths of the energy splittings in the water hexamer prism: doublet width 2|h_a|, triplet width 4|h_g| and a further triplet width 2|h_E + h_C| for the ground state and the 13 lowest-frequency excited vibrational modes. See the text for details.

The 13 lowest frequency normal modes of vibration of the hexamer prism are depicted in the ESI† in Fig. S1. Normal modes below ≈270 cm⁻¹, 15 in count, are the shape deformation modes and torsional modes of ‘free’ hydrogens 3, 7 and 11. Torsions are significantly mixed into modes 10, 12 and 13–15. The region between ≈270–1000 cm⁻¹ contains librational modes and has been accessed in a previous experiment,⁴¹ but is outside the convergence domain of the present theory. Intramolecular monomer modes lie above 1600 cm⁻¹. Torsional motion of monomer D in the excited mode 10 is seen to significantly enhance the splitting due to AD and ÃD/D̃A mechanisms. The excitation of the longitudinal mode, the one with the largest projection on the MAP (U_et) at ε, usually results in the enhancement of the associated splitting.⁶³ For the AD mechanism, the longitudinal mode is mode 3 and its excitation results in the decreased |h_a|. This is the result of the cancelation of two terms in the brackets in the 2nd line of eqn (7). Large F term (1st term in the bracket) is the usual cause of the enhancement, but, in this case, mode 3 is diabatically connected to mode 4 at the transition state and this results with a large U-term (2nd term in the bracket) of opposite sign, due to the alignment of U(S_Σ) with the eigenvector of Ā. Other MAPs have no distinct longitudinal mode, and we thus observe no large TS enhancements for the excited low-frequency modes.

In the excited modes 3–8, TM element |h_E| > |h_g| (h_C is non-manifest). The TS pattern for the excited mode 3 is shown in Fig. 6. The appearance is similar to the GS TS pattern, but the underlying mechanisms of ÃD/D̃A and ‘rot E/C’ switch roles. From left to right, the sets of levels show the doublet splitting induced by h_a, further splittings induced by the monomer rotations h_E/C (with h_g set to zero) and finally the full splitting pattern, including h_g, consisting of 18 energy levels, now in a different order relative to the GS, as the symmetry labels suggest. The mechanisms included in the mid set of levels in Fig. 6 do not form a group. The fine splitting of triplets into triplets due to ÃD/D̃A in Fig. 6 is below the current experimental resolution.³⁹

Fig. 6 Tunneling splitting pattern of the excited vibrational mode 3 of the water hexamer prism on MB-pol PES.^72–74 The sets of levels, from left to right, are the energy splittings due to AD, AD + rot E + rot C, and AD + ÃD/D̃A + rot E + rot C. See the text for notation.

In the excited mode 12, the sizes of all widths in Fig. 4 are reversed. The uncertainties in h_a, h_g, h_E and h_C are 30%, 20%, 50% and 90%, respectively. Inside of these error ranges, TM elements do not change order (apart from possibly h_E and h_C). The TS pattern in this state, shown in Fig. 7, is highly speculative, due to the high uncertainties in the TM element sizes. However, in this example, we show how the TS patterns in higher vibrational states might appear. The left-most set of levels displays the doublet of doublets splitting due to the largest TM elements, h_E and h_C (with h_a and h_g set to zero), whereby h_C is non-negligible here. The mid set of levels includes the h_g contribution, in which both doublets split their branches into a doublet and a triplet, which overlap. In the full splitting pattern on the right side of the figure, each branch is further split into a doublet by h_a apart from levels at energies ±(h_E − h_C).

Fig. 7 Tunneling splitting pattern of the excited vibrational mode 12 of the water hexamer prism on MB-pol PES.^72–74 The sets of levels, from left to right, are the energy splittings due to rot E + rot C, rot E + rot C + ÃD/D̃A, and AD + ÃD/D̃A + rot E + rot C. See the text for notation.

3.1.6 Isotopic substitutions. We also calculated the TS patterns in the GS and the excited low-lying modes 1–5 for the isotopically substituted hexamer prism (H₂¹⁸O)₆. The results are summarised in the ESI† in Table S2. TM elements in the GS are reduced by 17% for AD, ÃD/D̃A and ‘rot C’ and by 42% for ‘rot E’. Monomer rotation thus plays a smaller role in this case. In the excited modes 1 and 2, a small decrease of TSs relative to ¹⁶O-prims is also present, while the TSs in modes 3–5 are very similar for both isotopologues.

For the fully deuterated hexamer-d₁₂ prism, (D₂O)₆, the TS pattern is unobservable. Results are given in the ESI† in Table S3. The low-lying normal modes of hexamer-d₁₂ align well with the normal modes of the hexamer-h₁₂. The AD doublet width in the GS is 7.7 × 10⁻⁸ cm⁻¹, while the ÃD/D̃A triplet width is 10 times smaller. The mechanisms ‘rot E/C’ are negligible for the appearance of the TS pattern for all states (with a 10 times smaller triplet width in the GS), apart from the excited mode 5, where the mechanism ÃD/D̃A and ‘rot E’ give approximately equal contributions. The sizes of the TM elements for AD and ÃD/D̃A upon excitation vary in a similar way to the hexamer-h₁₂.

3.2 Cage

The most stable isomer of the water hexamer is the cage shown in its equilibrium geometry in Fig. 8 (1), where atom labels are also introduced. Its minimum is 89.1 cm⁻¹ higher in energy than the minimum of the hexamer prism in Fig. 1, on MB-pol PES, but the harmonic ZPE is 91.4 cm⁻¹ smaller than that of the prism, making it slightly more stable. More accurate calculations of the ZPEs⁸⁵ confirm this. The other three minimum structures of the cage, shown in Fig. 8 (2–4), are linked to it by a torsional flip of the terminal OH bonds in the doubly-bonded monomers A and B.^36,67 Potential energies of cage structures i = 1–4 at their respective minima, V_min,i, differ by no more than 56.6 cm⁻¹, as seen in Table 3. The total local GS energies E_i, including the ZPE in harmonic approximation, bring them closer together, within 21.5 cm⁻¹.

Fig. 8 Equilibrium geometries of four non-degenerate minima of the water hexamer cage ordered in energy 1–4 on MB-pol PES.^72–74

Table 3 Relative energies of the water hexamer cage isomers in Fig. 8 (1–4) on MB-pol PES^72–74 in cm⁻¹. The first column refers to the potential energies relative to cage 1; the second column refers to the total energy differences, including harmonic zero-point energy; the third and fourth column list the coefficients of the ground-state wavefunction as a superposition of localized single-well cage 1–4 states, , obtained using M-WKB and using isomer fractions from ref. 85, c^MM_i

Cage	V min,i	E i	c i	c ^MM i
1	0	0	0.559	0.574
2	12.8	3.63	0.480	0.537
3	31.3	11.7	0.519	0.454
4	56.6	21.5	0.433	0.418

---

3.2.1 Dynamics of flips. Single OH flip of monomer B connects cage isomers 1 ↔ 2 and 3 ↔ 4. The MAP for the former is shown in Fig. 9(a). The MAPs of the two flips are 168.2m_e^1/2a₀ and 168.6m_e^1/2a₀ long and have a slight asymmetry, with barriers at 50.6% and 51.3% along the pathway at 175 cm⁻¹ and 183 cm⁻¹ relative to cage 1 and cage 3 minima, respectively. The ZPE corrected barriers are −9.67 cm⁻¹ and −3.24 cm⁻¹ relative to cage 1 and 3, respectively. Flip of monomer A connects minima 1 ↔ 3 and 2 ↔ 4. The MAPs of these mechanisms are shorter, 111.6m_e^1/2a₀ and 111.3m_e^1/2a₀, have a more asymmetric potential profile along the MAP, with minima at 58.7% and 63.3% along the pathway, and lower barriers, 56.4 cm⁻¹ and 59.8 cm⁻¹, respectively. Corrected GS barrier heights, relative to the cage 1 and 2 GSs, are −57.9 cm⁻¹ and −52.7 cm⁻¹, respectively. Flips of monomer A and B are over-barrier motions and the cage GS is a superposition of four cage structures,^38,67,85 bringing further stabilization of its energy with respect to the hexamer prism.

Fig. 9 Minimum action paths of the rearrangement mechanisms in the water hexamer cage involving a doubly-bonded water monomer: (a) flip B, (b) bifurcation B̃ and (c) rotation of monomer B (rot B).

The localized GS wavefunctions ϕ_i of the four non-degenerate cage structures i may be used as a basis to construct a 4 × 4 Hamiltonian matrix, in an analogous way to the design of TM. Diagonal entries of the TM are the local GS energies E_i of cage i = 1–4 in Table 3. Off-diagonal TM elements h^ij_A/B quantify the interactions between the cage structures i and j and the letter in the subscript indicates the monomer involved in the motion. They are calculated using M-WKB and given in Table 4. Eigenvalues of TM are then the energies of the four states, delocalized over four cage structures, and form a quartet. The energy level shifts relative to the cage 1 GS are −115 cm⁻¹, −54.9 cm⁻¹, 73.8 cm⁻¹, and 133 cm⁻¹. It is known that the instanton theories overestimate the TM elements for the over-barrier states significantly.⁵⁸ The off-diagonal TM elements are similar in size to the diagonal energies E_i and also to the excitation energies of vibrational modes, which may consequently interact. The error is thus of similar size to the energy differences of vibrational states. It is therefore not possible to give meaningful estimates of the state energies using only the harmonic energies and M-WKB TM elements. Nevertheless, we list the coefficients c_i of the wavefunction, , in the lowest branch of the quartet, in Table 3. They are compared there with the coefficients c^MM_i, obtained as a square root of the isomer fractions of the cage, calculated in the DMC study of ref. 85 on the same PES. We use these coefficients to estimate the tunneling splittings in the lowest branch of the GS quartet due to bifurcations, below.

Table 4 Tunneling matrix elements of the water hexamer cage isomers in Fig. 8 (1–4) on MB-pol PES^72–74 in the ground state (GS) and the excited modes 1–15, in cm⁻¹. Labels are defined in the text

Mode	h B ^12	h B ^34	h A ^13	h A ^24
GS	−30.6	−28.8	−95.6	−92.2
	h B̃ ^12/10^−4	h B̃ ^34/10^−4	h à ^13/10^−5	h à ^24/10^−6
GS	−7.65	−8.09	−1.04	−6.16
1	−7.67	−8.11	−1.06	−6.18
2	7.70	8.11	−0.414	−3.08
3	8.07	8.51	1.08	6.39
4	−7.99	−8.42	−1.08	−5.74
5	−7.34	−7.74	−1.08	−7.13
6	−8.1	−5.5	−0.56	—
7	−9.0	3.8	3.0	0
8	−7.7	−8.9	−0.78	−5
9	−9.0	−9.7	−3	—
10	−13	−15	—	—
11	−7.0	5.0	−2	−6
12	−9.1	−11	−1	−12
13	−3.0	−4.6	−1	—
14	−70	60	—	−15
15	−3	−8	−100	150

---

3.2.2 Dynamics of bifurcations and symmetry analysis. Bifurcation mechanisms, which break and reform hydrogen bonds in the hexamer cage, were first identified in ref. 36 and 37 as rotations of monomers A and B in cage structure 1. Other mechanisms were also identified and found to be irrelevant to the TS pattern in a theoretical study of ref. 67. Subsequently, ref. 38 considered four cage isomers linked by flips between the structures, as described above, and also by bifurcations, that change structures and switch the positions of hydrogens 1 and 2 in monomer A, or 3 and 4 in monomer B, in the process. The MAP for bifurcation mechanism B̃, which links cage structures 1 and 2 and also switches the hydrogen 3 by 4 in the hydrogen bond, is shown in Fig. 9(b). The rotation of monomer B, shown in Fig. 9(c), can thus be viewed as a flip B from cage 1 to cage 2, followed by a bifurcation B̃ back from cage 2 to cage 1.

The switching of hydrogens (1 2) and (3 4) generates a 4-element molecular symmetry group of the cage that is isomorphic to the point group C_2v. Its character table is given in Table 3 of ref. 37. Here, we follow ref. 38 and use the C_2v irrep names. Each cage structure has four additional versions generated by the permutations of hydrogens. The localized GS wavefunctions of the 16 minima are denoted as φ^k_i, where i = 1–4 and k = E, (1 2), (3 4) or (1 2)(3 4). Ref. 38 constructs a 16 × 16 TM as a matrix representation of Hamiltonian in this basis. The states φ^k_i have energies E_i. TM elements for flips h^ij_A/B in Table 4 link the like, i.e., same k, permutational isomers of cage i and j. TM elements for bifurcations, denoted h^ij_Ã/B̃, link different permutational isomers of cage i and j, specifically, those that are connected by the symmetry element (3 4) for h_B̃¹² and h_B̃³⁴, and by the symmetry element (1 2) for h_ù³ and h_ò⁴. There are thus four non-zero off-diagonal TM elements in each row/column of TM. Diagonalization of the TM produces a TS pattern in which the flips split the GS into a quartet and bifurcations further split each quartet branch into another quartet, whereby we assumed that |h^ij_A/B| ≫ |h^ij_Ã/B̃|.

In order to determine the TS pattern due to bifurcations in the lowest branch of the flip quartet, we note that the GS of a particular symmetry in C_2v(M) can be written as a linear combination of ϕ_i's of that symmetry with coefficients c_i from Table 3. Each ϕ_i of a particular symmetry is uniquely obtained as a symmetry-adapted linear combination of φ^k_i, since all irreps of the group are one-dimensional. The energies of the four states comprising the bifurcation quartet in the lowest flip quartet branch are obtained analytically as expectation values of the TM in symbolic form, whereby we set energies E_i and TM elements for the flips to zero. These have been used to determine c_i above. The resulting energies are listed in Table 5.

Table 5 Analytic expressions for the tunneling energy levels of the water hexamer cage in terms of the coefficients c_i in Table 3 and the bifurcation tunneling matrix elements in Table 4. Coefficient c_i describes the parent vibrational state as a linear combination of cage structures i = 1–4. Energy levels are labeled by the irreducible representations of the C_2v group. Notation is defined in the text

E(A₁) = 2c₁c₂h_B̃¹² + 2c₁c₃h_ù³ + 2c₂c₄h_ò⁴ + 2c₃c₄h_B̃³⁴

E(B₁) = 2c₁c₂h_B̃¹² − 2c₁c₃h_ù³ − 2c₂c₄h_ò⁴ + 2c₃c₄h_B̃³⁴

E(A₂) = −2c₁c₂h_B̃¹² + 2c₁c₃h_ù³ + 2c₂c₄h_ò⁴ − 2c₃c₄h_B̃³⁴

E(B₂) = −2c₁c₂h_B̃¹² − 2c₁c₃h_ù³ − 2c₂c₄h_ò⁴ − 2c₃c₄h_B̃³⁴

---

The remaining task is to determine the TM elements for bifurcations using M-WKB. The MAPs for bifurcation B̃ linking cages 1 ↔ 2 and 3 ↔ 4 are 234.0m_e^1/2a₀ and 233.5m_e^1/2a₀ long, slightly asymmetric with barriers at 50.9% and 51.0% along the pathway, having heights of 1041 cm⁻¹ and 1055 cm⁻¹. The nearby transition states are at 989.3 cm⁻¹ and 1006 cm⁻¹ and the corrected barrier heights are 815 cm⁻¹ and 809 cm⁻¹, relative to the GSs of cage 1 and 3, respectively. Potential profiles along the MAPs for the flip and the bifurcation of monomer B linking cages 1 and 2 are shown in Fig. 10 (red line). The MAPs for bifurcation à linking cages 1 ↔ 3 and 2 ↔ 4 are significantly longer, 301.8m_e^1/2a₀ and 306.8m_e^1/2a₀, with larger barriers, 1597 cm⁻¹ and 1679 cm⁻¹ at 49.7% and 49.9% along the pathway, respectively. Transition state energies are 1503 cm⁻¹ and 1586 cm⁻¹, with corrected barrier heights of 1228 cm⁻¹ and 1299 cm⁻¹ relative to the GS of cage 1 and 2, respectively. Potential profiles along the MAPs for the flip and the bifurcation of monomer A linking cage structures 1 and 3 are also shown in Fig. 10 (black). The inset in the figure compares the MAP for the rotation of monomer A with that of the flip followed by a bifurcation. The paths link the same versions, but the ‘rot A’ mechanism is shorter. Due to a corner cutting effect, the MAP for rotation does not quite reach the cage 3 minimum. The direct route of rotation may contribute some additional flux in the TM matrix, but we neglect it here because it would lead to the double-counting of contributions due to a large overlap with the two-step process of the flip followed by a bifurcation. Moreover, it is not possible to calculate the TM element for all rotations in cages 1–4 using M-WKB. The rotation of monomer A in cage 3, rotation of monomer B in cage 2 and rotation of A and B in cage 4, all proceed along MAPs that pass near the lower-energy minima. The present M-WKB formalism cannot be applied to the regions of space where the potential is lower than that at the end points of the path (see eqn (2)). We also note here that the multiplication of the TM element for rotation, calculated on a single MAP, by two, to take into account the monomer rotation in the opposite senses, corresponds to two two-step processes of a flip followed by a bifurcation and vice versa, in our treatment here.

Fig. 10 Potential energy curves along the minimum action paths (MAP) for the flip (dashed line) and the bifurcation (full line) of the doubly-bonded monomers A (black) and B (red), positioned at the opposite vertices of the water hexamer cage. The inset shows the potential energy curve along the combined MAP for flip (A) and bifurcation (Ã) (black) alongside the MAP for rotation (blue) of monomer A.

TM elements for bifurcations in the GS and several excited vibrational modes, numbered in order of the increasing frequency, are listed in Table 4. Again, convergence within ≈2% is achieved for the states up to the excited mode 5. The accuracy of the results drops for higher states and is similar to that reported for the hexamer prism above. The excited modes 1–15 lie in the frequency range from 43 cm⁻¹ to 267 cm⁻¹ and they are associated with the cage shape distortions and torsions of ‘free’ OH bonds. Librations correspond to mode 17 and higher. Significant torsional motion of monomer A is present in modes 7–9, while the torsion of monomer B is present in modes 14, 11 and 12. Longitudinal mode, the one with the largest projection on MAP tangent at ε, is either mode 1 or 2, but the excitation in those modes does not result in any significant increase of the associated TM element. The TM elements are largely unchanged in the excited modes 1–5. We observed a large enhancement of the TM elements in the excited modes 14 and 15, but no case in which the TM element for mechanism à is larger than that for mechanism B̃. We have also calculated the TM elements for the rotation of other water monomers in the cage and found that, in the GS and other low-lying vibrational states, they lie in the region 10⁻¹⁰–10⁻⁸ cm⁻¹ and can safely be neglected.

3.2.3 Tunneling splittings. The TS pattern in the lowest flip quartet branch is obtained by inserting coefficients c_i from Table 3 and TM elements from Table 4 in the expressions given in Table 5. The splittings are schematically presented in Fig. 11 (left panel). The GS is split into a doublet with a width of 1.55 × 10⁻³ cm⁻¹, whereby each branch is further split into a doublet with a width of 1.72 × 10⁻⁵ cm⁻¹. The larger splitting is due to mechanism B̃ and the smaller due to mechanism Ã. We recalculated the widths using coefficients c^MM_i estimated from the isomer populations in the DMC study of ref. 85, and obtain practically unchanged results; the B̃ splitting is 1.56 × 10⁻³ cm⁻¹ and the à splitting is 1.64 × 10⁻⁵ cm⁻¹. If we only consider the cage 1 isomer and neglect the presence of other isomers in Fig. 8, as in ref. 36 and 37, the recalculated TS pattern, using the TM elements for the mechanisms ‘rot B’ and ‘rot A’, results in a B̃ width of 2|h_rot B| = 3.48 × 10⁻³ cm⁻¹ and an à width of 2|h_rot A| = 8.12 × 10⁻⁵ cm⁻¹.

Fig. 11 A qualitative diagram of the tunneling splitting pattern of the water hexamer cage in the ground state, obtained using theory, and in the excited state at 83 cm⁻¹, determined in the experiment of ref. 36.

In a terahertz laser spectroscopy experiment, ref. 36 and 37 interpreted the measured TS pattern in the hexamer cage at the band origin of 83 cm⁻¹ as an equally spaced triplet, having the width of 1.28 × 10⁻⁴ cm⁻¹. The TS pattern was rationalized in terms of two bifurcation mechanisms, ‘rot A’ and ‘rot B’, having equal TM elements and thereby causing an accidental degeneracy of the A₂ and B₁ energy levels, as shown schematically in Fig. 11 (right panel). Our calculations do not support the presumption that the two mechanisms give equal contributions either in the GS or the low-lying excited states. A large difference in the barriers along the two rearrangement paths, shown in Fig. 10, underlies our results, while the accuracy of the PES is validated by the correct structures it predicts. Ref. 37 presents detailed arguments in support of the accidental degeneracy, including absence of Stark shifts, the correct predicted intensity pattern of transitions and the rotational spacings. Coefficients c_i in the vibrationally excited states cannot be predicted using harmonic approximation and M-WKB TM elements, but it is evident that the excited states will make a highly delocalized mixture of localized vibrational states of different cage isomers. The coefficients c_i, with |c_i| ≤ 1 and , are not easily tuned to reproduce the TS pattern in Fig. 11 (right panel). They would need to satisfy c₁c₂ ≈ −c₃c₄ to bring about the partial cancellation of h_B̃¹² and h_B̃³⁴ terms (assuming h_B̃¹² ≈ h_B̃³⁴), in order to reduce the overall width of the TS pattern to ≈10⁻⁴ cm⁻¹ as observed in experiment. With this requirement, the bifurcation mechanisms à also interfere destructively in the TS pattern (assuming they are of the same sign) and, using values in Table 4, cannot quite reach the magnitude of 3.2 × 10⁻⁵ cm⁻¹, as observed. The c_i's in the relevant excited state could possibly be obtained in a separate calculation, e.g., using a reduced-dimensional model, in view of the fact that the dynamics of flips, which determines them, is predominantly localized at the two opposite vertices of the cage. The GS TS pattern that we report here is not very sensitive to the exact knowledge of c_i's and depends mainly on our estimates of the TM elements for bifurcations. Bifurcations are tunneling motions through large barriers and, in this regime, M-WKB is expected to deliver quantitative results, as exact quantum calculations in the hexamer prism⁵⁵ and the vinyl radical⁶⁵ show.

In another experiment,⁴⁰ the TS pattern of a hexamer cage was measured in a librational band at 525 cm⁻¹. TM elements for the mechanisms à and B̃ were estimated at 0.171 cm⁻¹ and 0.0892 cm⁻¹. This region of spectrum is outside the limits of M-WKB, but large enhancements for modes that are parallel to the rearrangement path is not surprising. We note, however, that the experiment implies |h_rot A| > |h_rot B|, which is certainly possible, although the barrier for mechanism à is substantially higher. Further investigations are therefore needed in order to reconcile theoretical results with the experimental findings.

3.3 Book

Hexamer book in its equilibrium geometry is shown in Fig. 12 (1). Atom labels are again introduced to distinguish different permutational isomers. The potential energy of the book isomer 1 in its minimum is 369 cm⁻¹ higher than that of the hexamer prism in Fig. 1, and 280 cm⁻¹ higher than that of the cage 1 in Fig. 8 (1). Its harmonic ZPE is 296 cm⁻¹ smaller than in the prism and 204 cm⁻¹ smaller than in cage 1. The GS energy of the book is thus 73.4 cm⁻¹ above the prism and 75.7 cm⁻¹ above the cage 1.

Fig. 12 Equilibrium geometries of four non-degenerate minima of the water hexamer book ordered in energy 1–4 on MB-pol PES.^72–74

Following the steps we have taken to determine the TS pattern in the cage, we looked for other minimum structures by flipping of ‘free’ H atoms of the book and found five relevant minima; four are shown in Fig. 12 and one in the ESI† in Fig. S2. The potential energies V_min,i of book structures i = 2–5 and their total local GS energies E_i, which include the harmonic ZPE contribution, relative to book 1, are listed in Table 6.

Table 6 Relative energies of the water hexamer book isomers in Fig. 12 (1–4) on MB-pol PES^72–74 in cm⁻¹. The first column refers to the potential energies relative to book 1; the second column refers to the total energy differences, including harmonic zero-point energy; the third and fourth columns list the coefficients of the ground-state wavefunction as a superposition of localized single-well book 1–4 states, , obtained using M-WKB using four (c_i), three (c^three_i) and two (c^two_i) book structures in the model. See the text for details

Book	V min,i	E i	c i	c ^three i	c ^two i
1	0	0	0.900	0.905	0.905
2	19.3	29.0	0.430	0.422	0.425
3	73.0	94.9	0.0581	0.0535
4	349	271	0.0291
5	200	106

---

3.3.1 Dynamics of flips. We then determined the MAPs that connect minima 1–5 to reveal the mechanisms. A simultaneous flip of monomers A and B, denoted AB, connects book structures 1 ↔ 2 and 3 ↔ 4. The geometries along the MAP for the double flip AB connecting book 1 ↔ 2 are shown superposed in Fig. 13(a). The MAP length is 256.8m_e^1/2a₀ and the potential along the MAP, shown in Fig. 14 (black line), has a barrier of 293 cm⁻¹ relative to book 1. The MAP for the double flip AB between book structures 3 ↔ 4 has a barrier of 370 cm⁻¹ relative to book 3 and the length of 276.1m_e^1/2a₀. Analogously, a double flip CD connects book structures 1 ↔ 3 and 2 ↔ 4 with MAPs with lengths of 349.7m_e^1/2a₀ and 302.9m_e^1/2a₀, and barriers of 359 cm⁻¹ and 491 cm⁻¹, relative to book 1 and 2, respectively. Book structure 5 connects by a flip of monomer D to book 1 with a MAP barrier of 1.65 cm⁻¹ at the distance of S_Σ = 18.32m_e^1/2a₀ from the minimum. The frequency of the book-5 normal mode, that is parallel to the MAP near the minimum, is 33.7 cm⁻¹. This kind of shallow well cannot accommodate a localized ‘bound’ state that can serve as a basis function in the TM approach and cannot be adequately treated using M-WKB. We thus use the localized GS wavefunctions of book 1–4 as basis functions ϕ_i in the TM model. This does not mean that the region around the minimum of book 5 is excluded from the treatment. The wavefunction along the paths that pass near it is constructed via M-WKB equations. Book structure 5 also connects via a flip of monomer C to book 3, with a barrier of 59.9 cm⁻¹. All paths have an asymmetric potential profile along the MAP and a derivative discontinuity at the barrier top, since they connect minima at different energies. Action is calculated using eqn (3), with V_min set at different values on the two sides of the dividing surface Σ.

Fig. 13 Minimum action paths of the rearrangement mechanisms in the water hexamer book, involving water monomers A and B (see Fig. 12): (a) flip AB, (b) bifurcation ÃB, (c) bifurcation B̃A, (d) rotation of monomer A (rot A), and (e) rotation of monomer B (rot B).

Fig. 14 Potential energy curves along the minimum action paths for the mechanisms of a double-flip AB, a bifurcation accompanied by a flip ÃB and B̃A, and rotations of monomers A (rot A) and B (rot B), all shown in Fig. 13. See the text for notation.

A 4 × 4 TM is constructed using the local GS energies E_i from Table 6 on diagonal and TM elements h^ij_AB/CD calculated via M-WKB in Table 7 at positions ij. There are two non-zero off-diagonal elements per row/column for the mechanisms of a double-flip AB and CD. Energies are obtained as −8.75 cm⁻¹, 36.3 cm⁻¹, 93.5 cm⁻¹ and 274 cm⁻¹. The interaction between the states of different book isomers brings a further stabilization of its energy, but does not change the energy ordering with the cage or the prism. TM elements for the flips are similar in size to the energy differences between the excited vibrational states (five normal modes have frequencies below 100 cm⁻¹) with which they may consequently interact. Due to the low vibrational frequencies of the book modes, the spread of the wavefunction is significant. Rearrangements start off from minima with a visually noticeable skeletal deformation over a flat region of the potential. Although the anharmonicity along the path is captured using M-WKB, it can still significantly affect the accuracy in other directions due to a large wavefunction spread. Therefore, we again concentrate on the lowest branch of the flip quartet to study its splitting pattern due to bifurcations below. The eigenvector of the TM that corresponds to its lowest eigenvalue gives us the coefficients c_i, listed in Table 6, of the delocalized GS wavefunction (of the lowest quartet branch), . We also list the coefficients in the three- and two-structure models c^three/two_i in Table 6, which only include book isomers 1–3 and 1–2, respectively, for later use.

Table 7 Tunneling matrix elements of the water hexamer book (H₂O)₆ on MB-pol PES^72–74 for rearrangement paths AB, ÃB, B̃A, CD, C̃D, D̃C, and rot A–D, connecting book structures 1–4, in the vibrational ground state, given in cm⁻¹. See the text for notation

h AB ^12	−17.5	h rot A ^11	−2.66 × 10^−4
h ÃB ^12	−1.09 × 10^−3	h rot B ^11	−5.22 × 10^−5
h B̃A ^12	−1.35 × 10^−4	h rot C ^11	−1.84 × 10^−5
h CD ^13	−6.11	h rot D ^11	−6.92 × 10^−6
h C̃D ^13	−5.81 × 10^−5	h rot A ^22	−5.02 × 10^−4
h D̃C ^13	−2.56 × 10^−5	h rot B ^22	−1.26 × 10^−4
h CD ^24	−16.5	h rot C ^22	−5.50 × 10^−6
h C̃D ^24	−2.11 × 10^−4	h rot D ^22	−1.53 × 10^−6
h D̃C ^24	−6.03 × 10^−5	h rot A ^33	−1.79 × 10^−5
h AB ^34	−18.0	h rot B ^33	−2.94 × 10^−5
h ÃB ^34	−4.79 × 10^−4	h rot C ^33	−7.12 × 10^−5
h B̃A ^34	−2.67 × 10^−4	h rot D ^33	−3.10 × 10^−5

---

3.3.2 Dynamics of bifurcations and symmetry analysis. Book structures linked by a double flip AB/CD can also be linked by bifurcations. The bifurcating monomer rotates around the axis perpendicular to the monomer plane to replace the hydrogen in the hydrogen bond by its ‘free’ hydrogen. Thereby, the in-bond hydrogen emerges on the opposite side of the book plane as ‘free’ hydrogen. This motion is accompanied by a flip of a monomer on the same side of the book. The mechanisms ÃB and B̃A, whereby the bifurcating monomer is capped by a tilde in our notation, are shown in Fig. 13(b) and (c), respectively. They link the same book structures as a double-flip AB, but the hydrogen atoms on the bifurcating monomer switch places in the final structure relative to the version reached by a double flip. For each of the four double flips discussed above, there are thus additionally two bifurcations mechanisms, consisting of a bifurcation accompanied by a simultaneous flip. These mechanisms are analogous to the mechanisms A1 (C̃A) and A2 (C̃B) in water trimer¹⁸ and the mechanisms ÃE and B̃C in water pentamer,⁵⁷ which link equivalent minima.

The mechanisms ÃB, B̃A, C̃D and D̃C involve a permutation of hydrogens (1 2), (3 4), (5 6) and (7 8), respectively. These four elements generate a 16-element commutative molecular symmetry group of the book, G₁₆. Each element is in a class of its own and we constructed its character table in the ESI† in Table S4 and named its irreps A₁ − A₁₆. Symmetries of tunneling states are obtained as . The internal wavefunction of (H₂O)₆ is of Γ_int = A₁₆ symmetry and of (D₂O)₆ is of Γ_int = A₁ symmetry. Nuclear-spin statistical weights of vibrational tunneling states, that govern the intensities of transitions, are given in the ESI† in Table S5 for the rotational ground state. No states have zero statistical weight.

In a two-structure model (book 1 and 2 in Fig. 12) discussed below, only mechanisms ÃB and B̃A are admitted in the TM model. The molecular symmetry group consists of four elements and is isomorphic to the C_2v point group. We choose the correspondence between its elements as (3 4) ≡ C₂, (1 2) ≡ σ_v(xz) and (1 2)(3 4) ≡ σ_v(yz) to define the irreps. The tunneling states are then obtained as Γ_tun = A₁ ⊕ A₂ ⊕ B₁ ⊕ B₂ and the internal wavefunction of (H₂O)₆ is of Γ_int = B₂ symmetry and of (D₂O)₆ is of Γ_int = A₁ symmetry. Nuclear-spin statistical weights of the rotation-less tunneling states are Γ_vib = 256A₁ ⊕ 768A₂ ⊕ 768B₁ ⊕ 2304B₂ for (H₂O)₆ and Γ_vib = 236 196A₁ ⊕ 118 098A₂ ⊕ 118 098B₁ ⊕ 59 049B₂ for (D₂O)₆.

Each book structure i = 1–4 therefore has 16 permutational isomers, labeled by k. TM is a 64 × 64 matrix in the φ^k_i basis set, where i labels structures and k labels group elements. TM elements for AB, ÃB, B̃A, CD, C̃D and D̃C, each linking two pairs of book structures, are inserted at appropriate positions in TM, following the procedure described for the cage above. We also consider internal monomer rotations below, linking like i structures, with its TM element denoted hⁱⁱ_rot X, where X = A–D, and i = 1–4. The eigenvalues of TM give the vibrational spectrum, whereby the lowest 16 states constitute the lowest flip quartet branch. Alternatively, we can derive the analytical expressions for the energy levels, as we have done for the cage in Table 5. We construct a state of a particular symmetry species in G₁₆ as , where χ_k are characters of the particular irrep (given in the ESI† in Table S4) and normalize it. We then combine the states ϕ_i using coefficients c_i from Table 6 and calculate the expectation value of TM in symbolic form, thereby setting the local GS energies E_i and the TM elements for flips h^ij_AB/CD to zero, for each symmetry species in G₁₆, to obtain analytical expressions for the energy shifts of vibrational states in the lowest quartet branch of the hexamer book. These are given in the ESI† in Table S6.

TM elements for all book mechanisms, which we determined using M-WKB, are listed in Table 7. We also calculated the TM elements for double bifurcations, e.g., ÃB̃, and for the rotations of the remaining two monomers (not labeled) and found them to be negligible, of the order 10⁻⁹–10⁻¹¹ cm⁻¹. The bifurcation mechanisms involving monomers ÖD̃ that start off from the book structure 1 have MAP barriers/MAP lengths of 951 cm⁻¹/317.0m_e^1/2a₀, 1249 cm⁻¹/359.0m_e^1/2a₀, 1195 cm⁻¹/419.9m_e^1/2a₀ and 1569 cm⁻¹/417.3m_e^1/2a₀, respectively. Potential energy profiles of ÃB and B̃A (1 ↔ 2) are shown in Fig. 14. MAP properties of bifurcation mechanisms involving monomers ÖD̃ that start off from the energetically higher book structures are similar, with path lengths within 15% and barrier heights within 20% of the corresponding monomer motions proceeding from book 1. Bifurcation MAPs of monomers B̃–D̃ become shorter in the higher structures, and all the barriers become higher (relative to the same monomer motion from the lower book structure). The sizes of the TM elements accordingly get larger in the order ÖD̃, which results in the clear separation of splitting widths in the TS pattern below. We also calculated the TM elements for internal monomer rotations, in Table 7, for book structures 1–3. The MAPs for monomer rotations in book 4 pass through regions where the potential energy falls below V_min,4 (apart from ‘rot D’), so they are not listed. Potential energy profiles for internal monomer rotations exhibit either two maxima or two inflexion points on one side of the dividing surface Σ; those for ‘rot A’ and ‘rot B’ in book 1 are shown in Fig. 14.

3.3.3 Ground-state tunneling splittings. The TS pattern of the hexamer book is schematically shown in Fig. 15 with symmetry labels of C_2v(M) and G₁₆, as appropriate. Inclusion of the bifurcation mechanism involving monomers ÖD̃, in that order, splits degenerate levels into equally-spaced doublets with widths w_A − w_D, as defined by Fig. 15. Table 8 lists the widths of the doublets obtained numerically, depending on which TM elements we include in the TM. In the first three columns, we include the TM elements for flips/bifurcations and exclude those for the monomer rotations from the model. The treatment is then analogous to the hexamer cage above. The first column lists the widths obtained by including all book structures 1–4 in a 64 × 64 TM. In the second column, we exclude book 4, which lies highest in energy, in a 48 × 48 TM model. In the third column, we also exclude book 3 and use a 32 × 32 TM. We note that the largest splittings w_A and w_B, due to the mechanisms involving Ã/B̃, are already converged in a two-structure (book 1 and 2) model. The TS pattern is a doublet of doublets, with widths of 1.7 × 10⁻³ cm⁻¹ and 2.1 × 10⁻⁴ cm⁻¹, respectively. It is qualitatively similar to the TS pattern obtained for the hexamer cage, but has an order of magnitude larger widths of the finer doublet splittings. The splittings due to bifurcations C̃/D̃ are an order of magnitude smaller than B̃. The inclusion of book 4 in the TM model brings an additional factor of ≈2 to their widths w_C/D.

Fig. 15 A qualitative diagram of the tunneling splitting pattern of the water hexamer book in its ground state. The diagram defines the widths w_A − w_D of the splittings due to the bifurcations of monomers A–D in Fig. 12. M-WKB estimates of the widths are given in Table 8. Levels are labeled using the C_2v and G₁₆ irreducible representations, as appropriate. See the text for details.

Table 8 Tunneling splitting widths w_A − w_D in the vibrational ground state of the water hexamer book, as defined by Fig. 15, in cm⁻¹. In the first set, columns 1–3, the tunneling matrix (TM) model includes all TM elements for flips and bifurcations listed in Table 7 involving book structures 1–4, 1–3, and 1–2, respectively. The second set, columns 1 and 2, additionally include the TM elements for rotations, rot A–D, in 4-structure and 2-structure TM models, respectively. See the text for details

w A	1.69 × 10^−3	1.66 × 10^−3	1.68 × 10^−3
w B	2.11 × 10^−4	2.06 × 10^−4	2.08 × 10^−4
w C	2.27 × 10^−5	1.12 × 10^−5	0
w D	8.37 × 10^−6	4.95 × 10^−6	0
w A	2.31 × 10^−3	2.30 × 10^−3
w B	3.42 × 10^−4	3.39 × 10^−4
w C	5.51 × 10^−5	0
w D	2.04 × 10^−5	0

---

It is not clear whether the exclusion of the TM elements for monomer rotations is justifiable since they do not clearly decompose into a combination of flips and a bifurcation, as they do in the hexamer cage. When the MAPs for rotations (using them as initial paths) are relaxed to MEPs, we find that ‘rot A’ decomposes into ÃB + B + A for all book structures 1–4. Similarly, ‘rot C’ decomposes into C̃D + D + C. On the other hand, the MEP for ‘rot D’ decomposes as D̃ + D in book structures 1–3, while ‘rot B’ does not decompose into flips/bifurcations. This suggests that only the TM elements for ‘rot B’ and ‘rot D’ should be included in TM, as the other mechanisms are already represented by the model, at least in part. Some additional flux through the dividing surface due to a direct mechanism of rotation, on top of that included by the two-step process of a flip and bifurcation, is expected, but cannot be quantified using the present approach. In the fourth column of Table 8, we list the widths obtained by including all TM elements for monomer rotations to expose the uncertainty introduced by the TM model (in addition to that of the M-WKB approximation). Widths w_C and w_D are more affected by inclusion of additional mechanisms (by a factor of ≈2.5), while w_A and w_B increase by 36% and 62%, respectively. In the fifth column, we list the corresponding widths in the two-structure (book 1 and 2) model. In general, the widths w_A − w_D are, to within 1%, given by the number in column 1 or column 4 in Table 8, depending on whether the particular monomer rotation A–D is excluded or not, respectively.

In summary, the splittings w_A and w_B are determined more reliably than the fine splittings due to bifurcations involving monomers C and D. Our estimates of w_C and w_D are exacerbated by the existence of book structure 5 and also by the fact that the TM elements that connect to the energetically higher book structures are larger in size. Our best estimate is that the TS pattern in the water hexamer book is a doublet of doublets. The larger width, due to Ã, is 1.7 × 10⁻³ cm⁻¹, where we disregarded the ‘rot A’ contribution, while the smaller splittings, due to B̃, are 3.4 × 10⁻⁴ cm⁻¹, where we included the ‘rot B’ mechanism. The finer splittings due to C̃ and D̃ are 2.3 × 10⁻⁵ cm⁻¹ and 2.0 × 10⁻⁵ cm⁻¹, where in the latter we included the rotation of monomer D. Finally, we note that the widths estimated using only the TM elements for monomer rotation in book 1 as w_X = 2|h_rot X¹¹|, for X = A–D, from Table 8, give results within a factor of ≈2–3 to our best estimates obtained above.

4 Conclusions

We applied the M-WKB method^62,63 in full dimensionality to calculate the TS patterns of energetically low-lying structural isomers of the water hexamer; in particular, the hexamer prism, the cage and the book. M-WKB is a semiclassical method, closely related to the instanton method, which gives equivalent results to it in the GS, but can also be used to calculate the TSs in the vibrationally excited states at little additional computational cost.

The GS TS pattern in the water hexamer prism has been determined experimentally and interpreted in terms of two rearrangement mechanisms, identified using instanton theory, in ref. 39. We calculated the TS patterns in the low-lying vibrationally excited states of the prism and identified an additional mechanism, the rotation of a double-donor monomer, that potentially plays a role in shaping of the TS patterns in the excited states. We find that there are no significant changes of the TS sizes in the excited modes 1 and 2 (numbered in order of their frequencies). In the excited modes 3–7, the TS sizes are reduced and the rotation of monomer C, in Fig. 1, competes with its TS size with the so-called geared double-flip mechanism of ref. 39. A significant enhancement (≈8×) of the overall TS width is found in the excited mode 10.

In the water hexamer cage, we determined the TS pattern in its GS, delocalized over four cage minima linked by torsional flips of the terminal OH bonds at the two opposite vertices of the cage, due to bifurcations. The TS pattern and the responsible mechanisms were first hypothesized in ref. 36 and 38. We find that the GS TS pattern is a doublet of doublets. The larger splitting is due to the bifurcation of monomer B, in Fig. 8, and is of the order of 10⁻³ cm⁻¹, while the smaller splitting, due to bifurcation of monomer A, is two orders of magnitude smaller. The TSs in the vibrationally excited states of the cage cannot be reliably estimated using the present method because the local harmonic energies, TM elements that link different cage structures, and the error associated with the M-WKB method used to determine them, are comparable in magnitude. Nevertheless, we find that, based on the TM element sizes for bifurcations in the excited states, it is difficult to reconcile the present findings with the experiment of ref. 36. The accidental degeneracy of the TM element sizes for the bifurcations of monomers A and B, used in the interpretation, is not reproduced here due to a large difference in the barrier heights along the associated rearrangement paths.

We also determined the TS pattern in the GS of the water hexamer book, delocalized over four minima linked by the double OH flips at four terminal water monomers, due to bifurcations. We found that the TS pattern is a doublet of doublets. The larger splitting is similar in size to that in the cage, while the smaller doublet splittings are an order of magnitude smaller. The responsible mechanisms are the monomer motions on one side of the book (monomers A and B in Fig. 12). The monomer motions on the other side of the book (monomer C and D in Fig. 12) cause further doublet splittings which are an order of magnitude smaller (≈10⁻⁵ cm⁻¹) and may play a role in the shaping of TS patterns in the vibrationally excited states.

The M-WKB method is expected to correctly predict the mechanisms responsible for the formation of TS patterns. In the systems with symmetry-related minima, such as the water hexamer prism, M-WKB results lie within a factor of two of the exact quantum methods on the same PES.⁵⁵ Previous studies found that it correctly predicts large enhancements of the TSs in the vibrationally excited states whenever they occur.⁶⁵ The study of the water hexamer cage and the book presents the first application of instanton theories to systems with symmetrically inequivalent minima. This brings an additional uncertainty in the results and prevents us from studying the TSs in the vibrationally excited states. Furthermore, low eigenfrequencies of the hexamer book imply a larger error associated with the anharmonicity of the potential energy in directions perpendicular to the rearrangement path. The associated error sizes are difficult to estimate, but the correct trends in the TS sizes and identifications of the responsible modes for the pattern formation are expected to hold. The use of perturbation theory⁸⁶ or higher-order WKB approaches⁸⁷ could potentially be employed in the future to estimate the errors associated with the anharmonicity in transverse directions relative to the rearrangement path. The minimum-action paths, determined here, could potentially facilitate in the development and application of sampling approaches, such as PIMD,⁵⁵ in more accurate treatments in the future.

Data availability

The data supporting this article have been included in part as part of the ESI.† A version of the code used to obtain the data in the article is available at the Digital Repository of the University of Zagreb at https://repozitorij.unizg.hr/islandora/object/pmf:12118.

Conflicts of interest

There are no conflicts to declare.

Appendix

Class structure and character table of the G₃₂ group of the water hexamer prism, which together define the irrep names, are given in Tables 9 and 10.

Table 9 Definitions of classes in the G₃₂ group of the water hexamer prism. The class name, number of elements in the class and a representative symmetry element of the class are given in columns 1–3, respectively

Class	Size	Element
1	1	E
2	2	(1 2)
3	2	(5 6)
4	2	(1 2)(5 6)
5	1	(1 2)(7 8)
6	2	(5 6)(7 8)
7	1	(5 6)(9 10)
8	2	(1 2)(5 6)(7 8)
9	2	(1 2)(5 6)(9 10)
10	1	(1 2)(5 6)(7 8)(9 10)
11	4	(A D)(B F)(C E)
		(1 8)(2 7)(3 11)(4 12)(5 9)(6 10)
12	4	(A D)(B F)(C E)
		(1 8 2 7)(3 11)(4 12)(5 9)(6 10)
13	4	(A D)(B F)(C E)
		(1 8)(2 7)(3 11)(4 12)(5 10 6 9)
14	4	(A D)(B F)(C E)
		(1 8 2 7)(3 11)(4 12)(5 10 6 9)

---

Table 10 Character table of the G₃₂ group of the water hexamer prism. Class names are defined in Table 9

	1	2	3	4	5	6	7	8	9	10	11	12	13	14
A 1a	1	1	1	1	1	1	1	1	1	1	1	1	1	1
A 1b	1	1	−1	−1	1	−1	1	−1	1	1	1	1	−1	−1
A 2a	1	−1	1	−1	1	−1	1	1	−1	1	−1	1	−1	1
A 2b	1	−1	−1	1	1	1	1	−1	−1	1	−1	1	1	−1
B 1a	1	−1	1	−1	1	−1	1	1	−1	1	1	−1	1	−1
B 1b	1	−1	−1	1	1	1	1	−1	−1	1	1	−1	−1	1
B 2a	1	1	1	1	1	1	1	1	1	1	−1	−1	−1	−1
B 2b	1	1	−1	−1	1	−1	1	−1	1	1	−1	−1	1	1
E 1a	2	0	2	0	−2	0	2	−2	0	−2	0	0	0	0
E 1b	2	0	−2	0	−2	0	2	2	0	−2	0	0	0	0
E 2	2	0	0	2	−2	−2	−2	0	0	2	0	0	0	0
E 3	2	0	0	−2	−2	2	−2	0	0	2	0	0	0	0
E 4	2	2	0	0	2	0	−2	0	−2	−2	0	0	0	0
E 5	2	−2	0	0	2	0	−2	0	2	−2	0	0	0	0

---

Acknowledgements

This work was fully supported by the Croatian Science Foundation Grant No. IP-2020-02-9932.

References

1. M. Ceriotti, J. Cuny, M. Parrinello and D. E. Manolopoulos, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 15591–15596.
2. M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie, A. Michaelides, M. A. Morales and T. E. Markland, Chem. Rev., 2016, 116, 7529–7550.
3. Y. Wang, X. Huang, B. C. Shepler, B. J. Braams and J. M. Bowman, J. Chem. Phys., 2011, 134, 094509.
4. K. Liu, J. D. Cruzan and R. J. Saykally, Science, 1996, 271, 929–932.
5. F. N. Keutsch and R. J. Saykally, Proc. Natl. Acad. Sci. U. S. A., 2001, 98, 10533–10540.
6. R. D. Suenram, G. T. Fraser and F. J. Lovas, J. Mol. Spectrosc., 1989, 138, 440.
7. G. T. Fraser, R. D. Suenram and L. H. Coudert, J. Chem. Phys., 1989, 90, 6077.
8. E. Zwart, J. J. ter Meulen, W. L. Meerts and L. H. Coudert, J. Mol. Spectrosc., 1991, 147, 27–39.
9. J. B. Paul, R. A. Provencal and R. J. Saykally, J. Phys. Chem. A, 1998, 102, 3279–3283.
10. T. R. Dyke, J. Chem. Phys., 1977, 66, 492.
11. R. Schwan, C. Qu, D. Mani, N. Pal, L. van der Meer, B. Redlich, C. Leforestier, J. M. Bowman, G. Schwaab and M. Havenith, Angew. Chem., 2019, 58, 13119–13126.
12. F. N. Keutsch, R. S. Fellers, M. R. Viant and R. J. Saykally, J. Chem. Phys., 2001, 114, 4005–4015.
13. M. R. Viant, J. D. Cruzan, D. D. Lucas, M. G. Brown, K. Liu and R. J. Saykally, J. Phys. Chem. A, 1997, 101, 9032–9041.
14. F. N. Keutsch, R. J. Saykally and D. J. Wales, J. Chem. Phys., 2002, 117, 8823–8835.
15. F. N. Keutsch, J. D. Cruzan and R. J. Saykally, Chem. Rev., 2003, 103, 2533–2578.
16. D. J. Wales, J. Am. Chem. Soc., 1993, 115, 11180–11190.
17. T. R. Walsh and D. J. Wales, J. Chem. Soc., Faraday Trans., 1996, 92, 2505–2517.
18. J. O. Richardson, S. C. Althorpe and D. J. Wales, J. Chem. Phys., 2011, 135, 124109.
19. J. D. Cruzan, L. B. Braly, K. Liu, M. G. Brown, J. G. Loeser and R. J. Saykally, Science, 1996, 271, 59–62.
20. J. D. Cruzan, M. R. Viant, M. G. Brown and R. J. Saykally, J. Phys. Chem. A, 1997, 101, 9022–9031.
21. D. J. Wales and T. R. Walsh, J. Chem. Phys., 1997, 106, 7193–7207.
22. M. G. Brown, F. N. Keutsch and R. J. Saykally, J. Chem. Phys., 1998, 109, 9645–9647.
23. H. A. Harker, M. R. Viant, F. N. Keutsch, E. A. Michael, R. P. McLaughlin and R. J. Saykally, J. Phys. Chem. A, 2005, 109, 6483–6497.
24. K. Liu, M. G. Brown, J. D. Cruzan and R. J. Saykally, J. Phys. Chem. A, 1997, 101, 9011–9021.
25. D. J. Wales and T. R. Walsh, J. Chem. Phys., 1996, 105, 6957.
26. M. T. Cvitaš and J. O. Richardson, Phys. Chem. Chem. Phys., 2019, 22, 1035–1044.
27. R. Schwan, C. Qu, D. Mani, N. Pal, G. Schwaab, J. M. Bowman, G. S. Tschumper and M. Havenith, Angew. Chem., 2020, 59, 11399–11407.
28. W. T. S. Cole, R. S. Fellers, M. R. Viant and R. J. Saykally, J. Chem. Phys., 2017, 146, 014306.
29. A. Rakshit, P. Bandyopadhyay, J. P. Heindel and S. S. Xantheas, J. Chem. Phys., 2019, 151, 214307.
30. D. M. Bates and G. S. Tschumper, J. Phys. Chem. A, 2009, 113, 3555–3559.
31. E. Miliordos, E. Aprà and S. S. Xantheas, J. Chem. Phys., 2013, 139, 114302.
32. K. Nauta and R. E. Miller, Science, 2000, 287, 293–295.
33. C. Pérez, M. T. Muckle, D. P. Zaleski, N. A. Seifert, B. Temelso, G. C. Shields, Z. Kisiel and B. H. Pate, Science, 2012, 336, 897–901.
34. Y. Wang, V. Babin, J. M. Bowman and F. Paesani, J. Am. Chem. Soc., 2012, 134, 11116–11119.
35. V. Babin and F. Paesani, Chem. Phys. Lett., 2013, 580, 1–8.
36. K. Liu, M. G. Brown, C. Carter, R. J. Saykally, J. K. Gregory and D. C. Clary, Nature, 1996, 381, 501–503.
37. K. Liu, M. G. Brown and R. J. Saykally, J. Phys. Chem. A, 1997, 101, 8995–9010.
38. D. J. Wales, in Theory of Atomic and Molecular Clusters: With a Glimpse at Experiments, ed. J. Jellinek, Springer-Verlag, Berlin, 1999, pp. 86–110.
39. 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, Science, 2016, 351, 1310–1313.
40. W. T. S. Cole, Ö. Yönder, A. A. Sheikh, R. S. Fellers, M. R. Viant, R. J. Saykally, J. D. Farrell and D. J. Wales, J. Phys. Chem. A, 2018, 122, 7421–7426.
41. W. T. S. Cole, J. D. Farrell, A. A. Sheikh, Ö. Yönder, R. S. Fellers, M. R. Viant, D. J. Wales and R. J. Saykally, J. Chem. Phys., 2018, 148, 094301.
42. C. Pérez, S. Lobsiger, N. A. Seifert, D. P. Zaleski, B. Temelso, G. C. Shields, Z. Kisiel and B. H. Pate, Chem. Phys. Lett., 2013, 571, 1–15.
43. J. O. Richardson, D. J. Wales, S. C. Althorpe, R. P. McLaughlin, M. R. Viant, O. Shih and R. J. Saykally, J. Phys. Chem. A, 2013, 117, 6960–6966.
44. C. Pérez, D. P. Zaleski, N. A. Seifert, B. Temelso, G. C. Shields, Z. Kisiel and B. H. Pate, Angew. Chem., 2014, 53, 14368–14372.
45. A. A. Reid, PhD thesis, University of Cambridge, 2014.
46. C. Leforestier, K. Szalewicz and A. van der Avoird, J. Chem. Phys., 2012, 137, 014305.
47. X.-G. Wang and T. Carrington Jr, J. Chem. Phys., 2018, 148, 074108.
48. D. Sabo, Z. Bačić, S. Graf and S. Leutwyler, Chem. Phys. Lett., 1996, 261, 318–328.
49. A. van der Avoird and K. Szalewicz, J. Chem. Phys., 2008, 128, 014302.
50. M. Losada and S. Leutwyler, J. Chem. Phys., 2003, 119, 304–312.
51. P. M. Felker, I. Simkó and Z. Bačić, J. Phys. Chem. A, 2024, 128, 8170–8189.
52. D. Blume and K. B. Whaley, J. Chem. Phys., 2000, 112, 2218–2226.
53. J. K. Gregory and D. C. Clary, J. Chem. Phys., 1995, 102, 7817.
54. J. K. Gregory and D. C. Clary, J. Chem. Phys., 1996, 105, 6626.
55. C. L. Vaillant, D. J. Wales and S. C. Althorpe, J. Phys. Chem. Lett., 2019, 10, 7300–7304.
56. Y.-C. Zhu, S. Yang, J.-X. Zeng, W. Fang, L. Jiang, D. H. Zhang and X.-Z. Li, J. Am. Chem. Soc., 2022, 144, 21356–21362.
57. M. Eraković and M. T. Cvitaš, Phys. Chem. Chem. Phys., 2021, 23, 4240–4254.
58. M. Eraković and M. T. Cvitaš, Phys. Chem. Chem. Phys., 2024, 26, 12965–12981.
59. G. V. Mil’nikov and H. Nakamura, J. Chem. Phys., 2001, 115, 6881–6897.
60. G. V. Mil’nikov and H. Nakamura, J. Chem. Phys., 2005, 122, 124311.
61. J. O. Richardson and S. C. Althorpe, J. Chem. Phys., 2011, 134, 054109.
62. M. Eraković, C. L. Vaillant and M. T. Cvitaš, J. Chem. Phys., 2020, 152, 084111.
63. M. Eraković and M. T. Cvitaš, J. Chem. Phys., 2020, 153, 134106.
64. M. T. Cvitaš and J. O. Richardson, in Molecular Spectroscopy and Quantum Dynamics, ed. R. Marquardt and M. Quack, Elsevier, 2020, ch. 10.
65. M. Eraković and M. T. Cvitaš, J. Chem. Phys., 2024, 160, 154112.
66. M. Eraković and M. T. Cvitaš, J. Chem. Theory Comput., 2022, 18, 2785–2802.
67. J. K. Gregory and D. C. Clary, J. Phys. Chem. A, 1997, 101, 6813–6819.
68. E. Mátyus, D. J. Wales and S. C. Althorpe, J. Chem. Phys., 2016, 144, 114108.
69. U. Góra, W. Cencek, R. Podeszwa, A. van der Avoird and K. Szalewicz, J. Chem. Phys., 2014, 140, 194101.
70. Y. Wang and J. M. Bowman, Chem. Phys. Lett., 2010, 491, 1–10.
71. Y. Wang, B. C. Shepler, B. J. Braams and J. M. Bowman, J. Chem. Phys., 2009, 131, 054511.
72. V. Babin, C. Leforestier and F. Paesani, J. Chem. Theory Comput., 2013, 9, 5395–5403.
73. V. Babin, G. R. Medders and F. Paesani, J. Chem. Theory Comput., 2014, 10, 1599–1607.
74. S. K. Reddy, S. C. Straight, P. Bajaj, C. Huy Pham, M. Riera, D. R. Moberg, M. A. Morales, C. Knight, A. W. Götz and F. Paesani, J. Chem. Phys., 2016, 145, 194504.
75. X. Zhu, M. Riera, E. F. Bull-Vulpe and F. Paesani, J. Chem. Theory Comput., 2023, 19, 3551–3566.
76. Q. Yu, C. Qu, P. L. Houston, R. Conte, A. Nandi and J. M. Bowman, J. Phys. Chem. Lett., 2022, 13, 5068–5074.
77. L. D. Landau and E. M. Lifshitz, Quantum Mechanics: Non-Relativistic Theory, Pergamon Press, Oxford, 2nd edn, 1965.
78. C. Herring, Rev. Mod. Phys., 1962, 34, 631–645.
79. A. I. Vainshtein, V. I. Zakharov, V. A. Novikov and M. A. Shifman, Sov. Phys.-Usp., 1982, 25, 195.
80. C. F. Karney, J. Mol. Graphics, 2007, 25, 595–604.
81. M. T. Cvitaš and S. C. Althorpe, J. Chem. Theory Comput., 2016, 12, 787–803.
82. M. T. Cvitaš, J. Chem. Theory Comput., 2018, 14, 1487–1500.
83. G. V. Mil’nikov, T. Ishida and H. Nakamura, J. Phys. Chem. A, 2006, 110, 5430–5435.
84. G. Mil’nikov, O. Kühn and H. Nakamura, J. Chem. Phys., 2005, 123, 074308.
85. J. D. Mallory and V. A. Mandelshtam, J. Chem. Phys., 2016, 145, 064308.
86. J. E. Lawrence, J. Dušek and J. O. Richardson, J. Chem. Phys., 2023, 159, 014111.
87. E. Pollak and J. Cao, J. Chem. Phys., 2022, 157, 074109.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5cp00155b

---