The hydrogen bond (HB) is a generally weak but ubiquitous intermolecular interaction. It is responsible, for example, for holding together the two strands of DNA molecules and many phases of water. As such, understanding the behaviors of HBs is lying at the core of many studies concerning the behaviors of water and biomolecules. The small mass of hydrogen (H) means that many of these problems are quantum mechanical in nature, ^{[ 1 – 7 ]} although all too often these nuclear quantum effects (NQEs) are neglected. These NQEs are essential for rationalizing many phenomena of te H-bonded systems, especially where quantum tunneling and zero-point motion are concerned. For example, starting from the 1950s, it has been known that replacing H by deuterium (D) in H-bonded molecular crystals with electronic structures unchanged can result in obvious variations of the lattice constants. ^{[ 8 ]} This is currently known as the Ubbelohde effect. ^{[ 9 ]} In enzyme reactions, the proton tunneling along HBs is referred to as “nature’s subway”. By using a neutron Compton scattering experiment, a picture concerning the impact of NQEs on the real space delocalization and vibrational properties of the proton in liquid water was also revealed. ^{[ 10 , 11 ]} Most recently, a cryogenic scanning tunneling microscope (STM) experiment even confirmed the isotope-dependent switching rate of water tetramer on salt between its two chirality states, which could only be rationalized by considering the NQEs. ^{[ 12 ]}

From the theoretical perspective, accurate computer simulations of these phenomena require an accurate quantum mechanical description of the finite-temperature nuclear fluctuations in real polyatomic systems. A natural choice is to construct a high-dimensional Schrödinger equation for the many-body entity of the nuclei. By solving the eigenstate wave-functions of this Schrödinger equation, one can observe not only the statistical but also the dynamical properties of the system under investigation related to the quantum nature of the nuclei. Over the last ∼ 20 years, this method has been extremely successful in describing gas-phase reactions. ^{[ 13 – 17 ]} However, because of the notorious scaling problem in mapping the high-dimensional potential energy surface (PES) and solving the Schrödinger equation, the application of this method has been seriously limited to systems with less than ∼ 6 atoms. For more realistic polyatomic systems and condensed matters, a practical method is highly desired. Thanks to the development of the path-integral representation of quantum mechanics starting from the late 1940s, ^{[ 18 – 21 ]} a scheme practical for large polyatomic systems and condensed matters in which not only the thermal effects but also the NQEs can be accounted for statistically has been systematically presented by Feynman and Hibbs. ^{[ 22 ]} Based on such a foundation, the molecular dynamics (MD) simulation technique was combined with this path-integral representation of quantum mechanics, and a series of path-integral molecular dynamics (PIMD) simulations were carried out in the 1980s, as reported by Chandler, Parrinello, and Berne et al. in Refs. [ 23 ]–[ 25 ], respectively. Analogous to these MD based simulations, the Monte–Carlo (MC) sampling technique can also be used. At roughly the same time, the properties of liquid helium, including its superfluidity, were systematically studied by Ceperley and Pollock using path-integral Monte–Carlo (PIMC) simulations. ^{[ 26 – 29 ]} In these PIMD/PIMC simulations, the NQEs are accounted for on the same footing as the thermal effects when the statistical properties are evaluated, as will be explained later in Section 2. Therefore, when comparing their results with the ones obtained in the MD simulations, in which the nuclei are treated classically, the NQEs can be addressed in a very clean manner.

These comparisons set up a rigorous framework for the NQEs to be investigated statistically, which is still used nowadays. However, in these early PIMD and PIMC simulations, researches often resorted to empirical potentials for descriptions of the interatomic interactions. When chemical reactions occur, e.g., a proton transfers along a HB, such a simple parametrization of the interatomic interactions may easily fail. To cope with these problems, in which the impact of the NQEs is usually more interesting, an effort toward combining the ab initio electronic structure calculations with the PIMD/PIMC sampling techniques has been made since the mid 1990s, first within the framework of Car–Parrinello (CP) MD (see e.g., Tuckerman, Marx, and Parrinello et al. ’s work in Refs. [ 1 ], [ 30 ], and [ 31 ]) and later by directly using Born–Oppenheimer (BO) MD/MC. ^{[ 6 , 7 , 32 – 34 ]} These methods allow for bond making and breaking events, as well as the thermal and quantum nuclear fluctuations to be addressed in a seamless manner based on the forces computed “on-the-fly” as the system evolves, therefore making PIMD/PIMC simulations for problems with complicated chemical environments possible. Currently, not only different functionals within the density-functional theory (DFT) ^{[ 35 – 37 ]} but also many traditional quantum chemistry methods such as MP2 ^{[ 38 , 39 ]} can be used in descriptions of the electronic structures and consequently the interatomic interactions. With such a wide range of choices for the electronic structures, one can safely rely on the results obtained from such simulations, if the accuracy of the interatomic interaction description and the sampling sufficiency are ensured. We note, however, that such an ensurance is far from being trivial. ^{[ 30 , 35 , 36 ]}

The above-mentioned recipe allows the NQEs to be described in a clean manner, if the molecular simulation time is long enough to address the issue of interest. A common problem for most molecular simulation techniques, however, is that the timescale one can afford is not long enough to allow the system to travel between different stable states, thereby hindering a direct evaluation of the free-energy differences between different meta-stable states. In such cases, to determine the relative stability of these different meta-stable states without resorting to the coarse-grained models, ^{[ 40 ]} the method one often chooses is to enhance the sampling efficiency either by imposing constraints or modifying the effective PES. ^{[ 41 – 44 ]} However, if it is the case that within the molecular simulation time, the system can stay within the local minima of the PES with all spatial configurations close to it explored, a comparison of the free energies associated with each local minimum directly without requiring the system to travel between different meta-stable states provides an alternative. In such cases, a method with which the free energy of each meta-stable state can be calculated is highly desired. When the influence of NQEs on these free energies needs to be quantified, a combination between the thermodynamic integration (TI) method and PIMD is a natural choice. ^{[ 45 , 46 ]} In the following, in addition to a general introduction to ab initio PIMD and its application to the studies of some H-bonded systems, a combination of the TI method with ab initio PIMD along with applications of this method to the free-energy calculations of some H-bonded systems will also be discussed.

The paper is organized as follows. An introduction to the ab initio PIMD method and its combination with TI, as well as the computational setups, are given in Section 2. In Section 3, we take a basic question in physics and chemistry, i.e., what is the impact of NQEs on the strength of HB, as an example to show how the ab initio PIMD method can be used to study such problems in H-bonded systems. As the second example, we quantitatively analyze the influence of NQEs on the free energy of some gas phase water clusters using a combination of ab initio PIMD and TI methods in Section 4. Conclusions and a short perspective are given in Section 5.

The first example of the ab initio PIMD simulations, in which the quantum nature of the nuclei plays an important role, concerns a fundamental question in physics and chemistry, i.e., what is the influence of NQEs on the strength of HBs? We start by using the structural properties, namely, the heavy-atom distances, as an indicator. Further investigations on the energetics will be discussed in Section 4. As mentioned earlier in this manuscript, it has been known for more than half a century that, in H-bonded molecular crystals, by replacing H with D, the strength of HBs and consequently the crystal structural parameters change. This phenomenon is known as the Ubbelohde effect. ^{[ 8 , 9 ]} The conventional Ubbelohde effect reflects an increase of the heavy-atom distance associated with the HB upon replacing H by D, which indicates that NQEs strengthen this intermolecular interaction. We note, however, that a negative Ubbelohde effect in which a decrease of the heavy-atom distance happens upon replacing H by D has also been observed in several materials. ^{[ 8 , 9 ]} From the perspective of molecular simulation, since the 1980s, when PIMD/PIMC becomes a conventional routine to investigate the NQEs in real polyatomic systems, a bunch of computer simulations have been carried out in a wide range of systems. A general conclusion is that the final effect is system-dependent. For example, in liquid HF, ab initio MD and PIMD simulations have shown that when NQEs are accounted for, the first peak of the F–F radial distribution function (RDF) sharpens and shifts to a shorter F–F distance, ^{[ 59 ]} indicating that NQEs strengthen the HB in this system. In contrast, a similar simulation for liquid water shows that the O–O RDF gets less structured when NQEs are included, indicating that NQEs weaken the HB in liquid water. ^{[ 2 ]} We note, however, that although this conclusion is probably right, it is the opposite of what has been reported in Ref. [ 5 ].

Besides this discussion concerning HBs in condensed phases, the impact of NQEs on the strength of HBs has also been widely discussed in gas-phase clusters. ^{[ 59 – 61 ]} For water clusters up to the hexamer, it was predicted that the NQEs weaken the HBs, as reflected by the change of the heavy-atom distances in Refs. [ 60 ] and [ 61 ]. For HF clusters, however, it was predicted that NQEs can strengthen or weaken the HBs depending on the size of the clusters. ^{[ 59 ]} For clusters smaller than tetramer, the NQEs weaken this intermolecular interaction, whereas for clusters larger than it, the strengthening of the HB is observed. In tetramer, the NQEs are negligible. With these system-dependent conclusions in mind, it is clear that a simple picture in which the influence of the NQEs on the strength of HBs and consequently the structural or energetic properties can be rationalized is highly desired. In Ref. [ 7 ], we have presented a simple picture which can be used to rationalize these different results. Ab initio PIMD simulations were used as the basic tool upon which the analysis was carried out. Here, using these simulations as an example, we will review this work in a way that how the ab initio MD and PIMD simulations are used to answer this question can be understood in a easy manner.

To be as unbiased as possible, a wide range of H-bonded systems have been chosen, including HF clusters (dimer to hexamer), H ₂ O clusters (dimer, pentamer, and octamer), charged, protonated, and hydroxylated water and ammonia clusters ( , , , and ), organic dimers (formic acid and formamide), and solids (HF, HCl, and squaric acid C ₄ H ₂ O ₄ ). For each system, both ab initio MD simulations, in which the nuclei are classical particles, and ab initio PIMD simulations, in which the nuclei are quantum particles, are carried out. By comparing the difference of the statistical structures obtained from these two simulations, one can analyze the influence of NQEs in a very clean manner. The quantities we focus on in characterizing the HBs include: (i) the heavy-atom distances (denoted by X – X , where X is either O, Cl, C, N, or F), which describe the intermolecular separations; (ii) the H-bond angles ( X –H··· X ), which characterize the HB bending (libration); and (iii) the X –H covalent bond lengths, which describe the stretching in the H-bond donor molecules. These quantities provide an indication of H-bond strength, as will be explained later in the discussion. As the main measure of H-bond strength, we use a standard estimator based on the computed red shift (softening) in the X –H stretching frequency of the HB donor molecule. We note that there is no perfect measure for the HB strength; ^{[ 62 ]} however, the red shift of the stretching frequency is a widely used measure [see, e.g. Refs. [ 63 ] and [ 64 ]). In Fig.  2(b) , it is shown that this estimator correlates well with the computed binding energy per HB in the neutral systems we study, which is defined as the difference between the total energy of the system and the sum over its unrelaxed components. ^{[ 62 ]} The larger the red shift of the stretching frequency (measured as the ratio of the X –H stretching frequency in the H-bonded cluster to that in the free monomer), the stronger the HB.

Fig. 2.

Figure Option View Download New Window

Fig. 2. (a) The differences between the shortest heavy-atom distances obtained from the PIMD and MD simulations , denoted by Δ ( X – X ), are plotted as a function of the HB strength index. Δ ( X – X ) characterizes the impact of the NQEs on the strength of the HBs. The HB strength is defined as the ratio of the X –H stretching frequency in the H-bonded system to that in the free monomer. (b) The correlation between this HB strength index and the binding energy per HB in the neutral systems. (c) A simplified schematic illustration of the isotope (Ubbelohde) effect. We suggest that regimes of positive, negligible, and negative Ubbelohde effects exist depending on the HB strength. For the HF clusters, labels 1–5 denote the HBs in the dimer to the hexamer. For the H 2 O clusters, labels 1, 2, 3a, and 3b refer to the HBs in the dimer, pentamer, and the long and short HBs in the octamer, respectively. For the charged clusters, labels 1–4 refer to , , , and . For the organic dimers, labels 1a, 1b, and 2 mean the red shifted and blue shifted HBs in the formamide and the red shifted HB in formic acid. For the solids, labels 1–3 refer to the HBs in HCl, HF, and squaric acid. The same labels also apply to Fig.  4 . [ 7 ]

With the definition of the above-mentioned quantities in mind, we first look at the results for the impact of the NQEs on the strength of HBs. Upon comparing these results for the various H-bonded systems, an interesting correlation can be observed between the H-bond strength and the change in intermolecular separations, as shown in Fig.  2(a) . Here, it can been seen that as the HB gets stronger, the heavy-atom separations in the PIMD simulations go from being longer than those in the MD simulations (positive Δ ( X – X )) to being shorter (negative Δ ( X – X )). In other words, the NQEs lead to longer HBs in weak H-bonded systems and shorter HBs in relatively strong H-bonded systems. We note that the HB strength increases upon going from small to large clusters and from water to HF. This trend is the main finding, and in the following we explain why it exists and discuss the implications it has for H-bonded systems in general.

To this end, we take the HF clusters (dimer, tetramer, and pentamer) as the guiding example, because upon increasing the cluster size, the HB strength increases, and the influence of the NQEs switches from a tendency to lengthen to a tendency to shorten the intermolecular separations (as seen in Ref. [ 65 ]). The results are summarized in Fig.  3 , where we plot the distance and angle distributions from MD and PIMD simulations for these three clusters separately. The left column shows the final results, where we can see that, in the dimer, the averaged F–F distance is increased by including the NQEs. In the tetramer, there is no difference between the averaged quantum and classical F–F distances. While in the pentamer, the F–F distance is clearly shortened by including the NQEs. The key to understanding this variation of the heavy-atom distances is realizing that there are also related differences between MD and PIMD in the covalent F–H bond lengths (center) and H-bond angles (right). Because of anharmonic quantum fluctuations, these two geometric properties also show systematic changes. The F–H bonds are longer in the quantum (PIMD) than in the classical (MD) simulations, and this elongation becomes more pronounced as the HBs get stronger. Besides this, the HBs are more bent in the quantum than in the classical simulations, and this bending generally becomes less pronounced as the HBs get stronger. In order to understand the influence of these variations in structure, analysis of various dimer configurations was performed. This analysis reveals that the covalent bond stretching increases the intermolecular interaction, whereas the HB bending decreases it. Taking the HF dimer as an example, a 0.04 Å increase in the F–H bond length of the donor leads to a 40 meV increase in interaction energy within the dimer; in contrast, a 21° reduction in the H-bond angle leads to a 16 meV decrease in interaction energy. In short, the F–F distance increases in the dimer as a result of a large decrease in the HB angle, but only a small increase in the covalent F-H bond length is noted upon including the NQEs. Whilst in the tetramer, the F–H stretching is sufficiently pronounced to compensate for the increase in HB bending, leaving the overall F–F distance unchanged; and in the pentamer, the F–F distance decreases because the F–H covalent bond stretching dominates over the HB bending.

Fig. 3.

Figure Option View Download New Window

Fig. 3. HF clusters as examples for detailed analysis of the QNEs. Distributions of the F–F distances (left), the F–H bond lengths (center), and the intermolecular bending (F–H···F angle, right) from MD (solid black lines) and PIMD (dashed red lines) for a selection of systems: the HF dimer (top), the HF tetramer (middle), and the HF pentamer (bottom). The MD and PIMD averages are shown in black and red vertical dashes, respectively. [ 7 ]

The above analysis provides a qualitative understanding of the trend observed. For a more rigorous examination of this picture and a quantitative description of this competition for all systems studied, we further calculate the projection ( X –H ^|| ) of the donor molecule’s covalent bond along the intermolecular axis (see the inset of Fig.  4 ). Since X –H ^|| increases upon intramolecular stretching but decreases upon intermolecular bending, this quantity itself allows the balance between stretching and bending to be evaluated to a certain extent. The influence of the NQEs is quantified by the ratio of the PIMD and MD projections, i.e., x = ( X – H ^|| ) ^PIMD /( X – H ^|| )MD. When this value is clearly greater than one, it indicates that when the NQEs are included, the main influence is on the stretching of the covalent bond, and when this value is clearly smaller than one, it indicates that when the NQEs are included, the main influence is on the bending of the HB. When one plots this ratio against the variations in intermolecular separations, y = Δ ( X – X ) (which we used to quantify the impact of the NQEs), a striking (almost linear) correlation is observed (Fig.  4 ). For all systems in which HB bending dominates ( x clearly smaller than 1), the heavy-atom distances are longer in PIMD than in MD ( y > 0). In cases where covalent bond stretching is dominant ( x clearly larger than 1), the heavy-atom distances are shorter in PIMD than in MD ( y < 0). With the increase of x , quantum fluctuations on the stretching mode become more dominant and the NQEs switch from weakening the hydrogen bonds to strengthening them. Thus, the overall influence of the NQEs on the HB interaction quantitatively comes down to this delicate interplay between covalent bond stretching and intermolecular bond bending. One notes that this explanation for the general case is consistent with what Manolopoulos and coworkers have elegantly shown for liquid water in Ref. [ 4 ].

Fig. 4.

Figure Option View Download New Window

Fig. 4. Quantification of the competition between the quantum fluctuations on the stretching and bending modes. The differences in average shortest heavy-atom distances between PIMD and MD simulations ( Δ ( X – X ), vertical axis) are plotted as a function of the ratio of the projection of the donor X –H covalent bond along the intermolecular axis from PIMD and MD simulations (horizontal axis). x larger (smaller) than 1 indicates a dominant contribution from the stretching (bending) mode when the NQEs are included. Negative values of Δ ( X – X ) indicate that NQEs decrease the intermolecular separation. An almost linear correlation between the two variables is observed. The inset illustrates the geometry used for projecting the donor covalent X –H bond onto the intermolecular axis. The curved red arrow represents the intermolecular bending and the straight blue arrow represents the intramolecular stretching. [ 7 ]

In Section 3, we have analyzed the influence of NQEs on the strength of HBs, taking the geometric changes as the indicator. A rigorous quantification of this influence, however, requires an explicit calculation of the contribution of the NQEs to the free energy of the H-bonded systems, using the method shown in Section 2. In Fig.  5 , we use two examples to show how such TI simulations for the calculation of Δ F [ n H ₂ O]/ n are carried out. The left panel corresponds to the results of water monomer and the right panel corresponds to those of the dimer. F ′( λ ) is as defined in Eq. ( 9 ), but divided by the number of water molecules in the simulation cell. Integration of this quantity from 0 to 1 gives us the final Δ F [ n H ₂ O]/ n . The temperature is 100 K. After integration, Δ F [ n H ₂ O]/ n equals 0.446 eV and 0.462 eV respectively in the monomer and dimer. To separate out the harmonic and anharmonic effects, we can perform a separate phonon frequency calculation for each system and integrate out the free-energy change due to NQEs at the harmonic limit through

Fig. 5.

	Figure Option View Download New Window
	Fig. 5. Thermodynamic integration based on ab initio PIMD sampling for the calculation of the free-energy change due to NQEs in the water monomer and dimer. F ′( λ ) is as defined in Eq. ( 9 ), but divided by the number of water molecules in the simulation cell. Integration of this quantity from 0 to 1 gives us the final Δ F [ n H 2 O]/ n . [ 57 ]

HB is generally weak, but ubiquitous and essential to life on earth. The small mass of hydrogen means that many properties of HBs are quantum mechanical in nature, although quite often their theoretical descriptions remain classical. Due to the development of computer simulation methods, this influence of NQEs on the structural and energetic properties of some hydrogen bonded systems has been systematically studied in recent years. Here we present a review of these works, by focussing on the explanation of the main methods behind the simulations, i.e., the ab initio PIMD method as well as its extension in combination with the thermodynamic integration method for the calculation of the free energies. As a practical example of these studies, ab initio PIMD simulations on the structural properties of the hydrogen bonded systems have shown that the NQEs have a tendency to strengthen the strong HBs and weaken the weak ones, which could be explained by a competition of the quantum fluctuations on the stretching and bending modes.