Hydration of NH + 4 in Water: Bifurcated Hydrogen Bonding Structures and Fast Rotational Dynamics
Jianqing Guo, Liying Zhou, Andrea Zen, Angelos Michaelides, Xifan Wu, Enge Wang, Limei Xu, Ji Chen
Hydration of NH +4 in Water: Bifurcated Hydrogen Bonding Structuresand Fast Rotational Dynamics Jianqing Guo ,
Liying Zhou,
Andrea Zen ,
Angelos Michaelides,
Xifan Wu, Enge Wang,
Limei Xu, † and Ji Chen ‡ International Center for Quantum Materials, Peking University, Beijing 100871, People ’ s Republic of China School of Physics, Peking University, Beijing 100871, People ’ s Republic of China Department of Physics and Astronomy, Thomas Young Centre and London Centre for Nanotechnology University College London,Gower Street, London WC1E 6BT, United Kingdom Department of Earth Sciences, University College London, Gower Street, London WC1E 6BT, United Kingdom Max Planck Institute for Solid State Research, Stuttgart 70569, Germany Department of Physics, Temple University, Philadelphia, Pennsylvania 19122, USA Collaborative Innovation Center of Quantum Matter, Beijing 100871, People ’ s Republic of China Songshan Lake Materials Lab, Institute of Physics, Chinese Academy of Sciences, Guangdong 523808, People ’ s Republic of China School of Physics, Liaoning University, Shenyang 110136, People ’ s Republic of China (Received 9 April 2020; accepted 4 August 2020; published 1 September 2020)Understanding the hydration and diffusion of ions in water at the molecular level is a topic of widespreadimportance. The ammonium ion (NH þ ) is an exemplar system that has received attention for decadesbecause of its complex hydration structure and relevance in industry. Here we report a study of thehydration and the rotational diffusion of NH þ in water using ab initio molecular dynamics simulations andquantum Monte Carlo calculations. We find that the hydration structure of NH þ features bifurcatedhydrogen bonds, which leads to a rotational mechanism involving the simultaneous switching of a pair ofbifurcated hydrogen bonds. The proposed hydration structure and rotational mechanism are supported byexisting experimental measurements, and they also help to rationalize the measured fast rotation of NH þ inwater. This study highlights how subtle changes in the electronic structure of hydrogen bonds impacts thehydration structure, which consequently affects the dynamics of ions and molecules in hydrogen bondedsystems. DOI: 10.1103/PhysRevLett.125.106001
The hydration and diffusion of ions and molecules inwater is one of the most fundamental processes in natureand in modern technology, having a direct impact onnucleation and crystallization, ion sieving, and aqueouschemical reactions, to name just a few examples [1 – þ ) in water. This systemhas a complex and debated hydration structure. In simu-lations the coordination number of NH þ is around five[6 – þ rotates rapidly in water [10,11]. The estimated rotationaldiffusion constant from nuclear magnetic resonance mea-surements is a few times larger than theoretical estimates. In 1999, Brug´e et al. proposed a mechanism involvingdiscontinuous rotational jumps associated with theexchange of water molecules [6]. Within this model, theNH þ ion forms four long-lived hydrogen bonds with water,and exchange occurs between the fifth water molecule andthe four water molecules bonded with NH þ in the firsthydration shell. In 2005, Intharathep et al. suggested thatNH þ undergoes free rotation due to a flexible hydrationstructure and a large coordination number [12]. However,neither of the suggested mechanisms are supported byaccurate investigations of the water-ammonium ion inter-actions and of the underlying potential energy surface.In order to correctly simulate the diffusion of ions inwater it is essential to use a theoretical method that isable to reliably describe hydrogen bonds and the energeticsinvolved [6,13,14]. In recent years, many breakthroughshave been reported in understanding the hydrogen bondingstructure of water and the dynamics of ions in water usingmolecular dynamics simulations based on density func-tional theory (DFT) [13,15 – Published by the American Physical Society under the terms ofthe Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution tothe author(s) and the published article ’ s title, journal citation,and DOI. Open access publication funded by the Max PlanckSociety. PHYSICAL REVIEW LETTERS = = = þ solution. Recently,quantum Monte Carlo methods, especially fixed nodediffusion Monte Carlo (DMC) calculations, have beenused in high-pressure ice, bulk water, molecular crystals,and other extended systems to provide highly accuratebenchmarks [20 – ab initio moleculardynamics (AIMD) simulations of an NH þ aqueous solu-tion using three different state of the art exchange corre-lation functionals. We then performed DMC calculationsto benchmark the accuracy of each functional and identifythe most reliable functional(s). We find that DMC calcu-lations and DFT, with the strongly constrained and appro-priately normed (SCAN) functional, correctly reproducethe energy ordering of the hydration structures, where themost favorable coordination consists of six water moleculesin the first hydration shell, which is in good agreementwith experimental indications. In the stable hydrationstructure, multiple water molecules appear at the so-calledbifurcated positions (a bifurcated position is a positionbetween two protons of NH þ ), forming bifurcated non-linear hydrogen bonds (a nonlinear hydrogen bond has alarge hydrogen bond angle). Such a bifurcated hydrationstructure directly leads to the fast rotation of NH þ in water,and the rotational diffusion constant calculated is in goodagreement with experimental measurements [11]. Wealso find that there is a relation between the rotationaldiffusion constant and the number of bifurcated hydrogenbonds predicted with different exchange correlationfunctionals.DFT-based AIMD simulations were carried out usingthe Q uantum ESPRESSO package [25,26] with the SCAN [27],the Perdew-Burke-Ernzerhof (PBE) [28], and the PBE þ vdW [29] functionals. The interactions between the valenceelectrons were treated with Hamann-Schlüter-Chiang-Vanderbilt pseudopotentials [30,31]. Each simulation wasperformed in the canonical ( NVT ) ensemble, and thetemperature was controlled by a Nos´e-Hoover thermostat[32]. Hydrogen atoms were replaced with deuterium inorder to reduce nuclear quantum effects and to maximizethe time step in the integration of the equations of motion.DMC calculations were performed using the
CASINO pack-age [33], with the size-consistent DMC algorithm ZSGMA[34]. The recently developed energy consistent correlatedelectron pseudopotentials [35] were used. To reduce oreliminate the finite size error (FSE), we used the modelperiodic Coulomb (MPC) interaction method [36 – þ vdW. Figure 1(a) shows the N-O radialdistribution function (RDF) G NO ð r Þ and the correspondingrunning coordination number (CN) n NO ð r Þ , where CN isdefined as the average number of water molecules in thefirst hydration shell. The first peak in G NO ð r Þ with SCAN ishigher and broader than those with PBE and PBE þ vdW,leading to a significantly larger CN of around 6.6, while the (a) SCANPBEPBE+vdW
Coordination number P e r ce n t a g e ( % ) SCANPBEPBE+vdW (b)
XC Functional A v e r a g e N u m b e r o f H B s HBGeneralized HB r (Å) G NO ( r ) n NO ( r ) (d) r (Å) G N ( r ) r r SCAN (300 K)PBE (300 K)EXP (300 K) (e) (f) (c) r (Å) G H N O ( r ) (c) SCANPBEPBE+vdW
FIG. 1. Hydration structures of NH þ in water obtainedfrom AIMD. (a) The N-O radial distribution function (solidline) and the corresponding running coordination number(dashed line). The running coordination number n NO ð R Þ ¼ πρ R R G NO ð r Þ r dr , where ρ is the density. The first minimum( r min ) of G NO ð r Þ is 3.64 Å, 3.38 Å, 3.55 Å for SCAN, PBE,PBE þ vdW, respectively. The coordination number (CN) is n ð r min Þ . The convergence of G NO ð r Þ is shown in Fig. S2(a).(b) The distribution of the coordination number of the firsthydration shell. (c) Calculated H N -O radial distribution func-tions. The first minimum of G H N O ð r Þ is 2.40 Å, 2.45 Å, 2.39 Å forSCAN, PBE, PBE þ vdW, respectively. (d) Average number ofNH þ -water hydrogen bonds (HBs) predicted by differentexchange correlation (XC) functionals. (e) The difference ofthe joint distribution in the first hydration shell between theSCAN and the PBE results as a function of r and β [defined in theinset of (d)]. Red color implies a larger value from SCAN. Thejoint distributions are shown in Fig. S3. (f) The N -water totaldistribution function G N ð r Þ . The blue solid line is from theneutron diffraction work of Hewish et al. [43]. The grey barsindicate the positions of the second maximum and the secondminimum of the experimental G N ð r Þ . PHYSICAL REVIEW LETTERS125,
FIG. 1. Hydration structures of NH þ in water obtainedfrom AIMD. (a) The N-O radial distribution function (solidline) and the corresponding running coordination number(dashed line). The running coordination number n NO ð R Þ ¼ πρ R R G NO ð r Þ r dr , where ρ is the density. The first minimum( r min ) of G NO ð r Þ is 3.64 Å, 3.38 Å, 3.55 Å for SCAN, PBE,PBE þ vdW, respectively. The coordination number (CN) is n ð r min Þ . The convergence of G NO ð r Þ is shown in Fig. S2(a).(b) The distribution of the coordination number of the firsthydration shell. (c) Calculated H N -O radial distribution func-tions. The first minimum of G H N O ð r Þ is 2.40 Å, 2.45 Å, 2.39 Å forSCAN, PBE, PBE þ vdW, respectively. (d) Average number ofNH þ -water hydrogen bonds (HBs) predicted by differentexchange correlation (XC) functionals. (e) The difference ofthe joint distribution in the first hydration shell between theSCAN and the PBE results as a function of r and β [defined in theinset of (d)]. Red color implies a larger value from SCAN. Thejoint distributions are shown in Fig. S3. (f) The N -water totaldistribution function G N ð r Þ . The blue solid line is from theneutron diffraction work of Hewish et al. [43]. The grey barsindicate the positions of the second maximum and the secondminimum of the experimental G N ð r Þ . PHYSICAL REVIEW LETTERS125, ¼ . ) and PBE þ vdW (CN ¼ . )are smaller. The distribution of CNs shown in Fig. 1(b)further demonstrates the differences of CN between SCANand the two other functionals. The most favorable co-ordination number of SCAN is six followed by seven andfive, while the CN distribution of PBE and PBE þ vdWpeak at four and five, respectively. As for the secondhydration shell, PBE and PBE þ vdW both have a clearsecond peak on G NO ð r Þ , whereas with SCAN the secondhydration shell mixes with the first shell. The flattening ofthe second peak predicted by SCAN has been observed inneutron diffraction and x-ray experiments [9,43].We now consider the number of hydrogen bonds formedby NH þ and the surrounding water molecules. Figure 1(c)plots the N N -O radial distribution functions from AIMD. Wedefine the standard hydrogen bonds with criteria reported forliquid water and aqueous solutions [44 – R NO < R NO C ; R H N O < R H N O C ; β < ° ; ð Þ where R NO C and R H N O C are the cut-off values for the N-O andH N -O distances obtained from the first minimum of thecorresponding radial distribution functions. The angle β isthe H N -N-O angle as shown in the inset of Fig. 1(d). We findthe number of NH þ -water hydrogen bonds predicted byPBE and PBE þ vdW are almost identical with an averagevalue of 3.50 and 3.48, respectively, suggesting a minorinfluence of vdW interactions on the NH þ -water hydrogenbonds. SCAN, however, predicts considerably weaker NH þ -water hydrogen bonds evidenced by a smaller average numberof 3.24. By ignoring the angular restraint in Eq. (1), i.e., onecan define the so-called generalized hydrogen bonds, includ-ing nonlinear, bifurcated hydrogen bonds that have β > °[47,48]. We find that SCAN predicts a larger number ofbifurcated NH þ -water hydrogen bonds.The hydration structure of NH þ in water can also bedemonstrated by the joint distributions of r and β , whichare plotted in Fig. S3 [41]. The peaks at β < ° and β ≈ ° represent the water molecules that form hydrogenbonds with NH þ and the increased intensity in between canbe attributed to the bifurcated water molecules that tend toform multiple nonlinear hydrogen bonds with the NH þ . Tohighlight the difference between different AIMD simula-tions, in Fig. 1(e) we plot the difference of the jointdistribution of SCAN and PBE. Because of the weakerNH þ -water hydrogen bonds and increased bifurcated watermolecules predicted by SCAN, the peaks at β < ° and β ≈ ° are reduced while the intensity in the middlebecomes stronger.Figure 1(f) further plots the N -water total distributionfunction [ G N ð r Þ ] obtained from our simulations and fromavailable neutron diffraction data. Experimental G N ð r Þ , asdefined in Eq. (2), is different from g NO ð r Þ in Fig. 1(a) andinvolves g NN and g NCl terms because experiments used . mol = L NH Cl water solutions [6,9]. G N ð r Þ ¼ . g NO ð r Þ þ . g NH ð r Þþ . g NCl ð r Þ þ . g NN ð r Þ − . : ð Þ Our simulated G N ð r Þ , however, contains only the first twoterms in Eq. (2). The difference in solution conditions isapparent from the difference between curves at large radii.Nevertheless, at small radii, the comparison indicates that theSCAN result is reasonable, in particular the position of thesecond maximum r max and the second minimum r min are inline with the experiment. Simulations with an inaccuratefunctional such as PBE do not correctly reproduce theseexperimental features. In addition, Figs. S4 and S5 [41]present vibrational spectra and density of states computedfrom our simulations, which may be compared in futureexperiments to confirm the proposed hydration structure.Overall, we find that the hydration structure predicted bySCAN is more consistent with experiments than thatobtained with the other two functionals, with the SCANresults characterized by weaker NH þ -water hydrogenbonds and an increased number of water molecules inthe first hydration shell. However, to confirm that theSCAN description of the hydration structure is indeedreliable and the agreement between theory and experimentis not merely fortuitous, it is necessary to examine theaccuracy of the underlying potential energy surface with thehelp of a higher level of theory. To this end we selectedfive distinct configurations with CN from four toeight to benchmark against DMC calculations, as shown inFigs. 2(a – e). Relative energies of these configurations arepresented in Fig. 2(f). According to DMC, the configura-tion with CN ¼ is the lowest energy state among theselected structures. The basin shape of the DMC predictedenergy curve is reproduced quite well by SCAN, indicatinga good description of the hydration structure and thehydrogen bonding network around NH þ . Indeed, theenergy ordering of the selected structures predicted byDMC and SCAN is qualitatively consistent with thedistribution of CNs shown in Fig. 1(b). Regarding theperformance of other DFT functionals, typical results areshown with PBE and PBE þ vdW. From this it can be seenthat they predict quite different energy ordering for thesame set of structures. Specifically, both PBE and PBE þ vdW tend to underestimate the stability of the CN ¼ configuration. Furthermore, SCAN correctly predicts thatthe configuration of CN ¼ is the lowest energy statewhile PBE and PBE þ vdW both overestimate the stabilityof the state of CN ¼ . In Fig. S6 [41], we show the resultscalculated using several other exchange correlation func-tionals and none of them predict a better energy ordering.Benchmarks on the gas phase NH þ -water clusters pre-sented in Fig. 2(g) further support our conclusions. It isworth noting that recent experimental measurements alsosuggested that hydrogen bonds of solvated NH þ are lessPHYSICAL REVIEW LETTERS125,
FIG. 1. Hydration structures of NH þ in water obtainedfrom AIMD. (a) The N-O radial distribution function (solidline) and the corresponding running coordination number(dashed line). The running coordination number n NO ð R Þ ¼ πρ R R G NO ð r Þ r dr , where ρ is the density. The first minimum( r min ) of G NO ð r Þ is 3.64 Å, 3.38 Å, 3.55 Å for SCAN, PBE,PBE þ vdW, respectively. The coordination number (CN) is n ð r min Þ . The convergence of G NO ð r Þ is shown in Fig. S2(a).(b) The distribution of the coordination number of the firsthydration shell. (c) Calculated H N -O radial distribution func-tions. The first minimum of G H N O ð r Þ is 2.40 Å, 2.45 Å, 2.39 Å forSCAN, PBE, PBE þ vdW, respectively. (d) Average number ofNH þ -water hydrogen bonds (HBs) predicted by differentexchange correlation (XC) functionals. (e) The difference ofthe joint distribution in the first hydration shell between theSCAN and the PBE results as a function of r and β [defined in theinset of (d)]. Red color implies a larger value from SCAN. Thejoint distributions are shown in Fig. S3. (f) The N -water totaldistribution function G N ð r Þ . The blue solid line is from theneutron diffraction work of Hewish et al. [43]. The grey barsindicate the positions of the second maximum and the secondminimum of the experimental G N ð r Þ . PHYSICAL REVIEW LETTERS125, ¼ . ) and PBE þ vdW (CN ¼ . )are smaller. The distribution of CNs shown in Fig. 1(b)further demonstrates the differences of CN between SCANand the two other functionals. The most favorable co-ordination number of SCAN is six followed by seven andfive, while the CN distribution of PBE and PBE þ vdWpeak at four and five, respectively. As for the secondhydration shell, PBE and PBE þ vdW both have a clearsecond peak on G NO ð r Þ , whereas with SCAN the secondhydration shell mixes with the first shell. The flattening ofthe second peak predicted by SCAN has been observed inneutron diffraction and x-ray experiments [9,43].We now consider the number of hydrogen bonds formedby NH þ and the surrounding water molecules. Figure 1(c)plots the N N -O radial distribution functions from AIMD. Wedefine the standard hydrogen bonds with criteria reported forliquid water and aqueous solutions [44 – R NO < R NO C ; R H N O < R H N O C ; β < ° ; ð Þ where R NO C and R H N O C are the cut-off values for the N-O andH N -O distances obtained from the first minimum of thecorresponding radial distribution functions. The angle β isthe H N -N-O angle as shown in the inset of Fig. 1(d). We findthe number of NH þ -water hydrogen bonds predicted byPBE and PBE þ vdW are almost identical with an averagevalue of 3.50 and 3.48, respectively, suggesting a minorinfluence of vdW interactions on the NH þ -water hydrogenbonds. SCAN, however, predicts considerably weaker NH þ -water hydrogen bonds evidenced by a smaller average numberof 3.24. By ignoring the angular restraint in Eq. (1), i.e., onecan define the so-called generalized hydrogen bonds, includ-ing nonlinear, bifurcated hydrogen bonds that have β > °[47,48]. We find that SCAN predicts a larger number ofbifurcated NH þ -water hydrogen bonds.The hydration structure of NH þ in water can also bedemonstrated by the joint distributions of r and β , whichare plotted in Fig. S3 [41]. The peaks at β < ° and β ≈ ° represent the water molecules that form hydrogenbonds with NH þ and the increased intensity in between canbe attributed to the bifurcated water molecules that tend toform multiple nonlinear hydrogen bonds with the NH þ . Tohighlight the difference between different AIMD simula-tions, in Fig. 1(e) we plot the difference of the jointdistribution of SCAN and PBE. Because of the weakerNH þ -water hydrogen bonds and increased bifurcated watermolecules predicted by SCAN, the peaks at β < ° and β ≈ ° are reduced while the intensity in the middlebecomes stronger.Figure 1(f) further plots the N -water total distributionfunction [ G N ð r Þ ] obtained from our simulations and fromavailable neutron diffraction data. Experimental G N ð r Þ , asdefined in Eq. (2), is different from g NO ð r Þ in Fig. 1(a) andinvolves g NN and g NCl terms because experiments used . mol = L NH Cl water solutions [6,9]. G N ð r Þ ¼ . g NO ð r Þ þ . g NH ð r Þþ . g NCl ð r Þ þ . g NN ð r Þ − . : ð Þ Our simulated G N ð r Þ , however, contains only the first twoterms in Eq. (2). The difference in solution conditions isapparent from the difference between curves at large radii.Nevertheless, at small radii, the comparison indicates that theSCAN result is reasonable, in particular the position of thesecond maximum r max and the second minimum r min are inline with the experiment. Simulations with an inaccuratefunctional such as PBE do not correctly reproduce theseexperimental features. In addition, Figs. S4 and S5 [41]present vibrational spectra and density of states computedfrom our simulations, which may be compared in futureexperiments to confirm the proposed hydration structure.Overall, we find that the hydration structure predicted bySCAN is more consistent with experiments than thatobtained with the other two functionals, with the SCANresults characterized by weaker NH þ -water hydrogenbonds and an increased number of water molecules inthe first hydration shell. However, to confirm that theSCAN description of the hydration structure is indeedreliable and the agreement between theory and experimentis not merely fortuitous, it is necessary to examine theaccuracy of the underlying potential energy surface with thehelp of a higher level of theory. To this end we selectedfive distinct configurations with CN from four toeight to benchmark against DMC calculations, as shown inFigs. 2(a – e). Relative energies of these configurations arepresented in Fig. 2(f). According to DMC, the configura-tion with CN ¼ is the lowest energy state among theselected structures. The basin shape of the DMC predictedenergy curve is reproduced quite well by SCAN, indicatinga good description of the hydration structure and thehydrogen bonding network around NH þ . Indeed, theenergy ordering of the selected structures predicted byDMC and SCAN is qualitatively consistent with thedistribution of CNs shown in Fig. 1(b). Regarding theperformance of other DFT functionals, typical results areshown with PBE and PBE þ vdW. From this it can be seenthat they predict quite different energy ordering for thesame set of structures. Specifically, both PBE and PBE þ vdW tend to underestimate the stability of the CN ¼ configuration. Furthermore, SCAN correctly predicts thatthe configuration of CN ¼ is the lowest energy statewhile PBE and PBE þ vdW both overestimate the stabilityof the state of CN ¼ . In Fig. S6 [41], we show the resultscalculated using several other exchange correlation func-tionals and none of them predict a better energy ordering.Benchmarks on the gas phase NH þ -water clusters pre-sented in Fig. 2(g) further support our conclusions. It isworth noting that recent experimental measurements alsosuggested that hydrogen bonds of solvated NH þ are lessPHYSICAL REVIEW LETTERS125, þ -waterinteraction is indeed weaker than previously thought, andthe relatively weak interaction leads a large coordinationnumber of NH þ and an increased number of nonlinearhydrogen bonds. Theoretically, the relatively weak hydro-gen bonds predicted by SCAN can be understood by thelower polarizability of the water molecules around NH þ (Fig. S5); a suggestion that is consistent with previousstudies in water [18,19].Having established that the hydration structure of NH þ is characterized by weak and bifurcated hydrogen bonds,and having found that the SCAN functional can reproducethis structure, we now discuss the rotation of NH þ inwater. To aid the quantification of the rotation, we define avector in the body-fixed frame and track the rotation of thevector as a function of time. The angle θ is defined as theangle between the initial vector and the vector after agiven time. Figure 3(a) shows a typical evolution of θ in apicosecond window when the rotation of NH þ occurs,which is described by a big change of θ of 80 –
90° in0.3 – – W to W arethe six water molecules in the first hydration shell of NH þ at T ¼ . ps. Initially two bifurcated water molecules( W and W ) form nonlinear hydrogen bonds with NH þ ,hence H and H form bifurcated hydrogen bonds with W =W and W =W , respectively. After the simultaneousswitching of the two bifurcated hydrogen bonds, N-H forms a stable hydrogen bond with W and N-H pointstowards W . Because of the bifurcated hydrogen bonds,such a mechanism does not involve complete breaking ofhydrogen bonds, hence facilitating the rotation. In addi-tion, the rotation of NH þ involves the rotation of two N-Hbonds instead of one, thus it is much more efficient whenthere is a pair of bifurcated hydrogen bonds. A statisticalanalysis of rotation times and the extent of the rotationalangle jumps is provided in Fig. S8, confirming that largerotational jumps dominate the overall rotational process.In contrast, with PBE and its hydration shell containingfewer bifurcated water molecules, the rotation involves thebreaking of hydrogen bonds, and the process is sloweddown (Fig. S9). CN = 6 CN = 4 CN = 5CN = 7 CN = 8
Coordination Number -5051015 R e l a ti v e E n e r gy ( k ca l/ m o l ) SCANPBEPBE+vdWDMC
Relative energy of DMC (kcal/mol) R e l a ti v e e n e r gy ( k ca l/ m o l ) SCANPBEPBE+vdWDMC (a)(f) (g)(b)(d) (e)(c)
FIG. 2. Benchmarks of different DFT functionals against DMCcalculations. (a) – (e) Snapshots of the selected configurations.Atoms in bright colors represent the first hydration shell of NH þ in water (blue: N; red: O; white: H). (f) Relative energiescalculated using different methods at different configurationswith CN from four to eight. The relative energy is defined as theenergy difference to the energy of the configuration with CN ¼ .(g) Relative energies of gas phase NH þ clusters calculated usingdifferent methods and plotted against the DMC results. Therelative energy is defined as the energy difference to the energy ofthe configuration with the lowest one. The gas phase NH þ clusters are randomly selected from the AIMD trajectories, whichconsist of NH þ and eight neighboring water molecules. H H H H W WWWWWWWWWWWW W W WWW W WWW W W H HHHH H H W WWWWWWWWW W WWW W W W W (b) T = 17.35 ps
17 17.2 17.4 17.6 17.8 18
Time (ps) () (a) H H H H H H H H (c) T = 17.70 ps FIG. 3. Rotation mechanism of NH þ in water. (a) Rotationangle θ (the angle change with respect to T ¼ . ) as a functionof time. Calculated θ is averaged by the vector NH (cid:2) ! þ NH (cid:2) ! ,NH (cid:2) ! þ NH (cid:2) ! , and NH (cid:2) ! þ NH (cid:2) ! . (b) – (c) Snapshots of a typicalrotational process. The blue sphere is N, red spheres are O, andwhite spheres are H. H to H belongs NH þ . Water molecules W to W represent the first hydration shell at T ¼ . ps. Rotationis associated with bifurcated hydrogen bonds formed withH and H . PHYSICAL REVIEW LETTERS125,
Time (ps) () (a) H H H H H H H H (c) T = 17.70 ps FIG. 3. Rotation mechanism of NH þ in water. (a) Rotationangle θ (the angle change with respect to T ¼ . ) as a functionof time. Calculated θ is averaged by the vector NH (cid:2) ! þ NH (cid:2) ! ,NH (cid:2) ! þ NH (cid:2) ! , and NH (cid:2) ! þ NH (cid:2) ! . (b) – (c) Snapshots of a typicalrotational process. The blue sphere is N, red spheres are O, andwhite spheres are H. H to H belongs NH þ . Water molecules W to W represent the first hydration shell at T ¼ . ps. Rotationis associated with bifurcated hydrogen bonds formed withH and H . PHYSICAL REVIEW LETTERS125, D R [7,12], which can bederived from the mean-square angular displacementsaccording to h θ ð t Þ i ¼ D R t: ð Þ The computed h θ ð t Þ i as a function of time t are shownin Fig. S10(a) [41]. From the linear part of h θ ð t Þ i curve,we can estimate the rotational diffusion constant D R which is a quarter of the slope. The estimate for D R usingSCAN is . (cid:2) . rad ps − at 363 K (the conver-gence of D R was shown in Fig. S11) while the values fromPBE ( . (cid:2) . ) and PBE þ vdW ( . (cid:2) . )simulations are significantly smaller. In Fig. 4(a) weplot D R as a function of temperature along with theexperimental values derived from the NMR measurements[7,10 – D R follows a good linear relationin the temperature regime studied, and a linear fit ishighlighted with the purple dashed line for NH þ inH O. Through the experimental data of ND þ in D O we draw a green dashed line parallel to the purple line, whichsuggests that our simulated D R with the SCAN functionalfalls in line with experiments. We note that to explicitlydiscuss the difference between the hydrogenated anddeuterated systems, nuclear quantum effects should beincluded using, e.g., path integral molecular dynamics inthe future [49,50]. Nevertheless, from the experimentalmeasurements it is evident that nuclear quantum effectshave a minor impact, of approximately . rad ps − , onthe rotational diffusion constant [11].To quantitatively demonstrate the correlation betweenhydration structure and rotational dynamics, we further plot D R at 363 K as a function of the average number ofbifurcated hydrogen bonds in Fig. 4(b). We find D R growsas the number of bifurcated hydrogen bonds increases,following a linear relation. Similar results are obtained atother temperatures, which are presented in Fig. S10 [41].The linear relation established offers strong evidence forour finding that the rotation involves the simultaneousswitching of bifurcated hydrogen bonds.To conclude, we have performed AIMD simulationsand DMC calculations of the NH þ aqueous solution. Wehave clarified the hydration structure of NH þ in water,finding that it is characterized by a coordination number ofapproximately six, with two pairs of bifurcated hydrogenbonds. Consequently, we find that this bonding natureleads to fast rotation of NH þ in water. A clear predictionfrom this study is that the rotation of NH þ may besignificantly suppressed in a different solvent where bifur-cated hydrogen bonds do not form or when it is confined toan extent that the hydration structure is partially broken. Inthe future, such studies are desirable to verify our con-clusions and confirm the diffusion mechanism suggestedhere. Last but not the least, this study brings attention to theimportance of achieving accurate electronic structures ofother aqueous and hydrogen bonded systems, and furtherimprovements in computational methods will graduallybring us towards an exact description of such complexprocesses.This work was supported by the National KeyR&D Program of China under Grant No. 2016YFA030091, the National Natural Science Foundation ofChina under Grants No. 11974024, No. 11525520, andNo. 11935002. A. Z. ’ s work is sponsored by the AirForce Office of Scientific Research, Air Force MaterialCommand, U.S. Air Force, under Grant No. FA9550-19-1-7007. X. W. ’ s work was conducted within theComputational Chemical Center: Chemistry in Solutionand at Interfaces funded by the DoE under Award No.DE-SC0019394. We are grateful for computationalresources provided by TianHe-1A and TianHe II super-computers, the High Performance Computing Platform ofPeking University, and the Platform for Data DrivenComputational Materials Discovery of the SongshanLake Materials Lab. Number of bifurcated H bonds D R (r a d p s - ) (b) SCAN (363 K)PBE (363 K)PBE+vdW (363 K)
260 280 300 320 340 360 380
Temperature (K) D R (r a d p s - ) (a) SCANPBEPBE+vdW EXP (NH in H O)EXP (ND in D O) FIG. 4. Diffusion of NH þ in water. (a) Rotational diffusionconstant D R from experiments and our simulations at differenttemperatures [11]. The purple dashed line is a linear fit of theexperimental rotational diffusion constant as a function oftemperature for H, and the green dashed line is a parallel lineto the purple dashed line crossing the only available experimentaldata for D as a guide to the eye. (b) D R calculated from oursimulations at 363 K as a function of the average number ofbifurcated hydrogen bonds. The gray line is a linear fit. PHYSICAL REVIEW LETTERS125,
Temperature (K) D R (r a d p s - ) (a) SCANPBEPBE+vdW EXP (NH in H O)EXP (ND in D O) FIG. 4. Diffusion of NH þ in water. (a) Rotational diffusionconstant D R from experiments and our simulations at differenttemperatures [11]. The purple dashed line is a linear fit of theexperimental rotational diffusion constant as a function oftemperature for H, and the green dashed line is a parallel lineto the purple dashed line crossing the only available experimentaldata for D as a guide to the eye. (b) D R calculated from oursimulations at 363 K as a function of the average number ofbifurcated hydrogen bonds. The gray line is a linear fit. PHYSICAL REVIEW LETTERS125, [email protected] † [email protected] ‡ [email protected][1] M. Smyth and J. Kohanoff, Phys. Rev. Lett. , 238108(2011).[2] G. C. Sosso, J. Chen, S. J. Cox, M. Fitzner, P. Pedevilla,A. Zen, and A. Michaelides, Chem. Rev. , 7078(2016).[3] J. Abraham, K. S. Vasu, C. D. Williams, K. Gopinadhan, Y.Su, C. T. Cherian, J. Dix, E. Prestat, S. J. Haigh, I. V.Grigorieva, P. Carbone, A. K. Geim, and R. R. Nair,Nat. Nanotechnol. , 546 (2017).[4] N. Agmon, H. J. Bakker, R. K. Campen, R. H. Henchman, P.Pohl, S. Roke, M. Thämer, and A. Hassanali, Chem. Rev. , 7642 (2016).[5] D. Pan and G. Galli, Nat. Commun. , 421 (2020).[6] F. Brug´e, M. Bernasconi, and M. Parrinello, J. Am. Chem.Soc. , 10883 (1999).[7] T.-M. Chang and L. X. Dang, J. Chem. Phys. , 8813(2003).[8] M. Ekimova, W. Quevedo, Ł . Szyc, M. Iannuzzi, P. Wernet,M. Odelius, and E. T. J. Nibbering, J. Am. Chem. Soc. ,12773 (2017).[9] G. Pálinkás, T. Radnai, G. I. Szász, and K. Heinzinger,J. Chem. Phys. , 3522 (1981).[10] C. L. Perrin and R. K. Gipe, Science , 1393 (1987).[11] C. L. Perrin and R. K. Gipe, J. Am. Chem. Soc. , 1088(1986).[12] P. Intharathep, A. Tongraar, and K. Sagarik, J. Comput.Chem. , 1329 (2005).[13] M. Chen, L. Zheng, B. Santra, H.-Y. Ko, R. A. DiStasio Jr,M. L. Klein, R. Car, and X. Wu, Nat. Chem. , 413 (2018).[14] Y. Bouchafra, A. Shee, F. R´eal, V. Vallet, and A. S. P.Gomes, Phys. Rev. Lett. , 266001 (2018).[15] W. Kohn, Rev. Mod. Phys. , 1253 (1999).[16] R. Car and M. Parrinello, Phys. Rev. Lett. , 2471 (1985).[17] M. J. Gillan, D. Alf`e, and A. Michaelides, J. Chem. Phys. , 130901 (2016).[18] M. Chen, H.-Y. Ko, R. C. Remsing, M. F. C. Andrade, B.Santra, Z. Sun, A. Selloni, R. Car, M. L. Klein, J. P. Perdew,and X. Wu, Proc. Natl. Acad. Sci. U.S.A. , 10846 (2017).[19] Y. Yao and Y. Kanai, J. Chem. Theory Comput. , 884(2018).[20] B. Santra, J. Klime š , D. Alf`e, A. Tkatchenko, B. Slater, A.Michaelides, R. Car, and M. Scheffler, Phys. Rev. Lett. ,185701 (2011).[21] G. Mazzola and S. Sorella, Phys. Rev. Lett. , 105701(2015).[22] A. Zen, J. G. Brandenburg, J. Klime š , A. Tkatchenko, D.Alf`e, and A. Michaelides, Proc. Natl. Acad. Sci. U.S.A. ,1724 (2018).[23] J. Chen, A. Zen, J. G. Brandenburg, D. Alf`e, andA. Michaelides, Phys. Rev. B , 220102 (2016).[24] K. G. Reeves, Y. Yao, and Y. Kanai, J. Chem. Phys. ,124705 (2016).[25] P. Giannozzi et al. , J. Phys. Condens. Matter , 395502(2009). [26] P. Giannozzi et al. , J. Phys. Condens. Matter , 465901(2017).[27] J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. , 036402 (2015).[28] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. , 3865 (1996).[29] A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. ,073005 (2009).[30] D. R. Hamann, M. Schlüter, and C. Chiang, Phys. Rev. Lett. , 1494 (1979).[31] D. Vanderbilt, Phys. Rev. B , 8412 (1985).[32] G. J. Martyna, M. L. Klein, and M. Tuckerman, J. Chem.Phys. , 2635 (1992).[33] R. J. Needs, M. D. Towler, N. D. Drummond, and P. L. Ríos,J. Phys. Condens. Matter , 023201 (2010).[34] A. Zen, S. Sorella, M. J. Gillan, A. Michaelides, and D.Alf`e, Phys. Rev. B , 241118 (2016).[35] J. R. Trail and R. J. Needs, J. Chem. Phys. , 204107(2017).[36] L. M. Fraser, W. M. C. Foulkes, G. Rajagopal, R. J. Needs,S. D. Kenny, and A. J. Williamson, Phys. Rev. B , 1814(1996).[37] A. J. Williamson, G. Rajagopal, R. J. Needs, L. M. Fraser,W. M. C. Foulkes, Y. Wang, and M.-Y. Chou, Phys. Rev. B , R4851 (1997).[38] P. R. C. Kent, R. Q. Hood, A. J. Williamson, R. J. Needs, W.M. C. Foulkes, and G. Rajagopal, Phys. Rev. B , 1917(1999).[39] J. G. Brandenburg, A. Zen, M. Fitzner, B. Ramberger,G. Kresse, T. Tsatsoulis, A. Grüneis, A. Michaelides, andD. Alf`e, J. Phys. Chem. Lett. , 358 (2019).[40] Y. S. Al-Hamdani, M. Rossi, D. Alf`e, T. Tsatsoulis, B.Ramberger, J. G. Brandenburg, A. Zen, G. Kresse, A.Grüneis, A. Tkatchenko, and A. Michaelides, J. Chem.Phys. , 044710 (2017).[41] See the Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevLett.125.106001 for moredetails of computational settings, convergence tests, addi-tional analyses, and extended data.[42] D. Alf`e and M. J. Gillan, Phys. Rev. B , 161101(2004).[43] N. Hewish and G. Neilson, Chem. Phys. Lett. , 425(1981).[44] A. Luzar and D. Chandler, J. Chem. Phys. , 8160(1993).[45] A. Luzar and D. Chandler, Phys. Rev. Lett. , 928(1996).[46] A. Chandra, Phys. Rev. Lett. , 768 (2000).[47] S. Chowdhuri and A. Chandra, Chem. Phys. Lett. , 79(2003).[48] S. K. Pattanayak and S. Chowdhuri, J. Mol. Liq. , 98(2013).[49] M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie,A. Michaelides, M. A. Morales, and T. E. Markland,Chem. Rev. , 7529 (2016).[50] W. Fang, J. Chen, Y. Feng, X.-Z. Li, and A. Michaelides,Int. Rev. Phys. Chem. , 35 (2019). PHYSICAL REVIEW LETTERS125,