Reconsidering Calcium Dehydration as the Rate-Determining Step in Calcium Mineral Growth.

The dehydration of cations is generally accepted as the rate-limiting step in many processes. Molecular dynamics (MD) can be used to investigate the dynamics of water molecules around cations, and two different methods exist to obtain trajectory-based water dehydration frequencies. Here, these two different post-processing methods (direct method versus survival function) have been implemented to obtain calcium dehydration frequencies from a series of trajectories obtained using a range of accepted force fields. None of the method combinations reproduced the commonly accepted experimental water exchange frequency of 10-8.2 s-1. Instead, our results suggest much faster water dynamics, comparable with more accurate ab initio MD simulations and with experimental values obtained using neutron scattering techniques. We obtained the best agreement using the survival function method to characterize the water dynamics, and we show that different method combinations significantly affect the outcome. Our work strongly suggests that the fast water exchange kinetics around the calcium ions is not rate-limiting for reactions involving dissolved/solvated calcium. Our results further suggest that, for alkali and most of the earth alkali metals, mechanistic rate laws for growth, dissolution, and adsorption, which are based on the principle of rate-limiting cation dehydration, need careful reconsideration.

Introduction

The dynamics of water molecules surrounding reacting cations is generally considered as rate-determining for the reactions involving these ions in aqueous solutions.[1−6] Almost 60 years ago, ultrasonic sound adsorption experiments performed by Eigen indicated a large range of water exchange frequencies around cations.[7,8] Years later, Nielsen found a correlation between the frequency of ions entering surface kink sites and the dehydration frequencies of the cations in this surface where the attachment frequency of the ions was approximately 10–3 times its dehydration frequency. Nielsen proposed that this difference was due to the ion diffusion from solution to a kink site, a process included in his overall attachment frequency. Nevertheless, as anions dehydrate faster, cations are considered to control the reaction rate.[9−11] These general findings have been incorporated into mechanistic models explaining crystal dissolution (e.g., Lasaga and Lüttge)[12,13] and crystal growth.[14−22] Essentially, in such models, it is assumed that the exchange frequency of water molecules between the first hydration shell around cation surface sites and bulk water controls the reaction rate.

The dynamics of water molecules can be investigated using molecular dynamics (MD) simulations, for example, of water molecules around dissolved ions,[3,23−25] at mineral–water interfaces,[26−29] and in amorphous materials.[30] Several studies have employed MD to study ion complexation and dynamics in water,[31−34] whereas the technique has also been used to analyze adsorption of small organics on mineral–water interfaces,[30,35] the structural changes in (water near) mineral interfaces,[36,37] and the influence of ion impurities on the kinetic stability of amorphous materials.[38]

In the case of calcite, combined experimental and computational work has shown that electrolyte ions impact water exchange frequencies and the crystal growth rate.[3] MD simulations also confirmed that the formation of the kink sites on the dissolving edge of the obtuse step of calcite is the rate-limiting step and the dissolution of this edge is thermodynamically favored[39] and the addition of the first units of calcium carbonate is enthalpically unfavorable.[40] Also, there is good agreement in the description of the structure of water at the calcite–water interface, as observed computationally,[41,42] using force-modulated atomic force microscopy,[43,44] small-angle X-ray reflectivity,[42] and neutron diffraction.[45]

Despite the promising agreement on the structure of the calcite–water interface obtained computationally and experimentally, large discrepancies remain in observations of water dynamics. Recent studies using ab initio MD (AIMD) calculations[46,47] show that exchange frequencies around dissolved calcium are almost 3 orders of magnitude higher than Eigen’s value. Comparably high frequencies have also been observed around dissolved calcium using incoherent quasi-elastic neutron scattering (IQUENS),[48] from which the authors concluded that the water exchange frequency around calcium ions is faster than 1010 s–1. Other MD simulations have indicated frequencies much closer to those measured by Eigen,[28] while classical MD simulations have resulted in water exchange frequencies around dissolved calcium and calcium in calcite surfaces[37] that are similar to those observed in ab initio MD and IQUENS. Besides, it is still debated if dehydration of the attaching species is rate-limiting in crystal growth[5] or that step-migration of carbonate ions is rate-limiting for kink nucleation.[49]

While an evaluation of the differences between the various experimentally observed frequencies is beyond the scope of the paper, the large range in computationally determined water exchange frequencies may be explained by differences between the force fields used in MD methods and/or the approach followed to obtain exchange frequencies from the MD output. Here, we present our evaluation of both potential areas for discrepancy. We have simulated identical systems of dissolved calcium using four sets of force fields, some of which have previously provided contradictory results. We next evaluated two generally accepted methods to determine water exchange frequencies. The results from our simulations and calculations provide new insights into the impact of the water dynamics around calcium on the reactivity of materials containing this cation. Our results dispute that calcium dehydration is the rate-limiting step in mineral formation and dissolution and suggest that other processes are rate-determining instead. We discuss implications for mechanistic growth models.

Computational Methods

Classical Molecular Dynamics

The MD simulations were conducted on a cubic simulation cell containing 809 water molecules and a single calcium ion. They were performed with DL_POLY 4.09,[50] which was modified to include a tail required for the SPC/fw-tail force field.[27] The Verlet leapfrog scheme and the Nosé–Hoover algorithm[51,52] were used to integrate the equations of motion with a time step of 0.2 fs and to maintain the temperature at 300 K with 0.5 ps relaxation time. The simulations were carried out in the isothermal-isobaric (constant NPT) ensemble using a Nosé–Hoover thermostat to maintain the average temperature and pressure at T = 300 K and P = 1 atm. The total simulation time for each of the systems was 144 ns including 200 ps of equilibration. An 8.0 Å cutoff was used for the van der Waals forces, except for the SPC/fw-tail force field in which the cutoff was defined at 9.0 Å with a tail from 6.0 Å.[27]

We tested four different force fields to analyze the water dynamics around the calcium cation: first, the shell-model (SM-lj) developed by de Leeuw and Parker and revised by Kerisit and Parker;[53,54] second, the simple point charge flexible water model (SPC/fw-tail).[55] Comparison between these first two is the focus of the main paper. However, in order to guarantee the robustness and comparability of these two water force fields, we also included an in-depth assessment of the dynamics and structures produced by SM-lj and SPC/Fw in comparison with two further force fields, that is, SWM4-NDP[56,57] and TIP3P,[58,59] which are also commonly used for simulations of calcium-containing systems or in electrolyte solutions. The full assessment of all four force fields can be found in the Supporting Information, which includes a study of the different structural parameters, including the dipole moment, and the three- and four-body order parameters. The dynamics of the system, such as diffusion coefficients and vibrational spectra, can also be found in the Supporting Information, as well as further details and a table with parameters for the different force fields used in this study.

Ab Initio Molecular Dynamics

AIMD simulations of one Ca2+ in 63 water molecules were conducted with the electronic structure code CP2K/Quickstep version 2.7.[60] CP2K implements the density functional theory (DFT) based on a hybrid Gaussian plane wave. The BLYP, PBE, and revPBE generalized gradient approximations for the exchange and correlation terms were used together with the general dispersion-correction DFT-D3.[61−65] Goedecker–Teter–Hutter pseudopotentials[66] were used to describe the core–valence interactions. All atomic species were represented using a double-ζ valence polarized basis set. The plane wave kinetic energy cutoff was set to 1000 Ry. k-sampling was restricted to the Γ point of the Brillouin zone. Simulations were carried out with a wave function optimization tolerance of 10–6 hartree that allows for 1.0 fs time steps with reasonable energy conservation. Periodic boundary conditions were applied throughout. Starting from the last configuration of 200 ps of classical MD simulation of 64 water molecules, one H2O molecule was replaced by one Ca2+, and we conducted 1 ns of classical MD simulations in the NPT ensemble (T = 300 K and P = 1 atm). The last configuration was taken as the starting point for AIMD simulations, which were carried out in the NVT ensemble using a Nosé–Hoover chain thermostat to maintain the average temperature at T = 300 K. The total simulation times for the different functionals were 500 ps for BLYP, PBE, and 40 ps revPBE. For comparison, AIMD simulations of pure bulk water were performed following the same approach.

Calculations of the Water Exchange Frequencies

To investigate the water dynamics around dissolved calcium ions, the mean residence time (MRT) of water molecules in the first hydration shell of the cation was determined, using both the “direct method”[46] and the Impey method.[67] An in-house Python code was developed to calculate the MRT using the direct method from the classical and ab initio MD trajectories (see the Supporting Information). The software tracks all the water molecules that leave or enter the first shell and, subsequently, stay outside or inside the first shell for a time longer than τDM = 0.5 ps.[46] In the direct method (DM), τDM = 0.5 ps is a generally accepted value, based on the average lifetime of a hydrogen bond between the solvent molecules.[3,26,46] The first hydration shell is taken to be within the distance from the Ca ion to the position of the first minimum in the Ca2+-Owater radial distribution function (RDF)[46] and is 3.5 Å for all force fields used here (see also Di Tommaso and de Leeuw[33] for the specific case of calcium).[47] The inverse of the average exchange frequency is taken to calculate the MRTwhere tsim is the total simulation time, CNav is the average coordination number for calcium (i.e., the number of water molecules around calcium obtained from integrating the RDF of calcium-oxygen), and Nexτ is the number of exchanges counted that lasted longer than 0.5 ps.

The survival function (SF) method (or the Impey method) determines the survival function for water molecules in the first shell. The first shell is defined in the same way as in the direct method, within 3.5 Å from the calcium ion. The SF method measures the residence time (RT) for each water molecule that enters the first shell. RT is recorded only when the leaving time τSF of a water molecule from the first shell meets the minimum time limit. τSF is advised to be set at 2 ps to guarantee that the water molecule truly moved to the bulk,[54,67] although a τSF of 1 ps has sometimes been used as well.[28] In this study, we adopt the initial proposal of the method and used τSF = 2 ps. After extracting all the residence times for all water molecules around calcium over the entire simulation period, the frequency (occurrences) of each RT is normalized by the total number of RT counted. When plotting the normalized frequencies against time, the survival function is constructed. This function gives the probability that a water molecule stays in the first shell as a function of time. A best fit exponential decay can then be plotted, and the MRT can be calculated by taking the integral of the function:

Results and Discussion

Four different force fields were assessed with respect to structural parameters, diffusion coefficients, and vibrational properties (see the Supporting Information for all details). Here, we focus on the results regarding water dynamics, as reflected by calculated exchange frequencies and mean residence times (MRTs) for water molecules in the first shell around the calcium cation, and the implications for the crystal growth rate.

Water Dynamics around Calcium Cations

The MRT values obtained using the direct method[46] and the survival function[67] are compared with other theoretical and experimental MRT evaluations in Table . There are large differences in MRT values obtained with different water force fields while using the same methodology to determine the residence time. For example, MRTSF from SM-lj (56 ps) is less than half the one obtained from MD simulations using the same Ca2+-Owater interaction potential but using the SPC/fw-buck water model (125 ps). The longest MRTSF is obtained based on the SPC/fw-tail (333 ps). The MRTSF value obtained from SM-lj output is the only value to agree with the experimental IQUENS results. The MRTDM follows a similar trend where the value obtained from SPC/fw-buck output reproduces the ab initio results.

Table 1

Mean Residence Times Calculated Using either the Direct Method (MRTDM) or the Survival Function (MRTSF)a

method	NH2O/Ca2+	MRTDM (ps)	CN	MRTSF (ps)	R2	ref
SM-lj	809	14.9	8.4	56	0.997	this work
SPC/fw-tail	809	116	7.0	333	0.999	this work
SPC/fw-buck	809	46	6.9	125	0.999	this work
BLYP-D3/CP2K	63	73.8	6.8			this work
revPBE-D3/CP2K	63	60.2	6.1			this work
CPMD	53	23.2	6			(47)
QM/MM-MD	199	42.5	7.8			(3)
SPC/fw-tail				200		(28)
SPC/E/TIP4P				158		(68)
SM-buck	2027	7.0				(37)
adsorption exp.		6310				[*][7,8]
IQUENS exp.		<100		<100		(48)

NH/Ca2+ is the number of water molecules per calcium in the simulation cell; the coordination number (CN) was used to normalize the exchange frequency for calculation of MRTDM, and the correlation coefficient R2 is for MRTSF. [*] Obtained using the exchange frequency of a molecule[4,5] using equation .

It was mentioned in earlier work that the higher number of water molecules in the first shell of calcium in SM-lj causes more exchanges and therefore shorter MRT, suggesting that the higher water density simulated with this force field is responsible for the higher coordination number (CN) of Ca2+.[28] We can add that the structure of the water around the cation in SM-lj is also more distorted with a CN = 8.4, whereas according to ab initio MD, the coordination geometry is octahedral (CN = 6[47]); the addition of extra molecules probably destabilizes this ideal geometry around Ca2+, contributing to the shorter MRTSF for SM-lj.

When comparing the two methodologies used to extract MRT from MD simulations, we note that the values obtained using the DM are consistently smaller than those produced by the SF, irrespective of the force field used in the MD simulations (Table ). Since the definition of the size of the first shell is the same for both DM and SF (3.5 Å, see Section ), the observed discrepancy can be explained by differences between the two methodologies to calculate (mean) residence times, as well as the different “time constraints” considered (τDM versus τSF).

In the direct method, an exchange event is counted when a water molecule enters or leaves the first shell around the cation, and this water molecule remains in/outside the first shell for at least 0.5 ps, a value for τDM that is commonly used[46] (see also Section ). The MRTDM is the mean of all residence times (RTs) observed, assuming that the RT distribution is a simple Gaussian distribution. However, the SF method shows that this distribution follows an exponential decay function (Figure a), suggesting that the Gaussian distribution is not a valid assumption. As a result, the MRTDM preferentially reflects fast exchanging molecules, resulting in shorter values of MRTDM.

Figure 1

(a) Probability distribution of different residence times (RT, black triangles) for first-shell water molecules obtained using SPC/fw-tail where the survival function is represented by the red and blue lines, which are two exponential fits. (b) Average distance between a coordinated water molecule and the calcium cation as a function of the residence time of that water molecule.

In contrast, in the SF method, the MRTSF takes into account the actual distribution of normalized RT frequencies. Such a distribution is shown in Figure a, based on MD simulations using SPC/fw-tail (SF based on results from SM-lj are shown in Figure S4). The MRTSF was calculated by taking the integral of the survival function (i.e., the integral of the fit to these graphs). In our survival functions for all RT probability distributions, a combination of two exponential functions gives the best fit to explain the data (e.g., Figure S6). This fit indicates that two underlying distributions are responsible for the observed MRTSF, one distribution corresponding to short MRTSF (<20 ps) and a second group of long MRTSF (>20 ps). A similar biexponential distribution of RTs for the calcium solvation shell was previously reported by De La Pierre et al.,[28] who related the short RTs to water librations and excluded those in their MRTSF calculations. Additionally, this bimodal MRTSF distribution may be due to differences in the average Ca2+-Owater distance. When the average distance is plotted against the RT (Figure b), one can see that, while generally all long RTs are found for water molecules at a distance of ≤2.40 Å, short RTs are found for water molecules in the full range of distances, from 2.50 Å to shorter than 2.40 Å. This suggests that another property affects RT besides the Ca2+-Owater distance, for example, the H-bonding network. To investigate if this is the case, we determined the percentage of H-bonds between first-shell water molecules and those outside the first shell (referred to as bulk water molecules) for short versus long RT water molecules (Figure ), following the definition of a hydrogen bond proposed by Chandra.[69] We used an in-house Python script to count the number of hydrogen bonds in the trajectory file. The molecules having a short RT show a larger number of H-bonds with bulk water molecules than those with long RT, suggesting that the H-bond network may affect residence time as well. Furthermore, SPC/fw-tail shows a larger number of H-bonds with bulk water than SM-lj (Figure ; Figure S7a). This more pronounced H-bond network around the cation leads to a less disturbed water structure and could therefore explain the longer MRTSF obtained from SPC/fw simulations.

Figure 2

Relative time a water molecule has H-bonds with bulk water molecules, during the 20 ps of each residence time (RT); the two gaussians represent the two different populations (short and long RT H2O).

If we group the distance-to-cation data based on their RTs, the average distance of the short RTs is 2.97 Å, while the average distance of the long RTs is only 2.48 Å (averaged over the first 20 ps after entering the first shell or, in the case of short RT, water molecules that left the first shell within these 20 ps, as long as they remained).

Besides the H-bond network and Ca2+-Owater distance, it is worth considering the structure of the first coordination shell. The average coordination number of SPC/fw water is found to be 7.3 (Table ). For 96% of the configurations, the short RT is only observed when the CN is +1 above this average (Table ). As mentioned above, CN is larger for SM-lj (Table ), and the coordination distances are on average also larger (Figure ), leading to greater distortion of the ideal bipyramid (octahedral) coordination for calcium. Moreover, the width of the first coordination shell is also controlled by the van der Waals interaction when comparing the different interactions using the same water model (SPC/fw-buck and SPC/fw-tail, Figure ). Using the Buckingham potential and the related parameters leads to a wider first shell, and the first maximum in the RDF is at a slightly larger distance from the cation compared to the tailed Lennard-Jones interaction proposed by Raiteri et al.[27] Based on the results for both force fields, it may be suggested that longer Ca2+-Owater coordination distances reflect first-shell water molecules in excess of the ideal for Ca2+ (CN = 6) that will be more liable to exchange than in the rigid model. Therefore, force fields yielding higher CN for calcium expectedly give rise to shorter average MRT.

Table 2

Coordination Number of Calcium in Relation to Residence Time of the Water Molecules in the First Shell Using SPC/fw-tail

	first shell configuration
CN	only long RT-water	1 ≥ short RT-water	all
6	1%		1%
7	75%	1%	70%
8	25%	96%	29%
9	0%	3%	0.16%

Figure 3

Radial distribution function reflecting Ca2+-Owater coordination distances for the different models.

Both the direct and survival function methods have advantages and disadvantages. The direct method results in shorter mean residence times because the fast exchanges are over-represented. Nevertheless, the relative differences in MRTDM we observe between the force fields tested (Table ) are very similar to those obtained with MRTSF. Moreover, the DM needs shorter simulation times than SF to obtain statistically meaningful output, and it can therefore also be used to investigate shorter trajectories generated by ab initio and QM/MM MD methods. However, due to the definition of a water exchange event in the DM method, it is likely to be less accurate. Nonetheless, this method can provide insight into the relative variations in local dynamics, such as a comparative assessment of the reactivity of (i.e., water dynamics around) different ions in aqueous solution or of the change in the reactivity of an ion in aqueous solutions with changing composition.[3]

In contrast, the survival function needs simulation times long enough to capture enough long residence times for the probability distribution function. Depending on the system investigated, this may be hundreds of nanoseconds, even microseconds. However, the SF results facilitate detailed differentiation of the MRT distribution in relation to CN, geometry, and H-bond network, but most importantly, the SF approach closely resembles the data post-processing, in terms of fitting exponential decay functions to determine MRT, used in experiments of Eigen and reviewed by Petrucci.[70,78] It is also similar to the theory behind the interpretation of IQUENS signals[48] to obtain MRT from neutron scattering data. As such, the SF results are in principle more comparable to available experimental data for 2 to 3 m calcium dissolved in water.

When we compare the MRT’s extracted from our MD results with those reported from experiments (Table ), neither methodology (irrespective of the force fields adopted) is able to reproduce the estimation of Eigen.[78] The MRT values obtained here are closer to IQUENS measurements,[48] ab initio,[47] and QM/MM MD simulations.[46] The different computer simulation methods still generate a range of water exchange frequencies, and comparability is hampered (i) because the SF cannot be used on the output generated by computationally demanding AIMD and QM/MM MD methods, which are generally limited to few hundreds of ps, and (ii) partly because of the intrinsic limitations of the classical MD force fields, such as the different ways of expressing the structure of the water with the different interatomic potentials.

The simulations conducted in this study considered a calcium ion in pure liquid water. However, Di Tommaso et al.[3] have shown that the frequency of water exchange in the first hydration shell of Ca2+ is reduced in the presence of halide ions; this effect is observed not only when the halide anions are directly coordinated to calcium but also when the alkali and halide ions are placed at or outside the second coordination shell of Ca2+. As natural aqueous solutions are far from pure water but rich in ions, in such environments, the water exchange dynamics around Ca2+ is likely to be influenced by the formation of contact and solvent shared ion pairs. An extended model, also considering the effect of solution speciation on the residence times of water around Ca2+, will therefore be part of future studies.

To summarize, residence times are affected by Ca-Owater distance and H-bonding network and therefore depend on which force field is used. The direct method can be used on MD as well as AIMD and QM/MM-MD output to obtain relative differences in dynamics. Based on the current comparison, we recommend the SPC/fw-buck force field to determine MRTDM since it most closely reproduces the ab initio results in MRTDM. However, the MRTSF based on this force field is too long compared to experimental IQUENS data where the SM-lj provides the best agreement. Furthermore, if at all feasible, we suggest that the use of the SF method to extract the MRTSF from the MD simulation output obtains the most accurate (absolute) value of the mean residence time.

Comparison with Water Exchange Kinetics around Carbonate Ions

Our atomistic simulations agree with recent experiments, supporting exchange frequency values for calcium that are 2 orders of magnitude larger than previously thought (Eigen). Consequently, the water exchange kinetics around Ca2+ approach has expected values around anions, for example, the carbonate ion (>1010 s–1).[37] To test this hypothesis, we ran a 13 ns simulation using SPC/FW and carbonate (details in the Supporting Information) to determine the MRT around a carbonate ion. We obtained an MRTSF of 8.4 ps (R2 0.999) and MRTDM of 9.2, which are both a factor of 5–40 shorter than for water residence times around calcium obtained with the same force field. As such, the difference in water exchange kinetics around the carbonate ion is not orders of magnitude faster than that around the calcium ion, although it remains noticeably faster. The slightly longer MRT around calcium is expected because the anion interactions are through hydrogen-bonds with a relatively shorter lifetime (∼2.5 ps)[71] than the metal-oxygen interactions as estimated by ab initio MD.[47] Furthermore, the free energy of hydration for calcium is larger than that for carbonate (−1447 ± 2 and −1301 ± 2 kJ mol–1, respectively),[72] indicating a stronger interaction between calcium and water than between carbonate and water.

Implications for Crystal Growth Rate

In order to discuss the impact of a smaller difference in water dynamics around calcium and carbonate ions than has been suggested thus far, we relate to crystal growth via ion-by-ion attachment. While there is strong evidence of ion-pair attachment during growth, for example, by the attachment of CaCO3(aq) pairs during calcite growth,[5,49,73,74] such a growth mechanism cannot explain the dependency of the growth rate on the solution stoichiometry at a constant degree of supersaturation observed in bulk[75,76] and microscopically.[17−19] According to thermodynamic equilibrium calculations, the concentration of the CaCO3(aq) ion pair is constant at a constant degree of supersaturation, irrespective of the solution stoichiometry, and the growth step velocity and bulk growth rate dependency on solution stoichiometry are therefore assumed to be caused by an ion-by-ion growth mechanism.[16,20−22]

In the ion-by-ion growth mechanism, the ion attachment frequencies are controlled by several, potentially simultaneous processes: cation and anion diffusion to lattice sites (surface anions and cations, respectively), dehydration of solvated ions, and surface lattice sites. Various experiments investigating calcite[10−14] surface growth strongly suggest that these processes are different for the structurally distinct growth steps at the calcite surface.[17−21] While obtuse step edges propagate through direct attachment of ions delivered via solution diffusion, the acute steps grow by attachment of ions delivered through surface diffusion. Moreover, Hong and Teng[18] have proposed, based on their atomic force microscopy observations, that obtuse step growth is limited by calcium dehydration and acute step propagation is limited by surface adsorption, reorientation, and rearrangement of carbonate ions. We would like also to indicate that interpretation of the kinetics of crystal growth, deconvolved from experimental data, using molecular simulations computer models is very challenging because the mechanisms controlling these processes (nucleation and growth) can be significantly more complex than those considered in simplified computer models.

One of the first direct relationships between crystal growth kinetics and ion dehydration was presented by Nielsen.[11] He obtained the values for the cation dehydration frequency νdh related to various minerals from Petrucci et al.[8] and obtained the correlation with ion integration frequency νin, as illustrated with blue dots and line in Figure . Expressing this correlation as an Eyring equation, it was found that the integration rate is approximately a thousand times smaller than the dehydration frequencywhere νD is the diffusion frequency, ΔGD≠ is the Gibbs free energy of diffusion, and the ratio between the diffusion rate and , with k and h being the Boltzmann and Planck’s constants, respectively, is .[11]

Figure 4

Logarithm of the cation integration frequency versus the logarithm of the cation dehydration frequency. Circles: as determined from parabolic function fitted to the dehydration frequency;[11] triangles based on IQUENS measurements;[48,77] cross based on QM/MM-MD simulations and DM;[3] gray square based on ab initio calculations (BLYP-D3 functional) and DM (this study); diamond from classical MD using SF (this study).

This correlation (Figure ) supported the notion that the cation dehydration is the limiting factor for growth (and dissolution). The cations retain water molecules longer than anions and are therefore rate-limiting.[3,28,37]

In general, the MRTs obtained with DM and SF support this notion, although the difference between cation and anion dehydration is not as extreme (orders of magnitude) as considered previously. Both cation and anion dehydration frequencies are in the order of, or slightly smaller than, the reported diffusion rate (νD) ∼1010 s–1.[11] However, adding our water exchange frequencies as well as the other values from Table to Figure reveals that the correlation previously reported by Nielsen[11] is not very strong, in particular in the case of cations with a dehydration frequency higher than 108 s–1 (Figure ).

Furthermore, if we use the updated values of dehydration frequencies, measured with the IQUENS method, and the values of the ion integration frequency νin, then the term (eq ) is not 10–3 for calcium anymore, but a factor of 10–5 is needed to correlate νdh with νin. This indicates that another process is affecting νin. Furthermore, the dehydration frequency derived from ab initio MD (Table ) is smaller than the diffusion frequency (νD), which suggests that the diffusion of the cation into the lattice site from a neighboring position is the rate-limiting step, in agreement with various experimental works.[17,18,78] Additionally, Petrucci[8] proposed based among others on data from Eigen and Wilkens[7] that, for some ions, including Ca2+, dehydration was not the limiting factor. Petrucci concluded that substitution of water by any ligand is the rate-limiting step only if the ligand’s exchange rate is less than 107 s–1. The latter is undoubtedly the case for most cations with small radii (<100 pm) such as Mg2+, Ni2+, Co3+, Zn2+, and Fe2+. Petrucci[8] compared the substitution rates of different ligands; for Mg2+, which has a very low water exchange frequency compared to Ca2+, there is indeed no difference in the reaction rate with a range of ligands, indicating that the dehydration of magnesium is the rate-controlling step. In contrast, he observed different reaction rates between Ca2+ and the same selection of ligands, suggesting that calcium dehydration is not the rate-determining step. Alternatively, it may also suggest that calcium dehydration can be affected differently by different ligands, as has been shown for halide ions.[3]

To summarize, the water exchange kinetics of Ca2+ and CO32– are more comparable than previously thought, which is substantiated by maximum bulk growth rates and step propagation rates in near-stoichiometric solutions.[16,75] Furthermore, the MRTs reported here, combined with the experimental IQUENS values for the dehydration frequency of calcium,[48] and the ligand experiments[8] shed a different light on the most likely rate-limiting step for calcium carbonate crystal growth. Rather than cation dehydration, (surface) diffusion of the cation into the lattice site from a neighboring position is the rate-limiting step, which is in agreement with the work by Hong and Teng,[18] Bracco et al.,[17] and Stack et al.[78] Alternatively, our results may support the proposition of De La Pierre et al.[49] that the attachment of the carbonate ion to the kink site is the rate-determining step of the calcium carbonate crystal growth.

Conclusions

The mean residence time of water molecules around dissolved calcium ions, obtained from our classical MD simulations, does not agree with the experimental mean residence time values around calcium estimated by Eigen,[7] as used in Nielsen’s analysis of crystal growth[11] and adopted by many recent calcite growth rate studies. In contrast, our MRT values are in the same order of magnitude as more accurate ab initio MD, whereas QM/MM-MD agrees with the limit obtained from IQUENS experiments. This finding indicates that calcium falls in the group of cations with fast water exchanges and its dehydration is therefore not rate-limiting during calcite crystal growth. This study also revealed that, of the two methods we tested to calculate the MRT from MD simulation output, only the survival function method provides accurate absolute values and can distinguish two different groups of water molecules around calcium, revealing that ion-water distance and H-bonding networks affect MRT values. Additionally, our work suggests that the direct method can be used to obtain relative differences in MRT, rather than accurate values, and is most suitable to post-process output from MD simulations of larger systems and ab initio calculations, as it can be implemented on shorter production runs.