General approaches for directly calculating the temperature dependence of dynamical quantities from simulations at a single temperature are presented. The method is demonstrated for self-diffusion and OH reorientation in liquid water. For quantities which possess an activation energy, e.g., the diffusion coefficient and the reorientation time, the results from the direct calculation are in excellent agreement with those obtained from an Arrhenius plot. However, additional information is obtained, including the decomposition of the contributions to the activation energy. These results are discussed along with prospects for additional applications of the direct approach.
I. INTRODUCTION
Molecular dynamics (MD) simulations are a ubiquitous tool for understanding the dynamics of chemical systems. In particular, they can be directly connected with experimental measurements by, for example, using time correlation functions (TCFs) to obtain observables ranging from kinetic rate constants to vibrational spectra. Inevitably, the data generated in a MD simulation dwarf the information produced by such analyses. It is thus important to seek methods by which the trajectory data can be examined to take advantage of the significant information content that is generally discarded.
Reaction rate constants, transport coefficients, and other important dynamical time scales for molecular systems are frequently obtained through the calculation of TCFs.1–3 For example, the rate constant for a chemical reaction can be determined from the flux-side TCF, Cfs(t), as4–7
where s defines the dividing surface between reactants () and products (), θ(x) is the Heaviside step function, and is the flux through the dividing surface with the velocity along s. The activation energy of the reaction,
where β = 1/kBT, is a quantity that is of significant interest due to the insight it provides into the effective barrier for the reaction. It is typically calculated by constructing an Arrhenius plot of versus 1/T based on measurements (or calculations) at multiple temperatures. In the context of Eq. (1), however, Ea is related to the temperature-dependence of a TCF such as Cfs(t). Indeed, it has been shown by Dellago and Bolhuis8 that this perspective can be used to calculate Ea directly from transition path sampling simulations at a single temperature,9–11 and a related approach has been demonstrated by Morita and co-workers for vibrational spectra.12,13 We have recently generalized this beyond transition path sampling and to other TCFs from which rate constants can be obtained, including those based on a quantum description.14
The measurement and calculation of activation energies by an Arrhenius analysis can raise significant issues due to the need to consider multiple temperatures. For example, many systems have a structure that is highly dependent on temperature, e.g., folded proteins, bilayer membranes or vesicles, and self-assembled structures. This can make the construction of an Arrhenius plot problematic by restricting the range of temperatures that can be considered. A similar issue arises even for bulk systems if one is interested in an activation energy at conditions that are near a phase transition. Yet it can be of significant interest to examine the activation energies of dynamical processes in the vicinity of these points of transformation.
It is also important to consider cases where quantities other than a rate constant are of interest, but for which an activation energy can still be measured and calculated. In this paper, we do just that by first showing a general expression for the temperature derivative of a TCF and then applying this result to two commonly considered dynamical properties: diffusion coefficients and reorientational time scales. The method is demonstrated by application to the relevant time correlation functions for these attributes in bulk liquid water. In addition, it is shown how deeper physical insight into these processes can be obtained than is possible from Arrhenius calculations.
II. TIME CORRELATION FUNCTIONS (TCFs) AND ACTIVATION ENERGIES
A. General expressions
A general time correlation function between two dynamical variables A and B in the canonical ensemble, , can be written more explicitly as
where Q is the partition function, H is the Hamiltonian, and Tr indicates an integration over phase space. In this correlation function expression, only Q and e−βH depend on the temperature. Thus, the derivative of CAB(t) with respect to β can be straightforwardly evaluated as
Noting the definition of the average energy, , this gives
where is the fluctuation in energy from its average value.
This shows that the temperature dependence of the TCF, and thus any associated time scales, can be determined by evaluating this simple time correlation function that is closely related to CAB(t) itself. Note that the interpretation of this result is intuitive as it relates the derivative to the (continuous) differences between CAB(t) when the system energy is initially greater than average () and when it is initially less than average (). In this sense, Eq. (5) is an analytical derivative instead of the numerical derivative that is obtained from an Arrhenius analysis.
As noted above, an analogous result has already been obtained in the context of reactive flux and other TCFs that can be directly related to rate constants.8,14 However, there are many contexts in which time scale or transport coefficients are calculated from a TCF and possess an activation energy, or even simply a dependence on temperature, that is of physical interest. In the following, we show how the general result in Eq. (5) can be used to evaluate such activation energies for the examples of self-diffusion and OH reorientation in liquid water.
B. Diffusion coefficients
We consider, as a first example, the application of this formulation to the activation energy of diffusion. The diffusion coefficient, D, can be calculated from the mean-squared-displacement (MSD) as
which becomes linear at long times such that
for motion in three dimensions. The activation energy associated with the diffusion coefficient is defined as
Then, taking the derivative of Eq. (7) with respect to β and dividing by D yields
Using Eqs. (5)–(7), this can be written as
where MSDH(t), defined by the relation above, is the MSD weighted by the energy fluctuation. Note that the limit of the ratio of the functions equals the ratio of the limits as long as both are well defined and the limit of the denominator is not zero; each condition is met here. Thus, this ratio of TCFs should approach a constant value at longer times that is equal to the diffusion coefficient activation energy.
C. Reorientation times
The same approach can also be applied to reorientation dynamics. Reorientation times, denoted , are typically calculated from the time decay of the -th order reorientational correlation function,
Here is the -th order Legendre polynomial and is a unit vector pointing along some molecular axis, e.g., the OH bond in a water molecule. The TCF, C2(t), is of particular interest as it can be directly measured by IR pump-probe anisotropy experiments,15–17 and an average rotational time defined as its integral can be obtained from NMR.16–20
In water and other hydrogen bonding (H-bonding) liquids, it has been shown that the C2(t) TCF is well described by a tri-exponential decay that distinguishes the time scales associated with reorientation due to inertial (τiner), librational (τlib), and H-bond making and breaking (τ2) dynamics.21 Taking this form for the general, -th order case, the TCF can be written as22
where α = iner, lib, and (corresponding to the three time scales), and Aα represent the amplitudes of the contributions of the three components to the overall reorientation dynamics. Applying the result in Eq. (5) to gives
where is the reorientational TCF weighted by the energy fluctuation. If the same derivative with respect to β is taken in Eq. (12), an additional expression for results,
Then, Eq. (13) can be used to calculate from MD simulations, while Eq. (14) gives the form to which it can be fit [constrained by Aα and τα obtained by fitting itself] to determine the temperature dependence of the amplitudes and time scales and, in some cases, the activation energy of the latter. In particular, an activation energy of one of the time scales can be obtained as
This again provides a simple method for calculating an activation energy from a single-temperature MD simulation.
III. COMPUTATIONAL METHODS
In principle, the energy fluctuation TCFs, MSDH(t) and , and the associated activation energies can be calculated from a single NVT MD simulation, as we have previously illustrated for H-bond exchange TCFs.14 However, there is naturally some effect due to the thermostat used to maintain the temperature that needs to be minimized. Here we adopt a different, nonequilibrium MD approach to illustrate the method that avoids any effect of the thermostat. Specifically, we sample initial conditions for short constant energy (NVE) trajectories from a long NVT MD simulation as illustrated in Fig. 1. Each NVE trajectory has an initial energy, sampled from the canonical ensemble, that defines δH(0) that weights the TCF giving MSD(t) or C2(t).
A schematic diagram of the nonequilibrium MD simulation approach in which NVE trajectories with different energies are initiated from a single NVT trajectory.
A schematic diagram of the nonequilibrium MD simulation approach in which NVE trajectories with different energies are initiated from a single NVT trajectory.
For the results presented here, 2000 NVE trajectories of 20 ps each were propagated starting from configurations and momenta sampled every 1 ps from a 2 ns NVT trajectory following a 0.1 ns equilibration. For each trajectory, MSD(t) and the reorientational TCF, C2(t), were calculated along with their energy fluctuation versions, MSDH(t) and C2,H(t). Uncertainties in the results are reported as 95% confidence intervals according to the Student’s t-distribution based on block averaging with 8 blocks of 250 trajectories each.
A fully-periodic cubic simulation cell with a side length of 21.725 311 Å was filled with 343 water molecules to give a density of 0.997 g/mol. The H2O molecules were modeled using the SPC/E force field.23 This model is a completely rigid, three-site model that treats intermolecular interactions as a combination of Lennard-Jones (LJ) and Coulombic interactions. In this model, point charges are placed on each atom, while only oxygen atoms are treated as LJ sites. The Lorentz-Berthelot mixing rules are used to calculate intermolecular interactions between unlike atom types.24,25
The MD simulations were performed using the Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).26,27 A simulation time step of 1.0 fs was used, with configurations in the NVE trajectories saved every 50 fs for the calculation of the correlation functions. The SHAKE algorithm was used to hold the OH bonds and H–O–H angles rigid,28 with a tolerance of 0.0001 that specifies the relative error in the iterative solution. Intermolecular interactions were treated with a spherical cutoff of 10.5 Å and long-range electrostatics were described with an Ewald summation with an accuracy parameter of 0.0001 (that specifies the root-mean-squared error of the per-atom forces relative to a reference system). Canonical (NVT) simulations were performed using a Nosé-Hoover thermostat,29,30 with a thermostat damping parameter of 100 fs.
For comparison, activation energies were also calculated in the usual way using the Arrhenius equation from NVT trajectories at T = 285, 298.15, 315, and 330 K. Each trajectory was propagated for 4.5 ns with the first 0.5 ns used for equilibration. The trajectories were split into ten 0.4 ns blocks for block averaging to obtain 95% confidence intervals from the Student’s t-distribution.
The results from the Arrhenius calculations are presented in Fig. 2. From linear fits of the Arrhenius plots, activation energies of and kcal/mol are calculated for D and τ2, respectively, in good agreement with previously reported values.33,34 These results will be used as a comparison for calculations using the energy fluctuation method described in Sec. II.
Arrhenius plots for (a) D and (b) τ2 are presented. Results from the MD simulations (filled black circles) are shown along with linear fits (dashed lines) based on the Arrhenius equation.
Arrhenius plots for (a) D and (b) τ2 are presented. Results from the MD simulations (filled black circles) are shown along with linear fits (dashed lines) based on the Arrhenius equation.
IV. RESULTS
In this section, we apply the approaches described above to directly calculate the full temperature-dependence of TCFs and determine the activation energies associated with relevant time scales. We consider two examples involving bulk liquid water: self-diffusion and OH-bond reorientation.
A. Diffusion coefficient activation energy
Using the nonequilibrium MD simulations described in Sec. III, we have calculated the mean-squared displacement, MSD(t), and the corresponding energy fluctuation TCF, MSDH(t) for the oxygen atom of water at 298.15 K. The results are presented in Fig. 3. A linear fit to MSD(t) at longer times (between 2 and 20 ps) gives the diffusion coefficient as cm2/s, in excellent agreement with reported values in the literature for the SPC/E water model.23,31,32 The time-dependence of MSDH(t) is generally similar to MSD(t) itself in that, following a short initial period, it is linear with time (with a slope of 5.3 kcal/mol Å2/ps).
(Bottom) The TCFs MSD(t) (black line) and MSDH(t) (red line) are plotted versus time. (Top) The ratio MSDH(t)/MSD(t) (red line) is plotted as a function of time. A fit of this ratio between t = 15–20 ps to a constant value is also shown (blue dashed line). Note: MSD(t) is in units of Å2/ps and MSDH(t) in kcal/mol Å2/ps.
(Bottom) The TCFs MSD(t) (black line) and MSDH(t) (red line) are plotted versus time. (Top) The ratio MSDH(t)/MSD(t) (red line) is plotted as a function of time. A fit of this ratio between t = 15–20 ps to a constant value is also shown (blue dashed line). Note: MSD(t) is in units of Å2/ps and MSDH(t) in kcal/mol Å2/ps.
As shown in Eq. (10), the activation energy for the diffusion coefficient can be calculated from the ratio MSDH(t)/MSD(t) at long times. This ratio is also shown in Fig. 3 as a function of time and does indeed reach a constant value for times longer than ps. The value of the activation energy was obtained by fitting the ratio to a constant value for t = 15–20 ps, yielding kcal/mol. This result is in excellent agreement with the value of kcal/mol obtained from the Arrhenius plot in Fig. 2. Note that an alternative method for calculating Ea,D would be to use the ratio of the slopes obtained from the linear fits to the respective correlation functions; this yields 3.52 kcal/mol when fit over the time range 15-20 ps.
The calculation of diffusion activation energies in this manner provides an effective alternative to the usual Arrhenius method. Trajectories are required at only a single temperature and thus no choice needs to be made of the conditions for the simulations at other temperatures, e.g., whether to keep the same density (as we have done in our Arrhenius calculations) or modify the density to correspond to the experimental or simulation model result for each value of T. Additionally, as this approach calculates the temperature dependence of the diffusion coefficient, it may be used for systems and conditions where the behavior is non-Arrhenius.
B. Reorientation time activation energy
The same nonequilibrium MD trajectories used to evaluate the diffusion coefficient were analyzed to calculate the reorientational correlation function, C2(t), and its temperature dependence via the energy fluctuation TCF, C2,H(t). These two TCFs are plotted as a function of time in Fig. 4 along with the tri-exponential fit to C2(t), Eq. (12), and the related fit to C2,H(t), Eq. (14). The former gives the three time scales for the reorientational dynamics as 25 fs, 0.49 ps, and 2.6 ps, corresponding to inertial, librational, and H-bonding breaking and making dynamics, respectively. These parameters are also used in the fit to C2,H(t) so that the fitting parameters are the derivatives with respect to β of each of the amplitudes and time scales. Note that the fits to both C2(t) and C2,H(t) are in excellent agreement with the calculated TCFs.
The reorientational TCF C2(t) (solid black line) is shown as a function of time along with a tri-exponential fit (dashed black line) to Eq. (12). The TCF including the energy fluctuation, C2,H(t), is also shown (solid red line) along with the fit (dashed blue line) to Eq. (14).
A key focus of the analysis of the temperature dependence of the time scales is on the activation energy associated with the reorientational time τ2. The fit to C2,H(t) gives this as kcal/mol. This agrees with the result of kcal/mol obtained from the Arrhenius plot, Fig. 2(b), as well as prior calculations that also yielded 3.5 kcal/mol.33,34 It is notable that the activation energy for OH reorientation is similar to that obtained for self-diffusion. This is indicative of the common molecular origin of the two in the breaking and making of H-bonds. Namely, the exchange (or “jump”) between two different H-bond acceptors that is required for OH to reorient is also the key molecular event for diffusion of a water molecule.
It is important to note that because the derivatives with respect to β are obtained from a simulation at a single temperature, they may or may not correspond to an activation energy. That is, it is not possible to determine from the derivative alone whether or not depends linearly on 1/T. Often one has some prior knowledge, such as in the case of D or τ2, that the property is activated. However, this is not true for the inertial and librational time scales, τiner and τlib. An indication that they may not obey the Arrhenius behavior is that our results give uncertainties that encompass zero activation energy: kcal/mol and kcal/mol. Indeed, plots of 1/τiner and 1/τlib versus 1/T show that these time scales do not exhibit the Arrhenius behavior.
Additional information is available in the form of the amplitude derivatives. The fit to C2,H(t) gives , , and . These indicate that as temperature increases (β decreases) the amplitudes of the inertial and librational components increase and those of the H-bond making and breaking component decrease. This is consistent with linear fits to the amplitudes obtained from the simulations at different temperatures, which give −0.08, −0.12, and 0.20 for Ainer, Alib, and A2, respectively, for T = 285–315 K. This is not the full temperature range we have simulated; however, the estimated derivatives change significantly (to 0.04, −0.12, and 0.08) when T = 330 K is included in the fitting. This is a further indication that in the energy fluctuation TCFs we are obtaining local derivatives that can have distinct quantitative and qualitative differences from that obtained from multiple-temperature simulations that can depend on the temperature range considered.
C. Energetic decomposition
A key advantage of the method proposed in this work is the ability to decompose an activation energy or, more generally, a derivative with respect to temperature, into individual contributions due to each component of the energy. Specifically, we can note that the fluctuation in the energy that appears in MSDH(t) and C2,H(t) can be written as
Here, δKE(0) and δV(0) are the fluctuations in the kinetic and potential energies, respectively. The second equality notes that the potential energy fluctuation can be further decomposed into the various types of interactions including this simplest example of the Lennard-Jones, δVLJ, and Coulombic, δVCoul, contributions to the water energy. Then, the activation energy can be likewise divided into such contributions as
where , , and are the components of the activation energy associated with fluctuations in the kinetic energy, LJ potential energy, and Coulombic potential energy, respectively. This approach builds on the pioneering studies of Tolman35 and Truhlar36 to elucidate the physical interpretation of an activation energy.
In the case of diffusion, the component of the activation energy associated with the Coulombic interactions is, for example, then given by
and similarly for and . Note that this is not the only way to obtain insight into the origins of the activation energy and a particularly simple choice. In general, the energy fluctuation, δH(0), can be divided up in any number of ways to gain insight into the nature of the activated process.
The resulting TCF for the Coulombic contribution to MSDH(t) given in Eq. (19) is presented as a function of time in Fig. 5 along with the LJ and kinetic energy results and the total MSDH(t). Each contribution to the activation energy can be obtained by fitting the constant value reached at longer times (t = 15–20 ps); this gives 1.1, −0.8, and 3.2 kcal/mol for the kinetic energy, Lennard-Jones energy, and Coulombic energy, respectively. Thus, it is clear from the data that the dominant contribution to the activation energy of diffusion is the Coulombic interactions between water molecules. This is expected given the central role of H-bonding in the mechanism of water diffusion. What is perhaps less obvious is the negative contribution from the Lennard-Jones interactions that are slightly more than that canceled by the kinetic energy component. These results point to the new insight that may be obtained by the present approach.
The contributions to the diffusion TCF MSDH(t)/MSD(t) associated with the kinetic energy (red line), Lennard-Jones potential energy (violet line), and Coulombic potential energy (blue line) are plotted versus time along with the total MSDH(t)/MSD(t) (black line).
The contributions to the diffusion TCF MSDH(t)/MSD(t) associated with the kinetic energy (red line), Lennard-Jones potential energy (violet line), and Coulombic potential energy (blue line) are plotted versus time along with the total MSDH(t)/MSD(t) (black line).
The same decomposition approach can be applied to the reorientational correlation function, C2(t). For example, C2,H(t) given in Eq. (13) can be written as the sum of
and the analogous contributions C2,KE(t) and C2,LJ(t). These three components are shown as a function of time, along with the total C2,H(t), in Fig. 6. As with the diffusion constant, the dominant contribution to the activation energy is Coulombic interactions. In fact, there is essentially complete cancellation of the kinetic energy and LJ contributions such that the Coulombic component is nearly equal to the total C2,H(t) for all times. This is again an indication of the central role of H-bond exchanges in OH reorientation in water, which has been extensively explored in the extended jump model of Laage and Hynes.33,37
The contributions to the reorientational TCF C2,H(t) associated with the kinetic energy (red line), Lennard-Jones potential energy (violet line), and Coulombic potential energy (blue line) are plotted versus time along with the total C2,H(t) (black line).
The contributions to the reorientational TCF C2,H(t) associated with the kinetic energy (red line), Lennard-Jones potential energy (violet line), and Coulombic potential energy (blue line) are plotted versus time along with the total C2,H(t) (black line).
V. CONCLUSION
This work demonstrates a general approach for evaluating the temperature-dependence of time correlation functions that can also yield the activation energy for transport coefficients or dynamical time scales. A key feature is that the activation energy is obtained from simulations at a single temperature. The method has been demonstrated for the self-diffusion and OH reorientation in bulk liquid water and gives activation energies in agreement with those obtained from standard Arrhenius calculations. The framework, however, is not limited to these examples and can be straightforwardly extended to other transport coefficients, dynamical time scales, or TCFs. Indeed, this approach gives the temperature derivative of a dynamical time scale at a given temperature and it does not require an assumption of the Arrhenius behavior.
This approach also provides additional insight into the origins of the activation energy. In particular, we have shown how the activation energy can be decomposed into components associated with the various contributions to the system energy, e.g., kinetic, Lennard-Jones, and Coulombic energies. For both diffusion and reorientation in bulk water, nearly the entire contribution to the activation energy arises from the Coulombic interactions, which is associated with the central role of hydrogen-bond dynamics in both processes. This kind of analysis should lead to a better understanding of the molecular-level interactions that influence the activation energy.
Because the activation energy calculations do not require simulations at multiple temperatures, the present method may be particularly useful in cases where changing the temperature is problematic. For example, biological or self-assembled systems, such as lipid bilayers or reverse micelles, can often display dramatic changes in the structure, e.g., protein unfolding, with relatively small temperature changes. Similarly, for systems near a phase transition, the range of temperatures for which an Arrhenius analysis can be used is strongly constrained. However, the approach presented here permits the calculation of activation energies or, more generally, the derivative of full time correlation functions with respect to temperature even in such cases.
ACKNOWLEDGMENTS
The authors of this work thank the NSF-EPA Networks for Sustainable Materials Design and Synthesis (NSMDS) under NSF Grant No. CHE-1339661 for funding. O.O.M. gratefully acknowledges support from a George E. Walrafen fellowship. The calculations were performed at the University of Kansas Center for Research Computing (CRC).
REFERENCES
See, e.g., Ref. 17.
It can be shown that the inertial dynamics give an initial Gaussian decay (see, e.g., D. A. McQuarrie, Statistical Mechanics, Harper-Collins, New York, 1976). However, it is more convenient to fit this component as an exponential decay like that associated with the librational and H-bond exchange dynamics. While this yields an approximate time scale for the inertial dynamics (which are not a focus of the present analysis), it does not affect the values obtained for τlib or τ2, as we have verified by different fitting approaches.
We have not accounted for finite-size effects32 in calculating the diffusion constants or the corresponding activation energies in any of the approaches used. Thus, while the comparisons are between consistent descriptions, the numerical results may be modified by such corrections.