|Home | About | Journals | Submit | Contact Us | Français|
Although it has long been recognized that multiple water molecules strongly associate with an extra proton in bulk water, some models and conceptual frameworks continue to utilize the classical hydronium ion (H3O+) as a fundamental building block. In this work, the nature of the hydronium ion in aqueous systems is examined using an ab initio energy decomposition analysis (EDA) that evaluates both the magnitude of and energetic stabilization due to charge transfer among H3O+ and the surrounding water molecules. The EDA is performed on structures extracted from dynamical bulk-phase simulations, and used to determine how frequently the pure hydronium ion, where the excess charge is primarily localized on H3O+, occurs under dynamic conditions. The answer is essentially never. The energetic stabilization of H3O+ due to charge delocalization to neighboring water molecules is found to be much larger (16 to 49 kcal/mol) than for other ions (even Li+) and to constitute a substantial portion (20% to 52%) of the complex's binding energy. The charge defect is also shown to have intrinsic dynamical asymmetry and to display dynamical signatures that can be related to features appearing in IR spectra.
Solvated (i.e., hydrated) protons play a central role in many areas of chemistry, biology, physiology and materials science. They act as essential reactants or products in acid-base chemistry, as the modulators of pH, as activators of protein folding or function via residue ionization, and as the source of energy in biological and inorganic energy transduction. Accordingly, the fundamental nature of the solvated proton has been studied extensively over the last century.1-4 It has long been accepted that when a bare proton is introduced to an aqueous environment it does not remain an independent charged nucleus (H+), but immediately associates with a water molecule to form H3O+. This entity, called a hydronium or oxonium ion, is traditionally defined as an oxygen atom strongly bound to three hydrogen atoms, possessing a single net positive charge and, in bulk water, strong hydrogen bonds to surrounding water molecules. Although this may be a sufficient definition for certain applications, a long history of theoretical, computational and experimental efforts have characterized the solvated proton in greater detail, and have revealed much more complexity in the mechanisms of proton solvation and transport. It is known, for example, that the proton diffuses via both vehicular and chemical diffusion (i.e., via the Grotthuss mechanism),1-4 that proton transfer from one oxygen to another is a complex process involving correlated hydrogen bond rearrangements in the 1st and 2nd solvation shells,5,6 that the excess charge often has structural and charge localization characteristics that have been described in terms of dynamically interchanging Zundel (H5O2+) and Eigen (H9O4+) cations evolving through intermediate structures,5-19 and that protons in confined environments, such as in ion channels or at the liquid-vapor interface,3 display qualitatively different behavior than in the bulk. Despite many such advances, the bare hydronium ion H3O+ still appears as a fundamental building block in a number of computer simulations20-29 and in most chemistry and biochemistry textbooks.
In this paper, the nature of the hydronium ion in bulk water is further examined by using established ab initio methods to probe both the magnitude and the energetic stabilization of charge transfer (i.e., the degree of delocalization or molecular covalency) in the hydrated proton. It is well accepted that there is always some degree of charge transfer involved in the hydrated proton complex, and that there is, by definition, charge transfer during a proton transfer event. The questions asked here are, how strong and distributed is the electronic charge delocalization? What are the magnitudes of both the charge transfer and the energetic stabilization due to charge transfer? Furthermore, can characterizing charge delocalization be used to better understand the nature of the hydrated proton?
The charge transfer properties are quantified herein with a recently developed ab initio energy decomposition analysis (EDA)30,31 that uses absolutely localized molecular orbitals (ALMO) to decompose interaction energies (e.g., for a collection of water molecules and one excess proton) into physically meaningful components, specifically geometric distortion, electrostatic, polarization, and charge transfer contributions. ALMO EDA is performed on many (H2O)nH+ geometries extracted from bulk-phase simulations containing hundreds of water molecules and one excess proton. The structures of the solvated species are then characterized based on the strength and nature of the interaction energies between the central cation (H3O+) and the surrounding water molecules. Within this characterization, the pure hydronium ion is defined as an H3O+ ion with a very localized charge distribution and with strong hydrogen bonds to, but minimal charge transfer from, the surrounding water molecules. The Eigen cation(H9O4+), is identified as an H3O+ ion having significant (further defined below) electronic charge delocalization among these water molecules, while the Zundel cation is labeled as any H5O2+ having dominant charge delocalization primarily between two water molecules.
Using these definitions, the pure localized hydronium ion is found to make a negligible contribution to the bulk-phase distribution of protonated structures. In agreement with previous studies,2-6,9-12,15-18,32 the most commonly occurring protonated complexes are best viewed as a continuum of structures between the “limiting” the Eigen (H9O4+) and Zundel (H5O2+) cations.
Comparing structures prevalent in the bulk to the minimum-energy structures of small clusters, it is found that the anisotropic nature of local interactions in aqueous media significantly destabilizes symmetric structures, resulting in very high occurrence of asymmetric structures, as also suggested in previous theoretical work.5,11,12,17,19,33 Therefore, the often used notion that the symmetric Eigen and Zundel cations, which are the minimum-energy configurations in gas-phase clusters, are also the “limiting” structures in bulk solvent is misleading. More accurate descriptions of the most probable structures are: 1) an asymmetric Eigen cation5,17 with two to four strong interactions to the 2nd solvation shell (relative to the central H3O+), and 2) an asymmetric Zundel cation (with two to four strong interactions to the 2nd-solvation shell) that forms transiently during a successful proton transfer (PT) event or repeatedly in what appear to be unsuccessful PT events. For the purposes of consistency, the water that shares the charge most strongly with H3O+ in the Zundel cation is denoted the 1st solvation shell in this manuscript, while the four waters interacting most strongly with H5O2+ are denoted the 2nd solvation shell. It is shown that the existence of the symmetric Zundel cation (i.e., one that shares the excess H equally and has equal charge transfer from all four 2nd-solvation shell waters) is rare. Much more common is the formation of an asymmetric Zundel cation5 with stronger 2nd shell interactions to one of the waters in the H5O2+ unit, resulting in a structure that resembles an asymmetric Eigen structure.17 As a result of this finding, this work supports previous suggestions that making a distinction between the Eigen and Zundel cations is difficult, if not impossible.10-12,16-19,34
In addition to monitoring structures characteristic of hydronium, Eigen, and Zundel species, we examine the extent of charge transfer from the 2nd and 3rd solvation shells around the central cation. We find that the charge transfer from the 2nd shell to the 1st is usually larger when the proton is present than it is in pure water in the absence of a proton. More importantly, we find that charge transfer from one solvation shell to the next is minimally altered by the presence or absence of other solvation shells. These results support the concept that charge transfer is a local phenomenon predominately dependent on the local atomic configuration, in contrast to the notion that charge transfer involves many-body interactions that cause the 2nd solvation shell to transfer more charge to the central H3O+ than the 1st shell.34 Moreover, our results demonstrate that placing a protonated water cluster in a polarizing medium does not cause significant charge localization, as would occur with a bare proton.
Finally, by monitoring the time evolution of the charge transfer during dynamical bulk-phase simulations, we show that periodic temporal signatures of the solvated proton's motion are reflected in charge transfer energies. The signature for the water molecule donating the most charge to the central cation during what we term Eigen dynamics is well correlated with the Eigen O-H stretching motion, which has a period of 12 fs and produces an experimentally established infrared (IR) spectroscopic peak near 2800 cm−1. The equivalent signature during what we define as Zundel dynamics produces peaks near 1800 and 2200 cm−1, which are both at higher frequencies than that of the bare Zundel (H5O2+) O-H stretching motion in gas-phase (~1000 cm−1),35 but fall in the range of experimentally observed IR peaks for the bulk phase.36 Perhaps the most interesting feature in the periodic signatures is a much slower feature found for the time evolution in the magnitude of the total charge transfer. This likely corresponds to either a symmetric bend along the principle axis through the central oxygen atom or a stretching motion of the 1st solvation shell.
Although other authors have discussed the magnitude of charge transfer in protonated clusters from population analysis or alternative approaches to EDA,34,37,38 this is, to the best of our knowledge, the first work to quantify both the magnitude of and the energetic stabilization due to charge transfer for the excess proton using structures representative of bulk-phase aqueous medium. The results clearly show that there is substantial charge transfer (i.e., delocalization of the electron density) among four to six water molecules, making the Eigen and Zundel cations significantly covalent in nature. They also show that the pure hydronium ion, with the vast majority of the excess charge localized on the central H3O+, is essentially never observed.
In the following section, the methods used in the bulk phase simulations, geometry optimizations, and energy decomposition analysis are described. In section three, the ALMO EDA results for the geometry-optimized isolated (i.e., in the absence of other water molecules) Eigen and Zundel cations are presented. Section four then turns to the bulk-phase results for cation clusters (H2O)nH+ (n = 4,6,10) whose geometries were extracted from ‘snap shots’ taken during dynamic bulk-phase simulations. Both ensemble-averaged and dynamical properties are considered. Closing remarks and implications are presented in section five.
To provide a baseline for the bulk-phase analysis, ALMO EDA was first performed on the geometry-optimized isolated Eigen cation (H2O)4+ and the Zundel cation with its 2nd solvation shell (H2O)6H+. These structures, shown in Figure 1, were optimized at the MP2/aug-cc-pVTZ39 and MP2/cc-pVTZ40 levels, respectively.41 All geometry optimizations were performed with the Gaussian03 program,42 and the Berney optimization algorithm with the tight convergence criteria. For comparison, the structures were also optimized with DFT using the hybrid B3LYP functional43 and the aug-cc-pVTZ basis set. The ALMO EDA is, thus far, limited to single-determinant wavefunctions, meaning Hartree Fock (HF) and density functional theory (DFT). Because the DFT and MP2 structures differ (averaged over all lengths and angles) by only 0.009 Å, 0.16°, and 0.72° in bonds, angles, and torsions, respectively, the B3LYP functional was used as a correlated method in the ALMO EDA analysis of protonated water clusters.44,45
Energy decomposition approaches are used to characterize intermolecular interactions by splitting the total interaction energy (i.e., binding energy) into physically meaningful components.30,31,46-48 In the present study, those components are geometric distortion EGD, frozen density (i.e. electrostatic) EFRZ, polarization EPOL, and charge transfer ECT.
In 2007 Head-Gordon and coworkers developed a new approach to EDA, which naturally and completely separates charge transfer and self consistently optimized polarization effects.30 Their method relies on absolutely localized orbitals (ALMOs), which had been previously introduced and used to increase the efficiency of SCF energy convergence,49 Each ALMO is expanded in terms of the atomic orbitals of only a single molecule in the system, as opposed to being delocalized over the entire system (e.g., as with canonical HF MOs). The resulting nonorthogonal (from one molecule to another) ALMOs require using a special SCF procedure to produce properly antisymmetrized many-electron wavefunctions out of such nonorthogonal components.49 Nevertheless, the resulting intermediate wavefunction is a variationally optimized electronic state (i.e., self-consistently polarized) that has, by design, prohibited charge transfer between molecules. This then allows one to keep separate the polarization and charge transfer contributions to the total interaction energy. ALMO EDA has been implemented in a beta version of the Q-Chem software package, which was used in all of the EDA calculations presented herein.50
The first energy component, EGD, is the energetic cost due to distorting the isolated molecules from their individual optimized structures to the structures they adopt in the molecular complex. This term can be calculated with standard SCF optimization and single-point techniques at any level of theory. The second term, EFRZ, is the energy due to bringing the separated, but distorted, monomers together to the complex geometry without further optimizing (i.e., relaxing) the MOs on either monomer,
where ΨX is the fully SCF optimized wavefunction for each isolated monomer × fixed in the geometry it has in the complex (i.e., with the monomer geometries and MOs ‘frozen’), and Ψ0 is the properly antisymmetrized wavefunction for the complex composed of the unrelaxed MOs of the monomers. The polarization term, EPOL, is the energy lowering due to relaxing each monomer's ALMOs in the field of the rest of the complex,
where ΨALMO is the intermediate determinant composed of SCF optimized ALMOs. Finally, the charge transfer, ECT, is the difference between the ΨALMO intermediate state's energy and that of the fully SCF optimized Ψ composed of delocalized MOs plus the basis set superposition error (BSSE),
Charge transfer and BSSE both arise from the delocalization of monomer MOs, but BSSE can be corrected for in the ALMO EDA method using the conventional counterpoise method. Moreover, BSSE decreases much more quickly with basis set enhancement than charge transfer, as will be shown below. A detailed discussion of the role of BSSE in ALMO EDA can be found elsewhere.30 The ALMO EDA approach can also separate the charge transfer into forward and reverse transfer,31 but this is not used in the present study as charge transfer is almost completely in one direction for protonated water clusters. Although the amount of physical charge transferred (measured in mill-electrons, mē) is reported, the energetic consequences of charge transfer are the focus in this work.
In all of the ALMO EDA analysis discussed herein, the clusters have been defined, unless otherwise mentioned, in terms of one monomer for the central cation (H3O+) and one monomer for each water molecule. In this way, the interactions between the central cation and each individual water molecule can be evaluated. However, only a single total value for all the surrounding waters in a given complex is presented for EGD, EFRZ, EPOL and EBIND. The charge transfer component, on the other hand, is broken down into each individual monomer-to-monomer contribution (excluding higher order terms as explained in reference 30).
Analyzing an entire bulk-phase system (i.e., including hundreds of water molecules) by applying ALMO EDA to all of the water molecules would be computationally intractable. Instead, EDA was performed on over 900 Eigen and Zundel conformations (as shown in Fig. 1) that were extracted from the bulk-phase simulations. The dynamical evolution of the bulk system was described using three approaches: multistate empirical valence bond molecular dynamics (3rd generation, MS-EVB3),51 Car-Parrinello molecular dynamics (CPMD)52 using the BLYP functional,53 and CPMD with the HCTH/120 functional.54 These simulations were used to generate a large ensemble of structures that likely arise in the bulk phase. Then, at least 150 structures for both Eigen and Zundel cations, from each type of simulation, were analyzed with ALMO EDA.
In the MS-EVB model3,15-18 used to generate some of the dynamical trajectories, proton solvation and transport is incorporated by describing the system as a linear combination of empirically motivated valence bond “states” (i.e., particular bonding topologies). A reactive potential energy surface is first defined as the lowest energy solution of a Hamiltonian matrix that contains the potential energy for each of the EVB states in the diagonal elements, and the coupling between states in the off-diagonal elements. Then the nuclei of the system are evolved classically along this reactive potential energy surface, resulting in the formation and breaking of chemical bonds and delocalization of the excess protonic charge defect, both of which are prerequisites for Grotthuss-like proton shuttling. More details on the MS-EVB method can be found elsewhere.3
At each time step in the MS-EVB simulations, the VB amplitudes, ci2, could, as in past work, be used to describe the structure of the protonated complex. The largest amplitude, cmax2, identifies the central cation. As described in reference 3, the Eigen cation has been defined as a structure with a cmax2 ≈ 0.6, meaning that ~ 60% of the excess charge is localized on the central cation and ~12% is distributed on each of the three surrounding water molecules. In contrast, the Zundel cation defined in MS-EVB has nearly equal charge sharing between two water molecules such that c12 ≈ c22 ≈ 0.45. The largest amplitude for a pure hydronium ion would have to be quite large (cmax2 ≥ 0.85). In this study, we also use geometric factors (described below) to distinguish the Eigen and Zundel cations, but we note that these identifications are consistent with the MS-EVB cmax2 definitions. An advantage of the geometry-based definitions is that they can also be used in the CPMD studies where ci2 values are not computed.
In the MS-EVB3 simulations, a cubic box (Lbox = 18.6 Å) holding 216 SPC/Fw55 water molecules and one excess proton was first equilibrated for 600 ps using the NVT ensemble (T = 298 K) to establish the proper density. The system was then equilibrated in the NVE ensemble for 50ps before starting production runs of 1 ns using a time step of 1 fs. Three such simulations provided more than enough candidate structures to use in the ab initio ALMO EDA described herein.
The CPMD simulations were performed with a plane-wave-basis DFT (cutoff = 80 Ry) in version 3.11 of the publicly available code.52,56 Troullier-Martins type pseudopotentials,57 were used to describe the ion-electron interactions, while the electronic interactions were described by either the HCTH/120 or BLYP exchange correlation functional. A fictitious electronic mass μ of 340 au was used, following the criteria of Schwegler et al.58 For both DFT functionals, the following procedure was applied. First, an excess proton was added to a pre-equilibrated system of 128 water molecules in a cubic box (Lbox = 15.6 Å). The pre-equilibration was carried out in classical MD using empirical force fields that were developed by force-matching CPMD simulations of liquid water using either the HCTH/120 or BLYP functionals.59 The system was then equilibrated under CPMD dynamical evolution for 8 ps before performing production runs of 100 ps in the NVE ensemble using a time step of 0.073 fs. The average temperature in the BLYP and HCTH/120 CPMD simulations were 301.1 and 295.8 K, respectively.
Once the MS-EVB or CPMD method had been used to obtain dynamical trajectories characteristic of the bulk-phase, ALMO EDA was preformed on over 150 structures from each of the three simulations (one MS-EVB and two CPMD). These structures were selected to maximize ensemble sampling. For example, over 150 structures were taken at times randomly distributed over the 192 ps CPMD/HCTH simulation. To compare how the EDA results varied from one simulation to the other, each ensemble was first analyzed at the HF EDA level. The CPMD simulation ensembles were also analyzed at the DFT EDA level using the same functional as that used in the CPMD simulation. The average energy components are later reported for each ensemble (i.e., MS-EVB, CPMD/BLYP, or CPMD/HCTH) and EDA method (HF or DFT).
In order to see if the presence of additional solvation beyond the clusters (i.e., those extracted from bulk-phase simulations and then studied herein) would qualitatively change the ALMO EDA charge transfer predictions, an additional approach to estimating the magnitude of charge transfer was utilized. The Natural Population Analysis (NPA)60 method was used to calculate the formal charges for each monomer in a molecular complex both with and without a polarized continuum model (PCM)61 representation of the surrounding solvation. The magnitude of charge transfer was then inferred, as is commonly done, from the changes in formal charges.
The 900 structures extracted from the MS-EVB and CPMD dynamical simulations were first categorized as Eigen-like or Zundel-like based on the value of the following asymmetry coordinate,
Referring to Figure 1, RO0-H* is the O-H bond length in the central cation (O0) and RH*--O1 is the distance between that hydrogen and one of the 1st-shell oxygen atoms (O1x,O1y, or O1z) in the Eigen cation. By definition, O1x has the smallest δ value, O1y the second smallest, and O1z the largest (δO1x < δO1y < δO1z).62 Thus, δ will refer to δO1x throughout the rest of the manuscript. Note that a perfectly symmetric Zundel cation would have δ = 0, whereas the perfectly symmetric Eigen cation (as in Fig. 1) would have δ = 0.53. In this study, where such perfectly symmetric geometries seldom occur in the dynamical simulations, the Zundel-like and Eigen-like structures are defined as those with δ ≤ 0.1 and δ ≥ 0.2 Å, respectively. It should be emphasized that there is no unambiguous cut-off between an Eigen–like and Zundel-like δ value. In fact, the δ values obtained from the dynamical simulations populate a left-skewed Gaussian distribution, shown in Figure 3 that spans δ = 0 to δ = 0.6 and peaks at δ ≈ 0.35 Å. Nevertheless, these definitions have been used in previous CPMD studies5,9,12,63 and were recently supported by a time-dependent method for identifying Zundel and Eigen cations.5 In the time-dependent approach, the identities of the central oxygen (O0) and the 1st-shell water molecule with the strongest interaction to the central cation (O1x) are tracked. The O0 and O1x water molecules are referred to as the ‘special pair’. In the same study, the dynamic behavior of the Eigen cation was described as the ‘special pair dance’ because of the switching (~ every 40 fs) of the ‘special pair’ partner (O1x) between the three water molecules in the Eigen's 1st solvation shell. Therefore, distorted (asymmetric) Eigen cations persist for substantial periods (~1−2 ps) during which O0 does not change identities, but O1x alternates among three different oxygen atom identities. Zundel cations, on the other hand, display dynamical interchange of the O0 and O1x identities, within the ‘special pair’, as the central hydrogen atom oscillates between the two equally sharing water molecules. The pure hydronium ion is difficult to identify with geometric or dynamic properties. Instead, given its dependence on the high degree of charge localization, the identification of the pure hydronium ion is herein based on the ALMO EDA results.
The central water (O0) is identified as that with the closest three oxygen atoms, which is generally consistent with finding the structure with the largest amplitude (ci2) in the MS-EVB3 simulations (~95% of the time). The 2nd solvation shell is defined as the six hydrogen bond-accepting water molecules that have oxygen atoms closest to the 1st-shell hydrogen atoms for the Eigen cation, and as the four hydrogen bond accepting water molecules closest to O0 and O1x for the Zundel cation. The three hydrogen bond-donating water molecules in the 2nd solvation shell of the Eigen cation are identified as those with hydrogen atoms closest to the 1st-shell oxygen atoms.
The Eigen and Zundel cations involve ten and six water molecules, respectively, when including their 2nd solvation shells. However, characterizing the dynamic nature (i.e., changes in magnitude and location) of charge transfer in a continuous, and therefore meaningful, manner requires analyzing clusters with the same number of water molecules. Therefore, clusters containing six water molecules were extracted from the MS-EVB and CPMD simulations over continuous time periods and analyzed with ALMO EDA. Although this approach neglects four waters in the 2nd solvation shell of the Eigen cation, it is computationally feasible whereas ALMO EDA of ten-water clusters is quite expensive. In order to verify that leaving these four water molecules out of calculations does not alter the qualitative results presented below, one time series of ten-water molecule clusters was also analyzed.
The ALMO EDA results for the geometry optimized Eigen and Zundel cations are presented in Tables 1 and and2,2, respectively. There are four points worth noting before discussing the implications of the ALMO EDA results. First, each of the energetic components appears to converge as the basis sets becomes locally complete. The largest changes with increasing basis set size are an increase (more favorable) in the polarization, EPOL, and decrease (less favorable) in the charge transfer, ECT, terms. This trend is inherent to the ALMO EDA method and is caused by increasing flexibility for electronic polarization to account for the total molecular interaction with increasing molecular basis size. This trend is limited, however, by the fact that the problem becomes ill posed and the method no longer valid when the monomer basis sets become large enough to be linearly dependent. In the present case, the quadruple-zeta basis set was found to be linearly dependent when Cartesian d orbitals were employed, indicating that the values presented in Table 1 are basis set saturated. The decrease in basis set superposition error (BSSE) as the basis set size grows additionally supports basis set saturation.
Second, in agreement with other work,30 the charge transfer energies (see Tables 1 and and2)2) from the DFT EDA are larger than those from the HF analysis. Given that DFT is known to underestimate the HOMO-LUMO gap and thus to often overestimate charge transfer, the DFT values can be considered upper bound to the values that would be obtained with ALMO EDA including more accurate intermolecular electron correlation. Third, the geometric distortion energy is consistently a small, unfavorable contribution to the total complex binding energy. When calculated for the geometry-optimized structures at the MP2/aug-cc-pvtz level, EGD is 2.32 kcal/mol for the Eigen cation and 2.75 kcal/mol for the Zundel cation. Given that this contribution has a small influence on the other contributing energies, it will not be calculated for the bulk-phase solvation structures. Nor will it be further discussed, but assumed to add ~2−3 kcal/mol to each structure's total binding energy.
The ALMO EDA results for the Eigen cation (see Table 1) demonstrate two remarkable features. First, the charge transfer component is a significant percentage of the total binding energy (in both the HF and DFT analysis). It constitutes 20% to 33% of the total intermolecular binding energy (depending on basis set and electronic structure method), indicating that the hydronium-water interactions are largely covalent in nature. Second, the charge transfer energies are large (i.e., from 14 to 22 kcal/mol). This means it would take 14 to 22 kcal/mol to localize the charge defect onto the central cation (i.e., to form a classically localized hydronium ion). This large energetic cost is why we essentially never observe the pure hydronium ion in bulk-phase simulations.
This magnitude of charge transfer energy in the Eigen cation is much larger than it is for the water dimer (~1 kcal/mol at the average bulk phase hydrogen bond length of 1.9 Å) or for other cation–water interactions. In fact, the charge transfer energy from water to most other monatomic cations, with the exception of Li+, is negligible.30 For Li+ the charge transfer energy to a single water molecule is 2.6 kcal/mol, which is only 8% of a total binding energy dominated by the electrostatic interactions, EFRZ, as expected for an ionic interaction. As more water molecules bind to Li+, their respective ECT values become negligibly small.30 These comparisons indicate that the bare proton is unique in its ability to delocalize charge (facilitate charge transfer) in the solvent environment. They also show that the Eigen cation involves far too much charge transfer (what is often considered to be covalency) to be considered a pure hydronium ion surrounded by three water molecules. It should be viewed as H3O+(H2O)3 with the charge highly delocalized.
The ALMO EDA results for the Zundel cation (see Table 2) show two clear similarities to those for the Eigen cation: 1) DFT predicts larger ECT than HF in the EDA, and 2) the Zundel cation is very different from a solvated, pure hydronium ion. There are, however, several differences that could have been anticipated from the close proximity of O0 and O1x of the ‘special pair’ water molecules in the Zundel structure (see discussion in section 2.4). The electrostatic interactions are only half as favorable due increased overlap, while ECT and EPOL are twice as favorable. The relative contribution of ECT to EBIND increases from 19% to 23%. Finally, the symmetry of the Eigen cation, in which each water molecule transfers an equal amount of electron density (see bottom rows of Table 1), is absent in the Zundel cation. The majority of the charge in the Zundel is transferred from the ‘special pair’ partner water, O1x. In fact, O1x transfers ~65% of ECT (16.94 kcal/mol in the HF/aug-cc-pVTZ analysis) while the remaining four water molecules transfer only 8 to 9% each.
To characterize the dependence of the various energetic components on the distance between the hydronium and water molecules, a “solvation” coordinate (see Fig. 1) was defined. This coordinate involves symmetrically pulling the three water molecules in the Eigen cation away from or toward the central cation. At each point along this coordinate, the structures were geometry-optimized at the MP2/aug-cc-pVTZ and B3LYP/aug-cc-pVTZ levels, and analyzed with ALMO EDA. This coordinate is not meant to represent a solvation process, but as a device that allows us to characterize the energy components in various structures discussed later as functions of distance between the oxygen atoms (ROO). Again, good agreement was observed between the MP2 and B3LYP optimized structures at all values of ROO. The largest MP2-B3LYP structural difference occurred at the shortest distance between oxygen atoms (ROO = 2.20 Å), and was 0.24 Å, 4.25°, and 42.51° for the average bond, angles, and dihedrals, respectively.
The EDA energy components calculated at the HF/aug-cc-pVTZ level as a function of ROO are shown in Fig. 2. A nearly identical plot was obtained with the DFT EDA data. As expected, the electrostatic (EFRZ) energy is attractive at large ROO, but becomes repulsive and rises quickly as the waters are pushed toward the central cation and charge densities are forced to overlap. The polarization (EPOL) and charge transfer (ECT) energies are attractive at large ROO and become more so as the monomers are brought together (for ROO ≥ 2.25) with the polarization term dropping at a slightly faster rate. For ROO ≤ 2.25 Å, ECT increases rapidly, indicating that there is a limit to increasing charge transfer with decreasing inter-atomic distance. At long distances, the binding energy is solely determined by ΔEFRZ. These interaction components provide insight that will prove predictive in the analyses described below, and may offer an appealing new way to define and parameterize more physically meaningful empirical potentials.
Given the above baseline ALMO analysis of geometry-optimized (gas-phase) clusters, we now turn our attention to structures extracted from the bulk phase. The average quantities presented in this section are the result of analyzing, as previously described, hundreds of structures taken at random time points during the MSEVB and CPMD dynamics simulations. Among the many structures taken from the bulk phase, a pure hydronium ion (defined by ci2 > 0.82 in the MS-EVB simulations and by ECT < 10 kcal/mol in all simulations) was never observed. The skewed Gaussian distributions (Fig. 3) of δ values (defined in Eqn. 5) identified the Eigen cation (δ ≥ 0.2 Å) in 67%, 74%, and 66% of the structures from the MS-EVB3, CPMD/BLYP, and CPMD/HCTH simulations, respectively. Similarly, the Zundel cation (δ ≤ 0.1 Å) was represented by ~ 16%, 11%, and 15% of these simulation structures. However, the distributions of δ-values show no sign of bimodal character, suggesting there is no demarcation line between Eigen and Zundel character.
An important question to address is whether or not the charge transfer calculated by applying ALMO EDA to four- to ten-water clusters would be significantly altered if a more extended solvent environment (i.e., a polarizing medium) were included. Basic physical principals tell us that a polarizing medium 1) causes a delocalized charge to condense, and 2) dampens intermolecular electrostatic interactions. Both phenomena suggest that surrounding a system with a polarizing medium should, on average, decrease the charge transfer between the central cation and the surrounding water molecules. It is not clear, however, if the charge transfer, for a given conformation such as those we are taking from bulk phase simulations, would be significantly decreased due to a surrounding polarizable medium. To address this question, we used Natural Population Analysis (NPA),60 both including and without a polarizable continuum model (PCM). We employed the NPA method because including a polarizable medium is not possible in the present version of the ALMO EDA code. Population analysis methods (e.g., Löwdin, Mulliken, and Natural PA)60,64-66 assign one-electron basis functions to nuclei and thereby define the formal charges. From these charges, the amount of charge transferred from one atom to another can be inferred. Population analysis methods differ in the way they define the basis functions and deal with orbital overlap. Although it has been shown that NPA can overestimate effective charges,67 it is more consistent across basis sets and molecular conformations than Mulliken PA.67 Moreover, the important factor in this study is the change in the effective charges with the addition of a PCM.
NPA of the geometry-optimized Eigen cation at the MP2/aug-cc-pVTZ level assigns a partial charge of 0.068 to the three peripheral waters, meaning that 68 me− is transferred from each peripheral water to the hydronium. Surrounding the Eigen cation with a PCM and using a dielectric constant representative of water (ε = 80) decreased that value by less than 1% (Δ q = −0.48 mē). We also tested several asymmetric structures extracted from the bulk-phase simulations. The NPA-assigned partial charges for these peripheral water molecules varied from 0.05 up to 0.14 for the water furthest and closest to the central cation's hydrogen atoms, respectively. Including the PCM decreased these water partial charges by anywhere from Δ q = 0.05 mē (0.02%) to Δ q = 0.34 mē (0.64%). Interestingly, charge transfer for the waters closest to the hydronium (i.e., those donating the most charge) was not decreased but increased by 0.16 to 0.55 mē when the PCM was included. Collectively, these results demonstrate that surrounding a cluster conformation with a polarizable continuum has a negligible effect on the local redistribution of electron density (i.e., charge transfer) within Eigen and Zundel-like species. Thus, our analysis of the charge transfer in clusters containing the essential surrounding water molecules, but without a PCM, is indeed a reasonable representation of the bulk phase environment.
As detailed earlier, many structures from each simulation (MS-EVB3, CPMD/BLYP, and CPMD/HCTH) were analyzed with ALMO EDA at the HF/aug-cc-pVTZ level for comparison across simulation approaches. The conformations extracted from both CPMD simulations were additionally subjected to ALMO EDA with the DFT functional used in the dynamics simulation to check for method and/or DFT functional dependence. The ALMO EDA results, averaged over the ensemble of Eigen structures extracted from the dynamics simulations, are presented in Table 4. The first point that must be emphasized is that the charge transfer in these asymmetric structures is not less, but consistently more than that in the symmetric (i.e., geometry-optimized) structure (shown again in column 1 of Table 4 for comparison). This increase is a result of the ‘special pair’ partner water molecule (O1x in Fig. 1) being pulled more closely to O0 and transferring more charge (see bottom rows in Table 4) than it had in it's geometry-optimized position. The water molecule with the second smallest δ value, O1y, also contributes more on average, though just slightly, while the third, O1z, contributes less. The bulk-phase average distances (Table 3), which are in excellent agreement with previous work and recently reported radial distribution functions,5,68 show that both O1x and O1y are closer to O0 in the bulk phase average than they are in the optimized structure, while O1z is further away. In fact, the O0-O1x distance was found to be shorter than 2.55Å (the optimized distance) for more than 93% of an entire 1 ns MSEVB3 trajectory. We note that the increase in charge transfer for waters O1x and O1y (and decrease for O1z) could have been anticipated from the distance dependence of ECT shown in Figure 2.
Comparing across simulations, the HF EDA results are very similar. However, there is a slight increase in ECT and EPOL in the CPMD/HCTH simulation, suggesting that this ensemble had more asymmetric structures with closer water molecules. This is supported by comparing the ensemble averaged O1X:O1y:O1z charge transfer ratios, which were 1.7:1.3:1, 2.4:1.1:1, and 2.4:1.8:1 for the MSEVB, CPMD/BLYP, and CPMD/HCTH ensembles, respectively. As was the case for the optimized structures, DFT-based EDA predicts somewhat more favorable ECT and EPOL components and less favorable EFRZ contributions than HF-based EDA. The differences in energy components (especially for HF and HCTH/120) are also very similar to what they were for the optimized structures. This indicates that HF and DFT perform similarly for the optimized and bulk phase structures. Moreover, the same trends are present in the DFT and HF EDA, suggesting that the trends are qualitatively correct.
The averaged ALMO EDA results for the bulk-phase Zundel cation are shown in Table 5. As with the Eigen cation, going from the geometry-optimized Zundel structure to the bulk-phase structures makes EBIND slightly lower (less favorable) and ECT larger. Therefore, the relative charge transfer contribution to the binding energy is also larger (36.4% for MSEVB3/HF). The increase in ECT is over 7% for the Zundel cation, larger than it was for the Eigen cation, jumping from 32.16 kcal/mol to 35.26 kcal/mol for the average HF EDA of the MS-EVB simulation. Also similar to the Eigen cation, EFRZ is less favorable and EPOL is more favorable in the bulk-phase averages. These results again suggest that, in the asymmetric bulk-phase ensemble, the water molecules are being pulled in more closely to the central H3O+ than they are in the geometry-optimized structure. Of course, the total binding energy, EBIND, is larger for the Zundel clusters than for the Eigen because there are six water molecules in the former as opposed to just four in the latter.
Comparing the MS-EVB, BLYP, and HCTH/120 simulations reveals the opposite trends from those demonstrated for the Eigen cation. Instead of being less favorable, EFRZ is more favorable in the BLYP and HCTH simulations than in the MS-EVB simulation. Similarly, EPOL and ECT are less, instead of more, favorable. These changes seem to be caused by weaker 2nd solvation shell interactions (i.e., O0-O1z, O1x-O2x1, O1x-O2x2). This is not the case, however, for O0-O1y. Again, this points to slightly more asymmetric structures in the CPMD simulations.
It has been previously suggested that charge transfer to the central cation occurs over multiple solvation shells and that the 2nd shell actually transfers more charge to the central H3O+ than the 1st shell does.34 Our results do not support this conclusion. Rather, charge transfer is found to be a predominantly local phenomenon (i.e., one solvation shell to the next) that is minimally altered by the surrounding environment. Analyzing the energetic components of a 20-water molecule cluster and sequentially removing shells verified that charge transfer between shells is independent. For example, although the 3rd shell transfers charge to the 2nd, removing the 3rd shell entirely decreases the charge transfer energy of the 2nd to the 1st and the 1st to the central cation by only 0.99 kcal/mol (1.4%) and 1.1 kcal/mol (3.2%), respectively. Similarly, the 2nd shell donates charge to the 1st, but removing the 2nd shell changes the charge transfer from the 1st to the central H3O+ by less than 7%. Thus, for a given conformation, ECT is primarily caused by the neighboring molecules. It should also be noted that the 2nd-to-1st shell charge transfer values were significantly larger than the 3rd-to-2nd shell values. The third shell's interactions are likely similar to those in bulk water, but the 2nd shell is more strongly involved in the total electron density defect.
To characterize the dynamic nature of charge delocalization, clusters containing six water molecules were extracted at every step (1 fs) during segments of the MS-EVB3 simulation and analyzed with ALMO EDA at the HF/aug-cc-pVTZ level. Six-water structures were chosen to maximize observable continuity as the protonated complex changes between Eigen-like and Zundel-like species. Analyzing ten-water clusters to capture the complete 2nd solvation shell of the Eigen cation over the entire time series was computationally expensive and only carried out for one time series. Analyzing the ten-water clusters verified that all six water molecules in the Eigen cation's 2nd solvation shell transfer approximately equal amounts of charge and result in average stabilization energies from 1.1 to 1.6 kcal/mol in the HF ALMO EDA. Moreover, it verified that the qualitative trends presented below are not altered from the absence of four waters in the six-water clusters.
This dynamical analysis revealed several interesting points that are completely absent in the ensemble averages discussed earlier. First, as one might predict from the distance dependence shown in Fig. 2 and the geometric averages shown in Table 3, the magnitude of charge transfer from each water molecule is predominantly determined by that water's hydrogen bond length to the central cation (RH--O1) and the symmetry between donor and acceptor oxygen atoms (i.e., δ in Eq. 5). This trend is shown in Figure 4 where the charge transfer energy from H2O1x to H3O0+ is well correlated with both RH--O1x and δO0-O1x. However, neither correlation is perfect, indicating that other factors, such as 2nd shell hydrogen bonds, intermolecular angles, or other 1st shell interactions, are also likely to influence charge transfer.
As described in section 2.4, the Eigen and the Zundel cations are defined according to their δO0-O1x values (Eqn. 5) or by tracking the O0 and O1x identities. To give the reader a better feel for these two approaches, Figure 5 shows both the δO0-O1x and O0, O1x identity analyses over 30 ps of the MS-EVB3 simulation. There is a strong correlation between more dense regions having δO0-O1x < 0.1 (i.e., Zundel as defined by δ) and those where the O0 identity is changing (Zundel). There are also, however, numerous occasions where δ1 is < 0.1 in regions where O0 does not change, suggesting that highly asymmetric structures occur frequently without acquiring Zundel-like dynamics.
Repeating temporal patterns can be seen in the charge transfer magnitudes for both Eigen- and Zundel- dominated sections of the trajectory shown in Figure 5. For example, a closer inspection of the segment near 30,000 fs in Figure 5 (Figure 6B) shows that the net charge transfer energy from the Eigen cation's 1st solvation shell and the net δ structural parameter track one another as they oscillate smoothly over time. Such oscillations might be expected for any single water-hydronium interaction, but not necessarily for the net (total) 1st-shell ECT. The δO0-O1x and O0/O1x identities shown in Figure 6A verify that this segment is representative of the Eigen cation. Even though δO0-O1x is < 0.1 several times in the last 30 fs, it takes another 50 fs before O0-O1x exchange occurs (i.e. reaching Zundel-like dynamics). Indeed, it was found that in almost every Eigen-like simulation segment, δO0-O1x drops below 0.1 multiple times. This suggests that the Eigen cation can approach Zundel-like structures without transitioning to Zundel-like dynamics (defined by rapid changes in the O0 and O1x identities). In these instances, the Eigen cation is clearly very asymmetric. This emphasizes how challenging it is to draw a clear distinction between the Eigen and Zundel cations based on structure alone.
The simulation segment shown in Figure 7 (from the time interval near 12,000 fs in Fig. 5) is representative of the Zundel cation according to both the δO0-O1x and O0/O1x identity criteria (Fig. 7A). Similar to the Eigen dominated-period, the 1st solvation shell net δ values track the corresponding net ECT quite well (results not shown). Collectively, these Eigen and Zundel results indicate that geometric properties (RO1-H and δO0-O1), which can be easily tracked over the entire simulation, are representative of the relative charge transfer energy from each water molecule to its hydrogen bond donor and thus can serve as useful and convenient descriptors of the charge transfer.
One interesting result for the Zundel cation is that the four water molecules strongly associated to the ‘special pair’ do not transfer charge symmetrically when δO0-O1x approaches zero. This supports recent structural analysis (see Fig. 5 in reference 5) that there can be significant asymmetry during an O0-O1x exchange. In Figure 7B the special pair bonds (RO0-H and RH--O1x) are plotted along with the net 2nd shell hydrogen bonds lengths for each water (i.e., O1y + O1z for O0 and O2x1 + O2x2 for O1x). Again, these hydrogen bond lengths are representative of the strength of charge transfer between each pair of water molecules. The distances are shown instead of the charge transfer energies in order to demonstrate the special pair oscillation, which would have been lost by examining a single O0-O1x charge transfer magnitude. When the bottom two lines in Figure 7B cross, the ‘special pair’ water molecules are equidistant from the central hydrogen atom (i.e., δO0-O1x = 0). There are several such δO0-O1x= 0 crossing points where the upper two lines do not cross, meaning that ECT from the 2nd solvation shell can be quite different for the O0 and O1x sides even when δO0-O1x= 0. There is additional asymmetry in the cluster due to differences between the two charge donors on the respective sides (e.g., different ECT from O1y and O1z). These results support the notion that some PT events, whether transient hops or successful transfers, can be described in terms of a very asymmetric Eigen cation on the donor side evolving into an asymmetric Eigen cation on the acceptor side.
It should be emphasized that the Eigen- and Zundel-dominated segments of the trajectory shown in Figures 6 and and77 are representative of patterns that are repeated throughout the simulations. They are not, however, meant to represent statistically determined macroscopic properties.
It is clear that the heterogeneous bulk-phase environment introduces a great deal of asymmetry in the electron density distribution around the central cation. Given that the optimized gas-phase structures are symmetric, it could be postulated that the asymmetry is caused by the fluctuations in the surrounding heterogeneous environment and that the complex is always tending toward its minimum energy, a symmetric conformation. However, the results shown in Figure 8 suggest this is not the case. This plot shows that the stability of the protonated complex, EBIND, is almost perfectly anti-correlated with the symmetry around the central cation (as defined by δO1y − δO1x + δO1z − δO1x). Therefore, asymmetry appears to have a stabilizing effect even in the interaction energy. This stabilization will likely increase in a free-energy analysis due to the entropic freedom of asymmetric structures compared to the configurational limitations of symmetric conformations.
As seen in Figures 4--8,8, oscillatory signatures emerge within the time evolution of various charge-transfer quantities, especially between the special pair (O0-O1x) and the net charge transfer from the 1st solvation shell (as defined for the Eigen cation). The frequencies of these oscillations can be measured through the entire simulation by tracking, for example, the geometric properties shown in Figure 9 (δO0-O1x and the sum of δ values for O1x, O1y, and O1z). Although, these charge transfer oscillations and their resulting frequency distributions cannot be directly related to the intensities and shapes of IR spectroscopic peaks, they are, we believe, worth discussing here to help characterize the dynamics of the delocalized charge defect.
The δO0-O1x frequency distributions shown in Figure 10 were taken from a 1 ns MS-EVB3 simulation by creating a histogram of the inverse of the temporal oscillation recurrence. During each oscillation, the structure was attributed to being undergoing either Zundel or Eigen dynamics depending on whether the mean δ value was ≤ 0.1 or ≥ 0.2, respectively. Given that the O0-H bond and H--O1x hydrogen bond lengths are anticorrelated and result in almost identical frequency distributions (results not shown), it is not surprising to see a large peak at ~2800 cm−1 (~12 fs recurrence time), a frequency range that has been associated with the Eigen OH stretch in previous theoretical and experimental work.36 The Zundel frequencies are seen to be lower than the Eigen's, peaking at 2200 cm−1 (~15 fs), but they are much higher than for the O-H..O stretch in the gas-phase (geometry-optimized) Zundel cation (1000 cm−1; ~33 fs).35 An increase in the frequency of this stretching motion was experimentally observed in larger protonated water clusters35 and was recently discussed by Buch et al..19
Perhaps the most interesting oscillations are the slower modes that are apparent in both ECT from O1x to O0 and net ECT from the 1st solvation shell. Figure 9 shows both the fast and slow oscillations in δO0-O1x and in the net 1st shell δ values. The slower motions produce a wide range of oscillation frequencies from 200 to 1000 cm−1 and appear to be correlated with a stretching motion of the 1st solvation shell around the central cation and with an symmetric bending motion along the principle axis of the central cation's oxygen atom.
This work explores the role of charge transfer in the structure and dynamics of the hydrated proton with ab initio energy decomposition analysis, and has come to the following conclusions:
Collectively, these results suggest that the degree of electron density delocalization and effective covalency in the hydrated proton complex are not consistent with the existence of a purely classical hydronium ion (H3O+), but are in agreement with previous simulations2,3,5,6,9-11,15-19,34 that have pointed to an entity of added complexity. Therefore, the description of acid-base equilibria (particularly its statistical mechanical formulation) and computer simulation models formulated in terms of classical hydronium ions will benefit from revision.
JMJS was supported by the National Institutes of Health under a Ruth L. Kirschstein National Service Award (F32 GM078905). The authors thank Dr. Sergey Izvekov for providing the results of the CPMD simulations, Brian W. Hopkins and Gregory A. Voth for helpful discussions, as well as Yihan Shao and Zhengting Gan for providing the beta version of QChem and the associated assistance.