Ab Initio Molecular Dynamics Simulation of the Structure and Proton

Robin L. Hayes , Stephen J. Paddison , and Mark E. Tuckerman. The Journal of ..... Ira Kurtz , Jeffrey Kraut , Vahram Ornekian , Minhtri K. Nguyen. Am...
0 downloads 0 Views 297KB Size
Download PDF
3712

J. Phys. Chem. B 2006, 110, 3712-3720

Ab Initio Molecular Dynamics Simulation of the Structure and Proton Transport Dynamics of Methanol-Water Solutions† Joseph A. Morrone Department of Chemistry, Princeton UniVersity, Princeton, New Jersey 08544

Kiryn E. Haslinger‡ Department of Chemistry, New York UniVersity, New York, New York 10003

Mark E. Tuckerman* Department of Chemistry and Courant Institute of Mathematical Sciences, New York UniVersity, New York, New York 10003 ReceiVed: September 22, 2005; In Final Form: December 21, 2005

Ab initio molecular dynamics simulations are employed to study the structural and proton transport properties of methanol-water mixtures. Structural characteristics analyzed at two different methanol mole fractions (XM ) 0.25 and XM ) 0.5) reveal enhanced structuring of water as the methanol mole fraction increases in agreement with recent neutron diffraction experiments. The simulations reveal the existence of separate hydrogen-bonded water and methanol networks, also in agreement with the neutron diffraction data. The addition of a single proton to the XM ) 0.5 mixture leads to an anomalous structural or Grotthuss-type diffusion mechanism of the charge defect in which water-to-water, methanol-to-water, and water-to-methanol proton transfer reactions play the dominant role with methanol-to-methanol transfers being much less significant. Unlike in bulk water, where coordination number fluctuations drive the proton transport process, suppression of the coordination number of waters in the first solvation shell of the defect diminish the importance of coordination number fluctuations as a driving force in the structural diffusion process. The charge defect is found to reside preferentially at the interface between water and methanol networks. The length of the ab initio molecular dynamics run (∼120 ps), allowed the diffusion constant of the charge defect to be computed, yielding a value of D ) 4.2 × 10-5 cm2/s when deuterium masses are assigned to all protons in the system. The relation of this value to excess proton diffusion in bulk water is discussed. Finally, a kinetic theory is introduced to identify the relevant time scales in the proton transfer/transport process.

1. Introduction The structural and dynamical properties of methanol-water mixtures are remarkably rich and complex.1-19 The most recent neutron diffraction experiments16 performed on methanol-water solutions at varying methanol mole fractions demonstrate that in the region XM∈[0.27, 0.54], both methanol and water form separate percolating hydrogen-bonded networks. The amphiphilic character of methanol causes the molecules to organize into extended structures in which the hydrophobic groups cluster together away from aqueous regions. In the aforementioned mole fraction regime, there is, then, sufficient methanol to give rise to the separate bi-percolating liquid character of the solution. This phenomenon has also been predicted in theoretical studies7,9,17 and thermodynamic measurements.15 This emergent picture of methanol-water mixtures has challenged long-held assumptions that enhanced structuring of the molecules in alcohol/water mixtures exist in the form of ice-like water cages surrounding the hydrophobic regions of the alcohols.20 It also †

Part of the special issue “Michael L. Klein Festschrift”. * Corresponding author. ‡ Present address: Cold Spring Harbor Laboratory P.O. Box 100, 1 Bungtown Rd. Cold Spring Harbor, NY 11724.

challenged a more recently introduced scenario interpreted from X-ray scattering data21 in which water molecules serve to stabilize ring-like structures of 6 or 8 methanols.21 The picture that emerges from the most recent investigations7,9,16,17 indicate the dramatic extent to which our understanding of the complex microscopic behavior of methanol-water mixtures has evolved over the decades. Given the unusual structural properties of methanol-water solutions, it is then natural to ask how topological defects in the form of excess protons are transported through the solution. This question is particularly relevant given that methanol crossover in methanol-based polymer-electrolyte membrane (or “PEM”) fuel cells leads to a mixture of methanol and water within the membrane pores.22-24 Moreover, the study of how proton transport occurs through a partially aqueous/partially amphiphilic environment such as a methanol-water mixture leads to important insights into how protons are transported in biological systems such as in regions near lipid membranes.25-28 In a methanol-water solution, proton transport occurs at anomalously high rates. If the experimental picture of separate percolating liquids is valid, then two competing mechanisms might seem to drive proton transport: one pertinent to water and the other to methanol. Let us recall that transport of excess

10.1021/jp0554036 CCC: $33.50 © 2006 American Chemical Society Published on Web 02/04/2006

Molecular Dynamics Simulation of Methanol-Water Solutions protons in bulk water is driven by a structural diffusion process in which specific solvent fluctuations lead to an interconversion between the 3-fold coordinated “Eigen complex”, H3O+‚(H2O)3, and a shared-proton “Zundel cation” [H2O‚‚‚H‚‚‚OH2]+, via proton transfer.29-38 Specifically, a hydrogen-bond breaking event between first and second solvation shell members of H3O+ gives rise to a change in the coordination number of a first solvation shell water from 4 to 3.30,31,35,37,39 In this new coordination state, the water molecule is coordinated like hydronium and is thus “prepared” to accept the transferring proton from the H3O+ ion. We call this notion the presolVation concept.38,39 The time required for this coordination number change appears to be rate-limiting, although other possible scenarios have been suggested recently.40 Not unexpectedly, the aqueous phase mechanistic picture contrasts rather sharply with that in liquid methanol. Experiments performed on small protonated methanol clusters41 suggest that second solvation shell fluctuations similar to those that govern proton transport in water also drive proton transport in methanol, while ab initio molecular dynamics studies42 in the bulk phase suggest that excess protons migrate randomly along longer quasi-linear structures characteristic of the pure liquid43,44 but that the diffusion of these chains is driven by outer solvation shell fluctuations. However, since MeOH and MeOH2+ have the same coordination number, outer solvation shell fluctuations are not needed to effect a reduction in the coordination number as occurs in water.42 The primary defect structure formed by an excess proton in liquid methanol appears to be a relatively long quasi-linear chain along which the center of charge can migrate in an essentially random manner.42 If the excess proton localizes on one end of the chain, then new hydrogen bonds can form at this end while hydrogen-bond breaking can occur at the opposite end, releasing methanols into the bulk (see Figure 10 of ref 42). Both events are driven by local solvent fluctuations. In this paper, ab initio molecular dynamics (AIMD) is employed to investigate the structure of methanol-water mixtures at different methanol mole fractions and to study the mechanism of excess proton transport in a mixture at a mole fraction of 0.5. In the latter case, an AIMD simulation of over 100 ps in length is used to enhance the statistical sampling. The structure of the pure mixture at XM ) 0.25 and XM ) 0.5 is found to agree well with the experimentally reported structure.16 While the experiment reports a high probability for the formation of small water clusters in the mixture, there is a small population of larger clusters that cannot be reproduced in the AIMD simulations due to limitations of system size. In the protonated mixture, four possible proton-transfer reactions can occur:

H3O+ + H2O f H2O + H3O+

(1)

H3O+ + MeOH f MeOH2+ + H2O

(2)

MeOH2+ + H2O f MeOH + H3O+

(3)

MeOH2+ + MeOH f + MeOH + MeOH2+

(4)

The branching ratios for these reaction channels are estimated from the simulations, and it is found that the first three can occur with probabilities of approximately 36%, 32%, 27%, respectively, while the fourth, with a branching ratio of roughly 5%, does not contribute significantly to the proton transport dynamics. More importantly, it is found that the charge defect

J. Phys. Chem. B, Vol. 110, No. 8, 2006 3713 is located preferentially at the boundary region between water and methanol phases, a fact that is supported by the estimated branching ratios. This result can be understood by appealing to recent studies of hydronium near the water-vapor interface,45-50 which suggest that when there is a dielectric mismatch, such as occurs at the liquid-vapor interface, the hydronium is most likely to be found at the surface. Since the boundary region between the methanol and water clusters also constitutes a substantial dielectric mismatch, one might, therefore, expect the hydronium to prefer this boundary region. Finally, a kinetic theory34,39,51-56 is introduced to extract the relevant time scales for proton transport in the liquid mixture. Two distinct time scales are uncovered corresponding to proton rattling within a hydrogen bond and true structural displacement events. This paper is organized as follows: In Section 2, the details of the AIMD protocol are reported. In Section 3, results of the structural properties of pure methanol-water mixtures are reported. In Section 4, results elucidating the microscopic transport mechanism of an excess proton in an equimolar methanol-water mixture are reported. In this section, we also present our kinetic model and analyze the results of the model in terms of characteristic time scales in the system. Finally, conclusions are given in Section 5. 2. Computational Details The studies carried out are ab initio molecular dynamics (AIMD) simulations in which the finite-temperature classical nuclear dynamics of a system are generated using forces obtained from electronic structure calculations performed “on the fly” as the simulation proceeds and are based on the CarParrinello extended-Lagrangian adiabatic dynamics approach.57 The electronic structure is represented within the Kohn-Sham formulation of density functional theory with exchange and correlation approximated using the BLYP58,59 generalizedgradient functional. Despite the recent controversy surrounding the ability of the BLYP functional to describe aqueous systems,60 the latest AIMD simulations suggest that BLYP leads to only a very slight overstructuring of liquid water in the second solvation shell.61,62 Similar studies also suggest that the computed diffusion constant is substantially lower than the experimental value; however, these results pertain only to very short runs and have not been corrected for finite-size effects.63 Therefore, we do not draw any serious conclusions about diffusion within the BLYP description at this point. The fictitious mass parameter in the Car-Parrinello extended Lagrangian is 600 au. The simulations performed here employ a plane-wave basis set in conjunction with Vanderbilt ultrasoft pseudopotentials.64 The latter allow an energy cutoff of 25 Ry to be employed in the truncation of the basis set. This protocol was shown to give a reasonable description of the structure of pure liquid methanol42 and of proton transport dynamics in water.29,30 Although some loss of accuracy is expected within this approach, we believe that it is reliable for the largely qualitative conclusions being drawn in the present study. Simulations of pure methanol-water mixtures were performed at XM ) 0.25 (24 waters and 8 methanols in a periodic box of length 10.69 Å) and XM ) 0.5 (16 waters and 16 methanols in a periodic box of length 11.46 Å) in the canonical or NVT ensemble at 300 K by attaching a Nose´-Hoover chain thermostat65 to each degree of freedom in the system. Prior to this, the systems were equilibrated for 200 ps using the AMBER force field66,67 and then further equilibrated for 6-10 ps using the preceding AIMD protocol. Simulations of the pure systems of approximately 3050 ps were performed. For the protonated system, a mole

3714 J. Phys. Chem. B, Vol. 110, No. 8, 2006

Figure 1. Water-water radial distribution functions (RDFs) as a function of methanol mole fraction for methanol-water mixtures. The two left panels show the computed RDFs, while the right panels show the experimental RDFs from the neutron diffraction data of ref 16.

fraction of XM ) 0.5 was simulated with one excess proton. The system was equilibrated for 18 ps using a Nose´-Hoover chain on each degree of freedom. Following this, the system was allowed to evolve in the microcanonical or NVE ensemble for approximately 102 ps over which time the temperature remained between 280 and 330 K with an average of approximately 310 K. In total, the length of the simulation is 120 ps. A simulation of this length, one of the longest AIMD simulations performed to date on a system of this size, was deemed necessary to give a reasonable statistical sampling of the space of possible reactions as described in (1)-(4) and was only possible within the ultrasoft pseudopotential scheme. To reduce the importance of nuclear quantum effects, all hydrogens are assigned the mass of deuterium. Nevertheless, our nomenclature will be to use “H” instead of “D”. The equations of motion were integrated with a time step of 0.12 fs using the nonorthogonal orbital approach of Hutter, et al.68 All simulations were performed using the CPMD Version 3.5 code.69 3. Structural Properties of Pure Mixtures Evidence for the formation of water clusters in methanolwater mixtures as a function of the methanol mole fraction is presented in Figure 1, which shows the theoretical and

Morrone et al. experimental16 water-water radial distribution functions (RDFs) gOWOW(r) and gOWHW(r) for three methanol mole fractions. The figure shows that theoretical RDFs reproduce the experimental trend very closely, in particular, the increase in the OWOW peak height for XM > 0. The XM ) 0 RDFs are taken from the data of ref 35. The theoretical RDFs for XM > 0 show more structuring in the second solvation shell than do the experimental ones, although they are able to reproduce the slight shift to the right in the location of the second peak. The overstructuring is most likely a consequence of finite-size effects. Although figure Figure 4 of ref 16 suggests that smaller water clusters (roughly 3-10 molecules in the mole fraction range [0.27,0.54]) are most probable (probability > 0.01), there is, nevertheless, a small fraction of larger clusters that cannot be reproduced in the present AIMD simulations because of the small number of water molecules. Thus, the second peak in the present theoretical RDFs is due entirely to inter-cluster contributions, while that in the experimental RDFs should have both intra- and inter-cluster contributions. While the latter could cause the peak location to shift, the former give rise to a much broader peak than is seen in the AIMD simulations. The theoretical RDFs also show a slight shift to the right in the first peak for XM > 0, which is due to the presence of a small shoulder at small r that is a consequence of the use of the ultrasoft pseudopotential approach. The enhanced structuring of water with increasing methanol mole fraction is accompanied by a decrease in the methanol structure, as is illustrated in Figure 2. The figure shows the RDF gOMOM(r) for methanol oxygens as a function of the methanol mole fraction. The results for XM ) 1.0 were taken from ref 42. The figure shows the dramatic decrease in the methanolmethanol peak height as the mole fraction is decreased. Although no such picture appears in ref 16, the authors comment on the fact that the tendency for methanol to form hydrogen-bonded chains, as is the case for the pure liquid,42-44,70 is suppressed as the mole fraction is decreased, which supports the picture contained in the RDFs of Figure 2. To assess the affect of the ultrasoft pseudopotential scheme on the liquid structure, a short (8.0 ps) AIMD simulation of the XM ) 0.5 simulation using the Troullier-Martins norm-conserving pseudopotential scheme71 was carried out using a plane-wave cutoff of 70 Ry. The purpose of this run was to compute, roughly, the RDFs in order to

Figure 2. The computed methanol oxygen-oxygen RDFs as a function of methanol mole fraction. The XM ) 1.00 RDF is taken from ref 42.

Molecular Dynamics Simulation of Methanol-Water Solutions

J. Phys. Chem. B, Vol. 110, No. 8, 2006 3715

Figure 3. Total HH, OH, CC, and OO RDFs from the XM ) 0.5 mixture comparing the norm-conserving (green) and ultrasoft pseudopotentials (red) to each other and to the nautron diffraction data of ref 16 (black) and force field based calculations based on the AMBER force field66,67 (blue).

Figure 5. (a) RDFs and integrated coordination numbers in a pure equimolar methanol-water mixture and in pure water. In the methanolwater mixture, the OW-O RDFs are shown, i.e., water oxygen with respect to either water or methanol oxygens. (b) RDFs and integrated coordination numbers for the first solvation shell neighbor of the defect oxygen O* participating in the most active hydrogen bond in the protonated equimolar methanol-water mixture. Here, the RDF is computed only for those configurations in which the first solvationshell oxygen is a water oxygen. The RDFs are computed for two δ regimes, namely, δ < 0.1 Å and δ > 0.4 Å. Figure 4. Left: snapshot from the classical MD simulation of the XM ) 0.5 methanol-water mixture. Only the water and methanol oxygens are shown with the former shown in red and the latter shown in yellow. Right: same for the AIMD simulation. Note that the AIMD box has been replicated 4 times in each spatial dimension for purposes of visualization.

compare to those from the much longer ultrasoft pseudopotential run. In addition, to examine finite-size effects, a classical MD simulation using the AMBER force field66,67 was performed with the PINY•MD code72 using a system of 1024 methanols and 1024 waters. The simulation was run for 240 ps in the NVT ensemble using a Nose´-Hoover chain65 in each molecule. Figure 3 shows the total HH, OH, CC, and OO RDFs for the ultrasoft and norm-conserving pseudopotentials together with the experimental RDFs16 and force field calculations. Here H refers to water and alcohol hydrogens. The three computed RDFs agree well with each other (and differ somewhat from experiment) for all but the CC RDF, where the norm-conserving pseudopotential is better able to reproduce the experimental result. Differences in the OH and HH RDFs with experiment could be due to neglect of nuclear quantum effects rather than from finite-size effects or the chosen functional, although the latter two certainly do have some effect. Given the complexity of the equimolar methanol-water mixture, the agreement between the two AIMD simulations and their overall agreement with experiment suggests that the qualitative results extracted from the simulation is likely reasonable. The clustering of the water molecules for XM ) 0.5 can be clearly seen in snapshots taken from the AIMD simulations. Such a snapshot is shown in Figure 4 (left panel). In the figure, the simulation cell has been replicated 4 times in each spatial dimension in order to make the clustering visually clearer. As

the figure shows, the AIMD simulations can only reproduce relatively small cluster sizes. A typical snapshot from the forcefield based simulation is shown in Figure 4 (right panel). Here, it is clear that cluster sizes larger than approximately 10 can be captured in the classical simulations appear, but that these are less probable than the smaller clusters that are seen in the AIMD simulations. 4. Transport of an Excess Proton in an Equimolar Mixture 4.1. Mechanistic Picture. As described in Section 1, presolvation effects largely dictate the microscopic details of structural diffusion in hydrogen-bonded liquids. Before applying this concept to proton transfer in methanol-water mixtures however, we must first understand the solvation structure of the protonated defect. The different solvation structures of water, methanol, and methanol-water mixtures lead to different structural diffusion patterns in each case. Figure 5(a) contrasts the radial distribution functions and integrated coordination numbers of a water oxygen to all other oxygens in a 50-50 methanol-water mixture and a pure water system. As is evident from the figure, water in the methanolwater mixture is undercoordinated with respect to pure water. This effect is likely due to the presence of the larger methanol molecules in the water solvation shell, which sterically hinder the formation of four hydrogen bonds. The undercoordination of water in an equimolar methanol-water mixture will play an important role in determining the mechanism for proton transfer, as will be shown below. The structural diffusion process may be described by utilizing a reaction coordinate defined as follows: Let O* denote the oxygen that, at any time step, to which the excess proton is

3716 J. Phys. Chem. B, Vol. 110, No. 8, 2006

Morrone et al.

Figure 6. Histograms showing the breakdown of nearest neighbors of the defect oxygen O* (independent of whether the defect oxygen is a water or methanol oxygen). The two histograms correspond to the number of water and methanol neighbors, respectively. Each column of the histograms is normalized by dividing by its coordination number and the total number of configurations, so that the fractions on the y-axis sum to 1 over the two histograms.

assigned. Clearly, this will be an oxygen with three associated hydrogens for water and two associated hydrogens for methanol. For each oxygen partner Oi of the hydrogen bonds in which O* participates, we compute the displacement coordinate

δi ) |dO*Hi - dOiHi|

(5)

where Hi is the shared proton in the ith hydrogen bond, and dO*Hi and dOiHi are the two oxygen-hydrogen distances in the hydrogen bond. The hydrogen bond with the smallest value of δi is identified as the most active hydrogen bond, i.e., the hydrogen bond through which the next proton transfer event is most likely, and the value of the displacement coordinate is denoted simply δ. In Figure 5(b), the OO radial distribution functions and integrated coordination numbers are shown for the first solvation shell water oxygen participating in the most active hydrogen bond with the defect oxygen O* for both small (δ < 0.1 Å) and large (δ > 0.4 Å) values of the reaction coordinate. These values are considered to represent “active” and “inactive” configurations with respect to proton transfer. It can be seen that the active configurations are only slightly undercoordinated with respect to their inactive counterparts. Recalling the discussion in Section 1, the receiving water must have a coordination number similar to that of hydronium, i.e., 3-fold, to facilitate proton transfer. However, as can be seen from Figure 5(a), water in a methanol-water mixture is often already undercoordinated, thereby decreasing the importance of second solvation shell hydrogen bond breaking in this mechanism. Thus, the reaction can be predicted to be barrierless, and an explicit calculation of the free energy profile along δ (not shown) confirms this prediction. Furthermore, on average, it appears that the water molecules around the defect have slightly smaller coordination numbers than in the nonprotonated mixture. A larger number of methanols around the water would likely cause this further lowering of the coordination number. This observation may provide evidence that the defect lies in an area abundant with methanol, namely at the methanol-water interface. To investigate this hypothesis more thoroughly, we have estimated the branching ratios for the proton-transfer processes

Figure 7. A typical snapshot from the AIMD simulation of an excess proton in a equimolar methanol-water mixture showing a protontransfer event at the interface between water and methanol regions. The color coding is as follows: Oxygen (red), carbon (blue), and hydrogen (white). The yellow spheres represent oxygens involved in a proton-transfer complex. The graphics were generated using the VMD program.79

(1)-(4). This calculation was performed simply by identifying the fraction of configurations for which the most active hydrogen bond corresponded to a water-hydronium pair, a methanolhydronium pair, an MeOH2+-water pair, or an MeOH2+methanol pair. The branching ratios obtained for reactions (1), (2), (3), and (4) are 36%, 32%, 27%, and 5%, respectively. These branching ratios indicate that proton transfer events between water and methanol play an important role in the proton transport kinetics. Thus, although water-water transfers constitute roughly 36% of the proton transfer reactions that can occur, suggesting that the proton can occasionally penetrate into a cluster, the defect appears to be mostly in contact with both methanol and water, thus lying in the interface between the two types of clusters. To quantify this finding further, we compiled a histogram of the number of waters and methanols in the first solvation shell of the defect (both as H3O+ and as MeOH+),

Molecular Dynamics Simulation of Methanol-Water Solutions

J. Phys. Chem. B, Vol. 110, No. 8, 2006 3717

Figure 8. (a) Mean-square displacement (MSD) for the structural defect in a protonated equimolar methanol-water mixture. The computed MSD is shown together with a least-squares fit to the linear regime. (b) Same for water and methanol.

arranged according to coordination number. This is shown in Figure 6. Since MeOH2+ has a coordination number of approximately 2 while H3O+ has a coordination number very close to 3, the histograms clearly indicate that the defect must exist in an environment in which both methanol and water are present. Indeed, on average, roughly sixty-five percent of the neighbors are water and forty-five percent methanol, again demonstrating that the defect is, on the average, in contact simultaneously with both water and methanol. A snapshot showing a typical configuration of a proton-transfer complex at the interface between water and methanol regions is shown in Figure 7. The finding that the defect is likely to be found at the methanol-water interfacial region is, perhaps, not unexpected given the dielectric mismatch between the two phases. Indeed, it has been observed in recent sum-frequency vibrational spectroscopy experiments45-47 as well as in molecular dynamics simulations48,49 that hydronium at the air-water interface, where a dielectric mismatch also exists, is most often located in the interfacial region. 4.2. Diffusion Constant. Taking advantage of the relatively long AIMD run, it was possible to calculate the mean square displacement (MSD) of the structural defect to reasonably good accuracy. The MSD was calculated for a continuous dynamical (NVE) trajectory having a length of 90 ps and is shown in Figure 8(a). It is clear from the plot that the MSD plot has a domain of [0,4.3ps], which is a factor of 21 less than the length of the simulation, indicating that the statistics are reasonably good. The diffusion constant of the structural defect in the methanolwater mixture, calculated via the Einstein relation, is given by one-sixth of the slope of the linear regime of the mean square displacement curve. This value is found to be 4.2 × 10-5 cm2/ s. This number should be contrasted with the experimental value for structural diffusion of D+ in heavy water, i.e., (6.69 × 10-5 cm2/s).31,73,74

At this point, it is difficult to say whether the computed diffusion constant is greater or smaller than the diffusion constant for protons in bulk water. Although it is less than the experimental value, the proton diffusion constant in water obtained in the study of ref 39 is 2.3 × 10-5 cm2/s for the BLYP functional, and even smaller values have been reported elsewhere.75,76 However, these numbers are based on shorter runs and have larger associated error bars. The diffusion constant reported here for the methanol-water mixture, although based on a relatively long run by current standards of AIMD simulations, is somewhat influenced by the use of ultrasoft pseudopotentials and finite-size effects. As noted earlier, protontransfer reactions appear to be barrierless even in the absence of nuclear quantum effects, from which one might expect faster diffusion than in bulk water. At the same time, molecular dynamics simulations of proton diffusion near a phospholipid membrane28 suggest retardation of the proton diffusion as the proton penetrates into the polar headgroup region of the membrane. Similarly, in a methanol-water mixture, given the observed branching ratios and the prediction that the charge defect is preferentially located at the interfacial region between water and methanol, one would expect slower diffusion of the defect in the mixture than in pure water as obtained in the present study. It is interesting to note that the value 4.2 × 10-5 cm2/s for the diffusion constant of D+ in a deuterated methanolwater mixture is comparable to the “pure” Grotthuss component for D+ in heavy water, estimated by subtracting the value 1.86 × 10-5cm2/s for heavy-water self-diffusion from 6.69 × 10-5 cm2/s as suggested by Agmon.31 This suggests that the diffusion in the methanol-water mixture might well be dominated by Grotthuss-type hopping events near the methanol-water interface with the hydrodynamic or vehicular component suppressed. Indeed, this is consistent with idea that the diffusion of water clusters through the mixture is likely to be a slow process. To confirm this, we have computed the diffusion constants for

3718 J. Phys. Chem. B, Vol. 110, No. 8, 2006

Morrone et al.

water and methanol in the system. The mean-square displacement plots are shown in Figure 8(b). The slopes of these curves yield Dwater ) 1.7 × 10-7 cm2/s and Dmethanol ) 1.3 × 10-7 cm2/s, which confirms that the diffusion of the charge defect is dominated by the Grotthuss mechanism. Nevertheless, a more accurate assessment of the extent to which it is retarded in the mixture must await larger system sizes, to give a more realistic distribution of cluster sizes, and longer simulation times. 4.3. Kinetic Theory. In this section, we develop a simple kinetic theory for the diffusion of the charge defect using the technique for population correlation functions.34,39,51-56 Here, the kinetic theory is developed based on the dynamics of what we termed the “most active” hydrogen bond.34,35,38 The definition of the most active hydrogen bond is based on the δ coordinate introduced in Section 4.1. If a particular hydrogen bond is identified as the most active hydrogen bond, then it is denoted by the symbol A. We assume that the kinetics are governed by two dominant processes: (1) an equilibrium involving an interconversion between the most active hydrogen bond and a normal hydrogen bond denoted B; and (2) the conversion of the normal hydrogen bond B to a normal hydrogen bond C that, due to diffusion, cannot be converted back into the most active hydrogen bond. These processes can be expressed as simple reactions:

AhB BfC

(6)

Let kf, kb, and k2 be the rate constants associated with the processes A f B, B f A and B f C, respectively. The rate equations then become

d[A] ) - kf[A](t) + kb[B](t) dt d[B] ) - (k2 + kb)[B](t) + kf[A](t) dt

(7)

We start by assuming that a particular hydrogen bond is identified at t ) 0 as the most active hydrogen bond, which means that the initial conditions to the rate equations are [A](0) ) 1 and [B](0) ) 0. The rate equations can now be solved easily using, for example, Laplace transforms, yielding the solutions:

[A](t) ) 1 [(λ - kf + κ)e-(kf+κ-λ)t/2 + (λ + kf - κ)e-(k1+κ+λ)t/2] 2λ [B](t) )

kf -(kf+κ-λ)t/2 [e - e-(kf+κ+λ)t/2] λ

(8)

where κ ) k2 + kb and

λ ) x(kf - κ)2 + 4kfkb

(9)

If a large ensemble of trajectories were initiated with the aforementioned initial condition, then these solutions can be regarded as predicting the fraction of trajectories for which the particular hydrogen bond initially identified as most active would remain the most active after time t. We note that the solutions for [A](t) and [B](t) decay according to two distinct time scales. There is a fast decay governed by the exponential exp[- (kf + κ + λ)t/2] and a slow decay governed by the exponential exp[- (kf + κ - λ)t/2]. Thus,

by calculating the rate constants in the theory, it should be possible to identify the two time scales with specific microscopic motions that govern the decay. To calculate the rate constants, it is necessary to identify a microscopic quantity that represents the populations [A](t) and [B](t). Following ref 39, we represent [A](t) in terms of a population correlation function

[A](t) )

〈f (0)f (t)〉 〈f〉

(10)

where f (t) is a function defined as follows: f (t) ) 1 if the most active hydrogen bond has the same identity at time t as it had at time t ) 0 and f (t) ) 0 otherwise. Note that if we define a population function F(t) such that F(t) ) 1 if the most active hydrogen bond retains its identity continuously up to time t and 0 otherwise, then a correlation function defined as

C(t) )

〈f (0)F(t)〉 〈f〉

(11)

can be used to calculate the lifetime of the most active hydrogen bond from

τlifetime )

∫0∞ dt C(t)

(12)

since eq 11 excludes simple proton rattling events. The correlation function in eq 10 was also proposed by Vuilleumier and Borgis34 for studying the active hydrogen-bond dynamics of an excess proton in water, although the definition of f(t) employed here differs somewhat from that of Vuilleumier and Borgis. Note that the biexponential form in eq 8 emerges a priori from the assumed form of the kinetics rather than from an a posteriori fit to the correlation functions. In a similar manner, the function [B](t) can be represented in terms of a correlation function

[B](t) )

〈f (0)(1 - f (t))g(t)〉 〈f〉

(13)

where g(t) is defined such that g(t) ) 1 if the most active hydrogen bond has a different identity at time t than it had at time t ) 0 but is situated to return to its original location. The bond is assumed to be situated to return to its original location if (1) the most active hydrogen bond is located at time t in an adjacent hydrogen bond with O* retaining its identity; or (2) the most active hydrogen bond moves to a position such that O* moves to a neighboring site at time t but the original O* is presolvated in such a way as to accept a transferring proton and thus become O* again. For simplicity, the theory is based on “ideal” H3O+ and MeOH2+ coordination numbers, although the analysis presented above suggests small variations from the ideal values. Thus, the definition of g(t) is based on a coordination number of 3 at time t if this oxygen is a water oxygen, and a coordination number should be 2 at time t if it is a methanol oxygen. The resulting correlation functions are fit to the functional forms:

[A](t) ) a0e-a1t + (1 - a0)e-a2t [B](t) )

a0a1+ a2(1 - a0) -a2t [e - e-a1t] a1 - a2

(14)

via a least-squares procedure, and the rate constants are then

Molecular Dynamics Simulation of Methanol-Water Solutions calculated from the fitting parameters a0, a1, and a2 according to

kf ) a1a0 + a2(1 - a0) kb )

(a1 - a2)2a0(1 - a0)

k2 )

a0a1 + a2(1 - a1) a1a2 a1a0 + a2(1 - a0)

(15)

Performing this procedure, we obtain the values kf ) 20.94 ps-1, kb ) 22.58 ps-1 and k2 ) 0.86 ps-1. We first note that the rate constants kf and kb are, within statistical error, essentially equal as expected. Moreover, we note that two distinct time scales emerge from this analysis. These time scales are τf ) 1/kf ≈ 48 fs, τb ) 1/kb ≈44 fs, and τ2 ) 1/k2≈1.2 ps. The faster of these time scales can be identified as that of the “rattling” of the excess proton within a hydrogen bond that causes O* to change its identity temporarily but return to its original location, or it can be identified as the time scale for fluctuations of the hydrogen-bond network that cause the most active hydrogen bond to oscillate between a particular location and an adjacent one. The slower time scale corresponds with the average lifetime of a hydrogen bond and, therefore, represents the time scale for actual local rearrangements of the hydrogen-bond network. These rearrangements cause the most active hydrogen bond to diffuse away from its current location, giving rise to the structural diffusion process. 5. Conclusions Ab initio molecular dynamics simulations have been employed to study the structural and proton transport properties of methanol-water mixtures. Structural characteristics were analyzed at two different methanol mole fractions, specifically XM ) 0.25 and XM ) 0.5, at which values the mixture is observed experimentally to consist of two bi-percolating hydrogen-bonded networks. The present simulation studies reveal enhanced structuring of water as the methanol mole fraction increases in agreement with neutron diffraction experiments. The simulations also reveal the existence separate hydrogen-bonded water and methanol networks, also in agreement with the neutron diffraction data. The addition of a single proton to the XM ) 0.5 mixture leads to an anomalous structural or Grotthuss-type diffusion or mechanism of the charge defect thereby created in which waterto-water, methanol-to-water and water-to-methanol protontransfer reactions play the dominant role, according to the estimated branching ratios for reactions 1-4. Methanol-tomethanol transfers appear to play a negligible role in the proton transport kinetics. In contrast to proton transport in bulk water, where second solvation shell fluctuations drive the proton transport process by reducing the coordination number of a first solvation shell water from 4 to 3 according to the presolvation concept, it is found that suppression of the water coordination number of waters in the first solvation shell of the defect diminish the role of coordination number fluctuations in the equimolar methanol-water mixture. Moreover, the charge defect is found to reside preferentially at the interface between water and methanol networks. Finally, a kinetic theory was introduced in order to identify relevant time scales in the proton transfer/ transport process. The analysis based on this theory leads to two distinct time scales. The faster of these is identified with

J. Phys. Chem. B, Vol. 110, No. 8, 2006 3719 the “rattling” of the excess proton within a hydrogen bond that causes the oxygen center of the defect to change its identity temporarily but return to its original location, or it can be identified with fluctuations of the hydrogen-bond network that cause the most active hydrogen bond to oscillate between a particular location and an adjacent one. The slower time scale corresponds with the average lifetime of a hydrogen bond and, therefore, represents the time scale for actual local rearrangements of the hydrogen-bond network. The picture that emerges from the recent experimental and present theoretical investigations into methanol-water solutions reveals a rich and complex behavior for both the pure and protonated systems. The theoretical studies provide further validation of the presolvation concept and provide a unique test case for the kinetic theory introduced in Section 4.3. Future work will involve bridging longer length scales in order to accommodate larger cluster sizes, thereby providing further refinement of the qualitative picture presented here. Acknowledgment. This work was supported by NSF CHE0310107 and by the Camille and Henry Dreyfus Foundation, Inc. (TC-02-012). J.A.M. would like to acknowledge the Fannie and John Hertz Foundation for funding during the latter phase of his work on this publication. References and Notes (1) Finney, J.; Soper, A. Phys. ReV. Lett. 1993, 71, 4346. (2) Bako, I.; Palinkas, G.; Heinzinger, K. Z. Naturforsch. A 1994, 49, 967. (3) Skaf, M. S.; Ladanyi, B. M. J. Chem. Phys. 1995, 102, 6542. (4) Ludwig, R. Chem. Phys. 1996, 195, 329. (5) Ladanyi, B. M.; Skaf, M. S. J. Phys. Chem. 1996, 100, 1368. (6) Skaf, M. S.; Ladanyi, B. M. J. Phys. Chem. 1996, 100, 18258. (7) Laaksonen, A.; Kusalik, P. G.; Svishchev, I. M. J. Phys. Chem. A 1997, 101, 5910. (8) Dixit, S.; Poon, W. C. K.; Crain, J. J. Phys.: Condensed Matter 2000, 12, L323. (9) Marcus, Y. Phys. Chem. Chem. Phys. 1999, 1, 2975. (10) van Erp, T. S.; Meijer, E. J. Chem. Phys. Lett. 2001, 333, 290. (11) Noskov, S. Y.; Kiselev, M. G.; Kolker, A. M.; Rode, B. M. J. Mol. Liq. 2001, 91, 157. (12) Kvamme, B. Phys. Chem. Chem. Phys. 2002, 4, 942. (13) Dixit, S.; Soper, A.; Finney, J.; Crain, J. Europhys. Lett. 2002, 59, 377. (14) Dixit, S.; Crain, J.; Poon, W.; Finney, J.; Soper, A. Nature (London) 2002, 416, 829. (15) Sato, T.; Chiba, A.; Nozaki, R. J. Mol. Liq. 2002, 101, 00. (16) Dougan, L.; Bates, S. P.; Hargreaves, R.; Fox, J. P.; Crain, J.; Finney, J. L.; Reat, V.; Soper, A. K. J. Chem. Phys. 2004, 121, 6456. (17) Allison, S. K.; Fox, J. P.; Hargreaves, R.; Bates, S. P. Phys. ReV. B 2005, 71, 024201. (18) Chang, T. M.; Dang, L. X. J. Phys. Chem. B 2005, 109, 5759. (19) Dougan, L.; Hargreaves, R.; Bates, S. P.; Finney, J. L.; Reat, V.; Soper, A. K.; Crain, J. J. Chem. Phys. 2005, 122, 174514. (20) Frank, H.; Evans, M. J. Chem. Phys. 1945, 13, 507. (21) Guo, J. H.; Luo, Y.; Augustsson, A.; Kashtanov, S.; Rubensson, J.-E.; Shuh, D. K.; Argren, H.; Nordgren, J. Phys. ReV. Lett. 2003, 91, 157401. (22) Vishnyakov, A.; Neimark, A. V. J. Phys. Chem. B 2001, 105, 7830. (23) Kreuer, K. D. J. Membr. Sci. 2001, 185, 29. (24) Kreuer, K. D.; Paddison, S. J.; Spohr, E.; Schuster, M. Chem. ReV. 2004, 104, 4637. (25) Heberle, J. J.; Riesle, G.; Thiedermann, D. O.; Oesterhelt, D.; Dencher, N. A. Nature 1994, 370, 379. (26) Scherrer, P. U.; Alexiev, T.; Marti, H.; Gobind, K.; Heyn, M. P. Biochemistry 1994, 33, 13684. (27) Nachliel, E.; Gutman, M. FEBS Lett. 1996, 393, 221. (28) Smondyrev, A. M.; Voth, G. A. Biophys. J. 2002, 82, 1460. (29) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M. J. Chem. Phys. 1995, 103, 150. (30) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M. J. Phys. Chem. 1995, 99, 5749. (31) Agmon, N. Chem. Phys. Lett. 1995, 244, 456. (32) Vuilleumier, R.; Borgis, D. J. Mol. Struct. 1997, 437, 555. (33) Vuilleumier, R.; Borgis, D. J. Phys. Chem. B 1998, 102, 4261.

3720 J. Phys. Chem. B, Vol. 110, No. 8, 2006 (34) Vuilleumier, F.; Borgis, D. J. Chem. Phys. 1999, 111, 4251. (35) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. Nature 1999, 397. (36) Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 1999, 111, 9361. (37) Day, T. J. F.; Schmitt, U. W.; Voth, G. A. J. Am. Chem. Soc. 2000, 122, 12027. (38) Tuckerman, M. E.; Marx, D.; Parrinello, M. Nature 2002, 417. (39) Tuckerman, M. E.; Chandra, A.; Marx, D. Acc. Chem. Res. in press. (40) Lapid, H.; Agmon, N.; Petersen, M. K.; Voth, G. A. J. Chem. Phys. 2005, 122, 014506. (41) Chang, H. C.; Jiang, J. C.; Lin, S. H.; Lee, Y. T. J. Phys. Chem. A 1999, 103, 2941. (42) Morrone, J. A.; Tuckerman, M. E. J. Chem. Phys. 2002, 117. (43) Yamaguchi, T.; Hidaka, K.; Soper, A. K. Mol. Phys. 1999, 96, 1159. (44) Yamaguchi, T.; Hidaka, K.; Soper, A. K. Mol. Phys. 1999, 97, 603. (45) Raduge, C.; Pflumio, V.; Shen, Y. R. Chem. Phys. Lett. 1997, 274, 140. (46) Miranda, P. B.; Shen, Y. R. J. Phys. Chem. B 1999, 103, 3292. (47) Shultz, M. J.; Schnitzer, C.; Simonelli, D.; Baldelli, S. Int. ReV. Phys. Chem. 2002, 19, 123. (48) Petersen, M. K.; Iyengar, S. S.; Day, T. J. F.; Voth, G. A. J. Phys. Chem. B 2004, 108, 14804. (49) Iyengar, S. S.; Day, T. J. F.; Voth, G. A. Int. J. Mass. Spectrosc. 2005, 241, 197. (50) Lee, H. S.; Brancato, G.; Tuckerman, M. E. (to be submitted). (51) Stillinger, F. H. AdV. Chem. Phys. 1975, 31, 1. (52) Rappaport, D. C. Mol. Phys. 1983, 50, 1151. (53) Luzar, A.; Chandler, D. Nature 1996, 379, 55. (54) Chandra, A. Phys. ReV. Lett. 2000, 85, 768. (55) Balasubramanian, S.; Pal, S.; Bagchi, B. Phys. ReV. Lett. 2002, 89, 115505. (56) Xu, H.; Stern, H. A.; Berne, B. J. J. Phys. Chem. B 2002, 106, 2054. (57) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (58) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (59) C. Lee, W. Y.; Parr, R. C. Phys. ReV. B 1988, 37, 785. (60) The recent controversy surrounding the structure of pure BLYP water deserves a brief comment. The present results are in agreement with other recently published structure using the BLYP functional77,61 yet appear to be at odds with those of refs 78. There appears to a growing consensus that the unexpectedly high peak heights presented in ref 78 are not

Morrone et al. representative of either the BLYP functional or of the Car-Parrinello technique. (61) Mantz, Y. A.; Chen, B.; Martyna, G. J. Chem. Phys. Lett. 2005, 405, 294. (62) Lee, H. S.; Tuckerman, M. E. J. Phys. Chem. A (submitted). (63) Yeh, I.-C.; Hummer, G. J. Phys. Chem. B 2004, 108, 15873. (64) Vanderbilt, D. Phys. ReV. B 1990, 41, 7892. (65) Martyna, G.; Tuckerman, M.; Klein, M. J. Chem. Phys. 1992, 97, 2635. (66) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (67) Chelli, R.; Ciabatti, S.; Cardini, G.; Righini, R.; Procacci, P. J. Chem. Phys. 1999, 111, 4218. (68) Hutter, J.; Tuckerman, M. E.; Parrinello, M. J. Chem. Phys. 1995, 102, 859. (69) Hutter, J.; Alavi, A.; Deutsch, T.; Bernasconi, M.; Goedecker, S.; Marx, D.; Parrinello, M.; Tuckerman, M. E. CPMD Version 3.5, MaxPlanck Institut fuer Festkoerperforschung and IBM Zurich Research Laboratory, 1995-1999. (70) Morrone, J. A.; Tuckerman, M. E. Chem. Phys. Lett. 2003, 370. (71) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43, 1993. (72) Tuckerman, M. E.; Yarne, D. A.; Samuelson, S. O.; Hughes, A. L.; Martyna, G. J. Comput. Phys. Commun. 2000, 128, 333. (73) Halle, B.; Karlstrom, G. J. Chem. Soc., Faraday Trans. 2 1983, 70, 1031. (74) Halle, B.; Karlstrom, G. J. Chem. Soc., Faraday Trans. 2 1983, 70, 1047. (75) Izvekov, S.; Voth, G. A. J. Chem. Phys. 2005, 123, 044505. (76) At present, it appears that the low-proton diffusion constant in water reported in ref 75 could be due to differences in the equilibration protocol. The proton diffusion constant reported in ref 39 contained over 30 ps of NVT equilibration within the employed AIMD protocol before the generation of a production run for determining the diffusion constant. (77) Kuo, I. F. W.; Mundy, C. J.; McGrath, M. J.; Siepmann, J. I.; VandeVondele, J.; Sprik, M.; Hutter, J.; Chen, B.; Klein, M. L.; Mohamed, F.; m. Krack.; Parrinello, M. J. Phys. Chem. B 2004, 108, 12990. (78) Grossman, J. C.; Schwegler, E.; Draeger, E. W.; Gygi, F.; Galli, G. J. Chem. Phys. 2004, 120, 300. (79) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33.