Tunneling splittings observed in molecular rovibrational spectra are significant evidence for tunneling motion of hydrogen nuclei in water clusters. Accurate calculations of the splitting sizes from first principles require a combination of high-quality inter-atomic interactions and rigorous methods to treat the nuclei with quantum mechanics. Many theoretical efforts have been made in recent decades. This Perspective focuses on two path-integral based tunneling splitting methods whose computational cost scales well with the system size, namely, the ring-polymer instanton method and the path-integral molecular dynamics (PIMD) method. From a simple derivation, we show that the former is a semiclassical approximation to the latter, despite that the two methods are derived very differently. Currently, the PIMD method is considered to be an ideal route to rigorously compute the ground-state tunneling splitting, while the instanton method sacrifices some accuracy for a significantly smaller computational cost. An application scenario of such a quantitatively rigorous calculation is to test and calibrate the potential energy surfaces of molecular systems by spectroscopic accuracy. Recent progress in water clusters is reviewed, and the current challenges are discussed.

Water is one of the most ubiquitous and vital substances on Earth. It is widely involved in chemical reactions and the cycles of matter and heat. The seemingly simple hydrogen atoms and oxygen atoms form relatively strong intermolecular hydrogen bonds and a complex hydrogen bond (HB) network. Therefore, water has many abnormal properties and abundant phases.1–5 Conventionally, these properties and phases were understood using the ball-and-stick model, i.e., electrons were treated using quantum mechanics, which results in chemical bonds, and the nuclei were treated as classical particles, which form the skeleton of the structures. In this scenario, difficulties for theoretical descriptions of water arise from the need to explore the large HB network and get accurate HB interactions. In recent years, advances in computer simulation methods mean that quantum mechanical descriptions of the nuclei are also possible,6–11 opening a door for studies of new exciting quantum effects of water. Assisted by high-resolution experimental techniques, such nuclear quantum effects (NQEs) had been found to be highly relevant to a series of nontrivial properties.12–21 The tunneling motion of hydrogen atoms, which involves the breaking and reforming of hydrogen bonds, is a prominent example.21–24 However, the condensed nature of bulk, which results in broadening of its vibrational and rotational spectra, makes it impossible for one to obtain high-resolution experimental data to decipher details of these processes in its bulk phases.

Water clusters provide a good entry point to circumvent this problem25–27 as its small size allows us to capture the specific hydrogen movement patterns without being obscured by the broadening caused by a large amount of chaotic movement like in the bulk phase. The high-resolution rovibrational spectrum of gas-phase water clusters can give much more precise and detailed information, where the splittings of the rovibrational energy levels reveal different tunneling behaviors of hydrogen atoms.28–35 These splittings have been investigated for decades. Permutation-inversion symmetry is useful for assigning the splittings.36 At very low temperatures, the geometry of the clusters can only be rearranged by tunneling motions, which results in equivalent conformations in different wells being connected through comparatively low barriers. Splitting of the energy levels happens because of couplings of these degenerate states.37 A molecular symmetry group can be generated by these permutation-inversion operation elements corresponding to the geometry rearrangements. Then, the splitting pattern can be described by the irreducible representations of this group, and the splitting sizes are determined by the ease of tunneling motions.38 

With the splitting sizes calculated, one can compare them with the experimental observations and analyze the atomic details of the tunneling paths. The accuracy for the value of the splitting sizes depends on two factors, i.e., the quality of the potential energy surface (PES) and the reliability of the method with which tunneling is described.39 Assuming that a rigorous and practical method is available to address the latter factor, experimental values of the splitting sizes can then be used as new calibrations for the improvement of the water potential.40,41 This is not only essential for simulations of bulk water in a more general sense but also in the cases when NQEs are of interest.41–46 Therefore, it seems that we have reached at a point in time when both the experiment and the first-principles based PES are almost ready, but there is still a lack of a rigorous and practical method to describe the quantum mechanical tunneling of the nuclei. With such a method, one can first calculate the ground-state splitting to understand the tunneling behavior more deeply and then analyze the finer structure of the infrared spectrum. In turn, the experimental splitting sizes serve as calibrations to further improve the quality of the water potential.40,41 Positive feedback from the quality of the PES to the theoretical method in describing tunneling, and vice versa, will push theoretical descriptions of water systems to the final stage with dreamed accuracy.

Along this route, several methods have been developed to find the tunneling paths and to calculate the splitting sizes. To date, only dimer splittings can be calculated by directly solving the Schrödinger equation in full dimensions47 due to the unfavorable scaling of the cost with system size. Diffusion Monte Carlo (DMC) is a promising method48,49 but suffers from the fixed-node approximation especially in the case of larger clusters and high-excited levels. The Wentzel–Kramers–Brillouin (WKB) method combined with the group theory is convenient in describing the splitting pattern50–53 but lacks accuracy in predicting the splitting sizes due to the approximations made in the method and the a priori specified tunneling path. The ring-polymer instanton method can avoid the latter problem by optimizing a minimum-action path as the optimal tunneling path in full dimensionality.54–57 It is also applicable to large water clusters, thanks to the relatively inexpensive computational cost and favorable scaling with system size. Now, it is becoming a mainstream method and has been used to calculate the finer structure of the splitting spectrum.58,59 However, the steepest-descent approximation made in the instanton method makes it short of spectroscopic accuracy for systems with anharmonic soft modes.35,56 For example, previous studies have exhibited departures from the experimental value within a factor of two and little dependency on the different water potentials used.35,56 To cure this problem, a path-integral molecular dynamics (PIMD) method making use of density matrix elements has recently been developed by Althorpe and co-workers.60–63 Their pioneering work showed that it can give an accurate splitting size of the acceptor path in the water dimer compared to benchmark results and showed promising results for larger water clusters as well.63 Yet, there was still room for improvement, e.g., resolving the relatively large sampling errors and other issues. Our recent theoretical work has reproduced the experimental values of the torsional tunneling splitting in a water trimer, which further demonstrated the potential of the PIMD method in providing the most accurate tunneling splitting results to date for large water clusters.64 

This Perspective focuses on these two path-integral based methods, i.e., ring-polymer instanton and PIMD. In Sec. II, we review the two methods in detail and then in Sec. III show that a connection exists between them, i.e., the instanton is exactly a steepest-descent approximation of the PIMD method although they originate from the quantum partition function and density matrix element, respectively. In Sec. IV, we introduce the progress on accurate calculations of the tunneling splittings in water clusters from the dimer to hexamer one by one and discuss the challenges in deciphering more finer structures of the spectra. Section V concludes our Perspective.

Here, we briefly introduce the tunneling splitting derived in the formalism of the ring-polymer instanton by Richardson and Althorpe, first in a simple two-degenerate-well system55 and then in a multiple-well system.56 

A molecular system with two degenerate wells has its ground-state energy levels E0 split by the mixing of the two states due to tunneling. In the low temperature limit, the contribution of higher energy levels can be neglected. This leads to
limβQ(β)2Q0(β)eβ(E0Δ/2)+eβ(E0+Δ/2)2eβE0=coshβΔ2,
(1)
where Q(β) (2Q0(β)) is the partition function of the system with (without) the inclusion of tunneling, E0 ±Δ/2 are the energy levels split by the tunneling, and β is 1/kBT. This expression provides a link between the tunneling splitting and the partition functions of the system.
Q0(β) is easy to evaluate. Applying the steepest-descent approximation, it corresponds to the situation of a collapsed ring-polymer in the bottom of one of the wells. This results in a simple harmonic vibrational partition function. With βN = β/N, N being the number of beads, one has
Q0(β)k1βNωk=1βNN1detG0.
(2)
ωk2 are the eigenvalues of the mass-weighted Hessian of the collapsed ring-polymer, the elements of which are
(G0)ii=2δiiδii1δii+1(βN)2+ωs2δii.
(3)
The first term comes from the spring interaction between the adjacent beads. ωs is the harmonic frequency of the well. The formulas in this subsection are for one-dimensional systems but can be easily generalized to realistic systems.

Unlike Q0, Q corresponds to the partition function of the entire double-well region with the barrier included due to the existence of tunneling. Again, applying the method of steepest descent to Q(β), one gets a saddle point in the ring-polymer space. This results in the instanton, i.e., a closed trajectory in the imaginary time that connects the two wells over a barrier with the least action. In the low temperature limit, with infinitely long imaginary time βℏ, the instanton trajectory spends most of its time staying in the wells and passes the barrier occasionally. Such a pass is called a “kink.” Each kink process requires a finite time, which can be ignored compared to the infinite total time. As will be seen in Sec. III, this finite time is of fixed size and can be defined as an imaginary “tunneling time.” The most simple trajectory involves two kinks to form a closed orbit. If no kink occurs, it will fall back to the situation of Q0. In the ring-polymer formalism where a large number of N beads are used to uniformly discretize the total imaginary time βℏ, most beads would stay in the wells and just a few on the barrier describing kink processes. One example is shown in Fig. 1. A finite number of M beads is needed for a “one kink” process, marked with hollow circles. The instanton trajectory in the figure contains four kinks, and each kink is exactly the same. In the gap between the kinks, it will stay for any length of imaginary time, which is equivalent to an arbitrary number of N1, N2, N3, and N4 beads in well a or b, represented by the dashed lines. This leads to N1 + N2 + N3 + N4 + 4M = N, where NM.

FIG. 1.

A schematic diagram of an instanton trajectory composed of four kinks in the ring-polymer formalism. The red curve represents the potential of a double-well system. The blue circles represent the beads connected by springs in the ring-polymer, forming a discrete instanton trajectory. The process from a to b through the barrier is called a kink, indicated by the M hollow circles. In order to be able to be seen clearly, we draw the four kinks separately, while, in fact, they are the same and completely overlapped. Between the kinks, a large number of beads stay in the wells indicated by the black dotted lines.

FIG. 1.

A schematic diagram of an instanton trajectory composed of four kinks in the ring-polymer formalism. The red curve represents the potential of a double-well system. The blue circles represent the beads connected by springs in the ring-polymer, forming a discrete instanton trajectory. The process from a to b through the barrier is called a kink, indicated by the M hollow circles. In order to be able to be seen clearly, we draw the four kinks separately, while, in fact, they are the same and completely overlapped. Between the kinks, a large number of beads stay in the wells indicated by the black dotted lines.

Close modal
As the beads in the wells do not contribute to the action, the more the kinks in a trajectory, the greater the action. In the low temperature limit, contributions from such multiple-kink trajectories as the one shown in Fig. 1 are non-negligible. Enumerating all such trajectories, one obtains
Q(β)=n=0,even2Nnn!Qn(β),
(4)
where Qn(β) represents the contribution from a trajectory containing n kinks. The prefactor is the number of ways of arranging n kinks in a ring-polymer of N beads. Given that N is large, CNnNn/n!. A trajectory can start at well a or well b, which contributes a factor of two. Only even numbers of n are included in the sum since otherwise the trajectory would not be closed. In the β limit, the kink processes will be fully separated by the long stay in wells, which allows us to treat these kinks as isolated and equivalent. Meanwhile, Qn(β) are to be divided by Q0 as shown in Eq. (1); the contribution of the beads in wells shall be offset. At last, Qn(β) is only related to n (the number of kinks) and the contribution of one kink (donated as θ) through
QnQ0=θn.
(5)
This θ can be calculated by a linear polymer composed of finite M beads connecting the two wells from a to b. The potential energy for this linear polymer is in a familiar form,
UM(r)=12(βN)2i=1M1(ri+1ri)2+i=1MV(ri).
(6)
We wrote mri as ri to shorten the notation. A one-dimensional model with mass m is considered for simplicity. With these treatments, the action reads
Skink=βNUM(r̃).
(7)
A stationary-state search algorithm can be used to find the minimum point r̃ of this action. The frequency of the wave-like mode, which corresponds to the case of a cyclic permutation of all beads moving along the barrier simultaneously (as shown in Fig. 2), tends to be zero when T → 0. Each kink has such a mode since they are assumed isolated.
FIG. 2.

Red points outline the energy of beads in the kink of the anti-geared tunneling in a water hexamer prism. Many beads collapse in the bottom of the wells with the same zero energy. The wave-like packet composed of beads on the barrier can move freely since it corresponds to a zero-frequency mode.

FIG. 2.

Red points outline the energy of beads in the kink of the anti-geared tunneling in a water hexamer prism. Many beads collapse in the bottom of the wells with the same zero energy. The wave-like packet composed of beads on the barrier can move freely since it corresponds to a zero-frequency mode.

Close modal
When applying the method of steepest-descent to Qn, these zero-frequency modes should not be integrated out by the Gaussian integral. They donate a factor of (βNSkink)n/2 to the partition function, and the resulting Qn reads55 
Qn(β)=12πβN2N2(βNSkink)n2k2πβNηk2enSkink/.
(8)
The first factor (2πβN2)N/2 comes from the integration of the kinetic energy term. ηk2 are the eigenvalues of the mass-weighed Hessian matrix G of the instanton. The prime indicates that the n zero frequencies are excluded from the product and replaced by (βNSkink)n/2. This Qn also reads
Qn(β)=1βNNnSkink2πn21detGenSkink/.
(9)
The fluctuation term detG describes the local curvature perpendicular to the instanton path, which comes from the Gaussian integral. It is the embodiment of the steepest-descent approximation. Because the kinks are separated and isolated by a large number of beads, the fluctuations of the beads in the wells can be offset by G0 through
detGdetG0=Φn,
(10)
where
Φ=detJdetJ012,
(11)
J is the mass-weighed Hessian of the linear polymer and J0 is the equivalent Hessian for the non-tunneling system with M beads collapsed in a well.
Finally, one has
θ=βNΦSkink2πeSkink/.
(12)
From Eqs. (4) and (5), the ratio of the partition functions is
Q(β)2Q0(β)=n=0,evenNnn!θn=cosh(Nθ).
(13)
β is always required in the derivation. Comparing with Eq. (1), one can get the final expression for tunneling splitting of a double-well system,
Δ=2βNθ=2ΦSkink2πeSkink/.
(14)
Let’s now consider tunneling splitting in a molecular system with G degenerate wells. Similar to the double-well system, mixing of these G states allowed by tunneling makes the ground-state energy E0 split into a set of G levels {Eν}. Equation (1) generalizes to
limβQ(β)GQ0(β)=limβ1Gν=1Geβ(EνE0).
(15)
Again, utilizing the method of steepest descent, the full system partition function Q(β) in the low temperature limit can be evaluated by the sum of the contribution of various instanton trajectories with different numbers of kinks. However, this time, the orbit does not locate on a single barrier. Any two of the G wells are connected by a unique tunneling path unless the barrier is very high to allow tunneling to occur. These tunneling paths can be marked as different kinds of kinks. Contribution of a closed instanton trajectory with n kinks must include all possible arrangements starting from a well and returning to this well through n kinks.
In other words, all kinds of kinks have the same opportunity to take part in the trajectories. If one defines Qn,ν as the contribution of the closed trajectories of n kinks that start and end at well ν, the partition function can be written as
limβQ(β)=ν=1Gn=0Nnn!Qn,ν.
(16)
To evaluate Qn,ν, a concept of connected graph, adjacency matrix A, is introduced. A kink connects two wells through a barrier. The element Aλμ is the number of such tunneling paths that connect well λ and well μ directly. In general, Aλμ is equal to 1 or 0; however, there are some special cases where there are multiple symmetrical paths connecting λ and μ.56 An example of such path is shown in Fig. 3. The red circles represent the degenerate wells. They are connected by solid lines if tunneling paths exist between them; then, the corresponding Aλμ would be non-zero. κ=1GAλκAκμ provides the number of connections composed by two kinks that start at λ and end at μ via an intermediate κ. Similarly, (An)λμ contains all the arrangement of n kinks connecting λ and μ. With the help of Nn/n!, all the closed n-kink trajectories that start and end at the same well ν can be found by (An)νν. The example in Fig. 3 shows a six-kink trajectory of the instanton. (A6)νν can help to traverse all possible connections of the kinks.
FIG. 3.

A schematic diagram of an instanton trajectory composed of six kinks in the ring-polymer formalism. The red circles represent the degenerate wells in a multi-well system. If connected by a solid line, a feasible tunneling path exists between the two wells, and the matrix element Aλμ is non-zero.

FIG. 3.

A schematic diagram of an instanton trajectory composed of six kinks in the ring-polymer formalism. The red circles represent the degenerate wells in a multi-well system. If connected by a solid line, a feasible tunneling path exists between the two wells, and the matrix element Aλμ is non-zero.

Close modal
After finding all the trajectories, one has to consider their contribution to the partition function. As in Eq. (5), the contribution of a trajectory is the product of the contributions of all n kinks. Here, it is assumed that there is no coupling between different kinds of kinks. One can take wells λ and μ as a double-well system, and θλμ has exactly the same form as θ in Eq. (12). Then, the contribution of a trajectory is the product of corresponding θλμ. With the help of adjacency matrix A, one can enumerate all possible sequences to obtain Qn,ν. A Hückel-type tunneling matrix W is defined for this purpose,
Wλμ=Aλμhλμ,
(17)
where hλμ = −θλμ/βN. Similar to Eq. (5), one has
Qn,νQ0=[(βNW)n]νν.
(18)
Combining with Eq. (16), one further gets
limβQ(β)GQ0(β)=ν=1Gn=0Nnn!1G[(βNW)n]νν=1Gtr[eβW].
(19)
The eigenvalues of W are the splittings EνE0.

The multi-dimensional integral in the partition function can be solved analytically through semiclassical approximation,65–67 i.e., the method of steepest descent. The computational cost is normally low. Errors of harmonic approximation, however, exist. It is natural to think of using molecular dynamics sampling to avoid errors of this harmonic approximation.

We start an introduction of the PIMD method from a density matrix element formalism,68,69 following the derivation of Mátyus et al.60 For a molecular Hamiltonian, the coordinate representation of its density matrix has a clear definition. By inserting a complete set of eigenfunctions ψn, it has the form
ρr,r;β=reβĤr=nψn(r)ψn*reβEn.
(20)
Consider a symmetric double-well system first, and mark the symmetry operation corresponding to rearrangement through tunneling as P̂. The two lowest-level eigenfunctions are symmetric and antisymmetric, respectively, under P̂, i.e., ψ0(P̂r)=ψ0(r) and ψ1(P̂r)=ψ1(r). In the low temperature limit, the higher levels can be neglected, and a ratio between the two density matrix elements is
limβρ(r,P̂r;β)ρ(r,r;β)|ψ0(r)|2eβE0|ψ1(r)|2eβE1|ψ0(r)|2eβE0+|ψ1(r)|2eβE1=tanhΔ2[ββ̄(r)],
(21)
where
β̄(r)=2Δlnψ1(r)ψ0(r).
(22)
The β̄(r) was interpreted as an imaginary “tunneling time.” Δ = E1E0 is the tunneling splitting to be determined.
In the density matrix element defined above, r is an arbitrary geometry as the molecular wavefunctions ψn spread in the coordinate space. The splitting size Δ, however, is independent of r and β. This means that one particular geometry is enough to determine the splitting in Eq. (21). The reactant r = a is a convenient and practical choice as P̂a=b, where a and b are the reactant and product geometry, respectively. They are stationary points in the two wells. A concise relationship between density matrix elements and the splitting size exists,
ρ(a,b;β)ρ(a,a;β)=tanhΔ2(ββ̄).
(23)
For brevity, the ratio ρ(a, bβ)/ρ(a, aβ) is labeled as I(β). The splitting size Δ can be calculated from I(β) with at least two different β as β̄ is also an unknown quantity. The extension to a multi-well system is easy by making use of the permutation-inversion symmetry group.61 Next, we shall concentrate on the calculation of I(β) using PIMD.
In a standard discretized path-integral formalism, the element ρ(a, b;β) is transformed into a multi-dimensional integral in a classical linear polymer phase space,
ρ(a,b;β)Adp1dpMdr1drMeβMHa,bM,
(24)
with the effective classical linear polymer Hamiltonian
Ha,bM=j=1fi=1Mpi,j22μi,j+i=1M112mjωM2(ri,jri+1,j)2+j=1f12mjωM2(ajr1,j)2+(rM,jbj)2+i=1MV(ri)+(V(a)+V(b))/2.
(25)
M beads are inserted between fixed end points a and b, with the coordinates r1, …, rM, each bead being a replica of the molecular system with f dimensions. Different from the convention that N is used to label the number of time slices in path-integral simulations, we use M deliberately since in the instanton method M is the number of time slices for the linear polymer. As will be shown later, the equivalence between the steepest-descent approximation of the path for ρ(a, b;β) in Eq. (24) and the linear polymer for one kink process in the instanton method is crucial to demonstrate that the instanton method is a steepest-descent approximation of the PIMD method. The imaginary time β is cut into M + 1 parts, βM = β/(M + 1) and ωM = 1/(βM). The pi,js are introduced by the M × f inverse Gauss integrals. μi,j can be arbitrarily chosen while a judicious choice can improve efficiency of sampling here as proposed in Ref. 62. Their contribution to the prefactor A will later be offset by the same prefactor in ρ(a, a;β).
Then, I(β) can be written as
ρ(a,b;β)ρ(a,a;β)dp1dpMdr1drMeβMHa,bMdp1dpMdr1drMeβMHa,aM,
(26)
which tends to be exact when M. Performing molecular dynamics sampling directly on this formula is not practical because a and b are far apart. It is natural then to think of dividing the space between a and b into sufficiently small segments. This strategy is equivalent to thermodynamic integration. A free energy function is defined as
F(λ,βM)=1βMlnρ(a,c(λ);β).
(27)
λ here denotes a reaction coordinate, ranging [0, 1], and geometry c(λ) connects c(0) = a and c(1) = b along a smooth path. Customarily, we take the instanton as the smooth path as it provides sensible orientations of the degenerate conformations by finding a minimum-action pathway. Then, one has
ρ(a,b;β)ρ(a,a;β)=exp[βM(F(1,β)F(0,β))]=exp[βMΔF].
(28)
The free energy difference can be calculated by dragging the end point c(λ) from a to b,
ΔF=01dλF(λ,β)λ=01dλHa,c(λ)Nλλ,
(29)
where ⟨⋯⟩λ represents a thermodynamic ensemble average of the linear polymer with an end point fixed at c(λ). The Hamiltonian Ha,c(λ)M is obtained by replacing b in Eq. (25) by c(λ).
It is shown that c(λ) is the only variable in the Hamiltonian that explicitly includes λ; therefore,
Ha,c(λ)Mλ=j=1f12V(c)cj+mjωM2(cjrM,j)cjλ.
(30)
The first term integrates to be zero since V(a) = V(b). Then, one has
ΔF=ωM201dλc̄λT[crMλ],
(31)
where c̄j=mjcj is used to shorten the notation. This is the final expression. c(λ) depends on the specific path chosen from a to b. The integral of λ from 0 to 1 can be implemented by numerical methods, such as Gauss–Legendre quadrature. The expectation values of rM at several λs are needed to evaluate ΔF in Eq. (31), as determined by
rM,jλ=dp1dpMdr1drMrM,jeβMHa,c(λ)Mdp1dpMdr1drMeβMHa,c(λ)M.
(32)

In the case of tunneling splitting calculations, the relationship between instanton theory and the PIMD method may be unclear because the former is derived from the partition function in Eq. (1) while the latter starts from the density matrix element in Eq. (21). In this section, through a simple derivation, it is shown that the instanton method is also a semiclassical approximation of the PIMD method in the calculation of tunneling splitting.

If a semiclassical approximation is made to the density matrix element ρ(a, b;β), we can get a stationary point in the Mf-dimensional linear polymer space corresponding to a fixed-end linear polymer connecting a and b over the barrier that minimizes the Euclidean action. It will be exactly the same as the kink trajectory. This equivalence is crucial in establishing the connection between the instanton and PIMD methods. After integrating out the momentum p1, …, pM, one has
ρ(a,b;β)12πβM2f(M+1)/2j=1fmj(M+1)/2×dr1drMeβMUa,bM.
(33)
Here, a large M is assumed. The prefactor would be offset by the same term in ρ(a, a;β). Ua,bM can be obtained by excluding the kinetic energy terms in Eq. (25). It is numerically equivalent to the potential energy in the instanton method in Eq. (6) with the only difference that fixed ends are used here. The Hessian matrix determined by Ua,bM, therefore, equals that of the instanton method.
Diagonalizing the mass-weighted Hessian matrix, one gets a zero-frequency mode, which describes a wave-like motion of beads where one bead moves forward and the adjacent bead takes its place. Beads on the barrier constitute a wave packet, which moves freely along the kink as there are many beads collapsed in a and b as shown in Fig. 2. In the instanton method, such permutation motion takes N steps to complete a cycle (N is the number of time slices of the trajectory). In the density-matrix based PIMD simulation, however, the ends of the linear polymer are fixed at a and b. If we assume that there are l beads on the barrier under the requirement that the action of the linear polymer is unchanged, only Ml steps of such permutation are allowed. Therefore, integration over this zero-frequency mode contributes a factor of (Ml)βMSkink. This is the only difference with the instanton method, while other terms such as action Skink and fluctuation Φ are all the same. Finally, we label the semiclassical approximation of the ratio of the density matrix elements I(β) as ISC,
ISC=(Ml)βMSkinkk2πβMηk2k2πβMωk2eSkink/=(Ml)βMΦSkink2πeSkink/,
(34)
where ηk2s are the eigenvalues of the mass-weighed Hessian matrix J of the linear polymer that stands for ρ(a, b;β). ωk2s are those of J0 from ρ(a, a;β). The prime indicates that the zero-frequency mode is omitted from the product.
As beads are distributed at equal imaginary time intervals βM, we further define the imaginary time experienced by these l beads on the barrier as the imaginary tunneling time, β̄inst=lβM. Then,
ISC=(ββ̄inst)ΦSkink2πeSkink/.
(35)
Making use of Eq. (23) and that tanh(x) ≃ x when x is small (one can arrive at the same conclusion without making the tanh(x) ≃ x approximation), we have
ΔSC2ΦSkink2πeSkink/.
(36)
Comparing with the Δinst in Eq. (14), we can get
ΔSC=Δinst.
(37)
This means that by applying the semiclassical approximation to the density matrix element, we will get the same result as in the instanton method. From this, the point that the instanton method for tunneling splitting is a semiclassical approximation of the PIMD method is clear. In addition, we get a by-product β̄inst. If using the half-height width of the barrier to make a rough measure of β̄inst as in Fig. 2, we get values consistent with the PIMD method in the order of magnitude.

In the rovibrational spectra of gas-phase water clusters from microwave and far-infrared experiments, splittings of the rovibrational levels were detected, which revealed profoundly the wave nature and quantum dynamics of nuclei in water. These splittings, coming from the couplings between degenerate wells as introduced in Sec. I, revealed a more complex molecular symmetry36 than an ordinary system with a rigid structure. In the literature, the study of tunneling splitting was accompanied by but more challenging than the resolution of the structure of water clusters.25,70–72 From a structure in the ball-and-stick form, a complete permutation-inversion group36,73 could be constructed, which traverses all potential tunneling rearrangements. The tunneling motions that are important to the splittings can be mainly categorized as flipping (torsional tunneling), geared (anti-geared) interchange, and bifurcation tunneling.74 “Interchange” refers to when two water molecules in an H-bond change their roles as the donor and acceptor. Early calculations use ab initio electronic structures to estimate the barriers of different tunneling pathways, which were competent to identify the possible tunneling rearrangements responsible for the observed splitting. The estimations of splitting sizes were based on the WKB approximation or low-dimensional models.50,75,76 Consistent with intuition, breaking fewer H-bonds leads to a lower barrier and hence a larger splitting size. Flipping tunnelings do not break H-bonds and tend to produce the largest splittings. Although numerically rough, the splitting patterns could be explained successfully by the permutation-inversion group generated from the most possible rearrangements.

With the advancement in computing power, the instanton, as one of the semiclassical methods, has achieved more precise tunneling pathways and, hence, numerical results that are more convincing.56 The complex tunneling dynamics and splitting patterns in water clusters have been reviewed by Cvitaš and Richardson.74 The instanton not only represents the most appropriate path but also gives results accurate in the order of magnitude. The sources of errors are then gradually clarified. There is now an opportunity to build on the instanton to enable more rigorous calculations and decipher the finer structure of the vibration–rotation–tunneling (VRT) spectra. As one of the most successful attempts, the torsional tunneling splitting in water trimer was the first to be reproduced rigorously with an unprecedentedly high accuracy using the PIMD method.64 However, for other clusters, various problems still exist and block the way of accurate calculation, such as the rotation– and vibration–tunneling couplings. The current two theories, i.e., instanton and PIMD, focus on calculation of the ground-state tunneling splittings. Some attempts have been made to include these couplings as will be introduced later. We review recent progress and discuss the challenges by enumerating the water clusters from the dimer to hexamer.

A water dimer is the smallest and earliest studied water cluster.28,77 It has one hydrogen bond with the two water molecules as donor and acceptor. A full 12-dimensional wavefunction method has been performed recently where the complete molecular Hamiltonian (including permutation-inversion symmetry) was solved by a kind of discrete variable representation method, and the whole VRT levels can be obtained.47 Using the ring-polymer instanton, up to five different tunneling pathways have been found between eight degenerate configurations.56 The corresponding five permutation symmetry operations can be generated from two basic ones, i.e., the acceptor tunneling and geared tunneling. A combination of these two operations can generate a doublet-of-triplet splitting.

Among them, the acceptor tunneling contributes a large doublet splitting as it involves no breaking of hydrogen bonds and thus has a small barrier.56,73,77,78 This tunneling is special in the sense that it is constituted by a tunneling path with barrier and a barrierless rotation operation, as demonstrated in Fig. 4. The swap of the two hydrogens of the acceptor [Fig. 4(a)] has to be coupled with the overall rotation around the inertia spindle [Fig. 4(b)] in order to finish the permutation [Fig. 4(c)]. A recent work has calculated the acceptor tunneling splitting using the PIMD method with the thermodynamics integration along an instanton path and obtained a value in excellent agreement with the variation method (diffusion Monte Carlo method).63 

FIG. 4.

The instanton pathway (a) of the acceptor tunneling has to couple with an overall rotation (b) to constitute a pure (AB) permutation operation (c). If taking twice the instanton operation without the rotation, the dimer will not return to the original orientation. ABCD label the four distinct hydrogen atoms. (d) is a schematic of the energy barrier experienced by the (AB) permutation.

FIG. 4.

The instanton pathway (a) of the acceptor tunneling has to couple with an overall rotation (b) to constitute a pure (AB) permutation operation (c). If taking twice the instanton operation without the rotation, the dimer will not return to the original orientation. ABCD label the four distinct hydrogen atoms. (d) is a schematic of the energy barrier experienced by the (AB) permutation.

Close modal

We noticed that a symmetry factor of 2 is applied to the PIMD result in the previous work based on the reasoning that there are two equivalent instantons in the system.62 Yet, according to our experience, such symmetry factor may be inappropriate for the PIMD method since the ratio of density matrix elements does not depend on the path of TI. Hence, the splitting size could be erroneously magnified twice, depending on the ergodicity of the sampling in the PIMD simulations. Due to the neglect of the rotation–tunneling coupling in the current PIMD method, we suspect that the PIMD result might not agree spot on with the benchmark for the water dimer. An attempt has been made to evaluate this coupling.58 They proposed an extension to instanton theory, inspired by the derivation of Mátyus and Althorpe in a complete SO(3) symmetry,61 and qualitatively evaluated the acceptor splitting of the first several rotational-excited levels in the water dimer. More theoretical efforts on deciphering this finer structure are desired.

Quantitative reproduction of experimental trimer splittings has long been viewed as an important criterion for testing the quality of the 3-body interactions of a water potential.79–81 To achieve this, a rigorous full-dimensional treatment is required, and on top of this, the accuracy depends on the quality of the water potential. The torsional tunneling splitting in a water trimer is probably the first one in all clusters that have been calculated rigorously from first principles with spectroscopic accuracy.64 It marks significant progress of the water potential and demonstrates the high predictive power of the PIMD + TI method.

A water trimer has a ring structure with three hydrogen bonds in the oxygen plane and three unbonded hydrogens up or down the oxygen plane.29,82 Flipping of one of the unbonded hydrogens to the other side of the plane leads to an equivalent conformation while experiencing a very low barrier. Flipping the three unbonded hydrogens one by one results in six permutationally equivalent conformations [Fig. 5(a)]. The flip (or torsional) tunneling operation generates a C6 permutation-inversion symmetry group, which then gives a quartet splitting with a 1, 2, 2, 1 degeneracy. Due to the low barrier, the splitting size is large (up to tens of wavenumbers). The rotational excitations are several orders of magnitude smaller than the torsional splitting and display as a fine structure in the spectrum. Fortunately, little rotational–tunneling coupling and vibrational–tunneling coupling have been observed, therefore, the ground-state splitting could then be inferred from the transitions observed in the VRT spectra. Extensive and detailed experimental measurements gave complete level information and a precise splitting size.29,33,82–84

FIG. 5.

The torsional tunneling splitting in a water trimer. (a) shows the six degenerate conformations. hi (i = 1, 2, 3) indicate adjacent and nonadjacent interactions between the six wells. (b) and (c) show that the calculated splitting sizes are in good agreement with the experiments. The lowest level of the quartet is set as the zero point. This figure is reprinted from Ref. 64 with permission. Copyright 2000 American Chemical Society.

FIG. 5.

The torsional tunneling splitting in a water trimer. (a) shows the six degenerate conformations. hi (i = 1, 2, 3) indicate adjacent and nonadjacent interactions between the six wells. (b) and (c) show that the calculated splitting sizes are in good agreement with the experiments. The lowest level of the quartet is set as the zero point. This figure is reprinted from Ref. 64 with permission. Copyright 2000 American Chemical Society.

Close modal

From the theoretical perspective, however, with nine more dimensions than dimer, the water trimer is very challenging to be treated with a full 21-dimensional wavefunction method.85 Semiclassical methods on the other hand are not quantitative enough for the task. The PIMD method here plays an important role to enable a full-dimensional and rigorous calculation.

Our recent theoretical work accomplished a reproduction of the experimental values with an unprecedented accuracy of just one wavenumber from first principles64 [Figs. 5(b) and 5(c)]. The Hückel approximation (considering only adjacent interactions) was shown to be no longer valid, which means that sticking to the tunneling matrix in Sec. II B of the present Perspective would encounter inaccuracy in a quantitative calculation. To achieve high accuracy, the nonadjacent terms that represent a delocalization of the Wannier-style nuclear wave packets were included in the torsional Hamiltonian as performed in Ref. 64. These nonadjacent interactions contribute the shift of the split levels from an even 1:2:1 spacing. Such effect may also exist in other multi-well systems with low barriers, e.g., cyclic water pentamer.

The water tetramer is an ideal platform to test the four-body term in a many-body water potential. It has a ring structure like the water trimer and four unbonded hydrogens up or down the oxygen plane alternatively34 (Fig. 6). To construct an equivalent conformation, the four unbonded hydrogens must flip together to opposite sides of the plane. The tunneling motion is between only two symmetric wells, which is consistent with the doublet splitting observed in the VRT spectra.34,86 Apart from this, no other type of tunneling splitting has been observed. This doublet splitting is much smaller in size than the torsional one in the water trimer (2.26 GHz for (H2O)4 and 5.6 MHz for (D2O)4)34,86 due to the fact that although no hydrogen bonds are broken, all the four unbonded hydrogen atoms are involved in the tunneling. The tunneling mechanism is still undetermined due to complexity in the multi-proton transfer process.49,86 There are questions on whether the four hydrogen atoms go through a concerted tunneling pathway or a step-wise one and whether the intermediate state contributes or not. Fortunately, the PIMD + TI method should not be affected by uncertainties in the path since the result of TI is path-independent, and an ongoing work of our group is aimed to answer this question.

FIG. 6.

Schematic of the tunneling splitting in a water tetramer. The observed doublet (the blue arrows) in spectra may be a sum of the splitting of the ground state and a vibrational-excited state.

FIG. 6.

Schematic of the tunneling splitting in a water tetramer. The observed doublet (the blue arrows) in spectra may be a sum of the splitting of the ground state and a vibrational-excited state.

Close modal

Another problem of water tetramer is the possible vibration–tunneling couplings. Because the splittings were observed in low-energy vibrational bands, there has been a possibility that the observed doublets are the sum or difference of splittings of the ground state and a vibrational-excited state86 as shown in Fig. 6. The vibration–tunneling couplings were considered to be weak because the vibration was supposed as a ring-deformation motion of the bonded hydrogens, which would have little influence on the flip motion.86 To achieve the goal of taking spectra data as the calibrator, methods that can accurately calculate the tunneling splittings of vibrational-excited states are desired. Early theoretical attempts were made in model systems.87 For molecular systems, an effort has been made recently to build a semiclassical method based on the instanton from a Hamilton–Jacobi formalism.59 Besides this, deciding whether the doublets are ground-state splittings or the sum still requires more evidence.

The cyclic water pentamer was supposed to have a similar torsional tunneling motion like the water trimer.88,89 Flipping the five unbonded hydrogens in sequence could result in ten equivalent conformations that split the rovibrational state into six levels with 1, 2, 2, 2, 2, 1 degeneracy. In addition, the splitting size was considered to be large up to tens of wavenumber as the water trimer because of the low barrier of torsional motion. Several bands in the VRT spectra have been identified as transitions between six levels,90 but the complete information of all the splitting levels is still lacking to our knowledge. Theoretically, Cvitaš and Richardson computed the ground-state splittings using instanton theory91 with a tunneling matrix of 320 equivalent minima, including the torsional and bifurcation tunneling motions.88 This effort might be helpful in the continuing effort to assign the experimental spectrum. However, the nonadjacent interactions should be presented in the torsional tunneling splitting of the water pentamer as well,64 which were not included in the tunneling matrix of the instanton method and will impact the ratio of the spacings between the torsional splitting levels. Therefore, a prediction from a high-quality calculation using the PIMD + TI method is desirable, yet such calculation would require an ab initio water potential with highly accurate four-body interactions (and perhaps accurate 5-body interactions may also be required).

The water hexamer has several possible structures that were considered to have similar energy. During the identification of the structures, ticks of tunneling splitting have been found in the pure rotational spectra,35,71 showing a unique doublet-of-triplet splitting pattern. In the prism structure, the interchange tunneling pathways have been identified to be responsible for the splittings; the anti-geared mechanism breaks one H-bond, while the geared mechanism concertedly breaks two. Using ring-polymer instanton theory, Richardson et al. discovered that the geared mechanism plays a crucial role in producing the unique splitting pattern. The absolute size of the splittings predicted by instanton theory is larger than the experimental result due to the semi-quantitative nature of the method and the accuracy of the water potentials available at the time.35 Vaillant et al. made a first PIMD calculation and obtained splitting sizes closer to the experimental results with the same water potential.63 However, the sampling uncertainties remain large, and from our experience, fully converging the PIMD + TI for the hexamer is very challenging due to the relatively high barrier and the complexity of the tunneling pathway. Improving the sampling efficiency is a prior task for accurate calculations of the tunneling splittings in large-size water clusters. On the other hand, the experiment in Ref. 35 observed that the transitions between different pairs of rotational levels exhibit different splitting sizes. This indicates that rotation–tunneling coupling is non-negligible in the hexamer, and thus, the exact ground-state splitting size has not been precisely derived. The rotation–tunneling coupling may be the next key point to be considered in future theoretical calculations.

Before concluding, we discuss a key issue, which is the impact of the water potential energy surface on the tunneling splittings. Previous studies have shown that the use of different water potential energy surfaces yields quantitatively different results. Specifically, in the case of the water trimer, a comparison was made between the MB-pol and “Our-pol” water potential energy surfaces.64 PIMD simulations using the former yielded I1(β = 8000) = 0.58(2) and h1 = 26(2) cm−1 while using the latter produced I1(β = 8000) = 0.520(3) and h1 = 22.1(2) cm−1. The tunneling splitting value from “Our-pol” was much closer to the benchmark obtained from spectra data, which was 21.7 cm−1. In a study of the water hexamer prism, Richardson et al. compared the tunneling splittings obtained using the HBB2-pol92 and MB-pol water potential energy surfaces.35 They found that the results obtained using the two water potentials could differ by 15%–30%. This disparity implies that the accuracy of the potential energy surfaces can be effectively evaluated by calculating tunneling splittings, provided that a rigorous treatment of the nucleus is achieved. Tunneling splittings can be used as the most stringent benchmarks to evaluate the accuracy of the potential energy surfaces.

In this Perspective, we reviewed the instanton and the PIMD methods on calculating ground-state tunneling splittings in water clusters. A clear connection is established between the two methods through a simple derivation. Instanton theory provides a computationally efficient way to predict tunneling splitting for larger systems, yet the splitting size it predicts is semi-quantitative due to the steepest-descent approximation made. The PIMD method improves the accuracy of calculation results by correcting this error, especially for systems with low barriers and strongly anharmonic modes.

The ability of the PIMD method to calculate rigorously and accurately the ground-state tunneling splitting has been illustrated in a water trimer. However, for other water clusters, before one can compare with the spectra data, more efforts are required. Each rovibrational level of an isolated molecular system should be split by the tunneling motions with the splitting size reflecting the ease of tunneling. The transitions observed in spectra are between rotational- or vibrational-excited levels. If the system exhibits a strong rotation– or vibration–tunneling coupling, it would be hard to figure out the splitting size of the ground state experimentally. We listed the current issues including the vibration–tunneling coupling in a water tetramer and rotation–tunneling coupling in a water dimer and hexamer prism. Some recent progress was summarized. The target of taking the spectra data of these clusters as calibrators of an ab initio water potential is challenging, and more theoretical efforts on solving rigorously the tunneling splittings of rovibrational-excited states are required. Besides, for large-size clusters and complex tunneling motions, there is still room for improving the sampling efficiency to enhance reliability and reduce computational cost. With the support of the two path-integral based methods and further developments in the methodology, we believe that accurate calculation of tunneling splittings in water clusters is possible, and we hope this also will facilitate constructions of water potentials with spectroscopic accuracy in the near future.

We thank Dr. Christophe Vaillant and Professor Stuart Althorpe for their guiding help on the PIMD method. The authors are supported by the National Science Foundation of China under Grant Nos. 11934003, 12234001, and 22288201; the National Basic Research Programs of China under Grand Nos. 2021YFA1400503 and 2022YFA1403500; the Beijing Natural Science Foundation under Grant No. Z200004; and the Strategic Priority Research Program of the Chinese Academy of Sciences under Grant No. XDB33010400. The computational resources were provided by the supercomputer center at Peking University, China.

The authors have no conflicts to disclose.

Yu-Cheng Zhu: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Writing – original draft (lead). Shuo Yang: Data curation (equal); Software (equal). Jia-Xi Zeng: Methodology (equal); Validation (equal). Wei Fang: Funding acquisition (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). Ling Jiang: Writing – review & editing (equal). Dong H. Zhang: Funding acquisition (equal); Project administration (equal). Xin-Zheng Li: Funding acquisition (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (lead).

The data that support the findings of this study are available within the article.

1.
Water: A Comprehensive Treatise
, edited by
F.
Franks
(
Plenum
,
New York
,
1982
), Vol.
1
.
2.
F. H.
Stillinger
, “
Water revisited
,”
Science
209
,
451
(
1980
).
3.
O.
Mishima
and
H. E.
Stanley
, “
The relationship between liquid, supercooled and glassy water
,”
Nature
396
,
329
(
1998
).
4.
P.
Ball
and
F. H.
Stillinger
, “
H2O: A biography of water
,”
Nature
401
,
850
(
1999
).
5.
C. G.
Salzmann
, “
Advances in the experimental exploration of water’s phase diagram
,”
J. Chem. Phys.
150
,
060901
(
2019
).
6.
D.
Marx
and
M.
Parrinello
, “
Ab initio path integral molecular dynamics: Basic ideas
,”
J. Chem. Phys.
104
,
4077
(
1996
).
7.
T. E.
Markland
and
D. E.
Manolopoulos
, “
An efficient ring polymer contraction scheme for imaginary time path integral simulations
,”
J. Chem. Phys.
129
,
024105
(
2008
).
8.
M.
Ceriotti
,
M.
Parrinello
,
T. E.
Markland
, and
D. E.
Manolopoulos
, “
Efficient stochastic thermostatting of path integral molecular dynamics
,”
J. Chem. Phys.
133
,
124104
(
2010
).
9.
M.
Ceriotti
and
T. E.
Markland
, “
Efficient methods and practical guidelines for simulating isotope effects
,”
J. Chem. Phys.
138
,
014112
(
2013
).
10.
S.
Habershon
,
D. E.
Manolopoulos
,
T. E.
Markland
, and
T. F.
Miller
III
, “
Ring-polymer molecular dynamics: Quantum effects in chemical dynamics from classical trajectories in an extended phase space
,”
Annu. Rev. Phys. Chem.
64
,
387
(
2013
).
11.
A. A.
Hassanali
,
J.
Cuny
,
V.
Verdolino
, and
M.
Parrinello
, “
Aqueous solutions: State of the art in ab initio molecular dynamics
,”
Philos. Trans. R. Soc., A
372
,
20120482
(
2014
).
12.
B.
Chen
,
I.
Ivanov
,
M. L.
Klein
, and
M.
Parrinello
, “
Hydrogen bonding in water
,”
Phys. Rev. Lett.
91
,
215503
(
2003
).
13.
D.
Marx
, “
Proton transfer 200 years after von grotthuss: Insights from ab initio simulations
,”
ChemPhysChem
7
,
1848
(
2006
).
14.
D.
Marx
,
A.
Chandra
, and
M. E.
Tuckerman
, “
Aqueous basic solutions: Hydroxide solvation, structural diffusion, and comparison to the hydrated proton
,”
Chem. Rev.
110
,
2174
(
2010
).
15.
J. A.
Morrone
and
R.
Car
, “
Nuclear quantum effects in water
,”
Phys. Rev. Lett.
101
,
017801
(
2008
).
16.
F.
Paesani
and
G. A.
Voth
, “
The properties of water: Insights from quantum simulations
,”
J. Phys. Chem. B
113
,
5702
(
2009
).
17.
R. H.
McKenzie
,
C.
Bekker
,
B.
Athokpam
, and
S. G.
Ramesh
, “
Effect of quantum nuclear motion on hydrogen bonding
,”
J. Chem. Phys.
140
,
174508
(
2014
).
18.
X.-Z.
Li
,
B.
Walker
, and
A.
Michaelides
, “
Quantum nature of the hydrogen bond
,”
Proc. Natl. Acad. Sci. U. S. A.
108
,
6369
(
2011
).
19.
M.
Ceriotti
,
J.
Cuny
,
M.
Parrinello
, and
D. E.
Manolopoulos
, “
Nuclear quantum effects and hydrogen bond fluctuations in water
,”
Proc. Natl. Acad. Sci. U. S. A.
110
,
15591
(
2013
).
20.
M.
Ceriotti
,
W.
Fang
,
P. G.
Kusalik
,
R. H.
McKenzie
,
A.
Michaelides
,
M. A.
Morales
, and
T. E.
Markland
, “
Nuclear quantum effects in water and aqueous systems: Experiment, theory, and current challenges
,”
Chem. Rev.
116
,
7529
(
2016
).
21.
X.
Meng
,
J.
Guo
,
J.
Peng
,
J.
Chen
,
Z.
Wang
,
J.-R.
Shi
,
X.-Z.
Li
,
E.-G.
Wang
, and
Y.
Jiang
, “
Direct visualization of concerted proton tunnelling in a water nanocluster
,”
Nat. Phys.
11
,
235
(
2015
).
22.
L. E.
Bove
,
S.
Klotz
,
A.
Paciaroni
, and
F.
Sacchetti
, “
Anomalous proton dynamics in ice at low temperatures
,”
Phys. Rev. Lett.
103
,
165901
(
2009
).
23.
C.
Drechsel-Grau
and
D.
Marx
, “
Quantum simulation of collective proton tunneling in hexagonal ice crystals
,”
Phys. Rev. Lett.
112
,
148302
(
2014
).
24.
L.
Lin
,
J. A.
Morrone
, and
R.
Car
, “
Correlated tunneling in hydrogen bonds
,”
J. Stat. Phys.
145
,
365
(
2011
).
25.
K.
Liu
,
J. D.
Cruzan
, and
R. J.
Saykally
, “
Water clusters
,”
Science
271
,
929
(
1996
).
26.
F. N.
Keutsch
and
R. J.
Saykally
, “
Water clusters: Untangling the mysteries of the liquid, one molecule at a time
,”
Proc. Natl. Acad. Sci. U. S. A.
98
,
10533
(
2001
).
27.
R.
Ludwig
, “
Water: From clusters to the bulk
,”
Angew. Chem., Int. Ed.
40
,
1808
(
2001
).
28.
E.
Zwart
,
J. J.
Ter Meulen
,
W. L.
Meerts
, and
L. H.
Coudert
, “
The submillimeter rotation tunneling spectrum of the water dimer
,”
J. Mol. Spectrosc.
147
,
27
(
1991
).
29.
N.
Pugliano
and
R.
Saykally
, “
Measurement of quantum tunneling between chiral isomers of the cyclic water trimer
,”
Science
257
,
1937
(
1992
).
30.
F. N.
Keutsch
,
N.
Goldman
,
H. A.
Harker
,
C.
Leforestier
, and
R. J.
Saykally
, “
Complete characterization of the water dimer vibrational ground state and testing the VRT (ASP-W) III, SAPT-5st, and VRT (MCY-5f) surfaces
,”
Mol. Phys.
101
,
3477
(
2003
).
31.
F. N.
Keutsch
,
R. J.
Saykally
, and
D. J.
Wales
, “
Bifurcation tunneling dynamics in the water trimer
,”
J. Chem. Phys.
117
,
8823
(
2002
).
32.
M. G.
Brown
,
F. N.
Keutsch
, and
R. J.
Saykally
, “
The bifurcation rearrangement in cyclic water clusters: Breaking and making hydrogen bonds
,”
J. Chem. Phys.
109
,
9645
(
1998
).
33.
K.
Liu
,
M. J.
Elrod
,
J. G.
Loeser
,
J. D.
Cruzan
,
N.
Pugliano
,
M. G.
Brown
,
J.
Rzepiela
, and
R. J.
Saykally
, “
Far-IR vibration–rotation–tunnelling spectroscopy of the water trimer
,”
Faraday Discuss.
97
,
35
(
1994
).
34.
J. D.
Cruzan
,
L. B.
Braly
,
K.
Liu
,
M. G.
Brown
,
J. G.
Loeser
, and
R. J.
Saykally
, “
Quantifying hydrogen bond cooperativity in water: VRT spectroscopy of the water tetramer
,”
Science
271
,
59
(
1996
).
35.
J. O.
Richardson
,
C.
Pérez
,
S.
Lobsiger
,
A. A.
Reid
,
B.
Temelso
,
G. C.
Shields
,
Z.
Kisiel
,
D. J.
Wales
,
B. H.
Pate
, and
S. C.
Althorpe
, “
Concerted hydrogen-bond breaking by quantum tunneling in the water hexamer prism
,”
Science
351
,
1310
(
2016
).
36.
H. C.
Longuet-Higgins
, “
The symmetry groups of non-rigid molecules
,”
Mol. Phys.
6
,
445
(
1963
).
37.
V. A.
Benderskii
,
D. E.
Makarov
, and
C. A.
Wight
,
Chemical Dynamics at Low Temperatures
(
John Wiley and Sons
,
2009
), Vol.
188
.
38.
D. J.
Wales
, “
Some estimates of tunneling splittings in small clusters
,”
J. Am. Chem. Soc.
115
,
11191
(
1993
).
39.
R. S.
Fellers
,
L. B.
Braly
,
R. J.
Saykally
, and
C.
Leforestier
, “
Fully coupled six-dimensional calculations of the water dimer vibration-rotation-tunneling states with split Wigner pseudospectral approach. II. Improvements and tests of additional potentials
,”
J. Chem. Phys.
110
,
6306
(
1999
).
40.
N.
Goldman
,
R. S.
Fellers
,
M. G.
Brown
,
L. B.
Braly
,
C. J.
Keoshian
,
C.
Leforestier
, and
R. J.
Saykally
, “
Spectroscopic determination of the water dimer intermolecular potential-energy surface
,”
J. Chem. Phys.
116
,
10148
(
2002
).
41.
V.
Babin
,
G. R.
Medders
, and
F.
Paesani
, “
Development of a ‘first principles’ water potential with flexible monomers. II: Trimer potential energy surface, third virial coefficient, and small clusters
,”
J. Chem. Theory Comput.
10
,
1599
(
2014
).
42.
V.
Babin
,
C.
Leforestier
, and
F.
Paesani
, “
Development of a ‘first principles’ water potential with flexible monomers: Dimer potential energy surface, VRT spectrum, and second virial coefficient
,”
J. Chem. Theory Comput.
9
,
5395
(
2013
).
43.
R.
Bukowski
,
K.
Szalewicz
,
G. C.
Groenenboom
, and
A.
Van der Avoird
, “
Predictions of the properties of water from first principles
,”
Science
315
,
1249
(
2007
).
44.
X.
Huang
,
B. J.
Braams
,
J. M.
Bowman
,
R. E. A.
Kelly
,
J.
Tennyson
,
G. C.
Groenenboom
, and
A.
van der Avoird
, “
New ab initio potential energy surface and the vibration-rotation-tunneling levels of (H2O)2 and (D2O)2
,”
J. Chem. Phys.
128
,
034312
(
2008
).
45.
A.
Shank
,
Y.
Wang
,
A.
Kaledin
,
B. J.
Braams
, and
J. M.
Bowman
, “
Accurate ab initio and ‘hybrid’ potential energy surfaces, intramolecular vibrational energies, and classical IR spectrum of the water dimer
,”
J. Chem. Phys.
130
,
144314
(
2009
).
46.
L.-P.
Wang
,
T. J.
Martinez
, and
V. S.
Pande
, “
Building force fields: An automatic, systematic, and reproducible approach
,”
J. Phys. Chem. Lett.
5
,
1885
(
2014
).
47.
X.-G.
Wang
and
T.
Carrington
, Jr.
, “
Using monomer vibrational wavefunctions to compute numerically exact (12D) rovibrational levels of water dimer
,”
J. Chem. Phys.
148
,
074108
(
2018
).
48.
J. K.
Gregory
and
D. C.
Clary
, “
Calculations of the tunneling splittings in water dimer and trimer using diffusion Monte Carlo
,”
J. Chem. Phys.
102
,
7817
(
1995
).
49.
J. K.
Gregory
and
D. C.
Clary
, “
Tunneling dynamics in water tetramer and pentamer
,”
J. Chem. Phys.
105
,
6626
(
1996
).
50.
W. H.
Miller
, “
Periodic orbit description of tunneling in symmetric and asymmetric double-well potentials
,”
J. Phys. Chem.
83
,
960
(
1979
).
51.
D. J.
Wales
, “
Theoretical study of water trimer
,”
J. Am. Chem. Soc.
115
,
11180
(
1993
).
52.
Y.
Watanabe
,
T.
Taketsugu
, and
D. J.
Wales
, “
An ab initio study of tunneling splittings in the water dimer
,”
J. Chem. Phys.
120
,
5993
(
2004
).
53.
M.
Takahashi
,
Y.
Watanabe
,
T.
Taketsugu
, and
D. J.
Wales
, “
An ab initio study of tunneling splittings in the water trimer
,”
J. Chem. Phys.
123
,
044302
(
2005
).
54.
G. V.
Mil’nikov
and
H.
Nakamura
, “
Practical implementation of the instanton theory for the ground-state tunneling splitting
,”
J. Chem. Phys.
115
,
6881
(
2001
).
55.
J. O.
Richardson
and
S. C.
Althorpe
, “
Ring-polymer instanton method for calculating tunneling splittings
,”
J. Chem. Phys.
134
,
054109
(
2011
).
56.
J. O.
Richardson
,
S. C.
Althorpe
, and
D. J.
Wales
, “
Instanton calculations of tunneling splittings for water dimer and trimer
,”
J. Chem. Phys.
135
,
124109
(
2011
).
57.
J. O.
Richardson
, “
Perspective: Ring-polymer instanton theory
,”
J. Chem. Phys.
148
,
200901
(
2018
).
58.
C. L.
Vaillant
and
M. T.
Cvitaš
, “
Rotation-tunneling spectrum of the water dimer from instanton theory
,”
Phys. Chem. Chem. Phys.
20
,
26809
(
2018
).
59.
M.
Eraković
and
M. T.
Cvitaš
, “
Tunneling splittings of vibrationally excited states using general instanton paths
,”
J. Chem. Phys.
153
,
134106
(
2020
).
60.
E.
Mátyus
,
D. J.
Wales
, and
S. C.
Althorpe
, “
Quantum tunneling splittings from path-integral molecular dynamics
,”
J. Chem. Phys.
144
,
114108
(
2016
).
61.
E.
Mátyus
and
S. C.
Althorpe
, “
Calculating splittings between energy levels of different symmetry using path-integral methods
,”
J. Chem. Phys.
144
,
114109
(
2016
).
62.
C. L.
Vaillant
,
D. J.
Wales
, and
S. C.
Althorpe
, “
Tunneling splittings from path-integral molecular dynamics using a Langevin thermostat
,”
J. Chem. Phys.
148
,
234102
(
2018
).
63.
C. L.
Vaillant
,
D. J.
Wales
, and
S. C.
Althorpe
, “
Tunneling splittings in water clusters from path integral molecular dynamics
,”
J. Phys. Chem. Lett.
10
,
7300
(
2019
).
64.
Y.-C.
Zhu
,
S.
Yang
,
J.-X.
Zeng
,
W.
Fang
,
L.
Jiang
,
D. H.
Zhang
, and
X.-Z.
Li
, “
Torsional tunneling splitting in a water trimer
,”
J. Am. Chem. Soc.
144
,
21356
(
2022
).
65.
S.
Coleman
, “
The uses of instantons
,” in
The Whys of Subnuclear Physics
edited by
A.
Zichichi
(
Springer US
,
Boston, MA
,
1979
), p.
805
.
66.
H.
Kleinert
,
Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets
,
5th ed.
(
World Scientific
,
2009
).
67.
J. O.
Richardson
, “
Ring-polymer instanton theory
,”
Int. Rev. Phys. Chem.
37
,
171
(
2018
).
68.
D. M.
Ceperley
and
G.
Jacucci
, “
Calculation of exchange frequencies in bcc 3He with the path-integral Monte Carlo method
,”
Phys. Rev. Lett.
58
,
1648
(
1987
).
69.
C.
Alexandrou
and
J. W.
Negele
, “
Stochastic calculation of tunneling in systems with many degrees of freedom
,”
Phys. Rev. C
37
,
1513
(
1988
).
70.
J. K.
Gregory
and
D. C.
Clary
, “
Structure of water clusters. The contribution of many-body forces, monomer relaxation, and vibrational zero-point energy
,”
J. Phys. Chem.
100
,
18014
(
1996
).
71.
C.
Pérez
,
M. T.
Muckle
,
D. P.
Zaleski
,
N. A.
Seifert
,
B.
Temelso
,
G. C.
Shields
,
Z.
Kisiel
, and
B. H.
Pate
, “
Structures of cage, prism, and book isomers of water hexamer from broadband rotational spectroscopy
,”
Science
336
,
897
(
2012
).
72.
C.
Pérez
,
D. P.
Zaleski
,
N. A.
Seifert
,
B.
Temelso
,
G. C.
Shields
,
Z.
Kisiel
, and
B. H.
Pate
, “
Hydrogen bond cooperativity and the three-dimensional structures of water nonamers and decamers
,”
Angew. Chem., Int. Ed.
53
,
14368
(
2014
).
73.
T. R.
Dyke
, “
Group theoretical classification of the tunneling–rotational energy levels of water dimer
,”
J. Chem. Phys.
66
,
492
(
1977
).
74.
M. T.
Cvitaš
and
J. O.
Richardson
, “
Quantum dynamics in water clusters
,” in
Molecular Spectroscopy and Quantum Dynamics
edited by
R.
Marquardt
and
M.
Quack
(
Elsevier
,
2021
), Chap. 9, p.
301
.
75.
M.
Schütz
,
T.
Bürgi
,
S.
Leutwyler
, and
H. B.
Bürgi
, “
Fluxionality and low-lying transition structures of the water trimer
,”
J. Chem. Phys.
99
,
5228
(
1993
).
76.
A.
van der Avoird
,
E. H. T.
Olthof
, and
P. E. S.
Wormer
, “
Tunneling dynamics, symmetry, and far-infrared spectrum of the rotating water trimer. I. Hamiltonian and qualitative model
,”
J. Chem. Phys.
105
,
8034
(
1996
).
77.
L. H.
Coudert
and
J. T.
Hougen
, “
Tunneling splittings in the water dimer: Further development of the theory
,”
J. Mol. Spectrosc.
130
,
86
(
1988
).
78.
T.
Taketsugu
and
D. J.
Wales
, “
Theoretical study of rearrangements in water dimer and trimer
,”
Mol. Phys.
100
,
2793
(
2002
).
79.
G. C.
Groenenboom
,
E. M.
Mas
,
R.
Bukowski
,
K.
Szalewicz
,
P. E. S.
Wormer
, and
A.
van der Avoird
, “
Water pair and three-body potential of spectroscopic quality from ab initio calculations
,”
Phys. Rev. Lett.
84
,
4072
(
2000
).
80.
A.
van der Avoird
and
K.
Szalewicz
, “
Water trimer torsional spectrum from accurate ab initio and semiempirical potentials
,”
J. Chem. Phys.
128
,
014302
(
2008
).
81.
R.
Schwan
,
C.
Qu
,
D.
Mani
,
N.
Pal
,
G.
Schwaab
,
J. M.
Bowman
,
G. S.
Tschumper
, and
M.
Havenith
, “
Observation of the low-frequency spectrum of the water trimer as a sensitive test of the water-trimer potential and the dipole-moment surface
,”
Angew. Chem., Int. Ed.
59
,
11399
(
2020
).
82.
F. N.
Keutsch
,
J. D.
Cruzan
, and
R. J.
Saykally
, “
The water trimer
,”
Chem. Rev.
103
,
2533
(
2003
).
83.
K.
Liu
,
J. G.
Loeser
,
M. J.
Elrod
,
B. C.
Host
,
J. A.
Rzepiela
,
N.
Pugliano
, and
R. J.
Saykally
, “
Dynamics of structural rearrangements in the water trimer
,”
J. Am. Chem. Soc.
116
,
3507
(
1994
).
84.
M. R.
Viant
,
J. D.
Cruzan
,
D. D.
Lucas
,
M. G.
Brown
,
K.
Liu
, and
R. J.
Saykally
, “
Pseudorotation in water trimer isotopomers using terahertz laser spectroscopy
,”
J. Phys. Chem. A
101
,
9032
(
1997
).
85.
T.
Carrington
, “
Perspective: Computing (ro-)vibrational spectra of molecules with more than four atoms
,”
J. Chem. Phys.
146
,
120902
(
2017
).
86.
W.
Lin
,
J.-X.
Han
,
L. K.
Takahashi
,
H. A.
Harker
,
F. N.
Keutsch
, and
R. J.
Saykally
, “
Terahertz vibration-rotation-tunneling spectroscopy of the water tetramer-d8: Combined analysis of vibrational bands at 4.1 and 2.0 THz
,”
J. Chem. Phys.
128
,
094302
(
2008
).
87.
G. V.
Mil’nikov
and
H.
Nakamura
, “
Instanton theory for the tunneling splitting of low vibrationally excited states
,”
J. Chem. Phys.
122
,
124311
(
2005
).
88.
D. J.
Wales
and
T. R.
Walsh
, “
Theoretical study of the water pentamer
,”
J. Chem. Phys.
105
,
6957
(
1996
).
89.
K.
Liu
,
M. G.
Brown
,
J. D.
Cruzan
, and
R. J.
Saykally
, “
Vibration-rotation tunneling spectra of the water pentamer: Structure and dynamics
,”
Science
271
,
62
(
1996
).
90.
H. A.
Harker
,
M. R.
Viant
,
F. N.
Keutsch
,
E. A.
Michael
,
R. P.
McLaughlin
, and
R. J.
Saykally
, “
Water pentamer: Characterization of the torsional-puckering manifold by terahertz vrt spectroscopy
,”
J. Phys. Chem. A
109
,
6483
(
2005
).
91.
M. T.
Cvitaš
and
J. O.
Richardson
, “
Quantum tunnelling pathways of the water pentamer
,”
Phys. Chem. Chem. Phys.
22
,
1035
(
2020
).
92.
V.
Babin
,
G. R.
Medders
, and
F.
Paesani
, “
Toward a universal water model: First principles simulations from the dimer to the liquid phase
,”
J. Phys. Chem. Lett.
3
,
3765
(
2012
).