Hydration of NH$_4^+$ in Water: Bifurcated Hydrogen Bonding Structures and Fast Rotational Dynamics

Understanding the hydration and diffusion of ions in water at the molecular level is a topic of widespread importance. The ammonium ion (NH$_4^+$) is an exemplar system that has received attention for decades because of its complex hydration structure and relevance in industry. Here we report a study of the hydration and the rotational diffusion of NH$_4^+$ in water using ab initio molecular dynamics simulations and quantum Monte Carlo calculations. We find that the hydration structure of NH$_4^+$ features bifurcated hydrogen bonds, which leads to a rotational mechanism involving the simultaneous switching of a pair of bifurcated hydrogen bonds. The proposed hydration structure and rotational mechanism are supported by existing experimental measurements, and they also help to rationalize the measured fast rotation of NH$_4^+$ in water. This study highlights how subtle changes in the electronic structure of hydrogen bonds impacts the hydration structure, which consequently affects the dynamics of ions and molecules in hydrogen bonded systems.

The hydration and diffusion of ions and molecules in water is one of the most fundamental processes in nature and in modern technology, having a direct impact on nucleation and crystallization, ion sieving, and aqueous chemical reactions, to name just a few examples [1][2][3][4][5]. A prototypical example of ion solvation with hydrogen bonding is the ammonium ion (NH þ 4 ) in water. This system has a complex and debated hydration structure. In simulations the coordination number of NH þ 4 is around five [6][7][8], whereas experiments indicated that the coordination number is larger [9]. In addition, there is an interesting but unresolved experimental observation that NH þ 4 rotates rapidly in water [10,11]. The estimated rotational diffusion constant from nuclear magnetic resonance measurements is a few times larger than theoretical estimates.
In 1999, Brugé et al. proposed a mechanism involving discontinuous rotational jumps associated with the exchange of water molecules [6]. Within this model, the NH þ 4 ion forms four long-lived hydrogen bonds with water, and exchange occurs between the fifth water molecule and the four water molecules bonded with NH þ 4 in the first hydration shell. In 2005, Intharathep et al. suggested that NH þ 4 undergoes free rotation due to a flexible hydration structure and a large coordination number [12]. However, neither of the suggested mechanisms are supported by accurate investigations of the water-ammonium ion interactions and of the underlying potential energy surface.
In order to correctly simulate the diffusion of ions in water it is essential to use a theoretical method that is able to reliably describe hydrogen bonds and the energetics involved [6,13,14]. In recent years, many breakthroughs have been reported in understanding the hydrogen bonding structure of water and the dynamics of ions in water using molecular dynamics simulations based on density functional theory (DFT) [13,[15][16][17][18][19]. However, the outcome of DFT crucially depends on the choice of the exchange correlation functional [17]. Therefore, performing molecular dynamics simulations based on well validated exchange correlation functionals is desired. In order to validate DFT, high level theories should be applied to evaluate the interactions in the NH þ 4 solution. Recently, quantum Monte Carlo methods, especially fixed node diffusion Monte Carlo (DMC) calculations, have been used in high-pressure ice, bulk water, molecular crystals, and other extended systems to provide highly accurate benchmarks [20][21][22][23][24]. Once the questions over the bonding and the hydration structure are answered, one may be able to clarify the rotational mechanism.
In this study, we carried out ab initio molecular dynamics (AIMD) simulations of an NH þ 4 aqueous solution using three different state of the art exchange correlation functionals. We then performed DMC calculations to benchmark the accuracy of each functional and identify the most reliable functional(s). We find that DMC calculations and DFT, with the strongly constrained and appropriately normed (SCAN) functional, correctly reproduce the energy ordering of the hydration structures, where the most favorable coordination consists of six water molecules in the first hydration shell, which is in good agreement with experimental indications. In the stable hydration structure, multiple water molecules appear at the so-called bifurcated positions (a bifurcated position is a position between two protons of NH þ 4 ), forming bifurcated nonlinear hydrogen bonds (a nonlinear hydrogen bond has a large hydrogen bond angle). Such a bifurcated hydration structure directly leads to the fast rotation of NH þ 4 in water, and the rotational diffusion constant calculated is in good agreement with experimental measurements [11]. We also find that there is a relation between the rotational diffusion constant and the number of bifurcated hydrogen bonds predicted with different exchange correlation functionals.
DFT-based AIMD simulations were carried out using the Quantum ESPRESSO package [25,26] with the SCAN [27], the Perdew-Burke-Ernzerhof (PBE) [28], and the PBE þ vdW [29] functionals. The interactions between the valence electrons were treated with Hamann-Schlüter-Chiang-Vanderbilt pseudopotentials [30,31]. Each simulation was performed in the canonical (NVT) ensemble, and the temperature was controlled by a Nosé-Hoover thermostat [32]. Hydrogen atoms were replaced with deuterium in order to reduce nuclear quantum effects and to maximize the time step in the integration of the equations of motion. DMC calculations were performed using the CASINO package [33], with the size-consistent DMC algorithm ZSGMA [34]. The recently developed energy consistent correlated electron pseudopotentials [35] were used. To reduce or eliminate the finite size error (FSE), we used the model periodic Coulomb (MPC) interaction method [36][37][38]. In previous studies, the above setup has been validated for noncovalent interactions through extensive comparisons with converged coupled cluster calculations and experimental evaluations [22,39,40]. Further computational details can be found in the Supplemental Material [41], which includes additional Ref. [42].
We begin our study by discussing the hydration structures from AIMD simulations, as obtained using SCAN, PBE, and PBE þ vdW. Figure 1(a) shows the N-O radial distribution function (RDF) G NO ðrÞ and the corresponding running coordination number (CN) n NO ðrÞ, where CN is defined as the average number of water molecules in the first hydration shell. The first peak in G NO ðrÞ with SCAN is higher and broader than those with PBE and PBE þ vdW, leading to a significantly larger CN of around 6.6, while the  Fig. 1(b) further demonstrates the differences of CN between SCAN and the two other functionals. The most favorable coordination number of SCAN is six followed by seven and five, while the CN distribution of PBE and PBE þ vdW peak at four and five, respectively. As for the second hydration shell, PBE and PBE þ vdW both have a clear second peak on G NO ðrÞ, whereas with SCAN the second hydration shell mixes with the first shell. The flattening of the second peak predicted by SCAN has been observed in neutron diffraction and x-ray experiments [9,43].
We now consider the number of hydrogen bonds formed by NH þ 4 and the surrounding water molecules. Figure 1(c) plots the N N -O radial distribution functions from AIMD. We define the standard hydrogen bonds with criteria reported for liquid water and aqueous solutions [44][45][46], i.e., where R NO  Fig. 1(d). We find the number of NH þ 4 -water hydrogen bonds predicted by PBE and PBE þ vdW are almost identical with an average value of 3.50 and 3.48, respectively, suggesting a minor influence of vdW interactions on the NH þ 4 -water hydrogen bonds. SCAN, however, predicts considerably weaker NH þ 4water hydrogen bonds evidenced by a smaller average number of 3.24. By ignoring the angular restraint in Eq. (1), i.e., one can define the so-called generalized hydrogen bonds, including nonlinear, bifurcated hydrogen bonds that have β > 30° [ 47,48]. We find that SCAN predicts a larger number of bifurcated NH þ 4 -water hydrogen bonds. The hydration structure of NH þ 4 in water can also be demonstrated by the joint distributions of r and β, which are plotted in Fig. S3 [41]. The peaks at β < 30°and β ≈ 100°represent the water molecules that form hydrogen bonds with NH þ 4 and the increased intensity in between can be attributed to the bifurcated water molecules that tend to form multiple nonlinear hydrogen bonds with the NH þ 4 . To highlight the difference between different AIMD simulations, in Fig. 1(e) we plot the difference of the joint distribution of SCAN and PBE. Because of the weaker NH þ 4 -water hydrogen bonds and increased bifurcated water molecules predicted by SCAN, the peaks at β < 30°and β ≈ 100°are reduced while the intensity in the middle becomes stronger. Figure 1(f) further plots the N-water total distribution function [G N ðrÞ] obtained from our simulations and from available neutron diffraction data. Experimental G N ðrÞ, as defined in Eq. (2), is different from g NO ðrÞ in Fig. 1(a) and involves g NN and g NCl terms because experiments used 5.0 mol=L NH 4 Cl water solutions [6,9].
Our simulated G N ðrÞ, however, contains only the first two terms in Eq. (2). The difference in solution conditions is apparent from the difference between curves at large radii. Nevertheless, at small radii, the comparison indicates that the SCAN result is reasonable, in particular the position of the second maximum r max 2 and the second minimum r min 2 are in line with the experiment. Simulations with an inaccurate functional such as PBE do not correctly reproduce these experimental features. In addition, Figs. S4 and S5 [41] present vibrational spectra and density of states computed from our simulations, which may be compared in future experiments to confirm the proposed hydration structure.
Overall, we find that the hydration structure predicted by SCAN is more consistent with experiments than that obtained with the other two functionals, with the SCAN results characterized by weaker NH þ 4 -water hydrogen bonds and an increased number of water molecules in the first hydration shell. However, to confirm that the SCAN description of the hydration structure is indeed reliable and the agreement between theory and experiment is not merely fortuitous, it is necessary to examine the accuracy of the underlying potential energy surface with the help of a higher level of theory. To this end we selected five distinct configurations with CN from four to eight to benchmark against DMC calculations, as shown in Figs. 2(a-e). Relative energies of these configurations are presented in Fig. 2(f). According to DMC, the configuration with CN ¼ 6 is the lowest energy state among the selected structures. The basin shape of the DMC predicted energy curve is reproduced quite well by SCAN, indicating a good description of the hydration structure and the hydrogen bonding network around NH þ 4 . Indeed, the energy ordering of the selected structures predicted by DMC and SCAN is qualitatively consistent with the distribution of CNs shown in Fig. 1(b). Regarding the performance of other DFT functionals, typical results are shown with PBE and PBE þ vdW. From this it can be seen that they predict quite different energy ordering for the same set of structures. Specifically, both PBE and PBE þ vdW tend to underestimate the stability of the CN ¼ 6 configuration. Furthermore, SCAN correctly predicts that the configuration of CN ¼ 6 is the lowest energy state while PBE and PBE þ vdW both overestimate the stability of the state of CN ¼ 5. In Fig. S6 [41], we show the results calculated using several other exchange correlation functionals and none of them predict a better energy ordering. Benchmarks on the gas phase NH þ 4 -water clusters presented in Fig. 2(g) further support our conclusions. It is worth noting that recent experimental measurements also suggested that hydrogen bonds of solvated NH þ 4 are less PHYSICAL REVIEW LETTERS 125, 106001 (2020) strong than previous theoretical predictions [8]. Here our DMC and DFT calculations clarify that the NH þ 4 -water interaction is indeed weaker than previously thought, and the relatively weak interaction leads a large coordination number of NH þ 4 and an increased number of nonlinear hydrogen bonds. Theoretically, the relatively weak hydrogen bonds predicted by SCAN can be understood by the lower polarizability of the water molecules around NH þ 4 (Fig. S5); a suggestion that is consistent with previous studies in water [18,19].
Having established that the hydration structure of NH þ 4 is characterized by weak and bifurcated hydrogen bonds, and having found that the SCAN functional can reproduce this structure, we now discuss the rotation of NH þ 4 in water. To aid the quantification of the rotation, we define a vector in the body-fixed frame and track the rotation of the vector as a function of time. The angle θ is defined as the angle between the initial vector and the vector after a given time. Figure 3(a) shows a typical evolution of θ in a picosecond window when the rotation of NH þ 4 occurs, which is described by a big change of θ of 80-90°in 0.3-0.4 ps. The details of the rotation are described in Fig. S7 and Movie S1 of the Supplemental Material [41].

Figures 3(b)-3(c)
show two snapshots just before and after the jump to highlight the key step of the rotation, namely the simultaneous switching of two bifurcated hydrogen bonds. In Figs. 3(b) and 3(c), W 1 to W 6 are the six water molecules in the first hydration shell of NH þ 4 at T ¼ 17.0 ps. Initially two bifurcated water molecules (W 5 and W 6 ) form nonlinear hydrogen bonds with NH þ 4 , hence H 3 and H 4 form bifurcated hydrogen bonds with W 3 =W 5 and W 4 =W 6 , respectively. After the simultaneous switching of the two bifurcated hydrogen bonds, N-H 3 forms a stable hydrogen bond with W 5 and N-H 4 points towards W 6 . Because of the bifurcated hydrogen bonds, such a mechanism does not involve complete breaking of hydrogen bonds, hence facilitating the rotation. In addition, the rotation of NH þ 4 involves the rotation of two N-H bonds instead of one, thus it is much more efficient when there is a pair of bifurcated hydrogen bonds. A statistical analysis of rotation times and the extent of the rotational angle jumps is provided in Fig. S8, confirming that large rotational jumps dominate the overall rotational process. In contrast, with PBE and its hydration shell containing fewer bifurcated water molecules, the rotation involves the breaking of hydrogen bonds, and the process is slowed down (Fig. S9). (g) Relative energies of gas phase NH þ 4 clusters calculated using different methods and plotted against the DMC results. The relative energy is defined as the energy difference to the energy of the configuration with the lowest one. The gas phase NH þ The above rotational process is characterized by the rotational diffusion constant D R [7,12], which can be derived from the mean-square angular displacements according to The computed hθðtÞ 2 i as a function of time t are shown in Fig. S10(a) [41]. From the linear part of hθðtÞ 2 i curve, we can estimate the rotational diffusion constant D R which is a quarter of the slope. The estimate for D R using SCAN is 0.156 AE 0.009 rad 2 ps −1 at 363 K (the convergence of D R was shown in Fig. S11) while the values from PBE (0.050 AE 0.007) and PBE þ vdW (0.073 AE 0.007) simulations are significantly smaller. In Fig. 4(a) we plot D R as a function of temperature along with the experimental values derived from the NMR measurements [7,[10][11][12]. We find that D R follows a good linear relation in the temperature regime studied, and a linear fit is highlighted with the purple dashed line for NH þ 4 in H 2 O. Through the experimental data of ND þ 4 in D 2 O we draw a green dashed line parallel to the purple line, which suggests that our simulated D R with the SCAN functional falls in line with experiments. We note that to explicitly discuss the difference between the hydrogenated and deuterated systems, nuclear quantum effects should be included using, e.g., path integral molecular dynamics in the future [49,50]. Nevertheless, from the experimental measurements it is evident that nuclear quantum effects have a minor impact, of approximately 0.02 rad 2 ps −1 , on the rotational diffusion constant [11].
To quantitatively demonstrate the correlation between hydration structure and rotational dynamics, we further plot D R at 363 K as a function of the average number of bifurcated hydrogen bonds in Fig. 4(b). We find D R grows as the number of bifurcated hydrogen bonds increases, following a linear relation. Similar results are obtained at other temperatures, which are presented in Fig. S10 [41]. The linear relation established offers strong evidence for our finding that the rotation involves the simultaneous switching of bifurcated hydrogen bonds.
To conclude, we have performed AIMD simulations and DMC calculations of the NH þ 4 aqueous solution. We have clarified the hydration structure of NH þ 4 in water, finding that it is characterized by a coordination number of approximately six, with two pairs of bifurcated hydrogen bonds. Consequently, we find that this bonding nature leads to fast rotation of NH þ 4 in water. A clear prediction from this study is that the rotation of NH þ 4 may be significantly suppressed in a different solvent where bifurcated hydrogen bonds do not form or when it is confined to an extent that the hydration structure is partially broken. In the future, such studies are desirable to verify our conclusions and confirm the diffusion mechanism suggested here. Last but not the least, this study brings attention to the importance of achieving accurate electronic structures of other aqueous and hydrogen bonded systems, and further improvements in computational methods will gradually bring us towards an exact description of such complex processes.