An Improved Multistate Empirical Valence Bond Model for Aqueous

Nov 14, 2007 - third issue concerns the proton self-diffusion constant in liquid water, which is ..... reliable guesses was chosen for manual adjustme...
1 downloads 0 Views 1MB Size
Download PDF
J. Phys. Chem. B 2008, 112, 467-482

467

An Improved Multistate Empirical Valence Bond Model for Aqueous Proton Solvation and Transport† Yujie Wu,‡,§ Hanning Chen,‡ Feng Wang,⊥ Francesco Paesani, and Gregory A. Voth* Department of Chemistry and Center for Biophysical Modeling and Simulation, UniVersity of Utah, 315 South 1400 East, Room 2020, Salt Lake City, Utah 84112-0850 ReceiVed: August 18, 2007; In Final Form: September 25, 2007

A new multistate empirical valence bond model (MS-EVB3) is developed for proton solvation and transport in aqueous solutions. The new model and its quantum version (qMS-EVB3) are based on the MS-EVB2 model [Day et al., J. Chem. Phys. 2002, 117, 5839] and recently developed flexible water modelssthe SPC/ Fw model [Wu et al. J. Chem. Phys. 2006, 124, 24503] and the qSPC/Fw model [Paesani et al. J. Chem. Phys. 2006, 125, 184507]sfor classical and quantum simulations, respectively. Using ab initio data as benchmarks, the binding energies and optimized geometries calculated with the new model for protonated water clusters, as well as the potential energy surface for proton shuttling between water molecules in a cluster environment, are improved in comparison to the MS-EVB2 model. For aqueous solutions, classical and quantum molecular dynamics simulations with the MS-EVB3 model yield a more accurate description of the solvation structure and diffusive dynamics of the excess proton. New insight is also provided into the proton solvation and hopping dynamics in water, as well as the “amphiphilic” nature of the hydrated proton that has been predicted to give rise to its enhanced concentration at aqueous interfaces and an effectively lower pH of the air-water interface [Petersen et al. J. Phys. Chem. B 2004, 108, 14804].

1. Introduction The importance of proton solvation and transport (PS&T) as a fundamental process in acid-base chemistry was realized more than two centuries ago. More recently, along with advances in biochemistry and molecular biology, the importance of PS&T has been further established in biology, where various important processes (e.g., mitochondrial respiration and ATP biosynthesis1) have been found to critically depend on PS&T.2 Proton conducting channels or pathways have been identified in many biomolecules. For example, the F0 component of the ATP synthase is in essence a proton channel,3 while the transmembrane M2 protein of the influenza A virus forms a proton channel that is essential for the viral uncoating process.4 Indepth computational studies of PS&T in different environments are important to our understanding of the related chemical, biochemical, or biophysical processes.5,6 A wealth of scientific effort has been invested in studies of PS&T during the last century. In general, two important properties of PS&T have been noted: First, PS&T is strongly influenced by the Grotthuss proton shuttle mechanism, which causes the delocalization of the excess proton to more than one molecular species.5 Because of this delocalization, the protonic charge is also distributed across two or more molecules along the hydrogen-bond network. This property of the excess proton thus has a strong polarization effect arising from the translocation of the nucleus, often called the “Zundel polarization”.7,8 †

Part of the “James T. (Casey) Hynes Festschrift”. * Corresponding author. Fax: 801-581-4353. Phone: 801-581-7272. E-mail: [email protected]. ‡ These authors contributed equally to this work. § Current address: Schro ¨ dinger, LLC, 120 W. 45th Street, 29 Floor, New York, NY 10036. ⊥ Current address: Department of Chemistry, Boston University 590 Commonwealth Ave. Boston, MA 02215.

Second, the nuclear quantum effects may be non-negligible in PS&T, even at room temperatures, because of the light mass of the proton. Despite important recent progress, many important aspects of the proton delocalization process, for example, its solvation structure and dynamics, remain incomplete and somewhat elusive.6 This is true even for systems as simple as bulk water, not to mention much more complicated biomolecular systems.2,9 One of the most interesting features of PS&T is the anomalously high proton mobility in liquid water under ambient conditions.10,11 It is observed that proton mobility is roughly 5-7 times higher than other “standard” cations (e.g., K+) that have an ionic radius similar to the hydronium ion (H3O+), the smallest stable hydrated form of the excess proton in water. This property of excess protons has been attributed to the so-called Grotthuss mechanism.12,13 Although it is generally agreed that the solvating-water dynamics is critical in this mechanism,13,14 the ratelimiting step of the proton translocation process and related structural/chemical details remain to be clarified.9,15,16 Another example is the microscopic hydration structures of the excess proton. It is generally accepted that, in aqueous solutions, the proton shuttle mechanism involves fast interconversions between two distinct proton hydration structures, namely, the Eigen17 ([H3O‚(H2O)3]+) and Zundel18 ([H‚(H2O)2]+) cations, proposed 30-40 years ago. A debate still remains on the relative populations of the two forms in the equilibrium solution.19,20 For this problem, definitive experimental measurements remain a challenge (e.g., see ref 21), and the difficulty is especially severe for dilute solutions. The main difficulty for a theoretical atomistic description of PS&T arises from the aforementioned properties of PS&T, which may require that electronic polarization and nuclear quantum effects be considered and that the Zundel polarization be included in a way consistent with the proton shuttle mechanism.

10.1021/jp076658h CCC: $40.75 © 2008 American Chemical Society Published on Web 11/14/2007

468 J. Phys. Chem. B, Vol. 112, No. 2, 2008 One significant advance within the past decade on the theoretical/computational side has been the development of the multistate empirical valence bond (MS-EVB) method. The MSEVB method allows explicit PT reactions to be simulated both efficiently and accurately in molecular dynamics (MD) simulations (see refs 5 and 6 for recent reviews). In principle, this in turn allows for a statistically reliable and accurate characterization of the detailed solvation structure and dynamics of excess protons in different environments. In recent years, force fields22,23 based on the MS-EVB method have been applied to investigate a large number of different fundamental PS&T problems such as the energies of protonated water clusters,24,25 the Grotthuss mechanism,15,23,26 the vibrational spectrum of the hydrated proton,27 PS&T in water-filled hydrophobic tubes,28,29 PS&T in supercritical water,30 and hydrated proton structure and the resulting enhanced acidity of the water-air interface.31 Moreover, the MS-EVB2 model23 has been applied to studies of PS&T in a variety of biomolecular systems, such as the LS2 and the influenza A M2 proton channels,32-35 the aquaporin channels,36-38 lipidmembranes,39-41 andcytochromecoxidase.42-44 In these studies, the characterizations of detailed proton hydration structures, dynamics, and proton permeation free-energy profiles have provided valuable mechanistic insight into molecular mechanisms such as proton-channel gating,34 selectivity,35 and blockage.36-38 The theoretical basis for the MS-EVB model is an extension/ generalization22,45-47 of the original empirical valence bond (EVB) method pioneered by Warshel,48-50 which drew from the valence bond (VB) ideas of Coulson and Mulliken.51,52 The EVB method describes a simple bond-breaking/formation reaction through “states”, typically in the form of one reactant-like state and one product-like state, where the energies of the two states are described empirically, that is, usually in the form of molecular mechanics force fields. This approach leads, in principle, to accurate yet computationally efficient descriptions of chemical reactions that involve bonding topology changes. By contrast, in the MS-EVB approach, the conventional EVB theoretical framework is extended to include multiple states that may also change their identity “on the fly” with the dynamical evolution of the system. This latter feature is essential for a physical description of reactions with multiple possible pathways and multistep coupled reactions where the product of one step is the reactant of the next, as well as reactive processes that also involve an interplay of reaction with molecular diffusion. In the case of aqueous PS&T, this multistate EVB extension is crucial since, in the Grotthuss mechanism, proton-transfer reactions are coupled with molecular diffusion and hydrogenbonding network rearrangements, so multiple pathways are possible from the same starting structure, and these pathways change with time. The delocalization of the excess proton also leads to the formation of the Eigen and Zundel complex as well as many possible intermediate structures. The MS-EVB extension is therefore fundamental for simulating in a deterministic fashion the dynamics of the migration of the excess proton along a constantly rearranging and diffusing aqueous hydrogenbonding network, which is in turn essential for a physically correct characterization of the PS&T process. During the late 1990s, both the Voth22,47 and Borgis45,46,53 groups independently developed their multistate models. Both of these models were extensions of the original two-state model for aqueous proton transfer of Lobaugh and Voth published in 1996.54 Our initial MS-EVB model,22 referred to in this paper as the MS-EVB1 model, provided good agreement with highlevel ab initio results for the geometries and binding energies

Wu et al. of small protonated water clusters.22,47 Good agreement with experimental measurements of the spectroscopic features and density of vibrational states,27 and of the relative stabilities of the Eigen and Zundel cations22 were also achieved for liquid water. This model was subsequently improved by introducing two repulsive terms to replace the ad hoc charge-scaling scheme for the hydronium-water electrostatic interactions and, importantly, by refining the MS-EVB state-selection algorithm. The resulting second-generation model, MS-EVB2,23 showed similarly good properties with somewhat more accurate energies for small protonated water clusters. However, it was shown that the MS-EVB2 model followed “real” energy-conserving Newtonian dynamics with higher fidelity by significantly alleviating the energy drift resulting from changes of the MS-EVB state identities during the MD simulation. This critical feature of any deterministic (Newtonian or extended Newtonian dynamics) MD model cannot be emphasized enough because, without good energy conservation, one cannot be assured that the proper statistical distribution is being sampled by the MD trajectory. The MS-EVB2 model has also been shown to predict the activation energy associated with proton diffusion in good agreement with the experimental value.15 The present study presents the results of our continuing effort to increase the accuracy and reliability of the MS-EVB model for aqueous and biomolecular PS&T. In particular, this work addresses three important issues. The first is that, although the MS-EVB2 model improved several aspects of the original MSEVB model as described above, the underlying water potential was not improved.55 The second regards an unphysical small artifact in the hydronium oxygen-water oxygen radial distribution function (RDF) that is due to the repulsive terms used in the MS-EVB2 model (see Results and Discussion section). The third issue concerns the proton self-diffusion constant in liquid water, which is underestimated by more than 50% relative to experiment in the classical MD implementation of the MS-EVB2 model.23 The new MS-EVB3 model described in this paper therefore provides a more accurate description of the structural and dynamical properties of the excess proton in aqueous solutions. For example, the ratio between the proton and water diffusion constants computed using classical MD simulations is improved within the MS-EVB3 framework with respect to the previous MS-EVB2 model. An explicit treatment of the nuclear quantum effects is also shown to achieve an even more accurate description of the diffusive dynamics of the excess proton, which results in better agreement with the experimental data. Most importantly, the MS-EVB state-selection algorithm is further refined to provide significantly better total energy conservation in the MS-EVB MD simulations. The remaining sections of this paper are organized as follows: In section 2, the theoretical framework of the MSEVB model is first reviewed, and then the new state-selection algorithm, the parameter optimization procedure, and details of the MD simulations are described. In section 3, extensive characterization of the new modelssMS-EVB3 and its quantum version qMS-EVB3sare presented and discussed. Conclusions and future directions are then given in section 4. 2. Methods 2.1. The Theoretical Framework of the MS-EVB Model for a Hydrated Excess Proton. Because a hydrated excess proton is delocalized among several solvating water molecules, and thus is indistinguishable from water hydrogen atoms, a conventional picture of Lewis bond arrangement cannot be

The MS-EVB3 Model for Aqueous Proton Transport

J. Phys. Chem. B, Vol. 112, No. 2, 2008 469

established unambiguously. This rules out any straightforward application of conventional force field methods. The delocalized solvation structure is described in the MS-EVB approach via a linear combination of several basis states (|i〉) that correspond to distinct bonding arrangements (see below). The linear combination gives a state function (|Ψ〉) of the excess proton for every configuration of the system nuclei x˜ as N

|Ψ〉 )

∑i ci|i〉

(1)

where N is the total number of the basis states, and ci’s are the coefficients. The state function |Ψ〉 is obtained by constructing the Hamiltonian matrix [H, with its elements hij ) 〈i|H|j〉 usually described by an empirical force field (see below for details)] and subsequently by solving the eigenvalue-eigenfunction problem:

Hc ) E0c

(2)

where E0 ) E0(x˜ ) is the lowest eigenvalue, corresponding to the ground-state energy of the system, and c is the corresponding eigenvector (whose elements are ci, i ) 1, 2, ..., N). The force exerted on each atom for every nuclear configuration can then be calculated using the Hellmann-Feynman theorem:

〈

| |

Fi ) - Ψ 0

∂H ∂xi

〉

Ψ0 ) -

cmcn ∑ m,n

be rigorously measured by its EVB amplitude [ci2, with ci given as in eq 1]. In reality, the ci2 are unknown before the identification of the significant states, hence leading one to an algorithmic challenge. To address this challenge, a geometrybased selection algorithm was developed in the original MSEVB1 model and further improved in the MS-EVB2 model. Although the algorithm in the MS-EVB2 model has been shown to significantly improve the energy conservation,23 it has been found to miss states that contribute non-negligibly to H under some circumstances (see Results and Discussion section). Thus, the critical MS-EVB state-selection algorithm is further improved in this work. Details of the new state-selection algorithm will be described in section 2.2. The diagonal elements hii of H in eq 2 are given by the potential energy function of each basis state, for which the present model uses the following expression: NH2O

hii ) V Hintra + + 3O

∑k

NH2O

V Hintra,k + 2O

∑k

NH2O

V Hinter,k + + 3O ,H2O

2

(4)

where the VHintra + term is the intramolecular potential of the 3O hydronium ion, defined as 3

V Hintra + ) 3O

∂hmn

∑j DOH[1 - e-a

j 0 OH(ROH-ROH)

]2 + 1

(3) ∂xi

It is well accepted that proton translocation involves, in part, a number of consecutive or concerted proton-transfer reactions along the hydrogen-bonding network formed by the water molecules. In each of these reactions, a hydrogen atom is transferred from a hydronium ion to a water molecule, turning the latter into a new hydronium ion. A single such reaction can be described by two EVB states, that is, a reactant state and a product state.54 Each of the two states consists of a hydronium ion plus all of the remaining water molecules, and we do not distinguish a product state from a reactant state in the MS-EVB formalism. This is very important in describing consecutive reactions since the product state of one reaction is the reactant state of the next. The Grotthuss translocation process12,13 can be described by including states from all possible proton-transfer reactions into the MS-EVB Hamiltonian H. In principle, any bonding topology that gives the correct numbers of hydronium ion(s) and water molecules will be a valid MS-EVB state. This will result in 20 states for even the simplest gas-phase protonated water dimer (H5O2+), hence enumeration of a complete set of states is intractable even for a small system. Fortunately, most of these states have a negligible contribution to the ground state eigenvalue E0(x˜ ), so that they can be safely omitted when constructing H. For the protonated water dimer, for example, two states are normally sufficient. In bulk water under ambient conditions, it is found that, on average, 22 states are normally required for describing the solvation state of a single excess proton. However, up to 40 states must sometimes be allowed in MS-EVB MD simulations in order to reasonably conserve the total system energy. The proper identification of the significant EVB states is critical for an efficient and accurate determination of H on the fly in the MS-EVB simulation. Inclusion of unnecessary states will waste central processing unit (CPU) cycles, whereas missing significant states will result in a less accurate H and usually cause poor energy conservation. The significance of a state can

V Hinter,kk′ ∑ O k rc 0, is multiplied on the repulsive terms [r is ROOk for eq 7 or RHOk for eq 8]. The switching ranges (rs-rc) are 2.85-3.05 Å and rep rep and VHO , respectively. 2.3-2.5 Å for VOO k k The off-diagonal (coupling) elements hij (i * j) in eqs 2 and 3 are given by

hij ) (V ijconst + V ijex)‚A(ROO,q)

(10)

if the two states |i〉 and |j〉 share a common hydrogen atom such that the system going from state |i〉 to state |j〉 corresponds to a proton-transfer reaction, while the hij terms are zero otherwise. In eq 10, Vijconst is a constant, while Vijex is the electrostatic potential between the H5O2+ Zundel complex (formed by the hydronium and an acceptor water) and the remaining NH2O - 1 water molecules, given by 7 NH2O-1 3

V ijex

)

∑ ∑k ∑ m n k

qnHk2Oqex m (11)

Rmnk

H2O where qex are the exchange charges of the H5O2+ m and qnk adduct and the atomic charges of the water molecule k, respectively. The term A(ROO,q) in eq 10 is a geometrydependent scale factor:

A(ROO,q) ) exp(-γq2)‚{1 + P exp[-k(ROO - DOO)2]}‚ 1 {1 - tanh[β(ROO - R0OO)]} + 2

{

}

P′ exp[-R(ROO - r0OO)] (12) where ROO is the distance between the two oxygen atoms of the H5O2+ adduct, and q ) (rO + rO′)/2 - rH* is the asymmetric stretch coordinate of the adduct (rO and rO′ are position vectors of the two oxygen atoms, and rH* is that of the central hydrogen atom). This definition of q is different from that in the original MS-EVB2 formalism,22 since the original definition gives a cusp in the potential when the hydrogen passes through the bisecting plane between the two oxygen atoms, as noted previously.56 2.2. The MS-EVB3 State-Selection Algorithm. The improved MS-EVB3 state-selection algorithm is described as follows: (1) In order to decide which states should be included in the MS-EVB Hamiltonian H at each MD time step, the state having the largest contribution to H must first be identified. During an

Figure 1. Schematic diagrams of the (a) Type-I and (b) Type-II bifurcated hydrogen bonds between water molecules included in the MS-EVB3 state searching algorithm.

MD simulation, the piVot state is defined as the one carrying the largest EVB amplitude at the previous time step. For the first step, the pivot state is chosen as the “best” hydronium configuration, which is normally formed by the excess proton and its closest water molecule. (2) The pivot state is regarded as the reactant state, and the algorithm then searches for significant product states to which this reactant state can transform by donating one of its hydrogen atoms. All water molecules that are within a 2.5 Å radial distance from each hydronium hydrogen atom are first selected. Each of these water molecules can potentially form a new hydronium ion by accepting the transferring hydrogen atom, thus forming a distinct product state. (3) For each potential product state, it is then treated as a reactant state, and the selection procedure is repeated as described in step 2. (4) Steps 2 and 3 are iterated until all states within (and including) the third solvation shell of the initial pivot hydronium (i.e., the hydronium in the pivot state) are included. It is important to note that there are two main differences between the present state-selection algorithm and that of the MS-EVB2 model.23 First, if the hydronium hydrogen atom forms two or more (bifurcated) hydrogen bonds with a water molecule within the 2.5 A cutoff distance, all of the these water molecules are allowed to form a state, whereas, in the previous algorithm, only the “best” one (that with the shortest hydrogen bond length) was selected. By virtue of this treatment, the present algorithm takes into account the Type-I bifurcated hydrogen bonds schematically shown in Figure 1a. Second, a water oxygen atom can participate in more than one state as the hydronium oxygen atom, as long as the hydronium hydrogen atoms are not all identical, whereas, in the previous algorithm, it can participate in only one state. The present algorithm also takes into account the Type-II bifurcated hydrogen bonds as shown schematically in Figure 1b. Although bifurcated hydrogen bonds are relatively rare in liquid water, they seem to be essential in hydrogen-bond swapping. The corresponding effects have turned out to be important for a more accurate MS-EVB Hamiltonian and hence for better energy conservation in MD simulations. 2.3. Parametrization of the MS-EVB Model. Two recently developed flexible water potentials55,57 were used within the MS-EVB3 framework for classical and quantum simulations, respectively. As discussed in detail in ref 57, for quantum simulations, use of a water model specifically parametrized as such is necessary to avoid a double counting of some of the quantum effects. The use of these new and more accurate underlying water models required nearly all parameters in the MS-EVB3 force field to be reoptimized. The parameter optimization was performed to fit the following ab initio (MP2/ aug-cc-pVDZ) results: (1) the binding energies and geometries of protonated water dimer and tetramer, and (2) the PES of the proton shuttling within a tetramer water wire. Specifically, the

The MS-EVB3 Model for Aqueous Proton Transport

J. Phys. Chem. B, Vol. 112, No. 2, 2008 471

TABLE 1: Force Field Parameters of the MS-EVB3 and qMS-EVB3 Modelsa DOH aOH R0OH ka a0 OOw σOOw HOw σHOw +

qHO3O

+

qHH3O B b b′ d0OO

88.96 89.94 2.1 2.14 1.0 1.0 77.4868 77.4868 111.7269 111.7269 0.1238 0.1238 3.142 3.142 0.0025115 0.0031705 1.582746 1.582746 -0.5 -0.5 0.5 0.5 11.2600138 11.6552448 1.1 1.0949993 2.12 2.12 2.4 2.4

kcal/mol

C

Å-1

c

Å

d0OH

kcal/mol/rad2 Vijconst degree

qex O

kcal/mol

qex H

Å

b qex H*

kcal/mol

g

Å

P

e

k

e

DOO

kcal/mol

b

Å-1

R0OO

Å-2

P′

Å

a r0OO

4.5715736 5.1195154 2.1 2.0936705 1.0 1.0 -23.1871874 -23.4683913 -0.0895456 -0.0908305 0.0252683 0.0322038 0.0780180 0.0528458 1.8302895 1.7940356 0.2327260 0.1945086 9.562153 12.2003098 2.94 2.94 6.0179066 5.4823731 3.1 3.12 10.8831327 8.4529834 10.0380922 14.5873449 1.8136426 1.9525044

kcal/mol Å-1 Å kcal/mol e e e Å-2

Å-2 Å Å-1 Å Å-1 Å-1 Å

a For each parameter in the table, upper and lower values are for the MS-EVB3 and qMS-EVB3 models, respectively. For the parameters of the SPC/Fw and qSPC/Fw water potentials, see refs 55 and 57, respectively. b H* is the central hydrogen atom in the Zundel structure.

relative errors of the binding energies and geometries of the three clusters (as marked by an asterisk in Table 2 later) and the root-mean-square deviation of the PES were calculated and minimized simultaneously. The MS-EVB3 model requires 29 independent parameters, which makes the full nonlinear optimization a challenging task. To make the optimization more tractable, a subset of 16 parameters for which we could give reliable guesses was chosen for manual adjustment. The simplex method was then used to optimize the remaining 13 parameters. Tens of initial values were examined, and the parameter sets generated by the optimization procedure were further examined for dynamic and structural properties of the excess proton in bulk water by performing an MS-EVB MD simulation of a box of 216 water molecules with periodic boundary conditions. The parameter sets summarized in Table 1 yielded the best overall performance. 2.4. Computational Details. The ab initio calculations were performed using the Gaussian 03 package.58 Unless otherwise specified, the MS-EVB3 simulations of an excess proton in the bulk water phase were carried out under the following conditions: The system consisted of a cubic box of 216 SPC/Fw55 (for classical simulations) or 216 qSPC/Fw57 (for quantum simulations) water molecules plus an excess proton. The latter was added into the water system by replacing a randomly selected water molecule with a hydronium cation. The system was simulated with a constant box size that provided a density of 1.0 g/cm3 and at a temperature of 298.15 K. In the classical simulations to calculate statistical properties, the temperature was maintained by a Nose´-Hoover thermostat59 with a 0.5 ps

relaxation time, and the equations of motion were integrated using the leapfrog algorithm with a time step of 1 fs, and a total of 24 ns long trajectory was obtained from four independent 6-ns long simulations. For the calculations of dynamical properties, the MS-EVB3 simulations were carried out under the same conditions as described above, except without the Nose´-Hoover thermostat, that is, in the constant NVE ensemble. It should be noted that, although the mean temperature of each NVE trajectory deviated by ∼4 K on average from the target temperature at which the system had already been equilibrated, there was no appreciable persistent heating or cooling during the MS-EVB3 simulations without the thermostat. The quantum-dynamical simulations were performed using the centroid molecular dynamics (CMD) method,60-66 which can account for various quantum features in a many-body system. A normal-mode representation of path-integrals was employed, with the centroid position of each atom defined as the zero-frequency mode of the corresponding quasi-particle chain.67 The propagation of the centroid degrees of freedom was carried out in the constant NVE ensemble with the centroid force computed “on the fly” according to the adiabatic time scale separation scheme,63 and with Nose´-Hoover chains of four thermostats68 attached to each nonzero frequency normalmode. A discretization of the quantum partition function with P ) 24 quasi-particles was used, which was shown to be sufficient for obtaining converged results.54,57,69 An adiabaticity parameter γ ) 0.06 and a time step ∆t ) 0.01 fs ensured a sufficiently large separation in time between the motion of the centroid and the other nonzero frequency normal-modes, as well as a proper integration of the latter.63,70 The quantum results were computed by averaging over 10 independent trajectories of 100 ps each. The temperature in each trajectory was T ) 298 ( 4 K. Running multiple relatively short trajectories was preferred over performing only one longer simulation in order to guarantee a better canonical average over the initial centroid variables. Periodic boundary conditions were employed in both classical and quantum simulations. The LJ interactions were treated using the spherical cutoff method (cutoff radius ) 9 Å), and corrections due to the missing long-range dispersion interactions were made for the energy using the standard method.71 The electrostatic interactions were treated using the Ewald sum method71 with an error tolerance of 10-6. Both quantum and classical simulations were performed using the DL_EVB program,39 which was derived from DL_POLY version 2.13 software.72,73 3. Results and Discussion At the outset of this section, it should be noted that two MSEVB models, denoted MS-EVB3 and qMS-EVB3, were developed in this work for classical and quantum (centroid) molecular simulations, respectively. The underlying water model is SPC/ Fw (for MS-EVB3) or qSPC/Fw57 (for qMS-EVB3). The major content of this section is focused on the classical model, but the results from the quantum model will also be presented and discussed wherever quantum effects are significant and interesting. 3.1. Energies and geometries of protonated water clusters. Table 2 summarizes the results of the binding energies and geometries of different protonated water clusters, whose configurations are depicted in Figure 2 for both the MS-EVB3 and the qMS-EVB3 models. The binding energies were calculated as ∆E ) E(H3O+‚n(H2O)) - E(H3O+) - n × E(H2O), where E(H3O+‚n(H2O)) is the energy of the optimized cluster, while

472 J. Phys. Chem. B, Vol. 112, No. 2, 2008

Wu et al.

TABLE 2: Binding Energies (∆E) and Geometries of Different Protonated Water Clusters (See Figure 2 for the Corresponding Optimized Structures) clustera 2wat*

3wat 4wat1*

4wat2* 5wat1 5wat2 5wat3 5wat4 5wat5

propertyb R(OO) R(OH) ∆E R(OO) R(OH) ∆E R(OO) R(OH) ∆E R(OO) R′(OO) R(OH) ∆E ∆E ∆E ∆E ∆E ∆E

(Å) (Å) (kcal/mol) (Å) (Å) (kcal/mol) (Å) (Å) (kcal/mol) (Å) (Å) (Å) (kcal/mol) (kcal/mol) (kcal/mol) (kcal/mol) (kcal/mol) (kcal/mol)

MS-EVB3

qMS-EVB3

MS-EVB2c

MP2/aug-cc-pVDZd

2.39 1.20 -33.06 2.50 1.07 -56.22 2.51 1.07 -77.84 2.37 2.59 1.19 -72.77 -91.27 -92.07 -91.92 -92.45 -88.22

2.39 1.20 -32.45 2.50 1.07 -56.91 2.53 1.05 -79.64 2.37 2.59 1.19 -73.88 -93.93 -94.36 -94.35 -88.28 -89.76

2.40 1.20 -33.93 2.52 1.06 -58.43 2.56 1.04 -82.08 2.37 2.62 1.19 -75.42 -93.91 -96.48 -96.42 -88.12 -89.64

2.39 1.20 -34.4 2.50 1.03 -57.6 2.56 1.01 -77.4 2.38 2.61 1.19 -73.9 -92.60 -91.72 -91.74 -88.47 -88.60

a The clusters marked with the asterisk (*) were the training set of the parameter optimization. b R(OO), R’(OO), and R(OH) are as defined in Figure 1. c Results for 2wat, 3wat, and 4wat1 are from ref 23. d Results for 2wat, 3wat, 4wat1, and 4wat2 are from ref 88; for 5wat1, 5wat2, 5wat3, 5wat4, and 5wat5, the results are from ref 89.

E(H3O+) and E(H2O) are the energies of the isolated hydronium ion and water molecule at their optimized monomer geometries, respectively. Compared to the MS-EVB2 model, both new models significantly improve the agreement with the MP2 binding energy for the 4wat1, 5wat1, 5wat2, and 5wat3 clusters. The other properties of the new models are generally as good as those of the MS-EVB2 model. The largest relative deviation from the MP2 results is only 4% (for the binding energy of the 5wat4 cluster in the MS-EVB3 column). 3.2. PESs of Proton Transfer in Small Water Clusters. Because the dynamics of an excess proton in liquid water involves fast proton shuttling between two water molecules on the subpicosecond scale, an accurate description of the PES is essential for a reliable description of both the dynamic and solvation properties of the excess proton. To this end, the new models were explicitly parametrized in this work to reproduce the MP2 PES for proton transfer in a linear tetramer water configuration (shown in Figure 3a), which mimics the single-file water chain structure as commonly observed in biomolecular proton-channel environments. The PES was also checked for proton transfer in a dimer configuration (shown in Figure 3b). In both cases, the PES was calculated by scanning the position10 of the central hydrogen between the two immediate water oxygen atoms, the latter being held at separations (dOO) varying from 2.2 Å to 2.8 Å, while all other atoms are allowed to relax. For both the tetramer and dimer cases, the MS-EVB2 model is in very good agreement with the MP2 results, which is encouraging given that the model was not explicitly parametrized to reproduce the MP2 PESs. However, in comparison to the MS-EVB2 model, the overall curves produced by the MS-EVB3 model are significantly improved. This is especially true for the 0.5 < |q| < 0.6 regions. The most remarkable improvement achieved by the MSEVB3 model is for the tetramer case with dOO ) 2.8 Å: it gives nearly perfect agreement with the MP2 result, whereas the MS-EVB2 model overestimates the barrier by ∼5 kcal/mol. For the dimer case, the agreement of the MS-EVB3 model is again generally good. In the parametrization, only the tetramer PES was included in the training set because the dimer represents a singular case where the Vijex term in eq 10 vanishes.

3.3. Excess Proton Hydration Structures in Bulk Water. As mentioned in the Introduction, two distinct proton hydration structures in bulk water, the Eigen17 and Zundel18 cations, have long been proposed.74,75 The former refers to structures where the hydronium ion is symmetrically solvated by three water molecules, resembling the global-minimum configuration of the protonated water tetramer in the gas phase (4wat1 in Figure 2); the latter refers to structures where a hydrogen is symmetrically shared by two water molecules, resembling the protonated water dimer in the gas phase (2wat in Figure 2). Recent studies have indicated that the proton hydration structure actually undergoes constant and fast interconversion between the Eigen and Zundel forms.22,76 The proton hydration structure can be described in terms of the largest EVB amplitude c21 and second largest EVB amplitude c22, which measure the degree of delocalization of the excess proton. The Eigen cation is in the form of a monohydrate plus three first solvation-shell water molecules and corresponds to c21 ≈ 0.6 and c22 ≈ 0.13, whereas the more delocalized Zundel cation is in the form of a dihydrate and thus corresponds to c21 ≈ c22 ≈ 0.45. Figure 4a shows the probability distribution functions of c21 and c22 computed in the simulations of an excess proton in bulk water. The distributions of c21 and c22 peak at 0.61 and 0.15, respectively, indicating that the Eigen cation is the prevalent hydration species. The prominent shoulders at ∼0.45 and 0.40 for c21 and c22, respectively, indicate a smaller probability of the Zundel cation. The MS-EVB2 and MS-EVB3 models are consistent with each other in their qualitative prediction of the relative populations of Eigen and Zundel cations. Quantitatively, the peak of the c21 curve is left-shifted by 0.05 for the MS-EVB3 model, which is accompanied by a slight right-shift of the peak of the c22 curve and enhanced populations of both shoulders. This difference suggests that the excess proton is more delocalized, and the Zundel cation is slightly more abundant in the MSEVB3 model. The relative free energy of the Zundel and Eigen cations can be easily calculated as a function of c21 using the formula ∆F( c21) ) -kBT ln P(c21) - C, where the constant C reflects an arbitrary choice of the zero value of the free energy. The freeenergy profile (Figure 4b) shows a minimum at c21 ≈ 0.61,

The MS-EVB3 Model for Aqueous Proton Transport

J. Phys. Chem. B, Vol. 112, No. 2, 2008 473

Figure 2. Optimized structures of different protonated water clusters with their labels shown in the figure.

which can be attributed to the prevalent Eigen cation, and a plateau at c21 ≈ 0.45, which can be attributed to the Zundel cation. From c21 ) 0.61 to c21 ) 0.45, the free energy increases monotonically by ∼1 kcal/mol, while, for c21 > 0.65 and c21 < 0.42, the free energy rises abruptly. The prediction by both models is in good agreement with a theoretical estimate (0.6 kcal/mol) by Agmon.19 In the qMS-EVB3 model, it is interesting to note that the population of the most stable solvation state is reduced in favor of the less stable ones, including both the more delocalized and more localized states. This result can be effectively attributed to the quantum delocalization effect for the system nuclei. 3.4. Free Energy Profiles of the Proton-Transfer Reaction. The transition free-energy for proton transfer can be analyzed using c21 - c22 as a coordinate22 because the pivot state changes its identity when c21 - c22 ) 0. The free energy profile (Figure 5) for the MS-EVB3 model shows a maximum of ∼1.0 kcal/ mol, 0.1 kcal/mol lower than that of the MS-EVB2 model. The coincidence of the free-energy maximum for proton transfer with the free energy difference between the Zundel and Eigen cations supports the notion that proton transfer in water is associated

with a hydration structural change between the two cations. However, it must be stressed here that the actual reaction coordinate for the proton hopping in aqueous solution is much more complicated, involving a collective hydrogen bond rearrangement.15 The quantum free energy profile obtained from the CMD simulations displays a shape similar to that of the corresponding classical result. However, the free energy peak is ∼0.2 kcal/ mol lower, and the overall profile is somewhat shallower. These features are a direct manifestation of quantum effects, which derive from the inclusion of the zero-point energy contribution as well as the associated nuclear delocalization. The structural aspects of the proton-transfer process can be further characterized via a two-dimensional free energy profile (Figure 6a). This profile exhibits a butterfly-like pattern with two oval-shaped basins centered at roughly ((0.15, 2.4) and a maximum that is ∼2.0 kcal/mol higher at roughly (0, 2.36). The profile indicates that proton-transfer mainly occurs at a donor-acceptor oxygen distance of dO*-O ) 2.3-2.4 Å, with activation energies lower than 3.0 kcal/mol. It is worth noting that, for dO*-O < 2.4 Å, the free energy barrier for proton

474 J. Phys. Chem. B, Vol. 112, No. 2, 2008

Wu et al.

Figure 4. The probability distributions of (a) the largest (c21) and the second largest (c22) MS-EVB amplitudes, and (b) the free energy profile of the proton solvation structure in bulk water as a function of c21.

Figure 3. PESs for (a) a protonated linear tetramer and (b) a protonated water dimer. The numbers 2.2, 2.4, 2.6, and 2.8 denote different dOO values in units of Å.

transport is still positive because of the solvation effect, even though the transition state in free energy is a minimum on the PES (Figure 3). Similar results were also obtained for the qMSEVB3 model. The two-dimensional free energy profile for the MS-EVB2 model (Figure 6b) is quantitatively similar to that for the MS-EVB3 model, although the free-energy basins shrink slightly. 3.5. Radial Distribution Functions. The solvating water structure of an excess proton can be characterized by interatomic RDFs (Figure 7). The most accurate experimental RDFs to date were determined using the empirical potential structure refinement (EPSR) method from neutron diffraction data.21,77 However, it should be noted that the neutron diffraction data were collected only for highly concentrated HCl solutions (6.15 M; whereas the proton concentration in our simulations is only 0.26 M) and that the EPSR calculations were based on a classical nonshuttling hydronium model. In the following text, we will compare and discuss our RDFs against these experimental results, but the reader should bear in mind that any lack of agreement with those results may not necessarily indicate an inaccuracy in the MS-EVB model.

Figure 5. The free energy profile of the proton-transfer reaction as a function of the difference between the two largest EVB amplitudes (c21 - c22).

Figure 7a shows the RDFs of the water oxygen atoms around the pivot hydronium oxygen atom O*. The first peaks in the MS-EVB curves are at rO*-O ) 2.5 Å and is ∼4 in magnitude, in good agreement with the experimental result, although the latter is slightly sharper, higher, and shifted to the left. The higher HCl concentration in the experiments results in the replacement of one first-shell water molecule with a Cl- ion, which may in turn cause the first peak to become sharper.21 The classical hydronium model employed in the EPSR calcula-

The MS-EVB3 Model for Aqueous Proton Transport

Figure 6. The free energy landscape as a function of dO*-O and q (see text) for the (a) MS-EVB3 and (b) MS-EVB2 models. The values of grayscale bars are in units of kcal/mol.

tions might also result in a higher and left-shifted first peak of the experimental curve since such a model localizes the +1 charge on a single molecule rather than distributing it into a number of water molecules as in the MS-EVB simulations. Following the first peak in Figure 7a is a depression at rO*-O ) ∼2.9 Å. This depression is so wide and deep it appears that the first-shell water molecules around the hydronium cation are “isolated” from those that remained. The MS-EVB3 curve is in good agreement with respect to the location of the bottom of the depression, while the MS-EVB2 RDF almost reaches zero in this region. The behavior of the MS-EVB2 model is an artifact that will be analyzed in more detail in the next subsection. The experimental curve shows a slightly shallower and left-shifted depression, which may again be attributed to its classical hydronium model. The MS-EVB and experimental curves both show a pronounced shoulder at ∼3.2 Å in Figure 7a. As will be illustrated in spatial distribution functions (SDFs; see below), this corresponds to the water distribution in the lone pair region of the pivot hydronium molecule. The second major peak in the experimental and MS-EVB curves is located at ∼4.5 Å, which is roughly at the same position as it is in the RDF for pure liquid water. The water potential used by the MS-EVB3 model gives nearly perfect agreement with the experimental RDF for pure liquid water for the second peak,55 but its lack of explicit electronic polarizability may underestimate the ionwater interactions, which could result in a less-structured second solvation shell. For the same reason, the MS-EVB2 model and the EPSR calculation might also underestimate the second-shell water structure. Figure 7b shows the RDFs of the water hydrogen atoms around the pivot hydronium oxygen atom O*. The first peak is at rO*-H ) ∼3.1 Å and is composed of 10-12 water hydrogen

J. Phys. Chem. B, Vol. 112, No. 2, 2008 475 atoms, as indicated by the cumulative-number curves. As will be shown by the analysis in section 3.6, these hydrogen atoms are from both the proton-acceptor water molecules and those at the “lone pair” region of the pivot hydronium (nonacceptor water molecules). Following the first peak is a dip at rO*-H ∼ 3.6 Å. Such a dip is, however, not present in the RDF of water hydrogen atoms around a water oxygen atom (e.g., see ref 55), so it may be ascribed to the delocalization of the excess proton, which “hydroniumizes” the first-shell water molecules by distributing ∼10% of the protonic charge onto each of them. This will in turn cause a more ordered structure compared to that of pure water. The experimental curve is much less structured beyond the first peak, and, furthermore, the first dip is replaced with a ramp as seen in pure water’s RDF. However, these differences relative to the simulation results may be ascribed to the classical hydronium model used in the EPSR calculations: the classical model is unable to delocalize the protonic charge to the first-shell water molecules, and thus the features for the second shell in the experimental RDF look similar to those for pure water. Figure 7c shows the RDF of the water hydrogen atoms around a hydronium hydrogen atom H*. The first peak at rH*-H ) ∼2.1 Å corresponds to the two hydrogen atoms of the acceptor water molecule. The smaller magnitude of this peak in the experimental curve can be ascribed to the replacement of one acceptor water molecule by a Cl- ion. Beyond the first peak, the experimental curve is less structured than the other simulation results, which may again be due to the use of the classical hydronium model in the EPSR model that cannot delocalize the protonic charge to the second-shell water molecules, and thus it underestimates the hydrogen-bond interactions between the second- and third-shell water molecules. The MS-EVB2 curve shows a small bump at 4.4 Å, which could be an artifact of the potential since it is not observed in all of the other curves. Last, Figure 7d shows the RDF of the water oxygen atoms around a hydronium hydrogen atom H*. The first peak, corresponding to the acceptor water oxygen atom, is located at rH*-O ∼ 1.45 Å and is separated from the rest of the curve by a wide depression, implying a strong hydrogen bond between the pivot hydronium and its first-shell water molecules. The RDF results described above demonstrate that generally good agreement with the neutron diffraction results can be achieved with the MS-EVB3 model in the description of the proton solvation structure. The differences between the MSEVB3 model and the experimental curves may be caused by the use of the classical hydronium model in the EPSR calculations and the high acid concentration in the experiments. Compared with the MS-EVB2 model, the new model MS-EVB3 is substantially improved, most notably in the O*-O RDF(rO*-O), where an artificial “hole” at 2.9 Å has been eliminated. 3.6. Spatial Distributions, Orientations, and Hydrogen Bonding of the Solvating Water Molecules. To further elucidate the solvation structures, SDFs were calculated according to the method described by Svishchev and Kusalik.78 The SDFs were calculated with respect to the pivot hydronium with the principal frame as defined in Figure 8. In a way different from that in a previous study,31 the current choice of frame distinguishes the pivot hydrogen (the one shared by the two most significant EVB states, orange-colored in Figure 8) from the other hydronium hydrogen atoms in the SDF results here. The SDF of water oxygen atoms within rO*-O < 3.5 Å around the pivot hydronium is illustrated in Figure 9a for the MS-EVB3 model and Figure 9b for the MS-EVB2 model. Each panel

476 J. Phys. Chem. B, Vol. 112, No. 2, 2008

Wu et al.

Figure 7. Interatomic RDFs (solid and dashed lines) and cumulative number functions (dot-dashed lines with the same color codes as the RDFs). The “Expt” curves are the experimental results from ref 21. The inset of panel (a) is a magnified view of the first peak of RDF(rO*-O).

Figure 8. The principal frame for the SDFs illustrated in Figures 9 and 10. The origin is at the pivot hydronium oxygen atom. The z-axis is perpendicular to the plane determined by the three hydrogen atoms. The x-axis lies in the plane determined by the z-axis and the vector from the origin to the pivot hydrogen (in orange color). The point P′ is the projection of the point P on the xy-plane, illustrating the definitions of the θ and φ angles.

shows side views and top views of the hydronium. The hydrogen-bond-acceptor water molecules constitute the three main lobes in this function. Clearly, they are the main contributions to the first RDF(rO*-O) peak (Figure 7a). For the range of 2.75 Å < rO*-O < 3.05 Å, the contribution is mainly from the molecules from the lone pair side. The density of solvating water molecules from this side is much larger for the MS-EVB3 model than for MS-EVB2, which evidently explains the hole on MS-EVB2’s RDF(rO*-O). The underestimation of water density by the MS-EVB2 model in this region is a small

artifact that again originates from the diagonal oxygen-oxygen repulsive term (eq 2 in ref 23), which was introduced to correct an overbinding between the hydronium and acceptor molecules. However, this term also overestimates the repulsions between the hydronium and non-acceptor molecules at the lone pair region. This interaction is corrected in the present model by multiplying by a term that is a function of rH*-O (eq 7) so that the repulsion is smoothly damped to zero when rH*-O is longer than the hydrogen bond length. It should be noted that the distribution of the above-mentioned solvating water molecules in the hydronium lone pair region is quite smeared, covering a spherical surface ∼2 Å above the three hydrogen-bond-acceptor water molecules and also trailing down somewhat into the space between these acceptors. A small portion of the relative density is also distributed beneath the hydronium. These features together contribute to the O*-O RDF shoulder at rO*-O ) ∼3.2 Å in Figure 7a. Figure 10a depicts the average orientations of the solvating water molecules with respect to the pivot hydronium for the MS-EVB3 model. It is obtained by calculating the average angle B vector (from a water oxygen atom to (R) formed by the OO* the pivot hydronium oxygen atom) and the water dipole vector for spatial regions with SDF(rO*-O) > 3. For the three hydrogenbond-acceptor water molecules, this angle is mainly populated between 130° and 160°. Such a nonlinear orientation is favorable partly because it allows the water oxygen atom of the acceptor to form another hydrogen bond. The three hydrogen-bondacceptor water hydrogen atoms are mostly directed downward (data not shown), which forms a local field that favors a negatively charged atomic group coming beneath the hydronium, as observed in Figure 9. The population between 120° and 130° is slightly enhanced for the hydrogen-bond-acceptor water of the pivot hydrogen, which is likely because these acceptor water

The MS-EVB3 Model for Aqueous Proton Transport

Figure 9. SDF of the solvating water oxygen atoms around the pivot hydronium molecule for the (a) MS-EVB3 and (b) MS-EVB2 models. The distributions with SDF(rO*-O) > 3 are shown.

molecules have some hydronium character so that they somewhat favor forming a trigonal pyramid geometry such that a smaller value of R results. Figure 10b shows the spatial distribution of the β angle (formed by the B OO* vector and the one B O B H bond vector that gives the smaller angle value) for water molecules within 2.75 Å< rO*-O < 3.5 Å. The molecules on the top are mainly orientated with 90°< β