Second-Generation ReaxFF Water Force Field: Improvements in the

Article pubs.acs.org/JPCB

Second-Generation ReaxFF Water Force Field: Improvements in the Description of Water Density and OH-Anion Diffusion Weiwei Zhang and Adri C. T. van Duin* Department of Mechanical and Nuclear Engineering, Pennsylvania State University, University Park, Pennsylvania 16802, United States S Supporting Information *

ABSTRACT: Hydronium (H3O+) and hydroxide (OH−) ions have anomalously large diffusion constants in aqueous solutions due to their combination of vehicular and Grotthuss hopping diffusion mechanisms. An improvement of the ReaxFF reactive water force field on the basis of our first-generation water force field (water-2010) is presented to describe the proton transfer (PT) mechanisms of H3O+ and OH− in water. Molecular dynamics simulation studies with the water-2017 force field support the Eigen−Zundel−Eigen mechanism for PT in acidic aqueous solution and reproduce the hypercoordinated solvation structure of the OH− in a basic environment. In particular, it predicts the correct order of the diffusion constants of H2O, H3O+, and OH− and their values are in agreement with the experimental data. Another interesting observation is that the diffusion constants of H3O+ and OH− are close to each other at high concentration due to the strong correlation between OH− ions in basic aqueous solution. On the basis of our results, it is shown that ReaxFF provides a novel approach to study the complex acid−base chemical reactions in aqueous solution with any pH value. H3O+ and OH− are not simple mirror images of each other and they provide different proton transport mechanisms.19−22 Various authors have shown that OH− prefers to be hypercoordinated by the surrounding water molecules,6,21,23 and one of the mirror species (H3O2−) could only be found as a short-lived transient solvation complex.24 At the viewpoint of theoretical studies, the diffusion of H3O+ and OH− in aqueous solution involves both vehicular and Grotthuss hopping mechanisms, where the latter refers to oxygen−hydrogen bond cleavage and formation. Therefore, traditional classical molecular dynamics (MD) simulations are insufficient.25 In principle, the ab initio MD (AIMD) simulations should lead to a more accurate result. However, the exact form of the exchange-correlation (XC) functional is unknown, and the solvation structures and diffusion mechanisms of H3O+ and OH− strongly depend on the choice of XC functionals.21 Additionally, it is still impractical to handle the systems at the nanoscale level currently (both time and system size) due to the high computational cost. To deal with big size systems in a large time scale, a methodology is required to bridge the gap between quantum chemistry and traditional force field methods to study the PT in aqueous solutions. Several groups have carried out the research on this field.26 A relatively straightforward solution is the so-called quantum mechanics/molecular mechanics (QM/MM) technique, which performs classical simulations in which particular regions are

1. INTRODUCTION Water is essential to the existence of life on earth. It has many unique features such as high freezing point, smaller density in the solid than liquid phase, high surface tension, etc. All of the properties are associated with its extremely high density of intermolecular hydrogen bonds (HBs).1 Among the water properties, the diffusion of hydronium (H3O+) and hydroxide (OH−) ions in aqueous solution attracts much attention due to their anomalously large mobility.2−4 Such high-speed diffusion is generally ascribed to the contribution of the so-called Grotthuss hopping mechanism featured by proton transfer (PT) between the cation or anion and its neighboring water molecule rather than a simple vehicular diffusion.5,6 Since H3O+ and OH− transport plays an important role in chemistry and biological processes, ranging from simple acid−base chemical reactions to complex enzymatic reactions,7−10 an understanding of the elementary steps of PT and proton transport mechanisms in the HB network is thus of fundamental significance. The current accepted picture of H3O+ diffusion in aqueous solution is that the H3O+ involves interconversions between two distinct hydration structures, namely, the Eigen (H3O+ (H2O)3) and Zundel ([H2O···H···OH2]+) cations, based on both experimental and theoretical investigations.7,10−17 Unfortunately, the characteristics of OH− in aqueous solution are less clear because it is historically assumed to be simply a “mirror image” of the H3O+. The OH− diffusion with such a mirror concept was described in detail by Agmon.18 However, during the past decades, the experimental and theoretical research has demonstrated that the hydration structures of © 2017 American Chemical Society

Received: March 17, 2017 Revised: May 31, 2017 Published: June 1, 2017 6021

DOI: 10.1021/acs.jpcb.7b02548 J. Phys. Chem. B 2017, 121, 6021−6032

Article

The Journal of Physical Chemistry B

Figure 1. (a) The bond dissociation energies for the O−H bond and (b) strain energies for the H−O−H angle in the water molecule. (c−e) Binding energies for water dimer as a function of oxygen−oxygen distance with the water-2010 force field, the water-2017 force field, and QM calculations, respectively.

treated quantum mechanically.27,28 It allows large simulations based on the efficiency of empirical potentials. However, this strategy requires specific regions to be treated as either quantum or classical and it affects the flexibility of the model. Lill et al. introduced the stochastic proton hopping scheme (Qhop) that is capable of simulating single excess proton transport in aqueous solutions.29 Such hopping algorithms are often useful, but it is intractable to associate actual physical interactions with the underlying dynamics. Selvan et al. proposed a reactive molecular dynamics algorithm, which incorporates the structural diffusion of a proton into a classical MD simulation to study the effect of temperature and concentration on the proton transport.30 Recently, a scaled model potential (sOSS2) for liquid water provided a good description of the solvation structure, mechanism of PT, and accurate diffusion constants of H3O+ and OH−.31 It also predicts the temperature dependence of the diffusion constants and the related activation energies that are in excellent agreement with experimental results after the necessary corrections were made from finite concentrations to infinite dilution. Another popular option is the multistate empirical valence bond (MS-EVB) method developed by Voth and coworkers,14,32,33 which is done at the semiempirical level, via introducing the diabatic states where the proton is assigned to different water molecules. This approach was successfully applied for investigating PT in aqueous solution and polymer electrolyte membranes.10 However, although the SCI-MS-EVB model was recently presented to treat the multiple excess protons in aqueous solution,34 the computational cost still linearly increases with the number of H3O+ in the system. Besides, this method currently cannot handle complicated systems, where both PT and complex chemical reactions are involved. The ReaxFF reactive force-field method provides a suitable alternative to study the mechanism of proton transport, as well

as associated chemical reactions.35 It has been successfully applied for the investigation of complex chemical transformations involving high-energy materials, silion/silicon oxides, and polymer decompositions in the gas, solid phase and at the interface.36−39 Moreover, the first-generation ReaxFF water force field (water-2010) is able to describe the dissociation of water and structural migration of H3O+, and it also succeeds in reproducing the experimental diffusion constants for both water and H3O+ in aqueous solution.40 However, we also note that there are some issues for the water2010 force field. First, the density of pure bulk water predicted by NPT simulation is significantly smaller than experiment data. Furthermore, the calculated diffusion constants of water, hydronium, and hydroxide ions are in the order of H2O < H3O+ < OH−, which is inconsistent with experiment (the diffusion constant of OH− is smaller than that of H3O+ but much bigger than water2,4). The second-generation ReaxFF water description presented here (water-2017 force field) manages to address these issues, and we present our simulation results in this paper, including the description of the structural and dynamical properties of the H3O+ and OH− in aqueous solution, as well as water itself. This paper is organized as follows: In section 2, we describe the force field development, validation, and simulation details. The results and discussion are shown in section 3, while conclusions are summarized in section 4.

2. COMPUTATIONAL METHODS 2.1. ReaxFF Force Field Development. ReaxFF, as a bond order dependent force field, can simulate bond formation and bond dissociation during the MD calculation. The general expression is given by 6022


Article


Figure 2. (a) Chemical structures of [H2O−H−OH2]+ and [HO−H−OH]− clusters. (b and c) PT barriers predicted with the water-2010 and water-2017 force fields as a function of O−O distance in [H2O−H−OH2]+, respectively. (d−f) PT barriers with the water-2010 and water-2017 force fields and QM calculations in [HO−H−OH]−, respectively.

water-2010 force field may not have described the hydrogen bond and nonbonded interactions properly. Therefore, we mainly retrained the parameters related to such interactions, for example, the parameters of oxygen−oxygen and oxygen− hydrogen van der Waals and hydrogen bond interactions. To improve the prediction of density, the experimental data is taken as the reference. Meanwhile, the QM results were used to improve the hydrated structure description of OH− ion (the details will be discussed as follows). Of course, the changes of these parameters also affect the bond order of the two atom pairs, making it is also necessary to retrain the parameters of bond and angle. In principle, the aim of our reoptimization of the force field is allowing it to produce a reasonable OH− diffusion in aqueous solution. Meanwhile, the good description for the properties of water and H3O+ should remain after fitting. From our previous studies,40 it has been shown that the water-2010 force field predicted good results compared to the QM data at the X3LYP/6-311++G** level, including the bond dissociation (Figure 1a), angle strain curve of water (Figure 1b), as well as binding energies of water dimer with different symmetry (the water dimer symmetries are Cs, Ci, C2, and C2v, as shown in Figure 1c−e). Here, we provide a comparison of the water-2010 and water-2017 force fields with QM data, also shown in Figure 1; good agreement between the two force fields with QM indicates that both of them are able to describe the water single molecule and dimer well. We also find that the binding energies of the water dimer are more strongly overestimated for the water-2017 force field compared to the QM data (about 1.00 kcal/mol deviation). This may be due to the water density retraining, which leads to a slightly stronger intermolecular interaction compared with the result from the water-2010 force field.

Esystem = E bond + Eval + Etor + Eover + Eunder + E lp + Evdw + Ecoul + E HB

(1)

which includes terms related to bond, angle, torsion, overcoordination, undercoordination, and lone pair energies. The nonbonded contribution refers to the van der Waals and Coulomb interactions, as well as the hydrogen bond energy. The details of the ReaxFF potential are shown in van Duin and Chenoweth publications.35,41 In this work, the force field parameters were optimized against the quantum mechanical (QM) calculations and/or experimental data using the successive one-parameter search technique42 ⎛ (xi ,ref − xi ,ReaxFF) ⎞2 ⎟ Error = ∑ ⎜ σi ⎠ i ⎝ n

(2)

Here, xi,ref and xi,ReaxFF are the QM and/or experimental data and the ReaxFF calculated values, respectively. σi is the weight assigned to ith data point. In our training procedure, the energies and structures were fitted to the QM data while the density of water and the diffusion constants of H2O, H3O+, and OH− were trained against the experimental values. The ReaxFF force field parameters, as improved in this work, are available in the Supporting Information. The simulation with the new force field parameter can be done using the standalone ReaxFF code, ADF43 and LAMMPS44 software packages. For example, all of the present MD simulations were performed by the parallel ADF-2016 program. From our previous studies, we have already known the water-2010 force field reproduced good results of structural and dynamics properties of water molecule and H3O+ ion in aqueous solution.40,45 However, it failed to predict the water density and properties of OH− ion (i.e., the hydrated structure and the diffusion constant). We believe the 6023


Article


Table 1. Comparison of Relative Energies (Units of kcal/mol) for OH−(H2O)m (m = 3 or 4) Clusters with QM and the Water2010 and Water-2017 Force Field Calculations

a

Reference energies for three- and four-water clusters.

Figure 3. Simulated and experimental radial distribution functions g(r) and their CNs for (a) H−H, (b) H−O, and (c) O−O in bulk water (the RDF is shown as a solid line and the CN as a dashed line). (d) O−O RDF predicted with different water models.

field predicts that the energy increases with the O−O distance varying 2.4, 2.5, 2.3, 2.6, and 2.7 Å while the order changes to 2.5 < 2.4 < 2.6 < 2.3 < 2.7 Å O−O distance in water-2017; the latter order is consistent with QM data at the MP2/aug-ccpVDZ level, as shown in Figure 2f.22 Another remarkable difference is that the energy of the minimum with 2.4 Å of O− O distance is much lower than that with 2.6 Å predicted by the water-2010 force field, while water-2017 becomes better. Furthermore, we compared the relative energies of some typical OH−(H2O)m clusters, as summarized in Table 1. For three water clusters, the water-2017 force field predicts that the three-coordinated structure is more stable than the twocoordinated one, which is consistent with the QM results at the

To test the PT reactions, we considered both cationic and anionic systems and the chemical structures of the two cases are shown in Figure 2a. It is clear that the energy curves obtained by the water-2017 force field are very close to that predicted by water-2010 for the system [H2O−H−OH2]+ (Figure 2b and c). Similar to ref 40, shortening the O−O distance removes the PT barrier but increases the overall energy, while lengthening the distance increases the barrier. However, the water-2010 and water-2017 force fields predict different trends for the [HO− H−OH]− system, as shown in Figure 2d and e. One significant difference are the energy orders in various O−O distances at the high symmetry-point (i.e., the distance between the sharedhydrogen and the O−O center is zero). The water-2010 force 6024


Article


Figure 4. Radial distribution functions (solid lines) and their CNs (dashed line) between atoms of H3O+ and water molecules. For each panel, results from the water-2010 (blue) and water-2017 (red) force fields and AIMD (cyan) simulations are shown. The atoms labels are defined in the text.

MP2/6-311++G(3df,3pd) level.46 However, the water-2010 force field gives the incorrect order. For the four-water case, the water-2010 force field prediction shows that the cluster prefers to form a three-coordinated structure and the relative energy of the four-coordinated structure is overestimated compared to the QM data.46 Instead, the energy of the four-coordinated cluster is reduced significantly predicted by the water-2017 force field. We should mention that all of the energies were calculated in the gas phase at 0 K. Considering the thermal fluctuation and different phase’s effect, the coordination number of OH− predicted with the water-2017 force field might be bigger than the water-2010 force field at room temperature in the liquid water phase due to the relatively small energy difference between three- and four-coordinated structures. The details will be discussed below. 2.2. ReaxFF MD Simulation Details. In order to simulate the density of liquid water, a cubic box containing 216 water molecules was constructed and a 200 ps NPT/MD simulation was preformed to equilibrate the periodic system at room temperature. Then, the final configuration from this NPT simulation was used as a starting point for a 200 ps NVT ensemble MD simulation to get the dynamical properties of bulk water. After that, randomly selected water molecules were replaced by H3O+ or OH− to consider the acidic or basic conditions. In the present work, four concentrations of ions were considered, namely, H3O+(OH−)/215H2O, 5H3O+(5OH−)/211H2O, 10H3O+(10OH−)/206H2O, and 20H3O+(20OH−)/196H2O. After introducing the H3O+ or OH−, 200 ps NVT/MD simulations were performed and the last 120 ps results were used for the statistical analysis of the structural property and calculations of the diffusion constant of H2O, H3O+, and OH−. During our NVT MD simulations, the

time step was set to 0.25 fs and a Berendsen thermostat was used with a temperature damping constant of 100 fs.

3. RESULTS AND DISCUSSION 3.1. Force Field Validation for Water. Density of Water. The liquid phase of water was simulated at 298.15 K, and the equilibrium densities are 0.92 and 1.01 g/cm3 with the water2010 and water-2017 force fields, respectively. Clearly, the water-2017 reproduces the experimental density significantly better, well within 5% of experimental data. Structural Properties of Water. In order to investigate the structural properties of water, we calculated the radial distribution function (RDF) and the coordination number of atom within a given radius shell47 g (r ) = CN =

n(r ) ρ 4πr 2Δr

∫0

r′

ρ 4πr 2g (r ) dr

(3)

where g(r) and CN are the RDF and coordination number, respectively. n(r) represents the atom number within a distance r′ of the central atom, and ρ stands for the number density. The RDF and CN were obtained by averaging over the trajectory. Figure 3a−c compares the ReaxFF and experimental H−H, H−O, and O−O RDFs for water at room temperature (here the peak within intramolecular distance was not shown due to its too large intensity), as well as their CNs. As mentioned in our earlier work,40 the water-2010 force field gives a good reproduction of the experimental data.48 In comparison, there are a small left shift and increasing first peak intensity for the 6025


Article


Figure 5. SDF of the water molecules around H3O+ (OH−) with the water-2017 force field: Top panels a and b depict the SDF of the water molecules around H3O+ within 2.8 and 3.3 Å, respectively; bottom panels c and d show the SDF with respect to OH− within 2.9 and 3.2 Å, respectively.

2017 force fields, as well as the AIMD results predicted by the Becke−Lee−Yang−Parr (BLYP) density functional with the Gaussian TZV2P basis set and empirical dispersion corrections, which leads to an improved description of water properties compared to the functional without long-range corrections.51 It is found that the water-2017 force field predicts similar results to water-2010, indicating both force fields can describe the hydration structure of H3O+ well. In addition, compared to the QM result, all RDFs predicted with ReaxFF simulations show a general agreement in terms of the peak positions and heights. One deviation is the first peak of O*−O (Figure 4d, O* and O stand for oxygen of H3O+ or OH− and oxygen in water molecule, respectively; similarly, H* and H represent hydrogen of H3O+ or OH− and hydrogen in water molecule in the following). The ReaxFF prediction shifts outward and becomes sharper and higher, which indicates that H3O+ and its first shell solvation water molecules are more loosely aggregated compared to the structure predicted with AIMD. However, both ReaxFF force fields are able to predict the same coordination number (CN) as AIMD for the first solvation shell and the loose structure of the hydronium−water cluster aids in the probability of proton transfer, thereby giving a more reasonable diffusion constant (the details will be discussed in the next section). On the contrary, AIMD calculations showed a much smaller diffusion constant of H3O+ because of the lower proton transfer probability, as well as water self-diffusion.51 Another discrepancy is that there is a small peak for the RDF of O*−H at 2.2 Å, as shown in Figure 4c; the hydrogen comes from the water molecule, which donates a HB to the lone pair site of H3O+. It seems to be a general trend for charge defects in describing relatively weak interactions compared to the strong

RDF of O−O from the results of the water-2017 MD simulations. Accordingly, the height of the first peak of H−O also becomes slightly higher than that of experiment. This might be caused by our training for the density of bulk water; as a result, the interaction between water molecules becomes slight stronger. However, the width of the O−O peak predicted with the water-2017 force field is much closer to the experiment. Furthermore, the water-2017 results can also be comparable with other popular classical models,49 such as SPC/ E, TIP3P, and TIP4P-2005, as shown in Figure 3d. In addition, all of the calculated CNs with the water-2017 force field are also consistent with water-2010 force field prediction. Most importantly, water-2017 produces a better density of bulk water and a unified reasonable description of structural and dynamical properties of H2O, H3O+, and OH−, such as diffusion constants as discussed in the next sections. 3.2. Proton Transfer in Water. In order to study the structural and transport properties referring to H3O+ or OH−, the index of oxygen (O*) of the H3O+ or OH− should be monitored during the simulation and distinguished from the oxygen atoms of water because the index of O* may change due to the PT between water molecule and H3O+ or OH−. We evaluated the sequence of PT events as was done in previous works.31,45,50 The same index number of O* between two adjacent frames in the MD trajectory indicates the vehicular diffusion without any PT, while the Grotthuss hopping occurs if the index changes between the adjacent frames.45,50 Structural Properties of H3O+ and OH−. The hydronium− water and hydroxide−water interactions can be characterized by their RDFs. Figure 4 shows the hydronium-water RDFs of H*−H (Figure 4a), H*−O (Figure 4b), O*−H (Figure 4c), and O*−O (Figure 4d), predicted with water-2010 and water6026


Article


Figure 6. Radial distribution functions (solid lines) and their CNs (dashed line) between atoms of OH− and water molecules. For each panel, results from the water-2010 (blue) and water-2017 (red) force fields and AIMD (cyan) simulations are shown. The atoms labels are defined in the text.

Similar findings were also reported in other works.24 Another remarkable difference is the first peak of the O*−O RDF (Figure 6d), as well as its coordination number. The water2017 force field predicts a stronger correlation and larger coordination number (around 5 and 4 predicted with the water2017 and water-2010 force fields, respectively) compared to water-2010, indicating more water molecules are involved in the first solvation shell of OH−. Apart from the top water molecule, which accepts a HB from OH− as the SDF discussion in the following, the water-2017 force field can capture the hypercoordinated solvation structure of OH−, in which it forms four accepting HBs. We expect this difference should have an effect on the proton hopping rate and diffusion constant, because the higher coordination might result in better stability of the hydroxide−water cluster; this is true, as proved from the diffusion constant results as follows. Compared to the AIMD results with the BLYP functional,51 the height of the first peak predicted with the water-2017 force field is much higher, as well as the RDF of O*−H. Actually, we note that the ReaxFF intensities are much closer to the experimental data.53 It indicates AIMD underestimated the correlation between OH− and water molecules. In addition, it is found that the intensity of the first peak for H*−O (Figure 6b) predicted with ReaxFF calculations is stronger than the AIMD prediction. Again, such strong correlation arises from the overestimated description between the OH− and water molecule accepting HB (H2O··· HO−), as discussed in the H3O+ case similarly. To study details of the solvation structure, SDFs of the solvating water oxygen atoms around OH− were shown in Figure 5c and d within rO*−O < 2.9 Å and rO*−O < 3.2 Å, respectively. We find that the orientations of the water molecules with respect to the OH− are

HB interaction; similar results were also reported in other reactive force field models.22,51 To further elucidate the solvation structure, the spatial distribution functions (SDFs) predicted with the water-2017 force field were visualized in Figure 5 (the SDFs with the water2010 force field are similar, and they are not shown here). The SDFs of the solvating water oxygen atoms around H3O+ within rO*−O < 2.8 Å and rO*−O < 3.3 Å are illustrated in Figure 5a and b, respectively. Clearly, three regions of maximum probability are visible (Figure 5a) in correspondence with the three water molecules accepting hydrogen bonds. It demonstrates the three-coordination nature of H3O+ and formation of strong hydrogen bonds between H3O+ and first solvation shell water molecules. With increasing distance, a broad cloud is observed on top of the H3O+ (Figure 5b); it corresponds to the additional water molecule donating a hydrogen bond via relatively weak interaction binding to H3O+, as indicated by the shoulder in the RDF of O*−O around 3.2 Å (Figure 4d). Our findings are in good agreement with experimental observation.52 For the solvation shell of OH−, two differences are observed from the comparison of RDFs of H*−H (Figure 6a), H*−O (Figure 6b), O*−H (Figure 6c), and O*−O (Figure 6d), calculated with the water-2010 and water-2017 force fields. One is the RDF of O*−H (Figure 6c); there is a small shoulder predicted with water-2010 at about 1.3 Å which is not observed with the water-2017 force field. From analyzing the structures of the hydroxide−water cluster, we found the shoulder corresponds to the Zundel-like intermediate [HO···H···OH]−. It indicates this low-coordinated intermediate becomes a shortlived transient species in the water-2017 force field simulations. 6027


Article


Figure 7. Configurations showing the mechanism of proton transfer in the simulation of (a) one H3O+ in 215 water molecules and (b) one OH− in 215 water molecules. (i), (ii), and (iii) show the stabilized structure before proton transfer, intermediate, and completion of proton transfer, respectively.

Zundel−Eigen mechanism and shuttling feature were reported in the AIMD simulations and other discussions previously.10,13,31 The mechanism in basic aqueous solution was also detected and shown in Figure 7b. The dominant fourcoordinated (Figure 7b(i), exclude the top water accepting a hydrogen bond from OH−) structure, which is supported with various experimental measurements,23,54 becomes “active” when its coordination number decreases to three during the reformation of the hydrogen bond network. In this configuration, OH− is tetrahedrally coordinated by forming three accepting hydrogen bonds (Figure 7b(ii)), and the proton can transfer from water molecule to OH− more easily. Unlike the hydration structure of H3O+, the lifetime of such a 3fold-coordinated structure is much longer than Zundel-like anion [HO···H···OH]− and such a three-folded structure governs the proton transfer probability and diffusion constants.21 From our ReaxFF simulations, we demonstrate the OH− is not a simply “mirror image” of H3O+. The hypercoordination structure and mechanism of proton transfer for OH− ion are similar to those predicted by scaled OSS2 simulations.31 As mentioned above, the proton shuttles between two oxygen atoms during the proton transfer process. Two types of hopping events are involved; one is oscillatory shuttling, and the other is Grotthuss hopping. Accordingly, the proton forward hopping function can be predicted by33

almost the same for both distances. One expects an overestimated hydrogen bond is formed between OH− and its top water because of its large distribution at the small distance of rO*−O. Different from our simulation results, experiments showed that four neighboring water molecules are around the OH− to form the OH−(H2O)4 cluster at a short distance and the fifth water is weakly bonded to the hydroxide− four water complex at a long distance.53 However, the fifth water molecule accepting a HB is of secondary importance compared to an accurate description of the hydration structure of OH−, where proton transfer is involved.22 Also, it is found that the SDF within a distance of 3.2 Å is in agreement with experimental observations.53 Mechanism of Proton Transfer in Acidic and Basic Solutions. In the above discussions, the static properties of H3O+ and OH− have been presented, but the mechanism of proton transfer in solution is still unclear. To achieve dynamical information, analysis of the structures as a function of time is needed to illustrate the proton transfer mechanism in water. In acidic aqueous solution, H3O+ is 3-fold-coordinated and the center water is 4-fold-coordinated before proton transfer occurs (Figure 7a(i)). With the thermal fluctuation, one hydrogen bond between the central water and its solvating water molecule donating a hydrogen bond breaks, leading to the solvating water molecule far away from the central water. Such solvation adjustment pushes the proton in H3O+ toward the undercoordinated central water molecule; as a result, the Zundel structure is formed, as shown in Figure 7a(ii), and the shared proton shuttles between two oxygen atoms and the product depends on the result of solvent reorganization. If the solvent returns to its former configuration (Figure 7a(i)), proton transfer does not take place. Alternatively, if the solvent reorganizes into the structure shown in Figure 7a(iii), proton transfer would occur and a new H3O+ is formed. Such Eigen−

h(Δt ) = h(Δt − 1) + Δh(Δt )

(4)

where Δt is the time step and h(0) = 0. Δh(Δt) are 0, 1, and −1 for no proton transfer, hopping to a new water molecule, and hopping back to the last one, respectively. For the purpose of comparison, Figure 8 depicts the forward hopping function with simulation time for H3O+/215H2O and OH−/215H2O systems calculated from averaging over 50 trajectories with the 6028


Article


contribution to the diffusion constant of ions which has been discussed in our previous work.50 3.3. Diffusion Constant of H3O+ and OH− in Water. From the time-dependent trajectory of O*, the diffusion constant of H3O+ and OH− can be estimated from mean square displacement (MSD) of O*55 MSD(m) = ⟨|r(t ) − r |2 ⟩ =

1 n

n

∑ |r(m + i) − r(i)|2 (5)

i=1

where r represents the position of the particle in the unfolded trajectory, t is the time, m is denoted the maximum point number, n is the number for averaging, m + n stands for the total number of frames, and i is the step counter. With the timedependent MSD, the diffusion constants were then estimated using the Einstein relation

Figure 8. Forward hopping functions as a function of time predicted with the water-2010 and water-2017 force fields.

D= water-2010 and water-2017 force field simulations. It is overall linear, and this characteristic allows the calculation of the forward hopping rate as the slope of the curve. In the H3O+/ 215H2O system, the fitted proton forward hopping rates are 0.50 and 0.53 ps−1 with the water-2010 and water-2017 force field simulation results, respectively; both of them are very close to each other, and they should predict similar proton transfer rates. Similarly, the estimated hopping rates in OH−/215H2O are 0.51 and 0.35 ps−1 with the water-2010 and water-2017 force fields, respectively. It illustrates the forward hopping rate of OH− is close to that of H3O+ for the water-2010 force field simulations, while it descends significantly predicted with the water-2017 force field. A lower hopping rate induces a lower proton transfer rate; as a result, it might yield a small diffusion constant because the Grotthuss hopping has a significant

1 ⟨|r(t ) − r |2 ⟩ 6Nt

(6) +

−

Here N is the atom number of H2O, H3O , or OH . To evaluate the diffusion constant of H3O+ and OH− at all concerned concentrations, we averaged over 50 trajectories and plotted the time dependence of MSD in Figure 9. All of the curves have linear characteristics basically, whose slopes are related to the diffusion constants directly. It is obvious that the slopes decrease with increasing concentrations for both H3O+ and OH−. In the same concentration, the slope of OH− predicted with the water-2010 force field is bigger than that of H3O+ (as shown in Figure 9a and b, respectively), indicating the diffusion constant of OH− is higher than that of H3O+, which is inconsistent with experiment. In contrast, the water2017 force field predicts the correct order of the diffusion constant of the two ions (the slope of OH− is smaller than that

Figure 9. MSD for H3O+ (left) and OH− (right) with the water-2010 and water-2017 force fields at the various concentrations. 6029


Article

The Journal of Physical Chemistry B of H3O+, as shown in Figure 9c and d, respectively). For clarity, the diffusion constants were listed in Table 2; the experimental

result has also been reported in the classical MD simulations by Han et al.56 We expect the strong correlation among OH− ions may help to enlarge the PT distance in solution, and as a result, the diffusion constant of OH− decreases slowly as the concentration increases compared to that of H3O+. It should be mentioned that we do not consider the counterion effect on the diffusion of H3O+ and OH− in this work. The presence of the counterion should reduce the diffusion constants of them, because the counterion disrupts the way of diffusion for H3O+ and OH− ions, especially at the high ion concentration.30

Table 2. Diffusion Constants (Å2/ps) of H2O, H3O+, and OH− at 298.15 K H2O H3O+

OH−

216H2O H3O/215H2O 5H3O/211H2O 10H3O/206H2O 20H3O/196H2O OH/215H2O 5OH/211H2O 10OH/206H2O 20OH/196H2O

expt.a

water-2010

water-2017

0.24 0.93

0.25 0.94 0.73 0.59 0.38 1.03 0.99 0.83 0.70

0.27 1.04 0.86 0.70 0.44 0.78 0.66 0.59 0.44

0.56

4. CONCLUSIONS In this work, we developed a second generation ReaxFF water force field to describe the proton transfer (PT) processes of hydronium (H3O+) and hydroxide (OH−) ions in aqueous solution. The force field parametrization is based on our previous water force field (water-2010). From analyses and comparisons, it demonstrates that the present force field (water-2017) is improved in the following aspects: (1) The density of water is consistent with the experimental result due to our training against the water density. (2) The solvation structure of OH− is further improved, and the coordination of OH− is in good agreement with experiment. (3) Most importantly, the water-2017 force field predicts a smaller diffusion constant for OH− compared with the water-2010 force field and it reproduces the correct order of the diffusion constants of H2O, H3O+, and OH−, which are in the order of H2O < OH− < H3O+. By means of ReaxFF reactive molecular dynamics simulations, we showed that ReaxFF supports the Eigen−Zundel− Eigen PT mechanism for H3O+ in acidic aqueous solution and it can also reproduce a hypercoordination hydration structure for OH− in basic aqueous solution at room temperature. From our analysis of proton hopping rates, we observe that a fast proton transfer rate induces a larger number of diffusion constant of ions; therefore, we expect the crucial factor to reproduce the diffusion is correctly describing the hydration structures of the ions. In addition, the studies of diffusion constants as a function of ion concentration show that the diffusion constants of both H3O+ and OH− decrease with the increase of their concentrations, and the diffusion constants of H3O+ and OH− are close to each other at a high concentration due to the strong correlation between OH− ions. On the basis of our studies, we expect ReaxFF could be employed to study acid−base chemistry in the fields of material science, chemistry, and engineering at a large scale simulation, especially for the proton transfer at a wide range of pH values in aqueous solution.

The experimental diffusion constants for H3O+ and OH− ions at infinite dilution.4 a

data of the H3O+ and OH− are also included in this table, as well as the self-diffusion constant of water. Here, we should mention that the experimental diffusion constants of H3O+ and OH− ions were measured at infinite dilutions while our calculations were done at finite concentrations. For instance, the lowest ion concentration is 0.26 M in our simulations. Extrapolated to finite dilution with a concentration of 0.26 M, the experimental diffusion constants are 0.82 Å2/ps for H3O+ and 0.45 Å2/ps for the OH− ion,31 which are approximately 10% smaller than that in infinite dilutions. It is clear that the diffusion constants of water predicted with both force fields are in good agreement with the experimental value. However, the water-2010 and water-2017 force fields predict a different order of diffusion constant for H3O+ and OH−, and the water-2017 force field prediction agrees with the experimental order. This indicates that the water-2017 force field is better than water2010 and it is indeed able to predict the diffusion constants of H3O+ and OH− reasonably. Furthermore, it is interesting to find that the diffusion constants of H3O+ and OH− are close to each other at their highest concentrations. To understand this phenomenon, we compared their RDFs of O*−O* for the 20ion systems, as shown in Figure 10. It is surprising to find that there are two high peaks for 20OH−/196H2O but no peak for H3O+/196H2O. The findings indicate the correlation between OH− ions is much stronger than that of H3O+ ions. A similar

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b02548. The second generation ReaxFF reactive force field parameters for water that were used in the present work (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

+

Figure 10. Radial distribution functions for O*−O* in 20H3O / 196H2O (blue) and 20OH−/196H2O (red) systems predicted with the water-2017 force field.

ORCID

Weiwei Zhang: 0000-0002-5255-7340 6030


Article

The Journal of Physical Chemistry B Notes

(21) Tuckerman, M. E.; Chandra, A.; Marx, D. Structure and Dynamics of OH-(Aq). Acc. Chem. Res. 2006, 39, 151−158. (22) Ufimtsev, I. S.; Kalinichev, A. G.; Martinez, T. J.; Kirkpatrick, R. J. A Multistate Empirical Valence Bond Model for Solvation and Transport Simulations of OH− in Aqueous Solutions. Phys. Chem. Chem. Phys. 2009, 11, 9420−9430. (23) Megyes, T.; Bálint, S.; Grósz, T.; Radnai, T.; Bakó, I.; Sipos, P. The Structure of Aqueous Sodium Hydroxide Solutions: A Combined Solution X-Ray Diffraction and Simulation Study. J. Chem. Phys. 2008, 128, 044501. (24) Chen, B.; Park, J. M.; Ivanov, I.; Tabacchi, G.; Klein, M. L.; Parrinello, M. First-Principles Study of Aqueous Hydroxide Solutions. J. Am. Chem. Soc. 2002, 124, 8534−8535. (25) Devanathan, R.; Venkatnathan, A.; Rousseau, R.; Dupuis, M.; Frigato, T.; Gu, W.; Helms, V. Atomistic Simulation of Water Percolation and Proton Hopping in Nafion Fuel Cell Membrane. J. Phys. Chem. B 2010, 114, 13681−13690. (26) Stinson, J. Modelling Microscopic Clusters of Sulphuric Acid and Water Relevant to Atmospheric Nucleation. UCL (University College London), 2015. (27) Payaka, A.; Yotmanee, P.; Tongraar, A. Characteristics of the “Hypercoordination” of Hydroxide (OH−) in Water: A Comparative Study of HF/MM and B3LYP/MM MD Simulations. J. Mol. Liq. 2013, 188, 89−95. (28) Somtua, T.; Tongraar, A. Correlation Effects on the Structure and Dynamics of the H3O+ Hydrate: B3LYP/MM and MP2/MM MD Simulations. Phys. Chem. Chem. Phys. 2011, 13, 16190−16196. (29) Lill, M. A.; Helms, V. Molecular Dynamics Simulation of Proton Transport with Quantum Mechanically Derived Proton Hopping Rates (Q-hop MD). J. Chem. Phys. 2001, 115, 7993−8005. (30) Selvan, M. E.; Keffer, D. J.; Cui, S.; Paddison, S. J. A Reactive Molecular Dynamics Algorithm for Proton Transport in Aqueous Systems. J. Phys. Chem. C 2010, 114, 11965−11976. (31) Lee, S. H.; Rasaiah, J. C. Proton Transfer and the Mobilities of the H+ and OH− Ions from Studies of a Dissociating Model for Water. J. Chem. Phys. 2011, 135, 124505. (32) Day, T. J.; Soudackov, A. V.; Č uma, M.; Schmitt, U. W.; Voth, G. A. A Second Generation Multistate Empirical Valence Bond Model for Proton Transport in Aqueous Systems. J. Chem. Phys. 2002, 117, 5839−5849. (33) Wu, Y.; Chen, H.; Wang, F.; Paesani, F.; Voth, G. A. An Improved Multistate Empirical Valence Bbond Model for Aqueous Proton Solvation and Transport. J. Phys. Chem. B 2008, 112, 467−482. (34) Wang, F.; Voth, G. A. A Linear-Scaling Self-Consistent Generalization of the Multistate Empirical Valence Bond Method for Multiple Excess Protons in Aqueous Systems. J. Chem. Phys. 2005, 122, 144105. (35) van Duin, A. C.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (36) Chenoweth, K.; Cheung, S.; van Duin, A. C.; Goddard, W. A.; Kober, E. M. Simulations on the Thermal Decomposition of a Poly (Dimethylsiloxane) Polymer Using the ReaxFF Reactive Force Field. J. Am. Chem. Soc. 2005, 127, 7192−7202. (37) Russo, M. F.; Li, R.; Mench, M.; van Duin, A. C. Molecular Dynamic Simulation of Aluminum−Water Reactions Using the ReaxFF Reactive Force Field. Int. J. Hydrogen Energy 2011, 36, 5828−5835. (38) Strachan, A.; van Duin, A. C.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A., III Shock Waves in High-Energy Materials: The Initial Chemical Events in Nitramine RDX. Phys. Rev. Lett. 2003, 91, 098301. (39) van Duin, A. C.; Strachan, A.; Stewman, S.; Zhang, Q.; Xu, X.; Goddard, W. A. ReaxFFsio Reactive Force Field for Silicon and Silicon Oxide Systems. J. Phys. Chem. A 2003, 107, 3803−3811. (40) van Duin, A. C.; Zou, C.; Joshi, K.; Bryantsev, V.; Goddard, W. A.; Asthagiri, A.; Janik, M. A ReaxFF Reactive Force Field for Proton Transfer Reactions in Bulk Water and Its Applications to Heterogeneous Catalysis. Comput. Catal. 2013, 14, 223.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by a grant from the U.S. Army Research Laboratory through the Collaborative Research Alliance (CRA) for Multi Scale Multidisciplinary Modeling of Electronic Materials (MSME).

■

REFERENCES

(1) Ohmine, I.; Tanaka, H. Fluctuation, Relaxations and Hydration in Liquid Water. Hydrogen-Bond Rearrangement Dynamics. Chem. Rev. 1993, 93, 2545−2566. (2) Atkins, P.; de Paula, J. Physical Chemistry; Oxford University Press: Oxford, U.K., 2002. (3) Eisenberg, D. S.; Kauzmann, W. The Structure and Properties of Water; Clarendon Press: Oxford, U.K., 1969; Vol. 123. (4) Light, T. S.; Licht, S.; Bevilacqua, A. C.; Morash, K. R. The Fundamental Conductivity and Resistivity of Water. Electrochem. SolidState Lett. 2005, 8, E16−E19. (5) Agmon, N. The Grotthuss Mechanism. Chem. Phys. Lett. 1995, 244, 456−462. (6) Tuckerman, M. E.; Marx, D.; Parrinello, M. The Nature and Transport Mechanism of Hydrated Hydroxide Ions in Aqueous Solution. Nature 2002, 417, 925−929. (7) Eigen, M. Proton Transfer, Acid-Base Catalysis, and Enzymatic Hydrolysis. Part I: Elementary Processes. Angew. Chem., Int. Ed. Engl. 1964, 3, 1−19. (8) Marx, D. Proton Transfer 200 Years after Von Grotthuss: Insights from ab Initio Simulations. ChemPhysChem 2006, 7, 1848−1870. (9) Mohammed, O. F.; Pines, D.; Dreyer, J.; Pines, E.; Nibbering, E. T. Sequential Proton Transfer through Water Bridges in Acid-Base Reactions. Science 2005, 310, 83−86. (10) Voth, G. A. Computer Simulation of Proton Solvation and Transport in Aqueous and Biomolecular Systems. Acc. Chem. Res. 2006, 39, 143−150. (11) Ando, K.; Hynes, J. T. Molecular Mechanism of HCl Acid Ionization in Aater: ab Initio Potential Energy Surfaces and Monte Carlo Simulations. J. Phys. Chem. B 1997, 101, 10464−10478. (12) Asthagiri, D.; Pratt, L.; Kress, J. ab Initio Molecular Dynamics and Quasichemical Study of H+ (Aq). Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 6704−6708. (13) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. The Nature of the Hydrated Excess Proton in Water. Nature 1999, 397, 601−604. (14) Schmitt, U. W.; Voth, G. A. The Computer Simulation of Proton Transport in Water. J. Chem. Phys. 1999, 111, 9361−9381. (15) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M. Ab Initio Simulations of Water and Water Ions. J. Phys.: Condens. Matter 1994, 6, A93. (16) Vuilleumier, R.; Borgis, D. Quantum Dynamics of an Excess Proton in Water Using an Extended Empirical Valence-Bond Hamiltonian. J. Phys. Chem. B 1998, 102, 4261−4264. (17) Zundel, G. Hydrogen Bonds with Large Proton Polarizability and Proton Transfer Processes in Electrochemistry and Biology, in Advances in Chemical Physics, Volume 111; Prigogine, I., Rice, S. A., Eds; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 1999. DOI: 10.1002/ 9780470141700.ch1 (18) Agmon, N. Mechanism of Hydroxide Mobility. Chem. Phys. Lett. 2000, 319, 247−252. (19) Chen, C.; Huang, C.; Waluyo, I.; Nordlund, D.; Weng, T.-C.; Sokaras, D.; Weiss, T.; Bergmann, U.; Pettersson, L. G.; Nilsson, A. Solvation Structures of Protons and Hydroxide Ions in Water. J. Chem. Phys. 2013, 138, 154506. (20) Marx, D.; Chandra, A.; Tuckerman, M. E. Aqueous Basic Solutions: Hydroxide Solvation, Structural Diffusion, and Comparison to the Hydrated Proton. Chem. Rev. 2010, 110, 2174−2216. 6031


Article

The Journal of Physical Chemistry B (41) Chenoweth, K.; van Duin, A. C.; Goddard, W. A. ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 1040−1053. (42) van Duin, A. C.; Baas, J. M.; van de Graaf, B. Delft Molecular Mechanics: A New Approach to Hydrocarbon Force Fields. Inclusion of a Geometry-Dependent Charge Calculation. J. Chem. Soc., Faraday Trans. 1994, 90, 2881−2895. (43) Te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J.; Snijders, J. G.; Ziegler, T. Chemistry with ADF. J. Comput. Chem. 2001, 22, 931−967. (44) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (45) Achtyl, J. L.; Unocic, R. R.; Xu, L.; Cai, Y.; Raju, M.; Zhang, W.; Sacci, R. L.; Vlassiouk, I. V.; Fulvio, P. F.; Ganesh, P.; et al. Aqueous Proton Transfer across Single-Layer Graphene. Nat. Commun. 2015, 6, 6539. (46) Morita, M.; Takahashi, K. Theoretical Study on the Difference of OH Vibrational Spectra between OH−(H2O)3 and OH−(H2O)4. Phys. Chem. Chem. Phys. 2012, 14, 2797−2808. (47) Islam, M. M.; Ostadhossein, A.; Borodin, O.; Yeates, A. T.; Tipton, W. W.; Hennig, R. G.; Kumar, N.; van Duin, A. C. ReaxFF Molecular Dynamics Simulations on Lithiated Sulfur Cathode Materials. Phys. Chem. Chem. Phys. 2015, 17, 3383−3393. (48) Soper, A. K. The Radial Distribution Functions of Water and Ice from 220 to 673 K and at Pressures up to 400 MPa. Chem. Phys. 2000, 258, 121−137. (49) Lu, J.; Qiu, Y.; Baron, R.; Molinero, V. Coarse-Graining of TIP4P/2005, TIP4P-EW, SPC/E, and TIP3P to Monatomic Anisotropic Water Models Using Relative Entropy Minimization. J. Chem. Theory Comput. 2014, 10, 4104−4120. (50) Zhang, W.; van Duin, A. C. ReaxFF Reactive Molecular Dynamics Simulation of Functionalized Polyphenylene Oxide Anion Exchange Membrane. J. Phys. Chem. C 2015, 119, 27727−27736. (51) Knight, C.; Lindberg, G. E.; Voth, G. A. Multiscale Reactive Molecular Dynamics. J. Chem. Phys. 2012, 137, 22A525. (52) Botti, A.; Bruni, F.; Imberti, S.; Ricci, M.; Soper, A. Ions in Water: The Microscopic Structure of a Concentrated HCl Solution. J. Chem. Phys. 2004, 121, 7840−7848. (53) Botti, A.; Bruni, F.; Imberti, S.; Ricci, M.; Soper, A. Solvation of Hydroxyl Ions in Water. J. Chem. Phys. 2003, 119, 5001−5004. (54) Cappa, C. D.; Smith, J. D.; Messer, B. M.; Cohen, R. C.; Saykally, R. J. Nature of the Aqueous Hydroxide Ion Probed by X-Ray Absorption Spectroscopy. J. Phys. Chem. A 2007, 111, 4776−4785. (55) van Duin, A. C.; Merinov, B. V.; Han, S. S.; Dorso, C. O.; Goddard Iii, W. A. ReaxFF Reactive Force Field for the Y-Doped BaZrO3 Proton Conductor with Applications to Diffusion Rates for Multigranular Systems. J. Phys. Chem. A 2008, 112, 11414−11422. (56) Han, K. W.; Ko, K. H.; Abu-Hakmeh, K.; Bae, C.; Sohn, Y. J.; Jang, S. S. Molecular Dynamics Simulation Study of a PolysulfoneBased Anion Exchange Membrane in Comparison with the Proton Exchange Membrane. J. Phys. Chem. C 2014, 118, 12577−12587.

6032


Second-Generation ReaxFF Water Force Field: Improvements in the

Recommend Documents