## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Unraveling quantum mechanical effects in water using isotopic fractionation

Contributed by B. J. Berne, February 28, 2012 (sent for review February 2, 2012)

## Abstract

When two phases of water are at equilibrium, the ratio of hydrogen isotopes in each is slightly altered because of their different phase affinities. This isotopic fractionation process can be utilized to analyze water’s movement in the world’s climate. Here we show that equilibrium fractionation ratios, an entirely quantum mechanical property, also provide a sensitive probe to assess the magnitude of nuclear quantum fluctuations in water. By comparing the predictions of a series of water models, we show that those describing the OH chemical bond as rigid or harmonic greatly overpredict the magnitude of isotope fractionation. Models that account for anharmonicity in this coordinate are shown to provide much more accurate results because of their ability to give partial cancellation between inter- and intramolecular quantum effects. These results give evidence of the existence of competing quantum effects in water and allow us to identify how this cancellation varies across a wide-range of temperatures. In addition, this work demonstrates that simulation can provide accurate predictions and insights into hydrogen fractionation.

Water within Earth’s atmosphere is naturally composed of the stable hydrogen isotopes hydrogen (H) and deuterium (D). During cycles of evaporation, condensation, and precipitation, these isotopes naturally undergo partial separation due to their differing masses, thereby leading to different H/D ratios in the two phases. This process of fractionation has a number of fortuitous consequences, which are utilized in hydrology and geology. For instance, by comparing the ratio of H to D, one can estimate the origins of a water sample, the temperature at which it was formed, and the altitude at which precipitation occurred (1, 2). Equilibrium fractionation, where the two phases are allowed to equilibrate their H/D ratio, is entirely a consequence of the effects of quantum mechanical fluctuations on water’s hydrogen bond network. Quantum mechanical effects such as zero-point energy and tunneling are larger for H due to its lower mass.

Despite numerous studies, the extent to which quantum fluctuations affect water’s structure and dynamics remains a subject of considerable debate. It has long been appreciated that one effect of quantum fluctuations in water is the disruption of hydrogen bonding, leading to destructuring of the liquid and faster dynamics (3⇓⇓–6). However, more recent work has suggested that a competing quantum effect may exist in water (7, 8), namely that the quantum kinetic energy in the OH covalent bond allows it to stretch and form shorter and stronger hydrogen bonds, which partially cancels the disruptive effect. This hydrogen bond strengthening has only been recently appreciated, as many original studies drew their conclusions based on models with rigid or harmonic bonds, which are unable to describe this behavior. The degree of quantum effect cancellation depends sensitively on the anharmonicity of the OH stretch and the temperature. These parameters tune the balance between the lower frequency hydrogen bonding disruption, which will dominate at lower temperatures, and the higher frequency hydrogen bond strengthening effect, which will dominate at higher temperatures when rotations become essentially classical.

If such a large degree of cancellation existed at ambient temperature, it would be highly fortuitous both in terms of the biological effects of heavy water, which is only mildy toxic to humans (9), as well as the ability to use heavy solvents in two-dimensional-IR and NMR spectroscopies, where deuteration is assumed not to dramatically alter the structure or dynamics observed. However, the size of this cancellation remains elusive because empirical quantum models of water are typically fit to reproduce its properties when used in path integral simulations and the two ab initio path integral studies performed have not produced a consistent picture (7, 10). In addition, many of these simulation studies compare the properties of water to those of its classical counterpart, but classical water is physically unrealizable even at relatively high temperatures, because water still has significant quantum effects present in its vibrations.

In this paper, we use equilibrium fractionation ratios as a sensitive probe to assess the magnitude of quantum mechanical effects in water. Fractionation ratios can be directly related to quantum kinetic energy differences between H and D in liquid water and its vapor and can be calculated exactly for a given water potential energy model using path integral simulations. The large number of accurate experimental measurements of these ratios allows for sensitive comparisons of theory and experiment over a wide-range of temperatures (11). In the present work, we show what features are needed in a water model to accurately predict these ratios by decomposing the contributions to the free energy difference leading to fractionation. This analysis in turn leads to a simple explanation of the inversion of the fractionation ratios seen experimentally at high temperatures, where D is favored over H in the vapor phase (11).

## Calculating Fractionation Ratios

The liquid-vapor fractionation ratio, *α*_{l-v}, is defined as [1]where *x*_{Z} is the mole fraction of isotope *Z*, *l* denotes the liquid phase, and *v* denotes the vapor phase. In the second equality, Δ*A* is the Helmholtz free energy corresponding to the process [2]In this work we consider the dilute D limit which reflects the situation found in the Earth’s atmosphere where it is 6,000 times less common than H. In this limit, we consider the free energy of exchanging a single D atom in a vapor water molecule with an H atom in a liquid water molecule, with all other molecules being H_{2}O. The free energy difference can be calculated from the thermodynamic integration expression (12) [3]where and are the kinetic energy expectation values for a hydrogen isotope of mass *m*_{Z} in a water molecule HO*Z* in the liquid and vapor phases, respectively. The kinetic energy can be calculated exactly for a given potential energy model of water using a path integral molecular dynamics (PIMD) simulation. These simulations exploit the exact isomorphism between a system of quantum mechanical particles and that of a set of classical ring polymers in which the spread of a polymer is directly related to that quantum particle’s position uncertainty (13). The kinetic energy for the particle *Z* in the molecule HO*Z* can be calculated from these simulations using the centroid virial estimator (14, 15) [4]where *T* is the temperature and *k*_{B} is the Boltzmann constant. Here, is vector from the *k*th bead representing particle *Z* to the center of the ring polymer and is the force on that bead as shown in Fig. 1. The first term in this expression is the classical kinetic energy and is independent of the surrounding environment, thus it is identical in the vapor and liquid phase. The second term is the kinetic energy associated with confinement of a quantum particle. This confinement depends on the forces exerted on the particle by the surrounding molecules. Equilibrium fractionation is thus an entirely quantum mechanical phenomenon.

PIMD simulations were performed using 1,000 water molecules with *n* = 32 ring polymer beads used for the imaginary time discretization.* Previously described evolution and thermostating procedures were used (16). The computational cost of these calculations was reduced by using the ring polymer contraction technique with a cutoff of 5 Å, which for this system size leads to more than an order of magnitude speedup compared to a standard PIMD implementation (17, 18). The integral in Eq. **3** was performed using the midpoint rule with 11 masses *m*_{Z} evenly spaced between *m*_{H} and *m*_{D}. Calculations were performed at the experimental coexistence densities (19).

To model the interactions within and between water molecules, we used the q-SPCFw (5) and q-TIP4P/F (8) models. Both models are flexible, use point charges, and have a harmonic description of the bending mode. However, whereas q-SPCFw uses a purely harmonic description of the OH stretch [5]q-TIP4P/F contains anharmonicity by modeling the stretch as a Morse expansion truncated at fourth order [6]Here *r* is the distance between the oxygen and hydrogen atom and the parameters are given in refs. 5 and 8. The anharmonicity in the q-TIP4P/F model makes the observed quantum mechanical effects much smaller than previously predicted from harmonic or rigid models and gave rise to the idea of competing quantum effects in water (8).

Both models have previously been shown to accurately reproduce many of water’s properties in PIMD simulations of liquid water. Because of their simple potential form, such models are generally less transferable to other phases than more sophisticated polarizable or ab initio descriptions. However, partially adiabatic centroid molecular dynamics simulations (20) have shown that the anharmonic stretch, Eq. **6**, allows reasonable agreement to be obtained in the observed frequency shifts in the infrared spectrum in going from liquid to gaseous water as well as from pure light to pure heavy water (8, 21, 22). These models were chosen to assess the importance of anharmonicity in the OH stretch using a zeroth order description of liquid water. Hence, as we show below, they offer a straightforward way to assess competing quantum effects in water.

## Results and Discussion

Fig. 2 shows the fractionation factors calculated from our PIMD simulations compared to the experimental data of ref. 11. For consistency with the experimental data, we plot 10^{3} ln *α*_{l-v}, which is simply -10^{3}Δ*A*/*k*_{B}*T*. Because Δ*A*/*k*_{B}*T* is generally small, *e*^{-ΔA/kBT} ≃ 1 - Δ*A*/*k*_{B}*T*. Thus an experimental value at 280 K of 100 corresponds to 10% more D residing in the liquid than the vapor. Above 500 K, the experimental data shows a well-characterized region of inverse fractionation where D becomes more favored in the vapor than in the liquid, as shown in the inset of Fig. 2.

Turning to the simulated data, we observe that the harmonic q-SPC/Fw model overpredicts the magnitude of fractionation at 300 K by a factor of 3, and does not fall to that experimental value until the temperature is raised to 450 K. In contrast, the q-TIP4P/F model is in error by only 25% at the lowest temperature and approaches the experimental values more closely at higher temperatures. It also correctly shows inverse fractionation above 540 K. A previous study using the rigid SPC/E model and *ℏ*^{2} perturbation theory, found an H/D fractionation of 450 at 300 K, which is about five times higher than seen experimentally (23). However, it is not clear whether this difference was purely due to the use of a rigid model or whether the approximate *ℏ*^{2} perturbation technique used to obtain the fractionation ratios was also at fault.

The data represented in Fig. 2 shows that the q-TIP4P/F model is in much better agreement with the experiment than the q-SPCFw model, however it is not yet clear what aspect of the parameterization causes this improvement. To better understand what causes q-TIP4P/F to be better, we constructed two models, which we denote Aq-SPC/Fw and Hq-TIP4P/F. In the former, the q-SPC/Fw water model has its harmonic OH stretch replaced by a fourth order Morse expansion using the parameters of q-TIP4P/F; in the latter, the Morse potential of q-TIP4P/F is truncated at the harmonic term (see Eqs. **5** and **6**). The anharmonic variant of q-SPC/Fw gives results as good as q-TIP4P/F, whereas the harmonic version of q-TIP4P/F fares as poorly as q-SPCFw. In other words, the accurate prediction of the fractionation ratios in liquid water is tied to the anharmonicity in the OH direction and is rather insensitive to the other parameters. This finding sheds light on a previous study where a sophisticated rigid polarizable model gave identical predictions for H/D fractionation to the simple fixed-charge rigid SPC/E model, i.e., varying the intermolecular potential alone does not give the flexibility required to accurately reproduce the experimental fractionation ratios (23).

To determine the reason for the inversion observed in fractionation above 500 K in both experiment and the q-TIP4P/F model, we decompose the contributions to the fractionation ratio by noting that the quantum contribution to the kinetic energy in Eq. **4** is the dot product of two vectors. The overall kinetic energy is invariant to the coordinate system used to evaluate it, thus when the kinetic energy is calculated in the standard Cartesian basis all three components will average to the same number due to the isotropy of the liquid. To gain further insight, we instead use the internal coordinates of the water molecule and determine the contribution to 10^{3} ln *α*_{l-v} arising from the OH bond vector, a vector in the plane of the molecule, and the vector perpendicular to the molecular plane as shown in Fig. 1. This approach is similar to the one taken by Lin et al. in the different context of investigating the proton momentum distribution in ice (24). The results of this decomposition are shown in Table 1 for 300 K, where D is experimentally seen to favor the liquid; and Table 2 for 620 K, where D is experimentally seen to favor the vapor.

From Table 1, we see that all of the models have largely similar contributions from the two directions orthogonal to the OH bond and that both are positive, and therefore, favor the D excess in the liquid. The contribution perpendicular to the plane is noticeably larger than the contribution in the plane. As shown in Eq. **3**, the values depend on the change in the quantum kinetic energy between H and D in the liquid and the vapor, which in turn is determined by how much the H or D atom’s position uncertainty is restricted by interacting with the other water molecules in the liquid or vapor. Because, in the vapor, there is little confinement in the plane orthogonal to the water, a larger contribution from that direction is expected; in the liquid, other molecules are present which restrict expansion in that coordinate.

In all cases, the OH contribution is negative, indicating that there is less confinement in the position of the H atom in that direction in the liquid than in the vapor. This reduction in confinement occurs because the hydrogen atom participates in hydrogen bonding allowing the OH chemical bond to stretch more easily. However, comparing the anharmonic models (Aq-SPCFw and q-TIP4P/F) with their harmonic counterparts (q-SPCFw and Hq-TIP4P/F), we observe a 10-fold increase in the values arising from the OH contribution. This increase gives rise to a larger cancellation of the positive contributions from the two orthogonal vectors. It is this cancellation that leads to the much better agreement with the experimental data at 300 K.

We now turn to Table 2, which shows the contributions to the fractionation ratio at 620 K, a regime where experimentally D is preferred in the lighter phase. This contribution is closer to the classical limit where the fractionation would be zero, thus, as expected, each component is reduced in magnitude compared to the lower temperature data in Table 1. However, the relative decrease in each component varies. The contributions arising from the in-plane and out-of-plane contributions orthogonal to OH decrease by a factor of 7–8, whereas the OH contribution falls by a factor of only about 4. For the anharmonic models, the negative OH contribution outweighs the positive components in the other two directions leading to an inversion of the fractionation compared to that seen at 300 K in agreement with the experimental observation of -2. The reason the OH component falls off more slowly is that this direction is dominated by stretching of the OH chemical bond, which is a high frequency coordinate, so even at high temperatures, quantum mechanics plays a noticeable role. In contrast, the two directions orthogonal to the OH direction are lower frequency, and so approach the classical limit more rapidly as the temperature is increased. Thus, the reason the H/D fractionation is low around 600 K is not due to the fact that all contributions are individually low, but rather that they nearly exactly cancel at this point due to the different rates at which the components approach the classical limit.

Finally, we computed the fractionation ratio for the TTM3-F water model (25), which is known to have a very large cancellation of its quantum mechanical effects at 300 K (8, 26). This model was fit to ab initio calculations using a potential form incorporating anharmonic flexibility, geometry dependent charges, and polarizability on the oxygen site. As such, it represents the current gold standard of parameterized water models and has been used extensively in recent studies probing the effects of quantum mechanical fluctuations on water (21, 27, 28). To see if the large cancellation of quantum effects predecided by this model is consistent with experimental fractionation ratios, we calculated the value at 300 K, which yielded a value of 10^{3} ln *α*_{l-v} = -57 and that D is favored in the vapor at all temperatures, a result that is in qualitative and quantitative agreement with the experimental results. This model therefore overpredicts competing quantum effects in water and hence care should be taken concerning its predictions on the effects of quantum fluctuations on water’s structure and dynamics. However, the overly anharmonic OH coordinate is likely much less of an issue when this model is used in classical simulations because the zero point energy in the OH bond equates to an effective 2,000 K temperature rise in that coordinate. As such, the large extensions of the OH bond into anharmonic regions experienced in quantum simulations would be very rare events in a classical simulation. In addition, based on our discussion above, it is likely that a small reparamaterization of the OH bond anharmonicity or monomer dipole moment surface could correct this discrepancy.

## Conclusion

In conclusion, we have shown that including anharmonicity in the OH bond when modeling water is essential to obtain agreement with the experimentally observed H/D fractionation ratios and that these ratios provide an excellent method to assess the accuracy of the quantum effects predicted by models of water. It will therefore be interesting to see if fractionation ratios can be used to assess the accuracy of ab inito predictions of the effect of quantum fluctuations on water’s structure and dynamics (7, 10). Since it has recently been shown that the competition between quantum mechanical effects applies to other hydrogen bonded systems (29), it is likely that many of our conclusions will be relevant to understanding isotopic fractionation in these systems. Additionally, whereas we only considered equilibrium fractionation in this work, which can be calculated exactly for a given potential energy model using PIMD simulations, many water processes occurring in the world’s atmosphere are nonequilibrium ones. Although including the effects of quantum fluctuations on dynamics is a much more challenging feat, the recent development of efficient condensed phase quantum dynamics approaches (30⇓–32) should allow insights to be gained into kinetic fractionation processes. These directions will form the basis of future work.

## Acknowledgments

The authors gratefully thank Joseph Morrone and David Selassie for helpful comments and a critical reading of this manuscript. This research was supported by National Science Foundation Grant NSF-CHE-0910943 (to B.J.B.).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: bb8{at}columbia.edu.

Author contributions: T.E.M. and B.J.B. designed research; T.E.M. performed research; T.E.M. analyzed data; and T.E.M. and B.J.B. wrote the paper.

The authors declare no conflict of interest.

↵

^{*}The simulations were repeated for 64 beads per particle and the results were found to be essentially the same as those reported in this paper for 32 beads per particle. Except for a small shift at*T*= 280 and 300 K, the conclusions are the same.

## References

- ↵
- Hoefs J

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Fanourgakis GS,
- Xantheas SS

- ↵
- ↵
- ↵
- ↵
- Li XZ,
- Walker B,
- Michaelides A

- ↵
- ↵
- ↵

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Physical Sciences
- Chemistry