Reaction Rate Theory in Coordination Number Space: An Application

Article pubs.acs.org/JPCC

Reaction Rate Theory in Coordination Number Space: An Application to Ion Solvation Santanu Roy,* Marcel D. Baer, Christopher J. Mundy, and Gregory K. Schenter* Physical Science Division, Pacific Northwest National Laboratory, 902 Battelle Blvd., Richland, Washington 99352, United States S Supporting Information *

ABSTRACT: Understanding reaction mechanisms in many chemical and biological processes requires application of rare event theories. In these theories, an effective choice of a reaction coordinate to describe a reaction pathway is essential. To this end, we study ion solvation in water using molecular dynamics simulations and explore the utility of coordination number (n = number of water molecules in the first solvation shell) as the reaction coordinate. Here, we compute the potential of mean force (W(n)) using umbrella sampling, predicting multiple metastable n-states for both cations and anions. With increasing ionic size, we find these states become more stable and structured for cations when compared to anions. We have extended transition state theory (TST) to calculate transition rates between n states. TST overestimates the rate constant due to solvent-induced barrier recrossings that are not accounted for. We correct the TST rates by calculating transmission coefficients using the reactive flux method. This approach enables a new way of understanding rare events involving coordination complexes.

■

The equilibrium structure of ion solvation37 has been studied by many different techniques such as the extended X-ray fine structure,38 X-ray39/neutron40 diffraction, and nuclear magnetic resonance.40−42 All of the aforementioned experimental techniques can determine the coordination number. The precise equilibrium structure of the first solvation shell and the coordination number for the diverse set of ions have been shown to be extremely sensitive to the form of the short-range molecular ion−water interaction.37 Exploring the dynamical aspects of the coordination number reaction coordinate will allow one to examine transition mechanisms and rates from undercoordinated to overcoordinated states during the dynamics of ion solvation. However, identifying multiple metastable coordination states and determining their kinetics with experimental techniques can be very challenging. These techniques provide only average coordination numbers and say nothing quantitative about the fluctuations. To this end, computer simulation incorporating a rigorous theoretical framework for rate theory is necessary to resolve this problem. Recently, in order to examine the metal ion−carbonate paring in water, the PMF computation for the coordination number of the ions and the calculation of the equilibrium transition rate constants between the metastable coordination states have been carried out.43 However, TST formulation and the correction due to solvent dynamical effects have not been attempted. Herein, we demonstrate a complete picture of the

INTRODUCTION Reaction rate theory1 is the building block for the investigation of many fundamental scientific problems such as ion solvation dynamics in bulk solvent, liquid interfaces, and ion channels.2−6 In the condensed phase, the theory centers around the potential of mean force (PMF) as a function of a reaction coordinate.6−10 A system with multiple metastable states can be identified by a PMF with minima separated by dividing surfaces known as transition states. Rate theory is a framework that describes long-time behavior of these metastable states and their interconversion mechanisms. It has advanced significantly from the simple picture of transition state theory (TST)10−13 to more accurate Kramers’ and Pollak’s theory,14,15 Grote−Hynes’ stable state theory,9 and Bennett’s and Chandler’s theory of reactive flux (RF)16−19 that account for solvent dynamical effects on reaction rates. Using these theories, extensive work on solvent exchange dynamics around ions and solvent mediated ion-pair association/dissociation were carried out in the last several decades.6,18−35 One of the fundamental assumptions in constructing the PMF is the choice of reaction coordinate. In general, the reaction coordinate for rare events can be multidimensional. However, the most frequently used reaction coordinate is the distance between reacting species. An exception, of course, is the rate theory of the solvent energy gap used in Marcus theory.36 Herein, we extend the rate theory to investigate the rare events of solvent exchange by exploring many different metastable states of the ion−solvent coordination number. This reaction coordinate is defined to be the number of solvent molecules in the first solvation shell of an ion. © XXXX American Chemical Society

Received: January 14, 2016 Revised: March 28, 2016

A

DOI: 10.1021/acs.jpcc.6b00443 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C reaction rate theory that provides a rigorous theoretical framework for computing the rates of change in coordination number around an ion. We apply this theory to model cations and anions solvated in water at room temperature (298 K) using classical empirical interaction potentials in conjunction with molecular dynamics (MD) simulations. Through the systematic investigation of kinetics of these ion−water coordination complexes with respect to ionic size and charge, we are able to extract general principles about how the asymmetry of water is reflected in both the PMF and dynamics in coordinations space.

H=

= K (pn ) + K (pq ) + V (n , q)

3N

Zn =

k TST =

i

1 − (ri /r †)a 1 − (ri /r †)2a

1 C

=

∑ f (ri) i

∫ dx1 ... dxl dp1 ... d pl δ(ξ − n)exp(−βH)

Qr =

1 Qr

∫ dq dpqdn dpn δ(n − n†)Znpn θ(Znpn )exp(−βH)

∫R dq dpq dn dpn exp(−βH)

(8)

We demonstrate the derivation of Zn and Qr in the Supporting Information, where we utilize n as given in eq 2 because this form treats n as a continuous variable. Now, inserting the following expressions in eq 7 Zn =

(2)

Λ μ

(9)

and Qr =

2πμ Λβ

∫

dpq exp[−βK (pq )]

∫0

n†

dn exp(−βW (n)) (10)

the TST rate constant (see the Supporting Information for details) is k TST =

Λ 2πμβ

exp( −βW (n†)) †

n

∫0 dn exp(−βW (n))

(11)

Here, μ is the ion−water reduced mass: μ = mion + mwater−1, given that mion and mwater are the masses of the ion and a water molecule, respectively. We show in the Supporting Information that Λ has the following expression: −1

(3)

⎡ N−1 2μ Λ = ⎢ ∑ (f i′ )2 + ⎢⎣ m ion i

Here, C is the normalization constant and β = 1/kBT with kB denoting the Boltzmann constant; H (function of the l = 3N Cartesian coordinates x = x1,..., xl and their conjugate momenta p = p1,..., pl) consists of kinetic and potential energy terms. This probability density is related to the PMF as follows: W (n) = −kBT ln Ω(n)

(6)

(7)

This polynomial form is equivalent to the step function in eq 1 for a large value of a. The smoothness of this form can be controlled by varying a. We choose a = 24 to make sure that the polynomial form is close to the step function, but it maintains the smoothness. The number of metastable states in coordination space should be independent of the choice of a; however, their distributions can be narrower or broader depending on the values of a. Note that we use this smooth function for the TST formulation of coordination number. Now, in analogy to a generalized theory by Schenter, Garrett, and Truhlar,10 we formulate TST for n as follows. For a system of N particles (one ion and N − 1 water molecules) in a canonical ensemble at a temperature T, the probability density Ω(n) associated with a generalized coordinate ξ having a value of n can be expressed in terms of the system Hamiltonian H: Ω(n) =

2

where θ(x) is a Heaviside step function (defined earlier). Qr is the reactant (R) partition function:

Here, ri is the distance between the ion and the center of mass of the ith water molecule. θ(x) is a Heaviside step function: θ(x) = 1 for x > 0 and θ(x) = 0 for x < 0. Therefore, only the water molecules inside the first solvation shell contribute to eq 1. An alternative form for n can be expressed in terms of a smooth function that allows continuous transitions of water molecules between the first and second solvation shells:

∑

1 ⎛ ∂n ⎞ ⎜ ⎟ mk ⎝ ∂xk ⎠

mk denotes the mass associated with the Cartesian coordinate xk. The first two terms in eq 5 represent kinetic energy terms, while the last term is the potential energy term. Given that the transition state is n = n† between the reactant and product states along the reaction coordinate n, the TST rate constant can be expressed as10

(1)

i

n=

∑ 1

THEORY AND SIMULATION Transition State Theory. To define coordination number (n) mathematically, let us consider the PMF W(r) as a function of the ion−water distance, r. If r† is the position of the transition state between the reactant (first solvation shell) and the product (second solvation shell) on W(r), n is defined as44

∑ θ(r † − ri)

(5)

where

■

n=

1 1 Znpn2 + pqTZqpq + V (n , q) 2 2

N−2 N−1

∑ ∑ i=1 j=i+1

−1

⎤ (rixr jx + riyr jy + rizr jz)⎥ ⎥⎦ rri j

f i′ f ′j

= [Λauto + Λcross] (12)

where f i′ =

(4)

df (ri) dri

is the derivative of the function f(ri) (see eq 2

for f(ri)) with respect to the distance between the ion and the ith water molecule, as given in eq 13. f ′j represents the same for the jth water molecule. Λauto and Λcross denote the auto- and cross-terms, respectively. rxi , ryi , and rzi denote the x, y, and z components of the distance vector between the ion and ith water molecule, whereas rxj , ryj , and rzj denote the same between

Upon a coordinate transformation from the Cartesian coordinates (x, p) to our coordinates of interest ((n, q), (pn, pq)), where q and corresponding conjugate momenta pq have 3N − 1 components, respectively, it is possible to choose a basis q such that the Hamiltonian takes the following form with no cross term involving pn and pq:45 B


Article

The Journal of Physical Chemistry C the ion and the jth water molecule. ri and rj are the magnitudes of these respective vectors. a

f i′ =

⎡

r† ⎢ ⎣

ri

2a − 1

()

2

r†

−

a−1

ri

()

⎡ ⎢⎣1 −

−

r†

ri

3a − 1 ⎤

() r†

⎥⎦

2a ⎤ 2

() r†

ri

every 100 fs. Additional 5 ns MD simulations using the time step of 0.1 fs were performed to determine Λ, where snapshots were saved every 5 fs. Unbiased MD simulations will usually not allow access to all possible metastable coordination states because of the size of the barriers and the length of the simulation. Therefore, to compute W(n), it was necessary to utilize umbrella sampling56 in conjunction with MD to ensure converged sampling of the coordination number space over a particular range. We patched PLUMED57 with GROMACS to perform the umbrella sampling. Coordination number was biased to sample from n = 3.5 to n = 11.5 with 80 windows between them, making sure that distributions of coordination number in these windows are overlapping with each other. The force constant for the biasing potential was chosen to be 1000 kJ mol−1 nm−2. WHAM58 was applied to the time series of coordination number obtained in all simulation windows, which resulted in W(n). Every umbrella sampling simulation was performed for 5 ns in canonical ensemble at 298 K. For every reactive flux calculation (eq 14), 4000 initial configurations in the vicinity of n† were generated. Then MD simulations in microcanonical ensemble were performed using a MD time step of 0.1 fs in the forward time direction to a time length of 1 ps. Three different sets of initial velocities (generated using Maxwell−Boltzmann distribution at 298 K) were used. κRF (t) was obtained by averaging over these sets.

⎥⎦

(13)

Λ is determined from the summation over all possible distances; therefore, it can be treated as a constant (i.e., a distance-dependent function provides a constant when integrated over the entire distance space). However, in thermal equilibrium, because of the fluctuations of the locations of the ion and water molecules, Λ also fluctuates. One can obtain a time series of Λauto, Λcross, and Λ from an MD trajectory, which can be used to determine their distributions. These distributions allow one to determine most likely values of Λ (from the widths of the distributions) to be utilized in eq 11 for calculating TST rates with error bars. Reactive Flux Method. TST assumes that once a system reaches the transition state from the reactant state, it immediately goes to the product state. However, there is a certain probability of barrier recrossing due to nonequilibrium solvent effects through the dynamical coupling of the solvent to the reaction coordinate. The transmission coefficient (κ) is the dimensionless quantity that contains the fraction of total attempts that reach the product state starting from the transition state. This dynamical correction can be incorporated to modify the TST rate through the relation κkTST. To determine κ, we opt for the RF method as given in eq 14.17−19 κnRF(t ) =

■

RESULTS AND DISCUSSION To compute n and W(n), a priori knowledge of W(r) is required. We determined W(r) by computing the ion−water radial distribution function g(r): W(r) = −kBT ln(g(r)) (see Figure S1 ). Consistent with previous studies,6,30,32 the barrier (W(r†)) between the reactant and the product decreases as ionic size increases, namely, ion−water binding gets weaker, which may result in faster water exchange between the first and second solvation shells. We extracted the values of r† for all the ion−water systems to determine n and W(n). W(n) computed by employing umbrella sampling and WHAM56−58 are depicted in Figure 1. For simplicity, we present and discuss the salient features of W(n) for five ions (A, C, and D) in the ascending order of ionic size. A1, C1, and D1 are the smallest ions and have identical sizes. However, they have very different W(n), which can be attributed to their charge differences. For A1, the most likely coordination state is n = 6, where the less likely coordination state n = 5 is accessible within 5.0 kcal/mol. The lower coordination state n = 4 is observed for C1, in addition to and within 3 kcal of n = 5 and n = 6, where n = 5 is most probable. For D1, higher coordination states n = 7 and n = 8 are found, which are within 8 kcal of the most probable state n = 6. With the increase of ionic size, higher coordination states become accessible for all ions, along with barrier reduction between the states. Coordination states for cations remain more stable and structured than anions, as clearly indicated by W(n) for A5, C5, and D5. The PMFs computed for all the model ions with different sizes generate a map W(σ, n) as depicted in Figure 2 (top panel). W(n) for a real ion can be extracted from this map if the corresponding Lennard-Jones parameter σ is known. To examine the quality of the extracted W(n) for real ions, we considered Na+, Cl−, and Ca2+ as test cases. We found an excellent match between W(n) computed using umbrella sampling MD simulation and that extracted from W(σ, n) for

⟨vn(0) δ(n − n†) θ(n(t ) − n†)⟩ ⟨vn(0) δ(n − n†) θ(vn(0))⟩

(14)

Here, n(t) represents the time evolution of coordination number, given that the starting location is n = n† and the dn starting velocity is vn(0) = dt |t = 0 . θ(x) is a Heaviside step function (defined earlier). The RF method determines the exact value of κRF from MD trajectories, provided ample n statistics are accounted for. MD Simulations. To apply the rate theory, we performed MD simulations of 10 monovalent anions denoted by Ai; i = 1− 10, 10 monovalent cations denoted by Ci; i = 1−10, and 10 divalent cations denoted by Di; i = 1−10. Increasing i indicates the increasing of the ionic size. We chose the Lennard-Jones parameter σ to represent the ionic size that is varied from 0.2 to 0.65 nm with an increment of 0.05 nm to generate these model ions. The other Lennard-Jones parameter, ϵ, for all ions was obtained from the exponential form of ϵ(σ) used in the OPLS force field.46−49 The masses, m(σ), were also obtained in a similar fashion. A complete list of ion parameters can be found in the Supporting Information. All MD simulations were performed using GROMACS MD.50 Every ion−water system comprises a box of 880 TIP4P water molecules51 and the ion. Each system was equilibrated for 10 ns in isothermal−isobaric ensemble at 1 atm and 298 K using Berendsen’s method;52 then, 100 ns MD trajectories were obtained for production run in canonical ensemble using the Nose−Hoover thermostat.53,54 These MD trajectories were used to calculate g(r), W(r), and r†. Note that an MD time step of 2 fs and Lennard-Jones and Coulomb potential cutoffs of 1.0 nm were used. The long-range electrostatics were treated with the particle-mesh Ewald method.55 MD snapshots were saved C


Article

The Journal of Physical Chemistry C

the Na+ and Cl− cases, while the agreement was reasonably good for the Ca2+ case. Experiments40 and simulations59−61 reported a trend for the average coordination numbers of these ions: Coordination number increases in the order Na+, Cl−, and Ca2+. This is consistent with our finding that the most probable coordination numbers are 6, 7, and 8 for Na+, Cl−, and Ca2+, respectively. The solvation structures corresponding to these coordination states are presented in Figure 2 (bottom panel). These case studies establish that the map can be used as a library, where one can simply find W(n) for real ions and avoid extensive calculations. Prior to calculating the transition rates between the metastable states on W(n), we need to determine the distributions of Λ. It makes sense to calculate transition rates for those ions that have at least two metastable states separated by a nonvanishing barrier. Therefore, we determine Λ for the following relevant ions only: A1−A4, C1−C6, and D1−D8. Λ depends on the derivative f ′ (see eqs 12 and 13), which has the maximum at the transition state (r†) between the reactant and product states along the reaction coordinate r, i.e., at the boundary between the first and second solvation shells, as depicted in Figure 3a,b for A1 and A4. f ′ decays for water

Figure 1. Potential of mean force as a function of the coordination number of monovalent anions (A) and monovalent (C) and divalent cations (D).

Figure 3. PMF for A1 and A4 along the reaction coordinate r to represent the boundary (r†) between the first and second solvation shells, indicated by the vertical dotted-lines (a). The derivative df as dr

given in eq 13 to show that water molecules at and in the proximity of r† contributes the most to Λ (b). Distributions of Λ indicated by black for A1 and by red for A4, while the corresponding distributions for Λauto and Λcross are indicated by green and blue, respectively (c, d).

molecules getting away from the boundary, and it is vanishingly small near the equilibrium ion−water distance in the first and second solvation shells. Therefore, water molecules at and in the neighborhood of r† contribute the most to Λ. In Figure 3c,d, we present the distributions P(Λ) for Λauto, Λcross, and Λ for A1 and A4, while for other ions these are presented in the Supporting Information (Figures S2−S4). Λcross has vanishingly small values with a slight bias toward negative magnitudes, resulting in small cancelations of Λauto. This is because the derivatives f i′ and f j′ are localized at the molecules i and j; therefore, their overlap, i.e., f i′f j′ in the crossterm Λcross of eq 12, should be close to zero. It is interesting to observe that, for smaller ions Λ is distributed around a single peak, whereas for larger ions P(Λ) has a two-peak structure. For smaller ions, it is less likely to find water molecules near r† because of sharp and large barriers; therefore, most of the

Figure 2. Potential of mean force W(σ, n) which can be used as a map to obtain W(n) for a real ion with known σ (top panel); ten contours are places between 0 (blue) and 10 kcal (red). W(n) obtained using the map is compared with that obtained from umbrella sampling MD simulation for real ions (middle panel). The solvation structure of the most probable coordination state of Cl−, Na+, and Ca2+ (bottom panel).

D


Article

The Journal of Physical Chemistry C contributions to Λ are coming from the water molecules near the reactant and product states. Because these contributions are smaller, P(Λ) are peaked at smaller values with long decaying tails for higher Λ. For larger ions, in addition to water molecules near the reactant and product states, one can also find water molecules in the proximity of r† due to smaller barriers, which results in the two-peak structure of P(Λ). Note that A3, A4, C5, C6, and D5−D8 have such two-peak like structures in their P(Λ). Because Λ is a part of the prefactor in the TST rate expression (eq 11), it establishes the fact that changing a coordination state of an ion by moving a water molecule between the first and the second solvation shells is slower (indicated by small Λ) than by moving a water molecule from the boundary to these shells (indicated by large Λ). Now, for TST rate calculations, we choose the Λ values for which P(Λ) = 1/2, so that we can assign a range to TST rates due to the most likely range of Λ. For anions, TST time scales, τi→j = 1/kTST, where i and j denote two metastable states, are listed in Table 1 (these are for forward transitions; see the

Table 3. Transition Time Scale from the Coordination States i to the Coordination State j for Divalent Cations Using TST

τ5→6 (ps)

A1 A2 A3 A4

6−14 ≤1

τ6→7 (ps) 13−34 1−4 ≤1

τ4→5 (ps)

τ5→6 (ps)

33−103 4−16 1−3

75−234 14−51 6−13 2−6 1−2

1400−2572 2−4

τ8→9 (ps)

9−15 2−5

432−783 19−42 4−8 1−2

τ9→10 (ps)

1023−2279 39−87 5−10 2−4 1−2

τ10→11 (ps)

37−79 6−15 2−6 1−2

(15)

where ⟨...⟩ denotes the ensemble average and t0 is an arbitrary initial time. Pn,eq is the equilibrium survival probability. It is assumed that C(t) can be expressed in a single exponential form going to zero at t → ∞; then, the lifetime of a coordination state can be obtained by integrating C(t): τc =

∫0

∞

C(t )dt =

∫0

∞

exp( −t /τc)dt

(16)

It is obvious that τc ∼ τi→j provided i and j are the only two coordination states of an ion−water system. It is clear from the PMFs presented in Figure 1 that ions other than A1 have more than two coordination states, which may result in C(t) taking on a multiexponential form for them. A1 has the coordination states n = 5 and n = 6. C(t) calculated for the n = 5 from a 25 ns MD trajectory was well approximated by a single exponential form (eq 16), as depicted in Figure 4. The corresponding τc is ∼1.3 ps, which is slightly under-

Table 2. Transition Time Scale from the Coordination States i to the Coordination State j for Monovalent Cations Using TST C1 C2 C3 C4 C5 C6

D1 D2 D3 D4 D5 D6 D7 D8

τ7→8 (ps)

C(t ) = ⟨δPn(t + t0)δPn(t0)⟩/⟨δPn(t0)δPn(t0)⟩

Supporting Information for backward transitions). τ5→6 and τ6→7 are on picosecond, subpicosecond, and tens of picoseconds time scales, and the transitions between the states get faster with the increase of ionic size. For A3 and A4, subpicosecond time scales are obtained if the right peak on P(Λ) is used. τi→j for monovalent cations are on the time scale of picoseconds, tens of picoseconds, and hundreds of picoseconds (Table 2). For C6 and larger cations, transitions are possible

ion

τ6→7 (ps)

peak of P(Λ) can not be distinctly separated from the TST time scales using the right peak on P(Λ). Therefore, using just the left peak will suffice. To examine and compare TST time scales, we followed the method suggested by Impey, Madden, and McDonald (IMM).62 According to the IMM method, the residence time of a water molecule in the first solvation shell of an ion can be calculated using the survival probability of the water molecule in that shell. Analogously, we assign the survival probability of a coordination state of the ion−water system Pn to 1 when the system is found in the designated coordination state, otherwise Pn is 0. The normalized time correlation function for the fluctuation δPn(t) = Pn(t) − Pn,eq is as follows:

Table 1. Transition Time Scale from the Coordination State i to the Coordination State j for Anions Using TST ion

ion

τ6→7 (ps)

87−208 13−40 3−10 1−5

only between n = 6 and n = 7, and the corresponding time scales are picoseconds or subpicosecond. Similar to anions, subpicosecond time scales are found for C5 and C6 if the right peak on P(Λ) is used. For divalent cations (Table 3), τi→j are found to be on nanosecond (ns), picosecond, and tens and hundreds of picoseconds time scales. For D8 and larger divalent cations, only n = 10 and n = 11 participate in interconversion with time scales close to 1 ps. For D5−D8, picosecond and subpicosecond time scales are found if the right peak on P(Λ) is used. Here we learn that, for any large ion, solvent-exchange is so fast (picosecond/subpicosecond), the TST time scales using the left

Figure 4. Time correlation function (black) for the survival probability for the coordination state n = 5 of A1. The red line indicates a single exponential fit to this correlation function. E


Article

The Journal of Physical Chemistry C

small values indicate significant barrier recrossing in coordination number space. The corrected time scales are then τ5→6 = 60−140 ps for A1, τ4→5 = 550−1717 ps for C1, and τ6→7 = 5.4− 9.9 ns for D1, which are 10, 16.6, and 3.8 times slower than the TST rates, respectively. Previous studies of water exchange dynamics around monovalent anions (Cl−, Br−, and I−),30 monovalent cations (Na+ and Li+),21,29 and divalent cations (Be2+ and Mg2+),34 utilizing the ion−water distance as the reaction coordinate, reported similarly small transmission coefficients (κRF n ∼ 0.05−0.2) that support significant barrier recrossing. In Figure 1, we find that most of the coordination states are separated by small barriers. The TST time scales listed in Tables 1 and 2 along with our reactive flux calculations that determine similarly small transmission coefficients suggest that the dynamical time scale separation for most of the coordination states are tenuous for monovalent anions and cations. However, because of the slightly larger transmission coefficient and distinct TST time scales (ranges from picoseconds to nanoseconds, see Table 3), a divalent cation’s coordination states can be dynamically distinguishable and separable. Note that due to this weak time scale separation problem it will be experimentally challenging to identify and separate different coordination states in terms of their kinetics; as aforementioned, experiments may provide only information about average coordination states. The smallness of the transmission coefficients can also suggest that coordination number could be a suboptimal reaction coordinate. To understand this, the committor probability distribution is a tool that should be employed to examine the quality of this reaction coordinate.65,66 If a system is initially located at the dividing surface between two metastable states along a reaction coordinate, a short time evolution will lead the system to commit to one of these two states. Now, if the probability of this commitment has a distribution peaked around 1 , it will be suggestive that the 2 reaction coordinate is an optimal choice for the reaction of interest. The committor67,68 for coordination number can be calculated by determining the recrossing factor Pn for the transition from ni to nf with the dividing surface n† as follows:

estimated compared to the TST time scale (τ5→6). In the future, it will be interesting to examine the longtime behavior of multistate ions in coordination space by employing methods such as the Markov state model.63 As previously described, TST provides only a qualitative picture of the total rate, and the nonequilibrium solvent effects that determine the transmission coefficient are necessary to complete the understanding of this rare event. Using the RF method (eq 14), we determined κRF n (t) for the smallest three model ions to examine how rapidly solvent dynamics can affect τ5→6 for A1, τ4→5 for C1, and τ6→7 for D1. To examine how the choice of initial coordination numbers and their velocities at the transition state n† can affect κRF n (t), we consider two cases: The first one (case A) accounts for all configurations in the vicinity of n† and the second one (case B) accounts for the subensemble of these configurations on the right side of n† as depicted in Figure 5 (top panel). It turns out that, κRF n (t) at t ∼

Figure 5. Initial coordination numbers (n(0)) and their velocities (ṅ(0)) at the transition state n† (top panel, red-dashed lines indicate exact positions of n† and black dots indicate all configurations at the transition state) used in the RF calculations for the n = 5 → n = 6, n = 4 → n = 5, and n = 6 → n = 7 transitions for A1, C1, and D1, respectively. The black and red lines in the middle panel denote κRF n (t) calculated for case A and case B, respectively. The committor probability distribution for finding A1 at n = 5, C1 at n = 4, and D1 at n = 6 after 200 fs, if their initial locations are the transition state configurations indicated by all black dots in the top panel (bottom panel).

Pn =

∫ hn(t′; n† , Xsol) Ptherm(pn )dpn

(17)

Here, hn(t′;n†, Xsol) is assigned to 1 if the coordination state is found to be ni at time t′ = 200 fs; otherwise, hn(t′; n†, Xsol) = 0. Xsol denotes the solvent bath coordinate. The integration represents averaging over initial coordination number momenta pn distributed thermally over Ptherm(pn). We calculated the distributions of Pn for the n = 5 → n = 6, n = 4 → n = 5, and n = 6 → n = 7 transitions for A1, C1, and D1, respectively, utilizing the trajectories used for the RF calculations. Pn determined from every 40 trajectories were used to calculate the distribution P(Pn). As depicted in Figure 5 (bottom panel), the distribution P(Pn) is centered around 0.3, 0.6, and 0.5 for A1, C1, and D1, respectively. This is an indication of a good reaction coordinate, especially for the D1 case, where P(Pn) is a narrowed distribution peaked at 1 and the transmission coefficient is 2 2.5−4.5 times greater than that for the A1 and C1 cases. To some extent, this is consistent with our recent findings68 and the argument by Mullen, Shea, and Peters66 that it is both the

0 are affected because of the different choices of initial coordination numbers and their velocities. κRF n (t) should ideally be close to 1 at t ∼ 0;17−19 case B reproduces that to some extent, whereas case A fails. This suggests that the transition state should be optimized (by employing methods such as variational TST13 and inertial likelihood maximization64) to get the appropriate κRF n (t ∼ 0). However, we are interested in the RF time-dependent behavior of κRF n (t), and for longer times, κn (t) decays are identical for both case A and case B. This insensitivity of κ nRF (t) to the choice of the starting configurations for the RF calculations suggests that the transition state ensemble will not be realized by simple translations in coordination number space. κRF n (t) displays two components as shown in Figure 5 (middle panel): a fast rapid decay up to t ∼ 0.2 ps followed by a plateau (for A1 and C1) or an asymptotic slow decay (for D1). Upon averaging over the last 0.8 ps, the actual κRF n is found to be 0.1, 0.06, and 0.26 for A1, C1, and D1, respectively. These F


The Journal of Physical Chemistry C maximization of κRF in addition to P(Pn) ∼ satisfied for a good reaction coordinate.

1 2

■

ACKNOWLEDGMENTS We gratefully acknowledge Liem Dang and Panos Stinis for useful discussions. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC0205CH11231. S.R., C.J.M., and G.K.S. were supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences. M.D.B. was supported by MS3 (Materials Synthesis and Simulation Across Scales) Initiative, a Laboratory Directed Research and Development Program at Pacific Northwest National Laboratory (PNNL). PNNL is a multiprogram national laboratory operated by Battelle for the U.S. Department of Energy.

that needs to be

■

SUMMARY AND CONCLUSION We demonstrated that rate theory in coordination number space can be applied to ion solvation in water. By studying a series of model anions and cations, we generated maps to determine the PMF in coordination space for real ions with known size and charge. We found that for a given size and charge, the first solvation shell of anions is dramatically less structured than that of the corresponding cation. These results once again showcase the subtle asymmetries of aqueous solvation of cations and anions. In coordination space, the TST rates between different metastable states on the PMFs becomes faster with increasing ionic size. We identified a sensitivity of these rates to Λ (a term in the prefactor in the TST rate equation), which is unique for every ion, and it mostly depends on the contributions from the water molecules near the boundary between the first and second solvation shells. On the basis of the TST calculations and RF analysis, we found that, while we could distinctly separate the transition time scales between the coordination states for divalent cations, such dynamical separation is weak for the monovalent cations and anions. It will be interesting to investigate in the future how a more accurate water model that includes polarization effects6,69,70 or ab initio MD71 may influence our findings because these methods are expected to affect the PMF in coordination space. The essential factor that influences solvent exchange rates is the strong coupling of the reaction coordinate to solvent fluctuations as determined by the small κRF. Future studies will be focused on quantifying this strong coupling by examining the spectral density of the solvent friction and its coupling to the barrier frequency. This theoretical study opens a new direction for investigating novel reaction coordinates to ionsolvation. The goal of simultaneously maximizing the reactive 1 flux with a committor distribution (P(Pn) ∼ 2 ) will likely require an investigation of a multidimension reaction coordinate that includes both the ion−water distance in addition to the coordination number. This will allow us to unravel mechanisms and identify multidimensional pathways of ion solvation and ion-pairing processes relevant to a variety of collective phenomena ranging from nucleation to chemical transformation in complex biological environments.

■

■

REFERENCES

(1) Pollak, E.; Talkner, P. Reaction rate theory: What it was, where is it today, and where is it going? Chaos 2005, 15, 026116. (2) Abad, E.; Reingruber, J.; Sansom, M. S. P. On a novel rate theory for transport in narrow ion channels and its application to the study of flux optimization via geometric effects. J. Chem. Phys. 2009, 130, 085101. (3) Koneshan, S.; Rasaiah, J. C.; Lynden-Bell, R. M.; Lee, S. H. Solvent structure, dynamics and ion mobility in aqueous solutions at 25 °C. J. Phys. Chem. B 1998, 102, 4193. (4) Frank, S.; Schmickler, W. Ion transfer across liquid-liquid interfaces from transition-state theory and stochastic molecular dynamics simulations. J. Electroanal. Chem. 2006, 590, 138. (5) Kim, H. J.; Hynes, J. T. A theoretical model for SN1 ionic dissociation in solution. 2. Nonequilibrium solvation reaction path and reaction Rate. J. Am. Chem. Soc. 1992, 114, 10528. (6) Annapureddy, H. V. R.; Dang, L. X. Understanding the rates and molecular mechanism of water exchange around aqueous ions using molecular simulations. J. Phys. Chem. B 2014, 118, 8917. (7) Hinsen, K.; Roux, B. Potential of mean force and reaction rates for proton transfer in acetylacetone. J. Chem. Phys. 1997, 106, 3567. (8) Guárdia, E.; Rey, R.; Padró, J. A. Potential of mean force by constraint molecular dynamics: A sodium chloride ion-pair in water. Chem. Phys. 1991, 155, 187. (9) Grote, R. F.; Hynes, J. T. The stable states picture of chemical reactions. II. Rate constants for condensed and gas phase reaction models. J. Chem. Phys. 1980, 73, 2715. (10) Schenter, G.; Garrett, B. C.; Truhlar, D. G. Generalized transition state theory in terms of the potential of mean force. J. Chem. Phys. 2003, 119, 5828. (11) Wigner, E. Ü ber das überschreiten von potentialschwellen bei chemischen reaktionen. Z. Phys. Chem. Abt. B 1932, 19, 203. (12) Eyring, H. The activated complex in chemical reactions. J. Chem. Phys. 1935, 3, 107. (13) Truhlar, D. G.; Garrett, B. C. Variational transition state theory. Annu. Rev. Phys. Chem. 1984, 35, 159. (14) Kramers, H. Brownian motion in a field of force and the diffusion model of chemical reactions. Physica 1940, 7, 284. (15) Pollak, E. Theory of activated rate processes: A new derivation of Kramers’ expression. J. Chem. Phys. 1986, 85, 865. (16) Bennett, C. H. Algorithms for Chemical Computations; ACS Symposium Series; American Chemical Society: Washington, DC, 1977. (17) Chandler, D. Statistical mechanics of isomerization dynamics in liquids and the transition state approximation. J. Chem. Phys. 1978, 68, 2959. (18) Rey, R.; Guardia, E. Dynamic aspects of the Na+-Cl− ion-pair association in water. J. Phys. Chem. 1992, 96, 4712. (19) Ciccotti, G.; Ferrario, M.; Hynes, J. T.; Kapral, R. Dynamics of ion-pair interconversion in a polar-solvent. J. Chem. Phys. 1990, 93, 7137.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b00443. Derivation of TST formulation in coordination number space; force field parameters for MD simulations; and W(r), P(Λ), and backward transition rates between metastable states in coordination space (PDF)

■

Article

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest. G


Article

The Journal of Physical Chemistry C (20) Rey, R.; Hynes, J. T. Hydration shell exchange dynamics for Na+ in water. J. Phys.: Condens. Matter 1996, 8, 9411. (21) Rey, R.; Hynes, J. T. Hydration shell exchange kinetics: An MD study for Na+ (aq). J. Phys. Chem. 1996, 100, 5611. (22) Spångberg, D.; Rey, R.; Hynes, J. T.; Hermansson, K. Rate and mechanism for water exchange around Li+(aq) from MD simulations. J. Phys. Chem. B 2003, 107, 4470. (23) Laage, D.; Hynes, J. T. On the residence time for water in a solute hydration shell: Application to aqueous halide solutions. J. Phys. Chem. B 2008, 112, 7697. (24) Møller, K. B.; Rey, R.; Masia, M.; Hynes, J. T. On the coupling between molecular diffusion and solvation shell exchange. J. Chem. Phys. 2005, 122, 114508. (25) Hermansson, K.; Wojcik, M. Water exchange around Li+ and Na+ in LiCl (aq) and NaCl (aq) from MD Simulations. J. Phys. Chem. B 1998, 102, 6089. (26) Spangberg, D.; Wojcik, M.; Hermansson, K. Pressure dependence and activation volume for the water exchange mechanism in NaCl (aq) from MD simulations. Chem. Phys. Lett. 1997, 276, 114. (27) Buchner, R.; Capewell, S. G.; Hefter, G.; May, P. M. Ion-pair and solvent relaxation processes in aqueous Na2SO4 solutions. J. Phys. Chem. B 1999, 103, 1185. (28) Buchner, R.; Chen, T.; Hefter, G. Complexity in “simple” electrolyte solutions: Ion pairing in Mg2SO4 (aq). J. Phys. Chem. B 2004, 108, 2365. (29) Dang, L. X.; Annapureddy, H. V. R. Computational studies of water exchange around aqueous Li+ with polarizable potential models. J. Chem. Phys. 2013, 139, 084506. (30) Annapureddy, H. V. R.; Dang, L. X. Water exchange rates and molecular mechanism around aqueous halide ions. J. Phys. Chem. B 2014, 118, 7886. (31) Wilkins, D. M.; Manolopoulos, D. E.; Dang, L. X. Nuclear quantum effects in water exchange around lithium and fluoride ions. J. Chem. Phys. 2015, 142, 064509. (32) Roy, S.; Dang, L. X. Computer simulation of methanol exchange dynamics around cations and anions. J. Phys. Chem. B 2015, 120, 1440. (33) Roy, S.; Dang, L. X. Water exchange dynamics around H3O+ and OH− ions. Chem. Phys. Lett. 2015, 628, 30. (34) Dang, L. X. Computational studies of water-exchange rates around aqueous Mg2+ and Be2+. J. Phys. Chem. C 2014, 118, 29028. (35) Yonetani, Y. Distinct dissociation kinetics between ion pairs: Solvent-coordinate free-energy landscape analysis. J. Chem. Phys. 2015, 143, 044506. (36) Schenter, G. K.; Garrett, B. C.; Truhlar, D. G. The role of collective solvent coordinates and nonequilibrium solvation in chargetransfer reactions. J. Phys. Chem. B 2001, 105, 9672. (37) Marcus, Y. Effect of ions on the structure of water: Structure making and breaking. Chem. Rev. 2009, 109, 1346. (38) Fulton, J. L.; Heald, S. M.; Badyal, Y. S.; Simonson, J. M. Understanding the effects of concentration on the solvation structure of Ca2+ in aqueous solution. I: The perspective on local structure from EXAFS and XANES. J. Phys. Chem. A 2003, 107, 4688. (39) Lincoln, S. F. Solvent coordination numbers of metal ions in solution. Coord. Chem. Rev. 1971, 6, 309. (40) Enderby, J. E. Ion solvation via neutron scattering. Chem. Soc. Rev. 1995, 24, 159. (41) Fratiello, A.; Lee, R. E.; Nishida, V. M.; Schuster, R. E. Proton magnetic resonance coordination number study of Al(III), Be(II), Ga(III), In(III), and Mg(II) in water and aqueous solvent mixtures. J. Chem. Phys. 1968, 48, 3705. (42) Carof, A.; Salanne, M.; Charpentier, T.; Rotenberg, B. On the microscopic fluctuations driving the NMR relaxation of quadrupolar ions in water. J. Chem. Phys. 2015, 143, 194504. (43) Raiteri, P.; Demichelis, R.; Gale, J. D. Thermodynamically consistent force field for molecular dynamics simulations of alkalineearth carbonates and their aqueous speciation. J. Phys. Chem. C 2015, 119, 24447. (44) Sprik, M. Coordination numbers as reaction coordinates in constrained molecular dynamics. Faraday Discuss. 1998, 110, 437.

(45) Darve, E.; Pohorille, A. Calculating free energies using average force. J. Chem. Phys. 2001, 115, 9169. (46) Chandrasekhar, J.; Spellmeyer, D. C.; Jorgensen, W. L. Energy component analysis for dilute aqueous solutions of Li+, Na+, F−, and Cl− ions. J. Am. Chem. Soc. 1984, 106, 903. (47) Lybrand, T. P.; Ghosh, I.; McCammon, J. A. Hydration of chloride and bromide anions: Determination of relative free energy by computer simulation. J. Am. Chem. Soc. 1985, 107, 7793. (48) McDonald, N. A.; Duffy, E. M.; Jorgensen, W. L. Monte carlo investigations of selective anion complexation by a Bis(phenylurea) ptert-Butylcalix[4]arene. J. Am. Chem. Soc. 1998, 120, 5104. (49) Aaqvist, J. Ion-water interaction potentials derived from free energy perturbation simulations. J. Phys. Chem. 1990, 94, 8021. (50) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 4, 435. (51) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926. (52) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684. (53) Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 1984, 52, 255. (54) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695. (55) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577. (56) Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23, 187. (57) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A portable plugin for free-energy calculations with molecular dynamics. Comput. Phys. Commun. 2009, 180, 1961. (58) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The weighted histogram analysis method for freeenergy calculations on biomolecules. I. The method. J. Comput. Chem. 1992, 13, 1011. (59) Bounds, D. G. A molecular dynamics study of the structure of water around the ions Li+, Na+, K+, Ca+2, Ni+2 and Cl−. Mol. Phys. 1985, 54, 1335. (60) Koneshan, S.; Rasaiah, J. C.; Lynden-Bell, R. M.; Lee, S. H. Solvent structure, dynamics, and ion mobility in aqueous solutions at 25 °C. J. Phys. Chem. B 1998, 102, 4193. (61) Khalack, J. M.; Lyubartsev, A. P. Car-Parrinello molecular dynamics simulations of Na+-Cl− ion pair in liquid water. Condens. Matter Phys. 2004, 7, 683. (62) Impey, R. W.; Madden, P. A.; McDonald, I. R. Hydration and mobility of ions in solution. J. Phys. Chem. 1983, 87, 5071. (63) Bowman, G. R.; Pande, V. S.; Noé, F. An Introduction to Markov State Models and Their Application to Long Timescale Molecular Simulation; Springer: Dordrecht, The Netherlands, 2014. (64) Peters, B. Inertial likelihood maximization for reaction coordinates with high transmission coefficients. Chem. Phys. Lett. 2012, 554, 248. (65) Dellago, C.; Bolhuis, P. G.; Geissler, P. L. Transition path sampling. Advances in Chemical Physics 2002, 123, 1. (66) Mullen, R. G.; Shea, J. E.; Peters, B. Transmission coefficients, committors, and solvent coordinates in ion-pair dissociation. J. Chem. Theory Comput. 2014, 10, 659. (67) Truhlar, D. G.; Garrett, B. C. Multidimensional transition state theory and the validity of Grote-Hynes theory. J. Phys. Chem. B 2000, 104, 1069. (68) Pluhařová, E.; Baer, M. D.; Schenter, G. K.; Jungwirth, P.; Mundy, C. J. Dependence of the rate of LiF ion-pairing on the description of molecular interaction. J. Phys. Chem. B 2016, 120, 1749. H


Article

The Journal of Physical Chemistry C (69) Dang, L. X.; Rice, J. E.; Caldwell, J.; Kollman, P. A. Ion solvation in polarizable water: molecular dynamics simulation. J. Am. Chem. Soc. 1991, 113, 2481. (70) Dang, L. X.; Chang, T.-M. Molecular dynamics study of water clusters, liquids, and liquid-vapor interface of water with many-body potentials. J. Chem. Phys. 1997, 106, 8149. (71) Timko, J.; De Castro, A.; Kuyucak, S. Ab initio calculation of the potential of mean force for dissociation of aqueous Ca-Cl. J. Chem. Phys. 2011, 134, 204510.

I


Reaction Rate Theory in Coordination Number Space: An Application

Recommend Documents