Ion dissolution mechanism and kinetics at kink sites on NaCl surfaces

Significance Ion desolvation limits the rate of elementary ion-attachment steps during ionic crystal growth from aqueous solution. However, the desolvation mechanism and solvent coordinates are unknown, so crystal dissolution/growth rates remain difficult to compute. This work identifies solvent reaction coordinates for detachment/attachment of Cl− and Na+ from kink sites on NaCl crystals. We find that the larger Cl− limits the overall rate of NaCl dissolution in water because its size prevents water from accessing the kink site during detachment. Our findings show how ion and solvent characteristics influence the kinetics of NaCl growth and dissolution. This work also provides generalizable solvent coordinates and rate expressions for efficient calculations of crystal growth and dissolution rates. Desolvation barriers are present for solute–solvent exchange events, such as ligand binding to an enzyme active site, during protein folding, and at battery electrodes. For solution-grown crystals, desolvation at kink sites can be the rate-limiting step for growth. However, desolvation and the associated kinetic barriers are poorly understood. In this work, we use rare-event simulation techniques to investigate attachment/detachment events at kink sites of a NaCl crystal in water. We elucidate the desolvation mechanism and present an optimized reaction coordinate, which involves one solute collective variable and one solvent collective variable. The attachment/detachment pathways for Na+ and Cl− are qualitatively similar, with quantitative differences that we attribute to different ion sizes and solvent coordination. The attachment barriers primarily result from kink site desolvation, while detachment barriers largely result from breaking ion–crystal bonds. We compute ion detachment rates from kink sites and compare with results from an independent study. We anticipate that the reaction coordinate and desolvation mechanism identified in this work may be applicable to other alkali halides.

C rystal growth rates at mild supersaturations are determined by rates of solute attachment and detachment at kink sites (1,2), which are special surface sites that resemble half-lattice positions. The barrier for solute detachment is believed to arise mostly from breaking solute-crystal interactions (3)(4)(5), but it cannot entirely result from this; otherwise attachment would be a barrierless, diffusion-controlled process. Thus, the attachment barrier (and by microscopic reversibility, a part of the detachment barrier) must result from changes in solvation.
In this work, we investigate the attachment/detachment of Na + and Cl − at kink sites of a NaCl crystal in water. Accurate rates can be computed using only ion positions without a complete understanding of the attachment/detachment mechanism (15,22) [e.g., via trajectory-based rate calculations (24)], but these calculations provide only numerical rate estimates. They do not provide the kind of mechanistic insight that can be used to construct generalizable theories for other conditions and/or other ions. Therefore, we examine the role of desolvation by identifying important solvent collective variables (CVs) for attachment/detachment at kink sites (25). We elucidate the desolvation mechanism via transition path sampling, likelihood maximization, and committor tests, from which we find that kink attachment can be described with one solute collective variable and one solvent collective variable. By computing the respective free energy surfaces, we gain additional insight on similarities and differences between Na + and Cl − detachment. Finally, we compute ion detachment rates from kink sites and make comparisons with other studies.

Methods
Simulation Setup. We use SPC/E (extended simple point charge) (26) water model and the  force field for crystalline NaCl as well as Na + and Cl − ions in solution. As suggested in earlier studies, the Lennard-Jones interactions are cut off at 10Å and shifted to zero, water geometry is held rigid with SHAKE (28), and long-range electrostatics are included via the particle-particle particle-mesh Ewald method. According to this force field, the solubility of NaCl in water is 3.6 mol/kg water (29)(30)(31). A brief table of its other properties can be found in SI Appendix. All properties and trajectories were computed at 298 K in the NVT ensemble (i.e., constant number of particles, volume, and temperature) with a Langevin thermostat in LAMMPS (large-scale atomic/molecular massively parallel simulator) (32). Details on solution equilibration and NPT (constant number of particles, pressure, and temperature) simulations to determine appropriate box dimensions can also be found in SI Appendix.
The system consists of a NaCl crystal slab (854 formula units) with pure water (1,800 molecules) on both sides. The crystal is tilted [similar to that in Klimes et al. (18) and Kolafa (33)] such that each exposed face contains two edges, each with a single kink site, resulting in an ever-kinked surface in periodic space (Fig. 1A). See SI Appendix for details on the crystal rotation. This system design eliminates kink-kink interaction effects because an ion pair can be removed without altering the distance between kink sites. Additionally, the slab thickness (ca. 36Å) is sufficient to mitigate finite size effects, which can alter computed properties like solubility (34). Therefore, we expect this system to adequately mimic the bulk crystal in the slab interior without fixing atom locations in an inner crystal layer during simulations.

Reaction Coordinate Identification.
To identify important solvent reaction coordinates, the flexible-length (36) aimless shooting variety of transition Significance Ion desolvation limits the rate of elementary ion-attachment steps during ionic crystal growth from aqueous solution. However, the desolvation mechanism and solvent coordinates are unknown, so crystal dissolution/growth rates remain difficult to compute. This work identifies solvent reaction coordinates for detachment/attachment of Cl − and Na + from kink sites on NaCl crystals. We find that the larger Cl − limits the overall rate of NaCl dissolution in water because its size prevents water from accessing the kink site during detachment. Our findings show how ion and solvent characteristics influence the kinetics of NaCl growth and dissolution. This work also provides generalizable solvent coordinates and rate expressions for efficient calculations of crystal growth and dissolution rates.
Clin kink solvated Clkink B Na + in kink solvated Na + kink path sampling (TPS) is used to separately harvest an ensemble of transition paths for Na + and Cl − (NVT simulations with a Langevin thermostat). This variant of TPS has one adjustable parameter, δt (time between candidate shooting points), which was set to 30 fs and yielded acceptance rates of 41% for Na + and 36% for Cl − . As shown in Fig.  1, the ion-attached basin corresponds to small ion-kink separation distances. The ion-detached basin is characterized by larger ion-kink separations and a water molecule in the kink site. See SI Appendix for details regarding basin definitions and determination of the ion-kink separation distance. We use the inertial version of likelihood maximization (37) (iLmax) to rank the accuracy of an extensive list of trial reaction coordinates of the form q = c 1 r ion,k + c 2 q solvent − q ‡ , where r ion,k is the ion-kink separation distance and q solvent is a solvent CV. The adjustable parameters c 1 and c 2 determine the direction/gradient of the trial coordinate in (r ion,k , q solvent ) space, while q ‡ is the predicted value of c 1 r ion,k + c 2 q solvent at the optimal dividing surface location. The trial solvent CVs can be generally characterized as distances, local water densities, coordination numbers, joint coordinations, solvent structures, and solvation energy gaps. Parameters within the trial solvent CVs, e.g., cutoffs for computing coordination numbers, are also optimized to maximize reaction coordinate accuracy (38). These simple trial reaction coordinates help to identify a clear role for solvent degrees of freedom in the mechanism.
We chose r ion,k as the solute CV not only because it is commonly used, but also because it outperformed (via iLmax) alternative choices of the solute CV (i.e., coordination number, height above surface, etc.). r ion,k even outperforms pairs of other solute CVs, such as a coordination number and height above the surface (17). However, r ion,k alone is not a good reaction coordinate, because it results in a hysteresis for Cl − (although not for Na + ) as shown in SI Appendix. The hysteresis results from poor sampling of the water molecule entering and leaving the kink site.
In iLmax, the reaction probability is p RX = (1 + erf[q + c Vq ])/2, where c V is an adjustable coefficient andq is the reaction coordinate velocity, determined from a finite difference after a single NVE (constant number of particle, volume, and energy) time step from the shooting point. The reaction coordinates are scored relative to the baseline r ion,k coordinate, where ln L[q] is the log-likelihood for reaction coordinate q and the Bayesian information criteria (BIC) helps determine whether differences between reaction coordinates are significant [BIC = ln(N traj )/2 = 4.95 in this work; N traj is the number of TPS trajectories]. The baseline log-likelihoods from iLmax for Cl − and Na + are −12,402 and −13,204, respectively. Committor tests are performed at the dividing surface for Na + and Cl − detachment. Configurations within the transition-state ensemble are collected from umbrella-sampling simulations, with each configuration separated by at least 2 ps. Two hundred configurations are tested, each with 20 trajectories initiated with different Boltzmann distributed velocities (accounting for rigid water molecules), and integrated under NVE conditions until the trajectory enters the ion-attached or ion-detached basin. The committor, p B , for each configuration is the fraction of trajectories which enter the ion-detached (product) state.
Rate Calculations. Free energy calculations are performed via umbrella sampling (39) with parabolic bias potentials on both r ion,k and the solvent CV ρw (discussed in Reaction Coordinate). This technique consists of separate simulations (i.e., umbrella windows) constraining the system to different parts of the transition between the ion-attached and -detached states via an added bias potential. Within each window, data were collected from at least 3 ns of thermostated molecular dynamics simulation. The overlapping distributions from the collection of biased windows and the weighted histogram analysis method (40,41) yield the overall free energy surface. Note that free energies are normalized by k B T = β −1 in this work, where k B is Boltzmann's constant and T is absolute temperature.
Rate constants were obtained from generalized transition-state theory (42) using the optimized reaction coordinate where κ is the transmission coefficient that accounts for dynamical barrier recrossing events and k TST is the transition-state theory rate constant. Note that the overall rate constant is independent of q. The effective positive flux method (43) was used to calculate κ. Briefly, NVE trajectories are initiated from a Boltzmann distribution of configurations at the dividing surface. Trajectories are tested to determine whether they cross in the positive (P) direction, whether they are the first (F) crossing, and whether they are effective (E) at reaching the product (detached) state. A PFE trajectory has χ = 1 and otherwise χ = 0. The transmission coefficient is where ... ‡ indicates a conditional equilibrium average over the dividing surface. Additional details on the rate constant calculation are provided in SI Appendix.

Desolvation and Kink Attachment Mechanism
Reaction Coordinate. According to iLmax, the mechanism for solute attachment/detachment at kinks involves a concerted exchange of solute and solvent for the kink site. The optimal reaction coordinate for Cl − and Na + combines r ion,k (solute CV) and the local water density near the kink site (solvent CV), where the sum is over water molecules, σ is an effective probe radius (given in Table 1), rw is the location of a water molecule (best characterized by the oxygen atom), and r probe is the probe location. We tested different probe locations in the vicinity of the kink site; for Na + , the kink site is the best probe location. For Cl − , the optimal probe location is slightly offset from/deeper into the kink site, which we anticipate to more closely reflect the preferred oxygen location in the kink site (SI Appendix). These differences are likely a function of different water coordination environments and different sizes of the kink sites (due to different ion sizes). By first comparing pure-solute and pure-solvent reaction coordinates, we can assess the importance of the solute and solvent (recall ∆ ln L[r ion,k ] = 0). For each ion, the pure-solvent reaction coordinate outperforms the pure-solute reaction coordinate (∆ ln L[ρw] = 11 for Cl − and 34 for Na + ). The synergistic benefit of accounting for the solute and solvent results in an increase of ∆ ln L by over 100, as shown in Table 1 along with reaction coordinate coefficients. Consistent with the observed hysteresis for Cl − , but not for Na + in F (r ion,k ), the results in Table 1 show that the solvent CV is more important for Cl − than for Na + . Additionally, the significantly nonzero cV coefficients indicate that the dynamics of attachment/detachment of Na + and Cl − are partially inertial, like the previously reported inertial association of these ions with each other (44).

Committor Tests.
A preliminary analysis of committors on the iLmax dividing surface revealed higher pB values on one end of the dividing surface than on the other. We used the solvent CV identified by iLmax, but adjusted the coefficients c1 and c2 [i.e., tilted the dividing surface in (r ion,k , ρw) space] to correct the iLmax prediction (see SI Appendix for details). The committor distributions at the optimized dividing surface for kink attachment/detachment are shown in Fig. 2A for Cl − and peaked at pB = 1/2, with no contributions in the limits of 0 and 1. Therefore, our reaction coordinate is not ideal. However, it is a significant improvement on r ion,k (Fig. 2 C and D), indicating that ρw is mechanistically important. It also greatly improves upon the only previously reported committor distribution for a similar process, where few realizations near pB = 1/2 were observed (16). Upon ion detachment and solvation of the kink, the ion is in an adsorbed state near the kink site ( Fig. 1 C and E); i.e., the ion is not a fully hydrated "aqua ion" (45). Similar "docked" or partially hydrated ion configurations near CaCO3 step edges have been predicted for CO 2− 3 in water (9). We have computed the probability, p sol , that the docked ion tends to go into bulk solution rather than reattach to the kink site. We find that p sol = 1 for Cl − and p sol = 0.86 for Na + , so detachment is essentially complete for both ions once they reach the docked state. This suggests that diffusion to/from the surface (via surface or from bulk solution) would be a negligible resistance to ion attachment/detachment. Our results confirm that the rate-limiting step is exchange of the ion between the adsorbed and the kink site and that solvent CV is important for the rate-limiting step.
Free Energy Surfaces. We gain additional mechanistic insight on desolvation and kink attachment from the computed free energy surfaces, shown in Fig. 3A for Cl − and Fig. 3B for Na + (note they have different color scales). The surfaces are plotted in terms of ln ρw to clearly visualize the impact of ρw. The basins are labeled according to whether the ion or a water molecule occupies the kink site. A high free energy ridge on the free energy surface for Cl − separates the solvated and desolvated kink sites (ca. r Cl,k = 1.8Å, ln ρw = −5.5 in Fig. 3A). The ridge, which is not present on the free energy surface of Na + , explains why sampling along r ion,k gave a hysteresis for Cl − and not for Na + . Committor estimates for configurations in the predicted transition-state ensemble are grayscale points to represent the committor (pB = 0 is black and pB = 1 is white). The intermingled black and white points on the dividing surface shows that some third component of the reaction coordinate remains unknown. The dotted lines indicate the minimum free energy path (MFEP) between the attached and detached states, as discussed below.
Overall, the free energy surfaces show that the kink attachment/detachment processes are similar between the two ions. As the ion detaches (increasing r ion,k ), a void grows at the kink site (i.e., water cannot immediately enter the kink). Similarly, the ion attachment process is initiated by void growth at the kink, resulting in kink desolvation (decreasing ln ρw). The work required to grow these voids will be reflected in the free energy barriers for ion detachment and ion attachment/kink desolvation.
As equivalent ion-crystal bonds are broken during detachment, the ion-attached basins for Cl − and Na + appear quite similar. However, the free energy surfaces begin to differ near the dividing surface as solvent motion begins to become important. Additionally, the different dividing surface locations in r ion,k show that a larger void must be grown for Cl − detachment than for Na + detachment. We hypothesize that this occurs because the larger Cl − must move farther away before a water molecule can access the kink site.
The kink desolvation work (i.e., moving to more negative ln ρw in Fig. 3 at fixed r ion,k ) appears to be markedly different between the two ions, with the Cl − kink being much more difficult to desolvate. We anticipate that this occurs due to the different water coordination environments in the vacant kinks, which consist of a "nest" of three counterions. As water binds more strongly to Na + than to Cl − , we expect, and our results (Fig. 3) show, that water binds more strongly to the vacant Cl − kink.
By analyzing a sample of TPS trajectories, we find that multiple attachment/detachment pathways exist; i.e., water can vacate the kink in different directions. However, a projection into the (r ion,k , ln ρw) space gives a free energy surface that appears to have just one pathway, i.e., along the minimum free energy path. The MFEPs for Cl − and Na + indicate distinct pathways; however, the MFEPs do not yield dynamical insights. In SI Appendix, we show the evolution of select trajectories on the free energy surfaces. Along a typical detachment pathway, r Cl,k increases until a sudden "jump" in ρw as the water enters the void left behind. The rapid change in ρw suggests a nonadiabatic rate formulation like that of Schenter and coworkers (46) for future attempts to develop a general kink attachment/ detachment model. The two-step nature of the Cl − attachment/ detachment pathway is likely due to water being weakly coordinated to Cl − . Conversely, the pathways for Na + indicate a more simultaneous motion of ion and water, which occurs because Na + drags water molecules with it due to its strong coordination to oxygen.

Kink Detachment Rates
The free energy profiles for kink attachment/detachment in Fig. 4 are obtained by projecting the free energy surfaces from (r ion,k , ρw) onto q. The transition state is located at q = 0 with the ion-attached state at q < 0 and the ion-detached/solvated kink state at q > 0. We note that projecting the free energy surfaces solely onto r ion,k decreases the free energy maximum by 2 kB T for Cl − and 1 kB T for Na + compared with F ‡ = F (q = 0), as expected from variational transition-state theory (47). The barrier for Cl − detachment is about 3.5 kB T greater than that for Na + detachment, which we attribute to the previously discussed larger void that must be formed before water enters the kink site.
Detachment rates are computed via where ∆F ‡ is the free energy difference between the transition state and the ion-attached state where q is a unit length of q. The average reaction coordinate velocities at the dividing surface are |q| ‡ = 0.0024 fs −1 and 0.0026 fs −1 for Na + and Cl − , respectively. For more details on arriving at Eq. 4, see SI Appendix. We note that k − is formally a detachment rate constant. It is equivalent to the detachment rate because the detachment process is zeroth order in solute concentration (48). From effective positive flux we find that κNa = 0.10 ± 0.01 and κ Cl = 0.25 ± 0.02. These lower transmission coefficients may result from intrinsic friction in the dynamics (i.e., a rigid solvent structure) and/or reaction coordinate error. Alternative reaction coordinates may improve the committor distributions, but we anticipate little improvement upon κ because Fig. 4, Inset shows that even the ideal dividing surface (pB = 1/2) is crossed several times along typical detachment trajectories. Thus, the maximum possible transmission coefficient is less than unity.
We find that the Cl − detachment rate (k − Cl = 6.3 × 10 6 s −1 ) is over an order of magnitude slower than the Na + detachment rate (k − Na = 1.6 × 10 8 s −1 ). Thus, Cl − detachment limits dissolution kinetics. From the direct simulation of NaCl dissolution in water by Lanaro and Patey (49) (using identical force fields), we estimate the average detachment rate of ions from a kink site is 4.5 × 10 7 s −1 (deconvolution of rates for Na + and Cl − is not possible with the reported data). This is in good agreement with the (harmonic) average detachment rate in this work, 1.2 × 10 7 s −1 . Note that Lanaro and Patey (49) examined dissolution of an approximately spherical NaCl nanoparticle with a rougher surface than our stepped surface. Therefore, a somewhat faster detachment rate is expected from their study.
The dissolution of NaCl has been the focus of many studies (17)(18)(19)(49)(50)(51)(52)(53)(54)(55)(56)(57), in particular for thin water films that occur at moderate relative humidities. While some works focus on direct dissolution simulations, a few use ab initio molecular dynamics to compute detachment barriers for ions at edge corner sites (i.e., the corner of a 2D nucleus) (17,19). While Liu et al. (17) found a similar Cl − detachment barrier compared with our results (although the barriers for Na + differ), this most likely results from a fortuitous cancellation due to the different surface sites, force fields, and reaction coordinates.
Stack et al. (15) determined the Ba 2+ detachment rate from a barite step site was approximately 10 4 s −1 . Darkins (23) computed similar detachment rates of CO 2− 3 from two different kink sites on calcite. We anticipate that the three orders of magnitude difference from our computed rates largely results from different ion charges, which impacts the bonds that must be broken to detach. In other works, Stack et al. (3) and Stack and coworkers (10) inferred even slower kink detachment rates by fitting models to experimental measurements of CaCO3 and MgCO3 crystal growth. Dogan et al. (22) found kink detachment rates on the order of 10 5 s −1 via molecular simulations for aspirin and alpha-lactose monohydrate crystals. While elec-trostatics cannot explain these slower rates, the differences may result from slower rotational motion of these larger molecules and/or the more collective nature of breaking many atomatom interactions. It is not yet clear how different solute sizes, charges, shapes, and solvent characteristics impact attachment/ detachment rates.

Conclusion
Ionic crystal growth and dissolution rates near saturation in aqueous solution are limited by desolvation and attachment of solute growth units at kink sites. We have examined the ion attachment mechanism in NaCl crystal growth/dissolution with path sampling and reaction coordinate identification techniques. The mechanism involves a competition between solute and solvent for the kink site. Likelihood analyses and committor distributions show that the local density of water molecules in the kink site is a mechanistically important solvent CV. Free energy surfaces, transition-state theory rates, and transmission coefficient corrections provide detachment rates that agree with independent estimates from direct simulations of dissolution (49). While the absolute dissolution rates may be sensitive to forcefield details, mechanistic conclusions are expected to be more robust. For example, Cl − detachment limits the overall kinetics because, compared with Na + , the larger Cl − must move farther away from the kink site before water can solvate the kink. The optimized solvent CV identified in this work is similar to a solvent CV that was accurate for ion-pair dissociation across several alkali halides (44,58). Thus, the reaction coordinates identified here may also be useful for understanding other alkali halide crystals. More generally, with the strategy used in this work it should be possible to fully understand how ion sizes, charges, shapes, and interactions with the solvent influence the rates of attachment and detachment at kink sites in ionic crystal growth and dissolution.