J. Phys. Chem. 1996, 100, 5611-5615
5611
Hydration Shell Exchange Kinetics: An MD Study for Na+(aq) Rossend Rey* Departament de Fı´sica i Enginyeria Nuclear, UniVersitat Polite` cnica de Catalunya, Campus Nord B4-B5, Barcelona 08034, Spain
James T. Hynes Department of Chemistry and Biochemistry, UniVersity of Colorado, Boulder, Colorado 80309-0215 ReceiVed: NoVember 22, 1995; In Final Form: February 14, 1996X
A molecular dynamics simulation method for the first hydration shell exchange kinetics is described and is illustrated for Na+ in water. The exchange process is found to have a complex character. A calculated small transmission coefficient, κ ≈ 0.21, indicates that transition-state theory seriously overestimates the rate. GroteHynes theory is unsuccessful in accounting for the κ value. This is traced to trajectory mechanistic aspects for which the theory’s assumptions break down. Inclusion of anharmonicity via the theory of Voth provides an improved estimate.
Introduction The exchange of a water molecule in the first hydration shell of an ion is a phenomenon of long-standing interest.1 In one limit, a multiply charged ion and its first hydration shell form a chemical unit due to strong electrostatic interactions, and the slow water ligand exchange process can be of bimolecular (SN2) or unimolecular (SN1) character.2 Here the hydration shell exchange is a variant of familiar chemical reaction classes. In the other limit, where the ion is singly charged, the situation is much less clear: electrostatic ion-solvent interactions are weaker and exchange times are variously estimated in the few ps to tens of picoseconds range.1,3-5 Nonetheless, the process is of obvious interest in connection with ionic transport,6 as well as in reactions of the ion with other species, where the hydration shell must rearrange.7 In this letter, we present a molecular dynamics (MD) simulation method for the hydration shell exchange process, here applied to the aqueous Na+ ion, based on a unimolecular dissociation perspective. The method is similar to that used for contact ion pair to solvent-separated ion pair interconversion reactions8-10 and in general allows the practical short-time simulation of a slow activated process. The dissociation rate constant k for Na+(aq) itself is sufficiently large that k can also be determined by direct simulation3 for comparison. The outline of this letter is the following. We first describe the model, theoretical framework, and simulation techniques. We then discuss the hydration shell exchange process, and indicate its complex character. We also discuss the exchange rate constant and put its value into theoretical perspective. In particular, the present reaction provides the first example where Grote-Hynes11 theory fails significantly for a realistic molecular model MD simulation. The reason for this is discussed in the context of the reaction dynamics. Approximate inclusion of anharmonicity effects12,13 improves the rate estimate. Methodology The dissociation rate constant can be written as the product k ) κkTST, where kTST is the transition state theory (TST) rate constant, and κ is the transmission coefficient. Here we define the reaction coordinate as the distance r from the Na+ ion to X
Abstract published in AdVance ACS Abstracts, March 15, 1996.
0022-3654/96/20100-5611$12.00/0
the water molecule center of mass, viewing the process in a unimolecular dissociation perspective. Given that each ion is surrounded by some six waters, the potential of mean force (pmf) can be determined with acceptable statistics from the radial distribution function14 W(r) ) -β-1 ln(g(r)). Once the pmf is computed, kTST can be readily determined from9
kT ST )
x
kBT 2πµ
rq e-βW(r ) 2
q
∫0r dr r2e-βW(r)
)
q
x
kBT 2πµ
e-βWeff(r ) q
∫0r dr e-βW q
eff(r)
(1)
which defines the centrifugally averaged effective potential Weff(r). µ is the ion-water molecule pair reduced mass, and rq indicates the barrier top position. The transmission coefficient is determined from the plateau value of the normalized reactive flux,15 computed in the constrained reaction coordinate ensemble;16 for the present case:9,10
k(t) ) 〈r˘ (0) θ[r(t) - rq]〉c/〈r˘ (0) θ[r˘ (0)]〉c
(2)
To compute the averages (〈...〉c), we start from equilibrium configurations with the reaction coordinate constrained at the TS, which after the release of this constraint and the generation of velocities according to a thermal distribution, are followed in time. The present problem also allows a direct rate constant computation3,4 given the short-time scale for the escape from the first hydration shell. Given an initial equilibrium configuration, we compute the escape time for each hydration shell molecule, which should be well represented by an exponential with characteristic time k-1. A molecule is assumed in the first hydration shell when its center of mass distance to the ion is smaller than rq, and no longer if it exceeds that distance for longer than 2 ps (which from trajectory analysis is characteristic of long excursions out of the first hydration shell). We have performed simulations of a sodium ion plus 107 water molecules in a cubic box with standard periodic conditions, with a side length giving a solvent density of 1 g/cm3, and the temperature fixed at a mean value of 298 K. The water model is rigid SPC.17 The ion-water interaction potential function is that in ref 18. Long-range forces were computed by the Ewald summation method.19 A leapfrog integration © 1996 American Chemical Society
5612 J. Phys. Chem., Vol. 100, No. 14, 1996
Letters
Figure 1. Potential of mean force W(r) (solid line) and centrifugally averaged effective potential Weff(r) (dot-dashed line). Vertical dashed lines indicate the distances for which the friction memory function has been computed. The calculated equilibrium constant Keq for second shell versus first shell residence is calculated from the equation in ref 9a, with an outer cutoff of r ) 5.54 Å (the second local maximum of W(r)). The result is Keq ≡ exp(-∆Gexch/RT) ) 3, and ∆Gexch ) -0.65 kcal/mol.
algorithm with coupling to a thermal bath20 has been used, with a 1.5 fs time step. Three different types of simulations have been performed. In the first set, the system is followed with no constraints applied (excepting those keeping the water molecules rigid via the SHAKE algorithm21). The radial distribution function is accurately computed, and statistics are collected for the hydration molecule survival time (with the definition rq ) 3.18 Å) from six runs of 300 ps each (after 20 ps initial stabilization for each set). An exponential is fitted to each resulting decay function; this set of values is used in the estimation of the mean value error. A second set of simulations involves the computation of the reaction coordinate friction functions ζ(t;r), determined from the fluctuations of the force acting on the reaction coordinate: 22
ζ(t;r) ≡ (β/µ)〈δF(t;r) δF(0;r)〉
(3)
and providing the basic ingredient of rate theories (see below). They are computed by constraining the reaction coordinate at a specific r value and subsequent determination of eq 3. This function may depend on the chosen value of r,9,23,24 and has been computed for r ) 2.4, 2.7, 3.18, 3.7, and 4.4 Å (see Figure 1). For r ) 3.18 Å, the determined TS, five sets of 300 ps each have been used, and a single run of 300 ps for the other separations with initial stabilizations in each case.25 A last set of simulations was required to compute the transmission coefficient κ from eq 2. First, initial configurations with the distance constrained at r ) 3.18 Å were picked every 3 ps from the five independent sets used in the friction memory function computation for r ) 3.18 Å; a total of 500 initial configurations were generated. Each initial configuration was followed forward and backward in time for 2.1 ps, with a total of 1000 trajectories used in the κ determination. The error is estimated as the standard deviation of the mean, after trajectories are grouped in 10 sets, and the transmission coefficient separately computed for each set.26
Figure 2. Solid line curve: MD reactive flux obtained. Horizontal dot-dashed line: GH theory prediction. Horizontal thin solid line: Kramers theory prediction. Horizontal thick solid line: Voth theory prediction. Dashed line curve: numerical solution of the anharmonic GLE (see eq 7).]
TABLE 1: Obtained Transmission Coefficients: KMD, Molecular Dynamics Result; KKr, Kramers Theory; KGH, Grote-Hynes Theory; KVTST, (Exact) Variational TST Theory for Anharmonic Potential; KV, Voth Theory κMD
κKr
κGH
κVTST
κV
0.21 ( 0.05
0.07
0.35
0.12
0.26
Simulation Results and Analysis The calculated pmf versus Na+-H2 O separation is displayed in Figure 1, together with its centrifugally corrected value. There are several important features to note. The first is the barrier height of ≈4 kBT for the reaction in the forward direction, and the second is the rather low barrier (≈2 kBT ) in the reverse direction, from a rather weak and broad minimum. The TST dissociation rate constant for the reaction is estimated from eq 1 to be kTST) 0.14 ps-1. The transmission coefficient correcting this kTST result is determined from the reactive flux (eq 2 and Figure 2) to be κMD) 0.21 ( 0.05 (cf. Table 1).30 There is thus a significant departure from TST due to barrier recrossing.28 Indeed, this κ value is among the smallest observed in MD computer simulation.9,10,31 We now check that this overall forward rate constant k ) κkTST is indeed correct by direct simulation as detailed in the previous section. A decay time of 34 ( 3 ps is obtained,32 which should be compared with k-1) 34 ( 8 ps. Both results are consistent, thus confirming the small transmission coefficient and slightly narrowing the confidence interval for κ. Since the rate constant is a steady-state flux,22 it can be evaluated with, and is independent of, any coordinate choice, although the partitioning of k between a TST value and a transmission coefficient depends on the particular choice. The choice employed above is attuned to an SN1 unimolecular dissociation process, in which exchange would occur (from the most typical first coordination shell population number of 6) by escape of a water to the second shell (followed much later by an uncorrelated replacement event). Figure 3a corresponds to such an event, although there is a small amplitude recrossing of the TS. However, trajectories often are of SN 2 character (cf. Figure 3b), where the entrance of one water and exit of another from the first solvent shell are nearly simultaneous events. Unsuccessful departures of single water molecules
Letters
J. Phys. Chem., Vol. 100, No. 14, 1996 5613
Figure 4. Distance Na+-water molecule center of mass for a molecule leaving the first hydration shell. Horizontal dashed line indicates the position of rq. We have also calculated the distribution of maximum displacements (∆rm ) for trajectories that recross the barrier. This gives a mean maximum displacement of 0.3 Å. The distribution is peaked at r ) rq (∆rm ) 0) but also has a small peak near ∆r ) -1 Å (≈15% of the central peak height) and a tail that is still ≈3% of the central peak value at ∆rm ) + 1 Å.
TABLE 2: Initial Value of the Friction Memory Function (ζ(0;r)) and Correlation Time (τ(r) ≡ (1/ζ(0;r))∫0∞dt ζ(t; r)) dist (Å)
ζ(0;r) (ps-2 )
τ(r) (ps)
2.40 2.70 3.18 3.70 4.40
1135 1230 1600 1712 1693
0.094 0.008a 0.188 0.094 0.078
a
This low value should be attributed not to an abrupt variation of the memory friction behavior for this distance but to the effect of statistical noise at large times, which results in a less-reliable correlation time in this particular case.
is a fairly complex one of disperse character, and that significant recrossings of a simple TS such as rq should be expected. The possibility that the process is more simple in a more complex coordinate, e.g., one associated with the entire first solvent shell, is a topic for future research. We now consider if the MD transmission coefficient can be accounted for theoretically. Grote-Hynes (GH) theory, which has been uniformly successful in this connection for a very wide range of reaction types,9,10,31 gives κ as11
κGH ) [κGH + (1/ωb) ∫∞0 dt e-ωbκGHtζ(t)]-1
Figure 3. Na+-water center of mass distance for water molecules that cross the TS. (a, top) SN1 type; (b, middle) SN2 type; (c, bottom) mixed type. Horizontal dashed line indicates the position of rq.
without accompanying insertions are also evident in Figure 3bsthese are recrossings of the simple rq TS. From these and other trajectories which are neither “clean” SN 1-like nor “clean” SN 2-like or which exhibit both characteristics in the span of a few picoseconds (Figure 3c), it is clear that the exchange process
(4)
in which ωb is the frequency of the effective potential barrier (ωb ≈ 21.6 ps-1 ) and the time-dependent friction corresponds to r ) rq. The initial values and correlation times of ζ(t; r) functions for selected distances are tabulated in Table 2. The ζ(t; r) character is similar to that found in related systems:9,10 a fast falloff on a short time scale (≈0.1 ps), followed by a slow decay on a time scale approximately 1 order of magnitude larger. The computed κGH value is 0.35, noticeably above the MD result κMD) 0.21 (plausible slight variations in the determination of ωeq, e.g., the range over which the inverse parabola is fitted, do not significantly change this conclusion). For comparison, the Kramers theory transmission coefficient in which the frequency dependence of the friction in eq 4 is ignored, resulting
5614 J. Phys. Chem., Vol. 100, No. 14, 1996
Letters
in the appearance of the friction constant ζ ) ∫∞0 dt ζ(t), has a value κKr ) 0.07, a clear underestimate of κMD. In the Kramers view, the fully time integrated friction is operative in the barrier crossing, with very extensive recrossing. Comparison of ζ(t;rq) (not shown) and Figure 2 shows that the long-time friction component, which helps to establish a large value of ζ (cf. Table 2), does not in fact have time to contribute on the time scale determining κ. The overestimate of κMD by κGH is striking and anomalous. What is the origin of the failure? A clear indication of this is revealed in Figure 4, where it is seen that the water molecule can often make very large excursions from the first hydration shell, into the second hydration shell, and recross the transition state toward the reactant side before stabilizing on the product side. This is inconsistent with the basic assumptions of GH theory,11,22,31b where κ is supposed to be determined by recrossing events in the neighborhood of the barrier top, such that the mean potential is locally parabolic and there is no serious spatial dependence of the time-dependent friction. We have also tested the theory proposed by Voth,12 which includes some aspects of anharmonicity and spatially dependent friction in an “effective” GH approach. Here, the transmission coefficient is obtained applying the same GH formula but with an averaged barrier frequency ωb,eff and initial memory function ζeff(0) given by12 (in mass weighted coordinates)
x x
ωb,eff2 ) -
ζeff(0) )
d2Weff(q + qq) -βωb2q2/2 βωb2 e dq ∫ 2π dq2
(5)
βωb2 2 2 ∫dq ζ(0;q + qq)e-βωb q /2 2π
(6)
We have applied these equations using the (interpolated)initial values of the memory function in Table 2, the potential Weff(r) (which has rq ≈ 3.14 Å instead of 3.18 Å), and the reference memory function at r ) 3.18 Å (which should not differ noticeably from that at r ) 3.14 Å). The result κV ) 0.26 is on the upper boundary of the κMD confidence interval. Although the initial friction varies with r (Table 2), a negligible change in κ results from eq 6: the deviation of the initial friction from its rq value is approximately antisymmetric over the symmetric averaging range. Essentially the entire reduction from κGH arises from the anharmonicity factor in eq 5. The violation of the parabolic barrier assumption is also considered by variational theories33,34 in which anharmonicity is incorporated into the calculation of κ. The best possible estimate from a subset of these theories would correspond to the solution of the generalized langevin equation:34
µr¨ ) -
dWeff(r) - µ∫t0dτ ζ(τ;rq) r˘ (t - τ) + δF(t) (7) dr
We have numerically solved this equation with standard methods.27,35,36 The transmission coefficient, κVTST ≈ 0.12, does not compare well with κMD, so that anharmonicity, in the GLE context, is unable to account for the discrepancy. In addition, Figure 2 shows that a plateau in k(t) is not clearly reached in the correct time scale (and κ is not well defined), as opposed to the MD result. Thus, eq 7 is not a useful description of the full dynamics.37 Concluding Remarks The present MD study has illustrated, via the application to Na+ (aq), a method for the simulation of hydration shell
exchange dynamics. It has also revealed the first breakdown in GH theory in a realistic reaction MD simulation and indicated its basic source; a recently proposed theory for anharmonicity is clearly closer to the exact value. The MD method employed will be used in future studies to clarify the details of the molecular mechanism for the Na+ hydration shell exchange. Further, it can be used for multiply charged ion cases, where a direct simulation of the very slow exchange process is not feasible. Acknowledgment. R.R. gratefully acknowledges fellowship support by the Spanish Ministerio de Educacio´n y Ciencia and support from DGICYT Project PB93-0971-CO3. This work was supported in part by NSF Grants CHE88-07852 and CHE93-12267 and a grant of CRAY-YMP time from the Pittsburgh Supercomputer Center. References and Notes (1) Friedman, H. L. Chem. Scr. 1985, 25, 42. (2) Orgel, L. E. An Introduction to transition-metal chemistry: ligandfield theory; Wiley: London, 1966. (3) Impey, R. W.; Madden, P.; McDonald, I. R. J. Phys. Chem. 1983, 87, 5071. (4) Gua`rdia, E.; Padro´, J. A. J. Phys. Chem. 1990, 94, 6049. (5) Hermansson, K., personal communication. (6) Wolynes, P. G. Annu. ReV. Phys. Chem. 1980, 31, 345. (7) For an example for the hydronium cation, see e.g.: Ando, K.; Hynes, J. T. J. Mol. Liq. 1995, 64, 25. For the Na+ cation, solvation shell rearrangement in water clusters (or in low-polarity liquids) could prove of interest for excited- state reactions of sodium halides (see e.g.: Rose, T. S.; Rosker, M. J.; Zewail, A. H. J. Chem. Phys. 1989, 91, 7415. Benjamin, I.; Wilson, K. R. J. Phys. Chem. 1991, 95, 3541). (8) Karim, O. A.; McCammon, J. A. Chem. Phys. Lett. 1986, 132, 219. (9) (a) Ciccotti, G.; Ferrario, M.; Hynes, J. T.; Kapral, R. Chem. Phys. 1989, 129, 241. (b) J. Chem. Phys. 1990, 93, 7137. (10) Rey, R.; Gua`rdia, E. J. Phys. Chem. 1992, 96, 4712. (11) Grote, R. F.; Hynes, J. T. J. Chem. Phys. 1980, 74, 2715. (12) Voth, G. A. J. Chem. Phys. 1992, 97, 5908. (13) See the appendix in: Haynes, G. R.; Voth, G. A.; Pollak, E. J. Chem. Phys. 1994, 101, 7811. (14) Hill, T. L. Statistical Mechanics; McGraw-Hill: New York, 1956. (15) Chandler, D. J. Chem. Phys. 1978, 68, 2959. (16) Carter, E. A.; Cicotti, G.; Hynes, J. T.; Kapral, R. Chem. Phys. Lett. 1989, 156, 472. (17) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981. (18) Pettitt, B. M.; Rossky, P. J. J. Chem. Phys. 1986, 84, 5836. (19) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford: New York, 1989; Bader, J. S.; Chandler, D. J. Phys. Chem. 1992, 96, 6424. (20) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; di Nola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3683. (21) Ryckaert, J. P. Mol. Phys. 1985, 55, 549. (22) Hynes, J. T. In The Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC: Boca Raton, FL, 1985; Vol. IV. (23) Straub, J. E.; Borkovec, M.; Berne, B. J. J. Phys. Chem. 1987, 91, 4995. (24) Rey, R.; Gua`rdia, E.; Padro´, J. A. J. Chem. Phys. 1992, 97, 1343. (25) While the dynamics should be Newtonian in the computation of these functions, we have kept the thermal coupling in order to avoid temperature drift over the long simulation runs required. Nevertheless, the thermal coupling constant (0.1 ps) has been chosen so that the effects on correlation functions be negligible (this was originally20 checked for velocity autocorrelation functionsswhich can be determined with much better accuracysusing the same water model). Additional test runs (of 300 ps) with different coupling constants (larger or bigger by 1 order of magnitude) indicate that the small differences between memory frictions do not show any systematic trend in this range and therefore should be attributed to statistical fluctuations. (26) This is a more conservative method than the usual estimation from the global reactive flux function.27 (27) Rey, R. J. Chem. Phys. 1996, 104, 1966. (28) The vibrational energy relaxation time in the reactant well has been estimated (by the Landau-Teller approach,29 with the time-dependent friction at the bottom of the reactant well) to be ≈0.1 ps. Thus vibrational energy relaxation is sufficiently fast to play no role in reducing κ below unity.
Letters (29) Oxtoby, D. W. AdV. Chem. Phys. 1981, 47, 487. Whitnell, R. M.; Wilson, K. R.; Hynes, J. T. J. Phys. Chem. 1990, 94, 8625. Rey, R.; Hynes, J. T. J. Chem. Phys. 1996, 104, 2356. (30) Considerations similar to those in footnote 25 apply here. While it is feasible in this case to remove thermal coupling, we have maintained it in the results we report for consistency with the method employed in the calculation of the friction memory functions. In any case, the computation of κ (t) without thermal coupling yields κ ) 0.19 ( 0.04 (for 1400 trajectories), so that both procedures agree within statistical indeterminacy, thus supporting the idea that thermal coupling barely affects the results. (31) (a) Bergsma, J. P.; Reimers, J. R.; Wilson, K. R.; Hynes, J. T. J. Chem. Phys. 1986, 85, 5625. Bergsma, J. P.; Gertner, B. J.; Wilson, K. R.; Hynes, J. T. Ibid. 1987, 86, 1356. Gertner, B. J.; Bergsma, J. P.; Wilson, K. R.; Lee, S.; Hynes, J. T. Ibid. 1987, 86, 1377. (b) Gertner, B. J.; Wilson, K. R.; Hynes, J. T. Ibid. 1989, 90, 3537. (c) Keirstead, W. P.; Wilson, K. R.; Hynes, J. T. Ibid. 1991, 95, 5256. Zichi, D. A.; Ciccotti, G.; Hynes, J. T.; Ferrario, M. J. Phys. Chem. 1989, 93, 6261. Smith, B. B.; Staib, A.; Hynes, J. T. Chem. Phys. 1993, 176, 521. Staib, A.; Borgis, D.; Hynes, J. T. J. Chem. Phys. 1995, 102, 2487. (d) Straub, J. E.; Borkovec, M.; Berne, B. J. Ibid. 1988, 89, 4833. Zhu, S. B.; Lee, J.; Robinson, G. W. J. Phys. Chem. 1988, 92, 2401. Tucker, S.; Truhlar, D. J. Am. Chem. Soc. 1990, 112, 3347. Roux, B.; Karplus, M. J. Chem. Phys. 1991, 95, 4856. Benjamin,
J. Phys. Chem., Vol. 100, No. 14, 1996 5615 I.; Lee, L. L.; Li, Y. S.; Liu, A.; Wilson, K. R. Chem. Phys. 1991, 152, 1. Lee, L. L.; Li, Y. S.; Wilson, K. R. J. Chem. Phys. 1991, 95, 2458. (32) This is also in reasonable agreement with a recent estimate of ≈40 ps for Na+ in SPC/E water (which has slightly different charges than the present SPC water)(Gua`rdia, E., personal communication). (33) A different variational theory (Lee, S.; Hynes, J. T. J. Chem. Phys. 1988, 88, 6853, 6863) cannot presently be applied due to the lack of knowledge of an appropriate collective solvent coordinate for the current problem. It should also be noted that GH theory in its standard form is predicted (Mathis, J.; Hynes, J. T. J. Phys. Chem. 1994, 98, 5445) to break down for the t-BuI SN 1 ionization, and a variational approach is required. (34) Pollak, E.; Tucker, S. C.; Berne, B. J. Phys. ReV. Lett. 1990, 65, 1399. Tucker, S. C.; Pollak, E. J. Stat. Phys. 1992, 66, 975. (35) Berkowitz, M.; Morgan, J.; McCammon, J. A. J. Chem. Phys. 1983, 78, 3256. (36) Tucker, S. C.; Tuckerman, M. E.; Berne, B. J.; Pollak, E. J. Chem. Phys. 1991, 95, 5809. (37) See e.g. ref 31b for a detailed discussion of the theoretical limitations of the GLE approach to reaction rate problems in the more restricted region of the transition-state barrier.
JP953429Z
