VOLUME 102, NUMBER 29, JULY 16, 1998
© Copyright 1998 by the American Chemical Society
LETTERS Multistate Empirical Valence Bond Model for Proton Transport in Water Udo W. Schmitt and Gregory A. Voth* Department of Chemistry and Henry Eyring Center for Theoretical Chemistry, UniVersity of Utah, Salt Lake City, Utah 84112 ReceiVed: April 8, 1998; In Final Form: May 22, 1998
A multistate empirical valence bond (MS-EVB) model for describing proton transport in aqueous systems is presented. In this approach the electrostatic interaction of the solvent water molecules with an exchange charge distribution is explicitly included in the off-diagonal elements of the EVB Hamiltonian. The MSEVB model is parametrized to reproduce geometrical and energetic quantities of selected H3O+‚(H2O)N clusters. The obtained geometries, formation energies, and energy barriers are in excellent agreement with results from high-level ab initio calculations. It is furthermore applied in a classical molecular dynamics simulation of condensed-phase water with an excess proton in order to estimate the proton-transfer rate.
1. Introduction Proton-transfer (PT) reactions represent a unique class of processes of paramount importance in many different fields in physics and chemistry and for that reason have achieved extensive experimental and theoretical attention over several decades.1-5 From the theoretical point of view, the accurate modeling of PT remains a challenge for mainly two reasons. First, owing to the small mass of the transferring particle, quantum effects such as tunneling or zero-point energy can have a major influence on the microscopic mechanism and calculated observables. The second challenge in the modeling of PT is to find a reliable and accurate description of the potential energy surface (PES) describing the ongoing formation and breaking of the participating hydrogen bonds. Since in PT reactions the differences in energy, e.g., reaction energy or barrier heights, are only about a few kcal/mol and the rate constants are strongly sensitive to it, the underlying PES needs to be represented with sufficient precision in order to calculate various quantities reliably. One possible route to model the reactive PES is given by the empirical valence bond (EVB) model,6 which has been successfully used in numerical studies of various charge-transfer
reactions in solution or enzymes,6 proton transfer in gas-7 and condensed-phase systems,8 and to construct the PES in reactive scattering calculations.9 Lobaugh and Voth8 first introduced a two-state EVB model in order to describe the ground-state PES for the migration of a proton between two water molecules in the liquid state. While this approach was able to model the dynamical delocalization of the proton between two water molecules and, e.g., reproduced the IR absorption spectrum characteristics for strong acids in aqueous solutions, it does not allow for explicit proton hopping between several water molecules. There has recently been an effort10 to improve upon the Lobaugh-Voth two-state EVB model through a higher level of ab initio calculations as the basis of the empirical modeling. Since these calculations do not include the effects of solute polarization, it is unclear whether they improve upon the original model. Very recently, Vuilleumier and Borgis extended the Lobaugh-Voth approach by introducing a nonadditive valencebond force field11,12 to model the migration of an excess proton in water clusters and bulk water, which goes beyond a twostate EVB description. The model has been used to calculate vibrational spectroscopic properties of small protonated water clusters and the transfer of an excess proton in bulk-phase water,
S1089-5647(98)01813-6 CCC: $15.00 © 1998 American Chemical Society Published on Web 06/29/1998
5548 J. Phys. Chem. B, Vol. 102, No. 29, 1998
Letters
but it differs in an important way from the model described herein as will be described below. In this Letter we present a methodology to describe the PES for PT reactions based on a multistate empirical valence bond (MS-EVB) picture, which allows one to model the migration of a single proton in aqueous environments, e.g., protonated water clusters, proton wires, and aqueous bulk-phase solutions. This modeling is accomplished in a numerically efficient but accurate fashion, and so it enables the dynamical treatment of large aqueous systems and provides a basis for including quantum effects stemming from the nuclei. This approach can be regarded as a multistate extension of the two-state Lobaugh and Voth8 EVB model in order to allow for unrestricted PT along three-dimensional hydrogen-bonded networks present in aqueous systems. 2. Multistate EVB Approach The EVB methodology, pioneered by Mulliken,13 allows in general the construction of potential energy surfaces (PES) for species undergoing chemical reactions in an accurate and numerically efficient fashion. It is based on the assumption that the ground-state PES of a reactive chemical system can be obtained from the ground-state energy of an effective EVB Hamiltonian
H ˆ (Q) )
|i〉hij(Q)〈j| ∑ i,j
(1)
where the matrix elements hij(Q) are a function of the set of nuclear degrees of freedom Q of the molecular system under consideration only. The diagonal elements hii(Q) are mostly assumed to have a simple functional dependence on the nuclear degrees of freedom and so can be represented for example by established molecular mechanics force fields.14,15 The functional dependence of the off-diagonal elements of Q has to be chosen to reproduce the PES in the reaction region where the approximations employed in the standard force fields break down. The PES is obtained by solving the eigenvalue problem
H(Q)c ) E(Q)c
(2)
and by calculating the expectation value of H(Q)
E(Q) )
c0i c0j hij(Q) ∑∑ i,j
(3)
with respect to its ground-state vector c0. Owing to its underlying simplicity and resulting numerical efficiency, the EVB methodology has been used to model various chemical reactions in the gas and condensed phase, and it seems to be a valuable tool for describing reactive PES in condensed systems, where the large number of atoms involved in the chemically reactive event prohibit a high-level electronic structure calculation. Lobaugh and Voth8 employed the EVB method to model the ground-state PES for the migration of a proton between two water molecules by using a linear combination of two EVB states, where the excess proton is bound to one distinct water molecule in each state, respectively. The diagonal elements were expressed as the sum of intramolecular contributions of the distinct molecular species (H3O+ and H2O) plus an intermolecular nonbonded interaction term. The off-diagonal elements were taken to be a function of a specified subset of the nuclear coordinates, i.e., the distance between the two oxygens in the H5O+ 2 dimer.
Figure 1. Schematic picture of the four valence-bond states necessary to describe the H9O+ 4 complex in an MS-EVB approach.
To properly treat proton migration along a three-dimensional hydrogen-bond network in aqueous solutions, it is necessary to extend the aforementioned two-level EVB approach to a generalized MS-EVB model. This results from the equivalence of all three hydrogens atoms bonded to the hydronium oxygen, which should be allowed to migrate to the water molecules constituting the first solvation shell. For example, the states necessary to correctly treat the proton migration in the Eigen cation H9O+ 4 are shown in Figure 1. While the central hydronium ion (state |3〉) can act as a proton donor to three H2O molecules, there are in total four states to be included. In general, a multistate EVB description results in an N × N EVB Hamiltonian (eq 2), where again the expectation value of H with respect to the eigenvector corresponding with the lowest eigenvalue of H gives the electronic ground-state PES. While the MS-EVB approach presented in this Letter also extends the Lobaugh-Voth two-state EVB description through a MS-EVB model, the implementation and parametrization of the off-diagonal matrix elements hij(Q) entering into the MSEVB Hamiltonian differs considerably and importantly from the approach presented in ref 11 in that they depend explicitly on the solvent configurations. To our best knowledge, the applications of the EVB approach so far have neglected this dependence, which is a natural consequence of the quantummechanical nature of the VB states (see below). As pointed out by Kim and Hynes,16 a molecular system described by a VB electronic Hamiltonian that is immersed in a polarizable solvent interacts with the solvent molecules in two ways. The first is due to the charge distribution Fii(QS) corresponding to the diagonal VB states |i〉 of the solute, i.e.
Fii(QS) ) -eΨ*i (QS)‚Ψi(QS)
(4)
where QS is the set of nuclear degrees of freedom of the solute complex only and
Ψi(QS) ) 〈QS|i〉
(5)
is the coordinate representation of the VB state |i〉. The second interaction comes from an exchange (or transition) charge distribution Fij(QS), i.e.,
Letters
J. Phys. Chem. B, Vol. 102, No. 29, 1998 5549
Fij(QS) ) -eΨ*i (QS)‚Ψj(QS)
(6)
resulting from the quantum-mechanical wave character of the electronic degrees of freedom of the solute. When one treats the electronic degrees of freedom as being described by a purely classical charge distribution, the transition charge distribution Fij and the corresponding exchange electric field
∫
E(Q)S ) ∇ dQ′S
Fij(Q′S) |Q′S - QS|
(7)
are not present,16 and this is not consistent with the quantum mechanics of the system. In the MS-EVB model presented in this Letter, the interaction of the transition charge distribution Fij with the solvent is explicitly taken into account in describing the off-diagonal elements in the EVB Hamiltonian (eq 2). This adds additional degrees of freedom to the parameter set used in modeling the functional form of the individual matrix elements hij in eq 2. The additional terms in the Coulomb interaction between a solute molecule k and a solvent molecule l
VklC )
∫∫dQkdQl
Fk(Qk)Fl(Ql) |Qk - Ql|
TABLE 1: Parameters Used in the MS-EVB Model DOH aOH R0OH kR R0 O′ σO′ qO′ qH′ B b d0OO R
266.3 kcal mol-1 1.285 Å-1 0.98 Å 73.27 kcal mol-1 rad-2 116.0 deg 0.155 kcal mol-1 3.164 Å -0.5 e +0.5 e 2.295 kcal mol-1 3.50 Å-1 2.50 Å 15.0 Å
r0OO β R0OO DOO P k γ ζ ω0 Vijconst qexch O qexch H
1.92 Å 4.50 Å-1 3.14 Å 2.875 Å 0.27 11.5 Å-2 1.85 Å-2 1.5 rad-2 174.3° -32.925 kcal mol-1 -0.1334 e 0.05336 e
(8)
can perhaps be seen more explicitly from the expectation value of the charge distribution operator with respect to the MS-EVB ground state:
Fk(Qk) ) 〈Fˆ k(Qk)〉 )
cicj Fkij(Qk) ∑∑i c2i Fkii(Qk) + 2 ∑∑ i
