Molecular Simulation of Nucleation and Dissolution ofAlkali HalidesbyGabriele LanaroA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Chemistry)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)March 2017c© Gabriele Lanaro 2017AbstractThe process of crystal nucleation, despite being so fundamental and ubiqui-tous in industrial and natural processes, is still not fully understood because of itsstochastic nature, and the high spatial and temporal resolution needed to observe itthrough experiments. This thesis investigates several aspects of nucleation throughthe use of molecular dynamics, a computational technique that is able to simulatesystems up to ∼ 1012 atoms (as of today’s computational power).The projects in this thesis focus on the nucleation from aqueous solution ofalkali halide salts, with supplementary studies on the related processes of dissolutionin water, and crystallization from the melt.The mechanism of NaCl nucleation from solution is examined in Chapter3 by direct simulation. The NaCl supersaturated solution was found to containmany small ionic clusters that continuously form and disappear from solution untilone (or more) of them nucleates and grows irreversibly. An original method wasdeveloped to detect and follow clusters in time, producing results useful in thestudy of their characteristics and lifetimes. Most importantly, it was found thatthe lifetime of transient clusters is ∼ 1 ns, and that both the cluster lifetime andnucleation probability are significantly higher if the cluster is more geometricallyordered. The dissolution of NaCl crystals was also investigated. The process wasfound to happen in stages, is characterized by an activation barrier, and can beiidescribed by a simple rate law.The crystal nucleation of LiF from supersaturated solution was observed, inour simulations, only at high pressure and temperature. The growth rate for analready nucleated crystal was found to have a temperature dependence that followsthe Arrhenius law, and further evidence suggests that the reason for such behavioris the high activation energy required to dehydrate the ions.The crystallization from the melt of the Joung-Cheatham and Tosi-Fumimodels for lithium halides was also investigated. We found that, for the Tosi-Fumimodel, all lithium halides crystallize as wurtzite. For the Joung-Cheatham model,LiF and LiCl crystallize as rock salt, while LiBr and LiI crystallize as wurtzite.iiiPrefaceThe Chapters in this thesis are based on work that has been, or will be,published by G. Lanaro and G. N. Patey.A version of Chapter 3 has been published by G. Lanaro and G. N. Patey,The Birth of NaCl Nanocrystals: Insights from Molecular Simulation, The Journalof Physical Chemistry B 120, 34 (2016), and was inspired by previous work by D.Chakraborty and G. N. Patey, How Crystals Nucleate and Grow in Aqueous NaClSolution, The Journal of Physical Chemistry Letters 4, 4 (2013).Versions of Chapters 4 and 5 are currently being prepared for publication.Chapter 6 is based on a publication by G. Lanaro and G. N. Patey, MolecularDynamics Simulation of NaCl Dissolution, The Journal of Physical Chemistry B119, 11 (2015).In all projects my contribution consisted of performing all of the simulationsand data analysis, writing codes for visualization and analysis, formulating hypoth-esis, designing the projects along with my supervisor, and producing initial draftsfor the manuscripts. Editing, reviews and production of final drafts were carriedout together with my supervisor.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Phase Transitions and Nucleation . . . . . . . . . . . . . . . . . . . . 11.2 Classical Nucleation Theory . . . . . . . . . . . . . . . . . . . . . . . 31.3 Computer Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Models and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.3 Interaction Potentials . . . . . . . . . . . . . . . . . . . . . . . . . . 142.3.1 Van der Waals Forces . . . . . . . . . . . . . . . . . . . . . . 142.3.2 Electrostatic Forces . . . . . . . . . . . . . . . . . . . . . . . 162.3.3 Water and Ion Models . . . . . . . . . . . . . . . . . . . . . . 162.3.4 Infinite Systems . . . . . . . . . . . . . . . . . . . . . . . . . 172.4 System Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . 23v2.4.1 Liquid Systems . . . . . . . . . . . . . . . . . . . . . . . . . . 232.4.2 Crystals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.4.3 Crystal-liquid Interfaces . . . . . . . . . . . . . . . . . . . . . 242.5 Time Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.5.1 Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.6 Temperature and Pressure Control . . . . . . . . . . . . . . . . . . . 282.6.1 Velocity Rescaling Thermostat . . . . . . . . . . . . . . . . . 292.6.2 Berendsen Barostat . . . . . . . . . . . . . . . . . . . . . . . 312.6.3 Parrinello-Rahman Barostat . . . . . . . . . . . . . . . . . . 323 Nucleation of NaCl from solution . . . . . . . . . . . . . . . . . . . . 343.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.2 Models and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.2.1 Simulation Details . . . . . . . . . . . . . . . . . . . . . . . . 373.2.2 Cluster Detection . . . . . . . . . . . . . . . . . . . . . . . . 393.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 473.3.1 Cluster Lifetimes . . . . . . . . . . . . . . . . . . . . . . . . . 503.3.2 Crystallinity and Nucleation . . . . . . . . . . . . . . . . . . 553.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 674 Crystallization of Lithium Halides from Molecular Simulation. . 694.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 694.2 Models and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 714.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 754.3.1 Lattice Energies . . . . . . . . . . . . . . . . . . . . . . . . . 754.3.2 Infinitely Periodic System Simulations . . . . . . . . . . . . . 824.3.3 Finite Size Cluster Simulations . . . . . . . . . . . . . . . . . 834.3.4 Temperature Effects . . . . . . . . . . . . . . . . . . . . . . . 894.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 905 Crystallization of Lithium Fluoride in Aqueous solution from Molec-ular Simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 915.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91vi5.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 925.3 Models and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 935.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 975.4.1 Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 985.4.2 Growth Barrier . . . . . . . . . . . . . . . . . . . . . . . . . . 995.4.3 Solution Structure and Dynamics . . . . . . . . . . . . . . . . 1015.4.4 Solubility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1065.5 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 1096 Dissolution of NaCl nanocrystals . . . . . . . . . . . . . . . . . . . . 1116.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1116.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1126.3 Models and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 1146.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 1186.4.1 Stages of Dissolution . . . . . . . . . . . . . . . . . . . . . . . 1186.4.2 Concentration Effects . . . . . . . . . . . . . . . . . . . . . . 1326.4.3 Temperature Dependence and Activation Energy . . . . . . . 1346.5 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 1367 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1387.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1387.2 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144AppendicesA Distribution of Various Properties for Nucleated and Failed Clus-ters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157viiList of Tables3.1 Joung-Cheatham Lennard-Jones parameters for NaCl.1 . . . . . . . . 373.2 Summary of the simulations performed. The number of water molecules(Water), the number of ions (Ions), the mole fraction xNaCl, the runlength (Length), and the number of clusters that achieved nucleation(Nuclei) are indicated. At xNaCl = 0.3, multiple nucleations occurredresulting in a polycrystal. The five simulations at xNaCl = 0.22 listedat the bottom of the table were used for the statistical analysis givenin the text. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394.1 Lattice energies for the JC model for the wurtzite (Ewurtz) and rocksalt (Erock) structures. The difference Erock- Ewurtzis reported as ∆E.The NaCl lattice energies are added for comparison. . . . . . . . . . . 804.2 Lattice energies (TF model) for the wurtzite (Ewurtz) and rock salt(Erock) structures. The difference Erock- Ewurtzis reported as ∆E. TheNaCl lattice energies are added for comparison. . . . . . . . . . . . . 804.3 Final structures obtained by melting and freezing periodic crystals.The final temperature is 300 K, and the periodic box contained 1000ions. The results are consistent with the lattice energy calculations(Tables 4.1 and 4.2). . . . . . . . . . . . . . . . . . . . . . . . . . . . 824.4 Final structures obtained by melting and freezing finite size clus-ters using the JC and TF models, the final temperature is 300 K. Awurtzite-like structure corresponds to a structure that is not perfectlycrystalline but shows hexagonal motifs similar to that shown in theleft panel of Figure 4.5. . . . . . . . . . . . . . . . . . . . . . . . . . . 844.5 LiY lattice energies at 0 K and 300 K. The value Ew was obtainedby simulating a box to initialized from the perfect wurtzite structure,while Ew† is the energy of the crystal obtained from freezing of themolten salt. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 895.1 Joung-Cheatham Lennard-Jones parameters for LiF1 adapted for theSPC/E water model.2 . . . . . . . . . . . . . . . . . . . . . . . . . . 93viii5.2 Summary of the simulations performed to study LiF nucleation andgrowth. The temperature at which a simulation was conducted isindicated by T, the length of the simulation is indicated by Length,the concentration in mole fraction is abbreviated as m.f., and the timeto nucleation (ttn) is the time at which the first nucleus was observed.LiF and H2O indicate the number of molecules of each type used in thesimulation. The pressure values for the crystal growth simulations arereported as a range for some simulations because, when substantialcrystal growth occurs, the pressure drifts to higher values. . . . . . . 945.3 Diffusion constants for Li+ (D+) and F− (D−) at different temperatures. 995.4 Growth rates at different temperatures . . . . . . . . . . . . . . . . . 1015.5 Average number of clusters (nc) per time frame at different temper-atures, for LiF and NaCl supersaturated solutions (see Chapter 3 fordetails). nc is normalized by the number of ion pairs in the simulationsto allow comparisons between simulations with different numbers ofions. The concentration in mole fraction is reported in the m.f. column.1066.1 Joung-Cheatnam NaCl parameters1 for the Lennard-Jones potential. 1146.2 Summary of crystal shape, the number of ion pairs, the number ofwater molecules, and the temperatures used in the simulations. . . . . 1166.3 Goodness of fit parameters R2 obtained for different rate laws. Thebest fits are indicated in bold. . . . . . . . . . . . . . . . . . . . . . . 130ixList of Figures1.1 Kinetic model of CNT. Droplets can grow only by single moleculeattachments and shrink by single molecule detachments. The forwardrate a(n) = A(n)β is proportional to the surface area of the nucleus,A(n), and the rate of arrival of single molecules on the droplet, β. Thebackward rate b(n) = A(n)α is proportional to the rate of detachmentper unit area, α, and the surface area. . . . . . . . . . . . . . . . . . 42.1 The LJ potential with parameters σ = 1 and = 1. The σ correspondto the x-intercept of the potential, while corresponds to the well depth. 152.2 Positions of the charges in the models employed in our simulations.The SPC/E water model on the left represents the atoms as pointcharges, with a O-H distance of ∼ 0.1 nm and a HOH angle of ∼109.47 deg. The alkali and halide ions are modeled as single pointcharges. The short-range interaction parameters are not shown inthe picture. For the water model, the masses of the oxygen andhydrogen atoms were employed for the negative and positive chargedsites, respectively. Similarly, the mass for the corresponding elementwas chosen for the alkali halide models. . . . . . . . . . . . . . . . . . 172.3 Depiction of PBC. The central system is repeated in all directions.A particle exits from the central box on the right side (black arrow)and reenters from the opposite side. . . . . . . . . . . . . . . . . . . . 192.4 The short-range and long-range contribution to the electrostatic po-tential for κ = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.1 Density of q8 values in a sample crystal (∼ 2000 ions)(blue) and asupersaturated solution (xNaCl = 0.20) (green). The two distribu-tions show excellent separation. The vertical orange line indicatesthe selected threshold (q8 > 0.35). . . . . . . . . . . . . . . . . . . . . 41x3.2 Schematic diagram illustrating the DBSCAN algorithm for η = 4.Top left panel: an ion is chosen at random and ions within arecounted; the value is 4 so the ion is marked as core (red). Its neigh-bours are labeled as border (pink). Top right panel: the procedure isrepeated for one of the border ions which becomes a core ion. Bottompanel: all ions in the cluster have been processed, and no more ionsare reachable from the core ions. The isolated ion is labeled as noise(blue). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.3 Schematic diagram illustrating the cluster time resolution algorithm.At t = 0 the clusters are labeled employing the DBSCAN algorithm.A cluster is compared with clusters in the following frame, and aconnection is made if the clusters are sufficiently similar, judged bythe Jaccard index as described in the text. The procedure is repeatedrecursively for the following frames. By performing this procedurefor all clusters, it is possible to obtain a connectivity graph (bottom)where each connected component represents the time evolution of asingle cluster. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.4 Size (top left) and crystallinity (top right) histograms for clusters atfirst detection. The bottom panel shows the joint distribution of sizeand crystallinity. Note that very few clusters begin with both highsize and high crystallinity, suggesting that this feature develops asclusters evolve in time. The red points indicate the thirteen clustersthat achieve nucleation, and we note that these show no obviouspreference for any region of the joint distribution. . . . . . . . . . . . 493.5 Examples of clusters of different sizes and crystallinities (given inparenthesis below each cluster). Na+ and Cl− ions are colored pur-ple and green, respectively. At inception larger clusters tend to beelongated and of lower crystallinity (top left), whereas larger highcrystallinity clusters (not observed at early stages) tend to be morecompact (top right). For small clusters structural features associatedwith low (bottom left) and high crystallinity (bottom right) are notso apparent. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.6 Histograms of “failed” cluster lifetimes (left) and the overall survivalfunction (right) for all clusters detected in the simulations. . . . . . . 523.7 Survival curves for clusters of size 10 grouped by high (> 0.40) andlow (≤ 0.40) crystallinity. . . . . . . . . . . . . . . . . . . . . . . . . . 543.8 Ratio of median lifetimes for high versus low crystallinity clusters ofdifferent size. The high crystallinity clusters always have the highermedian lifetime. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55xi3.9 Trajectory shown in (size, crystallinity) space for a nucleated (blue)and a failed (red) cluster as they evolve in time. In early stages bothclusters oscillate in size up to ∼ 30 ions, but the crystallinity of thecluster that eventually achieves nucleation reaches higher values thanthe failed case. This suggests that it is a combined effect of both sizeand crystallinity that promotes nucleation. . . . . . . . . . . . . . . . 563.10 Comparison of growth and crystallinity profiles for nucleated (col-ored) and failed (gray) clusters. After a short period, the nucleatedcluster manages to reach quite a high crystallinity and maintain itssize. In contrast, the failed cluster, while maintaining its size, expe-riences a steady decrease in crystallinity. . . . . . . . . . . . . . . . . 583.11 Size (blue) and crystallinity (orange) profiles for a nucleated cluster.The fluctuation in crystallinity is fairly high, and, as the crystal in-creases in size, the crystallinity reaches a plateau at ∼ 0.47. Thesnapshots represent the cluster at 0 ns (left), 40 ns (center) and 80ns (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593.12 Crystallinity distributions (blue histogram) of failed clusters of dif-ferent size (number of ions). Clusters that achieved nucleation areindicated by single orange lines. Note that clusters that achieve nu-cleation come preferentially from the upper part of the crystallinitydistribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.13 Probabilities of high (q¯8 > 0.40) and low (q¯8 < 0.40) crystallini-ties achieving nucleation. The error bars (vertical black lines) repre-sent one standard deviation. The lower panel shows the results forsmall clusters on an expanded scale. Note that for all cluster sizesP (N |high) is always larger than P (N |low), and the high/low ratio isshown as a gray line in the upper panel. . . . . . . . . . . . . . . . . 623.14 Joint distribution of crystallinity and cluster energy (kJ mol−1) forclusters of size 6 (left) and size 10 (right). Values corresponding toclusters that eventually nucleated are displayed in orange. The cor-relation coefficient is negative and very low at size 6, and moderatelylow at size 10, indicating that crystallinity carries different informa-tion than energetic stability. . . . . . . . . . . . . . . . . . . . . . . . 64xii3.15 Running coordination numbers for: high q8 Na+ with all Cl− (blue);high q8 Cl− with all Na+ (green); and all Na+ and Cl− (red). In thefirst coordination shell (∼ 0.35 nm) both high q8 ions are surroundedby more counterions (3-4) than average (1). Also, high q8 Na+ aresurrounded by more counterions than high q8 Cl−, suggesting thaton average they occupy positions deeper within the clusters. Thedifference in the first shell coordination number decreases at longertimes, as larger crystals develop in the simulation, and the internalpreference of Na+ becomes less noticeable. . . . . . . . . . . . . . . . 653.16 ∆q8 as a function of cluster size. The shaded area indicates onestandard deviation. The average value is slightly positive indicatingthat Na+ ions tend to occupy ordered environments than Cl− ions. . . 664.1 Unit cells used for the rock salt (right) and wurtzite (left) crystalstructures. The configurations are dependent on a single cell pa-rameter a, all other parameters being constrained. In the rock saltstructure, a is the length of any edge, while for the wurtzite structurea is the length of any of the two short edges. . . . . . . . . . . . . . . 734.2 Radial distribution functions for an infinite LiCl crystal for the TFmodel (left panel) and JC model (right panel). The initial rock saltconfiguration and the liquid phase are colored in blue and red, respec-tively, and do not show noticeable differences between models. Therdf corresponding to the crystallized structures are displayed in greenand correspond to wurtzite on the left and rock salt on the right.The inset plots show the rdf for ideal wurtzite on the left and idealrock salt on the right (black lines) superimposed onto the simulationresult (green). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.3 TF and JC potentials for the pairs Li+-Cl− (top panel), Li+-Li+ (mid-dle panel) and Cl−-Cl−(bottom panel). The TF model has a shortercontact radius for the Cl− ion, and is more slowly varying functionat short range. Notice also how at very short range (∼ 0.1 nm),the TF model becomes attractive (but this does not influence thesimulations). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784.4 Lattice energies as a function of the cell length a for the JC (leftcolumn) and TF (right column) potentials. The solid red lines indi-cate the energy of the wurtzite structure, while the dashed red linesindicate the rock salt structure. The energy minima are annotatedbelow each curve. Blue and green lines represent the contributions ofthe positive and negative ion sites, respectively. . . . . . . . . . . . . 81xiii4.5 Spontaneously crystallized LiCl/TF (left panel) and LiCl/JC (rightpanel) clusters of 216 ions. For the cluster on the left, the structurehas hexagonal motifs characteristic of the wurtzite structure, but ismostly hollow inside. On the right, the cluster crystallizes into arock salt structure, without substantial surface deformations. In thisfigure, Li+ is represented in purple and Cl− is in green. . . . . . . . . 854.6 LiI/JC rdfs at 300 K for the perfect wurtzite crystal (red), and forthe interior of a finite size crystal of 1000 ions (blue). The peaks arewell defined for the perfect wurtzite crystal, while for the finite sizecrystal there is wide broadening and merging of adjacent peaks. . . . 864.7 Example of surface reconstruction for a LiI crystal (Li+ and I− ionsare colored in purple and grey, respectively). On the left, surfacecharges are present on the wurtzite crystal face (0001). After simu-lating the crystal in vacuum at 300 K, the resulting surface (on theright) rearranges to reduce the dipole moment. The dipole momentper ion is reduced from ∼ 6.3 D for the structure on the left, to ∼ 0.7D for the structure on the right. . . . . . . . . . . . . . . . . . . . . . 874.8 Spontaneously crystallized LiBr/JC cluster (1000 ions) using the LiBr/JCmodel (Li+ in purple, and Br− in brown). One notices both hexago-nal motifs, characteristic of the wurtzite structure, as well as squaremotifs characteristic of the rock salt structure. . . . . . . . . . . . . 885.1 Clustering based on neighbor distance. The dashed lines representdistances lower than a fixed threshold, and the colors represent dif-ferent clusters. Free particles form a cluster of their own. . . . . . . . 975.2 Schematic diagram of the free energy landscape at different temper-atures for the crystal growth process. If the driving force (which isthe difference in free energy of initial and final stages) increases andthe barrier is weakly affected by the increase in temperature, the ac-tivation barrier (which is the height of the curve between inital andfinal stages) of the process decreases. . . . . . . . . . . . . . . . . . . 985.3 Linear regression to determine the activation energy for the diffusionprocess for Li+ ions. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1005.4 Arrhenius fit to the growth rates . . . . . . . . . . . . . . . . . . . . 1015.5 Fraction of ions in connectivity based clusters of a certain size atdifferent temperatures. . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.6 Li+F− radial distribution functions at different temperatures. Thecontact peak, characteristic of ion pairs increases with temperature,while the second peak (solvent-separated pairs) decreases. . . . . . . . 103xiv5.7 Diagram of surviving pairs. Initially, the {1, 2, 3, 4} cluster contains6 unique pairs, after the cluster splits, 2 of the pairs are still associated.1045.8 Surviving pairs at a time delay of t ns for a range of different tem-peratures. As the temperature increases, there are more associatedpairs that decay at a faster rate. . . . . . . . . . . . . . . . . . . . . . 1055.9 Solution concentration profiles obtained by dissolving a crystal at 300K (blue), by growing a crystal at 500 K (red, top line), and dissolvinga crystal at 500 K (red, bottom line). The estimated solubilities are∼ 0.025 mole fraction at 300 K and ∼ 0.01 at 500 K. . . . . . . . . . 1085.10 Possible free energy landscape for LiF nucleation from solution. Afirst barrier is required to go from the free ions (solution) to the clus-ter aggregates. The cluster aggregates are metastable, continuouslyform and disappear from solution, and are more easily found at hightemperature. The second barrier reflects the fact that most clustersdo not nucleate but dissolve back into solution, depending on clusterrelated properties such as size, crystallinity, and surface tension. Inabsence of nucleation rate measurements, it is not possible to estimatethe relative magnitude of the two barriers, and two possible optionsare displayed as solid and dashed lines. . . . . . . . . . . . . . . . . . 1106.1 A representation of the crystals simulated. The number of NaCl ionpairs is given below each illustration. . . . . . . . . . . . . . . . . . . 1156.2 An example of the order parameter applied to a surface ion. The cen-tral atom (yellow) is surrounded by 14 neighbors (red). The 0.6 nmcutoff applied to the central atom is highlighted with a gray trans-parent sphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1186.3 Dissolution profiles for the cubical and spherical crystals displayed asnumber of ions in the crystal vs time. The starting and final sectionsof the profiles are detailed in the zoomed-in plots. . . . . . . . . . . . 1206.4 Snapshots corresponding to different points in the dissolution profileof the cubic crystal. After 20 ns, only ions located at edges andcorners are removed. The edges and corners gradually get consumed(100 ns) and the crystal becomes roughly spherical after about 185ns. After that point the shape doesn’t change until the final stage ofdissolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121xv6.5 Comparison of the dissolution profiles for different shapes expressed interms of fraction of crystal dissolved fcry(t), the derivatives dfcry/dt(on the order of −10−3 ns−1) are depicted in the inset. The colorcoding is: dark blue for the sphere; green for the cube; red for thetablet; light blue for the rod. The rod-shaped and tablet-shapedcrystals show a higher overall dissolution rate. . . . . . . . . . . . . . 1226.6 Snapshots along the rod-shaped crystal trajectory. Ions are detachedmainly from the ends of the rod. . . . . . . . . . . . . . . . . . . . . . 1226.7 Snapshots along the tablet-shaped crystal trajectory. The top andbottom surfaces are never attacked by water molecules. . . . . . . . . 1236.8 Three quantities obtained for the spherically shaped crystal as func-tions of the number of neighbor in the first coordination shell (within0.37 nm), Left panel: The relative probability of detachment (within1.6 ns) of a suface ion, given the number of neighbors. Upper rightpanel: The total number of ions on the surface at t = 0. Lower rightpanel: The total number of ions detached over 1.6 ns. All quantitiesare averaged over 50 time slices betwen 20 and 40 ns of the dissolutionsimulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1256.9 Radial ion density profile starting from the center of a spherical crys-tal at 300 ns. The dashed horizontal line in the inset indicates thedensity at saturation. The calculated radius of the crystal is ∼ 2 nm,and a concentration gradient (that could be related to a diffusionlayer) is not observed. . . . . . . . . . . . . . . . . . . . . . . . . . . 1316.10 Snapshots of the last stage of dissolution of the spherical crystal.When the crystal reaches ∼64 ions, the structure begins to disinte-grate. In the final snapshot, water quickly penetrates the structureand breaks apart the remaining crystal nucleus. . . . . . . . . . . . . 1326.11 Concentration effects on rate. Curves (a), (b), and (c) are for NaClmole fractions (calculated assuming complete dissolution of the entirecrystal) of 0.0254, 0.0163, and 0.0127. Note that curve (a) showsa sharp decrease in rate at ∼120 ns, while there is no substantialdifference between the profiles (b) and (c). . . . . . . . . . . . . . . . 1336.12 Detachment and reattachment events for the spherical crystal. Dur-ing the course of dissolution, the detachment rate systematically slowswhile the attachment rate remains substantially constant. . . . . . . . 1346.13 Fits to the cube root law for the spherical crystal at temperatures of300, 320 and 340 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135xvi6.14 Fits to the Arrhenius equation of the rate constants obtained for thespherical crystal at 300, 320, and 340 K. Error bars represent onestandard deviation of uncertainty. . . . . . . . . . . . . . . . . . . . . 135A.1 Distributions of various properties for failed clusters (blue histogram)of size 10. Clusters that achieved nucleation are indicated by singleorange lines. Note that the property values of the clusters that achievenucleation are found around the mean of the distribution. . . . . . . . 159A.2 Distributions of various properties for failed clusters (blue histogram)of size 30. Clusters that achieved nucleation are indicated by singleorange lines. At this size, some preference for low radius of gyration,low surface area and higher sphericity can be observed. . . . . . . . . 160A.3 Two dimensional representation of the algorithm used to estimate thevolume. The cluster is encased in a grid, and the grid points that liewithin the cluster are represented in red. . . . . . . . . . . . . . . . . 161xviiList of SymbolskB Boltzmann constantK Kinetic energyV Potential energyTr(A) Trace of matrix Abcc Body centered cubicCNT Classical Nucleation Theoryfcc Face centered cubicJC Joung-CheathamLJ Lennard-JonesPBC Periodic boundary conditionsPME Particle-mesh Ewaldrdf radial distribution functionTF Tosi-FumixviiiTo my parentsxixChapter 1Introduction1.1 Phase Transitions and NucleationPhysical matter can exist in different phases, and phase transitions are theprocesses that describe how matter can transform from one phase to another, after achange in one or more thermodynamic parameters. In equilibrium thermodynamics,the stable phase of a system at constant temperature T and pressure p is the onewith the lowest chemical potential µ. In general, as the temperature and pressure ofthe system are changed, the chemical potential of each phase will vary by differentamounts, defining the regions of stability in the phase diagram.While thermodynamics describes systems at equilibrium, it gives no infor-mation on the time evolution, nor on the atomistic details of phase transitions. Thedynamical aspects of phase transitions are, however, extremely important as thetime scales involved can be very large, and their understanding is crucial in manyapplications.A system that is stable to small fluctuations in thermodynamic parameters,but can transition to a more stable state given a large enough perturbation, is termed1metastable. Metastable states are due to the fact that the system is in a local freeenergy minimum, and to reach the global minimum it is necessary to overcome afree energy barrier. An extreme example is the graphite-diamond transition. At298 K and 1 bar, graphite is the stable phase of carbon, with a free energy thatis ∼ 2.90 kJ mol−1 lower than that of diamond.3 However, the time to convertdiamond to graphite at room temperature and pressure is extremely long, makingdiamond remarkably metastable.The dynamics of phase transitions are usually explained through a processtermed nucleation. In order to form a new phase, particles in the metastable phasespontaneously arrange to form the initial seed of the new phase, that will grow asthe transition progresses. The formation of the initial seed can happen directly inthe metastable phase (homogeneous nucleation), or on a surface, such as the wall ofthe container, or can be initiated by impurities (heterogeneous nucleation).Nucleation is at the core of many physical and biological processes. For ex-ample, some organisms in cold environments have developed strategies to preventfreezing of bodily fluids through the use of antifreeze proteins that inhibit ice nucle-ation.4,5 Another example are high clouds (clouds that occur between heights of 5 to13 km above the earth’s surface), where water droplets exist in a supercooled state,and freezing can be triggered6 at about -15 ◦C, catalyzed by other particles, suchas clay particles, carried into the atmosphere by the wind.7,8 Once water dropletsfreeze they can grow and coalesce into snowflakes or hail that precipitate to theground.9In the remainder of this Chapter, we will describe the classical kinetic theoryof nucleation, the experimental and theoretical attempts to validate it, and howcomputer simulations have proven to be a valuable tool in the study of nucleation.21.2 Classical Nucleation TheoryThe kinetic treatment of nucleation dates back to 1926, when Volmer andWeber10 postulated that the nucleation rate J (number of nuclei formed per unittime and volume) depends on the negative exponential of the free energy requiredto form a new nucleus. Later developments and refinements11–13 gave rise to thestandard treatment for nucleation, often called classical nucleation theory (CNT).In the original treatment of CNT, applied to a homogeneous gas-liquid tran-sition, the metastable vapor is characterized to be a mixture of single molecules andliquid droplets of various sizes. The system is assumed to be in an equilibrium14characterized by a droplet distribution N(n), that is proportional to the negativeexponential of the free energy ∆G needed to form a droplet of size nN(n) = N0 exp(−∆G(n)kBT), (1.1)where N0 is the number density of the nonassociated molecules, and kB is the Boltz-mann constant. The equilibrium distribution N(n) describes the state of the vaporbefore the phase transition, and to characterize the kinetics of the process it is neces-sary to consider the time evolution of the actual, instantaneous, droplet distributionf(n, t).To recover this dynamical quantity, CNT assumes that the interchange ofmolecules between the liquid droplets and the vapor happens by single moleculeattachment or detachment, each collision between a molecule and a droplet resultingin an attachment, and every attachment or detachment event is uncorrelated withthe previous ones (Figure 1.1). If this is true, the droplet formation process describedabove can reach a steady state and the nucleation rate, J , can be found solely from3Figure 1.1: Kinetic model of CNT. Droplets can grow only by single molecule attach-ments and shrink by single molecule detachments. The forward rate a(n) = A(n)βis proportional to the surface area of the nucleus, A(n), and the rate of arrival ofsingle molecules on the droplet, β. The backward rate b(n) = A(n)α is proportionalto the rate of detachment per unit area, α, and the surface area.time independent quantities, such that9J =1∞∑n=11βA(n)N(n), (1.2)where A(n) is the surface area of a droplet of size n, and β is the rate of arrival ofsingle molecules at the droplet per unit area.Other assumptions are necessary to simplify Equation (1.2). According toCNT, droplets are postulated as large, homogeneous, spherical, and incompressibleobjects. The inside of the droplet is made of bulk liquid, and the surface free energyof the cluster can be described by the product of a size-independent surface tension,γ, and the cluster surface area A(n). This set of assumptions is usually referred toas the capillary approximation.If the capillary approximation applies, the free energy of formation for adroplet of size n will depend on the free energy of formation of the new bulk phase,as well as on the free energy required to form a new surface,∆G(n) = −n∆µ+ γA(n), (1.3)4where ∆µ is the difference between the chemical potential of the liquid and thevapor, and γ is the surface tension of the liquid. Even though CNT was originallyformulated for vapor to liquid transitions, the theory can be applied to other tran-sitions, such as crystal nucleation, by using the chemical potentials and interfacialfree energies of the phases of interest.When the cluster is very small, the free energy required to create the newsurface will be higher than the free energy gained by the formation of the morestable phase, and as the size of the droplet grows, ∆G will go through a maximum,then decrease to the point of being negative. The size n∗ at which the maximum isfound, also called critical size, isn∗ =32pi3v′2γ3∆µ3, (1.4)and the maximum ∆G(n∗), also termed the nucleation barrier, is∆G(n∗) =16pi3v′2γ3∆µ2, (1.5)where v′ is the volume per particle in the condensed phase.Given the above assumptions, the sum in Equation (1.2) can be expressed asa solvable integral9 and the rate J is found to depend on ∆G(n∗) and a prefactorJ0J = J0 exp(−∆G(n∗)kBT). (1.6)5The factor J0 can be expressed asJ0 = ZβA(n∗)N0, (1.7)Z =√−∆G′′(n∗)2pikT, (1.8)where Z, also termed Zeldovich factor,9 physically represents the fact that there isa probability that critical nuclei can revert to vapor, while the term A(n)β is thefrequency of arrival of single particles on a nucleus of size n. The term ∆G′′(n∗) isthe second derivative with respect to n of ∆G, evaluated at n∗. The main parametersaffecting the nucleation rate, J , are the temperature, the difference in the chemicalpotential of the two phases, and the surface tension.The temperature dependence is expressed through the prefactor J0, the ex-ponential dependency on 1/kBT , as well as on the degree of undercooling that willaffect the chemical potential ∆µ. While lower temperatures tend to increase therate because of the higher degree of undercooling, the rate of arrival β will substan-tially decrease because of decreased mobility. The chemical potential dependencyis expressed through ∆G(n∗). As ∆µ increases, the barrier ∆G(n∗) will decrease,causing an increase in nucleation rate J . A higher surface tension, γ, leads to anincrease in nucleation barrier and a decrease in nucleation rate.Early experiments using expansion chambers,15–18 showed that classical nu-cleation theory predicts the correct dependence of nucleation rate on supersatu-ration, but fails to reproduce the temperature dependence, a fact that has beenattributed to an inadequate description of the prefactor J0.19Auer and Frenkel20 were able to characterize the crystal nucleation rate ofthe hard-sphere model using molecular simulation, and compared their results with6experiments on colloids.21,22 They showed that, albeit CNT adequately describesthe functional form of ∆G(n∗), the nucleation rates are more strongly dependent ontemperature in models than in experiments, suggesting that the problem is, again,in the description of the temperature dependent prefactor. Further experiments oncolloidal systems also showed issues in considering the surface tension as concentra-tion independent.23Very interesting developments on the microscopic mechanism of nucleationwere evidenced by experiments on CaCO3 crystal nucleation.24–26 In this system,crystal nucleation proceeds quite differently from the molecule-by-molecule growthprocess described by CNT. It was observed24 that in the saturated solution, amor-phous calcium carbonate clusters are first formed, these subsequently undergo tran-sitions to vaterite and only then, convert to the final structure calcite. This indicatesa nucleation process that happens in stages, similar to a mechanism that was firstpostulated by Ostwald27 and termed rule of stages.In contrast with what was found for CaCO3, the microscopic process de-scribed in CNT was shown to apply to the description of certain systems. For exam-ple, Sleutel et al., performed direct observation (through atomic force microscopy)of nucleation in two dimensions of the protein glucose isomerase and showed thatsubcritical nuclei have the same structure as the bulk crystalline phase.28 Measure-ments of the nucleation rate of NaCl crystals also seem to be in good agreement withestimates from CNT.29 Additionally, calculations of the critical size for glycine andNaCl aqueous solutions were found to be consistent with estimates of the averagecluster size.30Attempts to describe the first stages of nucleation have proven hard to per-form experimentally, because of the intrinsic small size and the short lifetimes of7most nuclei.31 Additionally, estimates of the homogeneous nucleation rates are sub-ject to significant uncertainties, and often complicated by unwanted heterogeneousnucleation.32Computer simulations are a great way to avoid many of the shortcomingsthat affect experiments. As discussed above, the foundation of CNT relies uponassumptions about the microscopic nature of the metastable and nascent phases,and through molecular simulation it is possible to observe nucleation in its veryearly stages for a variety of model systems and physical conditions.1.3 Computer SimulationsComputer simulations have long been used to study phase transitions of sim-ple atomic systems,33–36 and early studies of nucleation were performed on Lennard-Jones (LJ) liquids.37–40 It was found that, despite the fact that the stable solidstructure for the LJ system is face centered cubic (fcc), the LJ liquid nucleates atmoderate supercooling with a body centered cubic (bcc) structure and, starting fromthe interior, gradually converts to the fcc structure, suggesting that the nucleationprocess, even for this very simple system, happens in stages.39Other early work includes simulations of hard spheres that were used to com-pare crystal nucleation rates with experimental data for colloids.20,41 More recently,thanks to a dramatic increase in parallel computational power, simulations are beingused to study homogeneous and heterogeneous nucleation molecular systems,42–45including aqueous solutions.46–50One of the most studied processes is perhaps the nucleation of ice in liquidwater.32 Homogeneous nucleation measurements for this system vary by orders of8magnitude between theory and experiments and also between different studies.32Attempts to study homogeneous nucleation through unbiased simulations have alsoproven to be problematic32 given the very long time scales involved in ice nucle-ation, but has been observed when an electric field is applied to the system.51–53Coarse grained simulations have suggested that ice nucleation involves both cubicand hexagonal variants of ice,42 also observed in unbiased simulations of heteroge-neous ice nucleation,54–56 and in experiments.57A considerable body of work, including some of the projects in this thesis,has been dedicated to the NaCl/water system. Early work on NaCl crystal nucle-ation46,47 used small systems to study the early stages, and showed that the Na+ iontend to be positioned near the center of the cluster aggregates present in NaCl so-lutions.47 The mechanism of cluster formation was further explored by Hassan,58,59showing that the supersaturated NaCl solution is characterized by a combinationfree ions and ionic clusters with relatively long lifetimes. Later work on larger sys-tems48–50 gave evidence of a two-stage nucleation process.49,50 The formation of theNaCl nucleus was characterized by a density fluctuation followed by an orderingtransition, where the cluster rearranges to form a more regular crystal structure,gradually dehydrating as the crystal grows.In order to accurately describe highly concentrated alkali halide solutions,Joung and Cheatham1 developed models tuned to reproduce solid properties, such aslattice energies and lattice constants, as well as solution properties (such as the freeenergy of solvation). Joung-Cheatham (JC) parameter sets have been extensivelyused to investigate solubility and crystal nucleation from solution.60–66The main focus of this thesis is the study of crystallization of alkali halidesfrom solution. In addition to using more recent models,1 we also employ advanced9data analysis techniques to observe and study the morphology and dynamical be-havior of ionic clusters in solution.In Chapter 2, we describe the molecular dynamics techniques and algorithmsused for the simulations performed in the subsequent chapters.In Chapter 3 we employ large-scale molecular simulations to investigate thenucleation of NaCl in water using the JC parameter set, and develop computationalmethodology to detect and follow in time a large number of pre- and post-nucleationclusters. Statistical analysis shows that a very large number of clusters are formed,and that the cluster lifetime and nucleation probability is strongly affected by thegeometrical arrangement of the ions contained in the cluster.In Chapter 4, we evaluate pair potentials for pure lithium halides by calcu-lating lattice energies, and by melting and freezing finite size clusters to assess theirpossible role in nucleation studies. The results indicate that for most lithium halidemodels, the stable structure is not rock salt but wurtzite, and that, in certain cases,the finite size structure can differ from that of the bulk.In Chapter 5, we examine the crystallization of LiF from solution. We per-form molecular simulations of supersaturated LiF solutions at temperatures rangingfrom 300 to 500 K and, while nucleation was observed only at a temperature of 500K, we were able to measure the growth rate of an already nucleated crystal at lowertemperatures. We found that the growth process is activated and that the barrierto growth is much larger than that of diffusion. We suggest that the barriers togrowth and nucleation are likely related to ion dehydration.In Chapter 6, we focus on the NaCl dissolution process. The study of disso-lution provides insights into the dynamics of ion attachment and detachment, and10into the stability limit of small ionic clusters in dilute solutions. By simulating thedissolution of NaCl nanocrystals of different sizes, shapes, and at different temper-atures, we find that the dissolution process is activated, and that the activationenergy is affected by the position of ions on the crystal surfaces.A global summary, further remarks, and future directions are presented inChapter 7.11Chapter 2Models and Methods2.1 OverviewMolecular dynamics is a simulation method widely used to study model sys-tems at the microscopic level. To perform a simulation, it is necessary to specifya system in terms of its initial state and interaction forces that apply. Numericalintegration of the equations of motion is then performed to obtain a trajectory. Inthis chapter we explain the main algorithms used to perform molecular dynamicssimulations, and how to include other effects such as periodic boundary conditions,constraints, temperature, and pressure controls.2.2 Molecular DynamicsMolecular dynamics is a computational technique used to simulate the evo-lution of a molecular system given the initial conditions and a description of theinteratomic forces. The result of a molecular simulation is a time trajectory that12can be used to calculate the thermodynamic, structural and dynamic properties ofthe system under study.Atoms in molecular systems are subject to intra and intermolecular inter-actions that originate from chemical bonds, electrostatic charges, short-range re-pulsions, and dispersion forces. In molecular simulations, these interactions areusually approximated employing simple potentials that can be used in the frame-work of classical mechanics. Even though molecular interactions are best describedby quantum mechanics, this approximation has proven quite successful in describingmany systems of physical and biological interest, especially since molecular dynamicssimulations are typically concerned with nuclear motion.The time evolution of a system, where the initial interaction site positionsri(t = 0), velocities r˙i(t = 0), and forces fi are specified, can be described byNewton’s equations of motionfi = mir¨i, (2.1)where mi is the mass associated with the interaction site, and r¨i is the acceleration.To obtain a trajectory ri(t), the equations of motion are integrated numerically.In this chapter we discuss the functional forms of the molecular interactionsused in the following chapters, and some of the techniques used to initialize theatomic positions in our systems of interest. We also describe how numerical integra-tion is performed, and how to impose further conditions such as periodic boundaries,bond constraints, temperature and pressure controls.132.3 Interaction PotentialsThe first step in performing molecular dynamics is the definition of the mod-els that represent the interactions between the atoms and molecules of interest. Thesimulations performed in this thesis involve water and alkali halide ions which aresubject to van der Waals and electrostatic forces. The chemical bonds between oxy-gen and hydrogen are not explicitly modeled, rather the bond distances and anglesare kept fixed by using constraints (See Section 2.5.1 for details.)2.3.1 Van der Waals ForcesVan der Waals forces are commonly expressed through the use of Lennard-Jones (LJ) type potentials. In this model the interaction between a pair of atoms ischaracterized by a strong repulsive core and a soft attractive term. The functionalform of the LJ potential requires the specification of two parameters that dependon the pair of atoms considered, σij and ij,uLJ(rij) = 4ij[(σijrij)12−(σijrij)6], (2.2)where rij = |rij| = |rj − ri| is the interatomic distance. The σij and ij parametersare usually termed the size and energy parameters, respectively, and they correspondto the x-intercept and well depth of the potential as shown in Figure 2.1.Quite often, σij and ij are obtained by combining parameters that dependon the individual particles σi, σj, i, j, effectively reducing the number of parametersneeded from 2(N2)to 2(2N), where N is the number of atom (interaction site) types.The combination rules depend on the specific parameter set and, in this14Figure 2.1: The LJ potential with parameters σ = 1 and = 1. The σ correspondto the x-intercept of the potential, while corresponds to the well depth.thesis, the Lorentz-Berthelot rules are used. According to these rules, σi and σj arecombined using the arithmetic mean, while the geometric mean is used to combinei and jσij =σi + σj2, (2.3)ij =√ij. (2.4)The LJ and other short-range potentials aimed at describing van der Waalsforces decay rapidly with distance, and are typically calculated by truncating theinteractions after a certain distance, also termed the cutoff radius.152.3.2 Electrostatic ForcesThe electrostatic interactions are modeled using Coulombic potentials, thatdepends on the interatomic distance, rij, and the partial charges, qi and qj,uC(rij) =14piε0qiqjrij, (2.5)where ε0 is the dielectric constant in vacuum. The electrostatic potential doesnot decay rapidly with distance, and cannot be truncated by imposing a cutoff.Strategies have been developed to deal with electrostatic interactions, as explainedin Section 2.3.4 below.2.3.3 Water and Ion ModelsThe modelling of atoms and molecules requires the parameterization of thepotential functions. Parameters can be obtained by fitting molecular propertiesobtained from ab initio calculations,67–69 or by fitting the parameters to reproduceknown experimental properties2,61 (e.g. hydrogen bond lengths, densities, localstructures, hydration energies).To model water, a simple approach is to use single point charges on theoxygen and hydrogen sites, and to represent the molecule as a sphere centered onthe oxygen atom. In the SPC/E (extended simple point charge) water model2 usedin our simulations, partial charges are positioned on the oxygen and hydrogen atomsaccording to the geometry shown in Figure 2.2. The masses for the three interactionsites correspond to the masses of the oxygen and hydrogen atoms. The hydrogenatoms lack LJ interaction terms (σH = 0), while the oxygen σO = 0.3166 nm. The16Figure 2.2: Positions of the charges in the models employed in our simulations. TheSPC/E water model on the left represents the atoms as point charges, with a O-Hdistance of ∼ 0.1 nm and a HOH angle of ∼ 109.47 deg. The alkali and halide ionsare modeled as single point charges. The short-range interaction parameters are notshown in the picture. For the water model, the masses of the oxygen and hydro-gen atoms were employed for the negative and positive charged sites, respectively.Similarly, the mass for the corresponding element was chosen for the alkali halidemodels.parameters are obtained to reproduce experimental properties such as the radialdistribution function, diffusion constant, and density at 300K.2The alkali halide ions used in our simulations are modeled as point chargesplus LJ interaction terms.2.3.4 Infinite SystemsCurrently, typical system sizes used in molecular simulation are on the or-der of ∼ 105 − 107 atoms with a record holding simulation of ∼ 1012 atoms.70As macroscopic systems are much larger (a mole of substance is ∼ 1023 particles)than molecular simulations can handle, special techniques are usually employed asdescribed below.A large system can be approximated as a periodic repetition of one of itsparts. Using this construction, the simulation can be performed on a small system17that interacts with identical copies that extend in all directions, a constructiontermed periodic boundary conditions (PBC).A depiction of a two dimensional cell under PBC is shown in Figure 2.3.The central cell is surrounded by identical copies and, as the system evolves intime, particles that exit from one side of the box are reinserted from the oppositeside. As there are an infinite number of copies of the central cell, to calculate thedistance between particles it is necessary to specify a particular periodic image. Inmolecular simulation, the image is chosen so that the distance between the atomsin consideration is minimum. This choice is termed the minimum image conventionand is implemented by calculating each component of the displacement vector rij ={xij, yij, zij} as followsxij =xj − xi − Lx if xj − xi > Lx/2,xj − xi + Lx if xj − xi < −Lx/2,xj − xi otherwise,(2.6)where Lx is the box length in the x direction. The calculation is performed similarlyfor the yij and zij components.The treatment of interactions is simple for the short-range (van der Waals)terms. As already mentioned, since they become negligible at a relatively shortdistance (e.g. ∼ 1 nm), they can be calculated by truncating all interactions at adistance larger than a given cutoff radius.For the long-range electrostatic interactions, the expression of the potential18Figure 2.3: Depiction of PBC. The central system is repeated in all directions. Aparticle exits from the central box on the right side (black arrow) and reenters fromthe opposite side.19energy felt by a charge in the central cell can be written asui =∞∑n=−∞∑jqiqj|rij + n| i 6= j, when n = {0}, (2.7)where n = {na, nb, nc} is the displacement of the periodic cell, and rij is the dis-placement vector between the ith and jth atoms in the central cell. Notice that theterm |rij + n| corresponds to the position of the charge in the periodic cell definedby n. The infinite series expressed in Equation (2.7) is conditionally convergent (i.e.its result depends on the way the terms are summed), and special methods havebeen developed for its evaluation.71One of the methods used to calculate electrostatic energies is the Ewaldsummation,72 that was developed to calculate Madelung constants for ionic crystals,and can be used to calculate the electrostatic interactions in a simulation subjectedto PBC. It was later shown71 that the evaluation of the sum of Equation (2.7) givesdifferent results depending on the medium surrounding the infinite system and, inthe case of the original Ewald summation,72 the medium is a perfect conductor (itsrelative dielectric constant r =∞.)In the Ewald summation method, the contribution of the Coulombic potentialis split into short-range and long-range parts by adding and subtracting an erf(κr)/rterm1r=1 + erf(κr)− erf(κr)r=erfc(κr)r+erf(κr)r, (2.8)erf(r) =2√pi∫ ∞re−z2dz, (2.9)erfc(r) = 1− erf(r). (2.10)20Figure 2.4: The short-range and long-range contribution to the electrostatic poten-tial for κ = 1.The first term on the right of Equation (2.8) decays rapidly with r (Figure 2.4)and, similar to the other short-range interactions, its contribution to the sum canbe calculated directly by using a cutoff, while the second term can be Fourier trans-formed and converges rapidly in reciprocal space.73 The parameter κ can be usedto modulate the convergence of the two terms as explained below.The expression of the Ewald summation for the potential energy felt by ancharged site i in a cell of volume V can be expressed as follows73ui = uri + uki + uselfi , (2.11)uri =∑n∑jqiqjerfc(κ|rij + n|)|rij + n| i 6= j when n = {0}, (2.12)uki =1piV∑k 6=0∑jqiqj(4pi2|k|2 ) exp(−|k|24κ2) cos(k · rij), (2.13)uselfi = −κpi1/2q2i . (2.14)The first term in Equation (2.11) uri is simply the direct sum, similar to Equation21(2.7), of the rapidly decaying part, and n is truncated to include charges up to acertain cutoff (often the truncation will only include the central cell.) The secondterm uki is the sum of the slowly varying part of the potential over the reciprocallattice vectors k, and uselfi is a correction needed to remove the interaction of aparticle with itself, introduced by the term uki . The parameter κ can be tuned tomodulate the convergence rate of either of the sums and is often set to 5/L, whereL is the box length. Given the value κ = 5/L, to achieve convergence in reciprocalspace, ∼ 200 reciprocal vectors k are required.73To retrieve the total electrostatic potential of the cell, it is sufficient to sumover all particles i in the following wayU =12∑i(uri + uki ) +∑iuselfi . (2.15)Additionally, if the system is assumed to be surrounded by a medium with a certaindielectric constant r, a correction that depends on the dipole of the central cellmust be added to the totalUdip =2pi(1 + 2r)V∣∣∣∣∑iqiri∣∣∣∣2. (2.16)Notice that, in the case of ionic systems (and in the original Ewald formulation),this term is neglected, corresponding to a system surrounded by a perfect conductor(r =∞).In terms of computational time, the Ewald summation method (when a con-stant cutoff on the short-range part is applied) scales as O(N2), where N is thenumber of charges. Most simulation programs, however, implement other variantssuch as particle-mesh Ewald74 (PME) that evaluates the reciprocal term (i.e., the22Fourier transform of the charge density) on a grid using the fast Fourier transformalgorithm.75 This reduces the time required to evaluate the electrostatic potentialto O(N log(N)).2.4 System InitializationGenerating the initial configuration is one of the most important aspect ofmolecular simulation, where the computational scientist sets up “experiments” tar-geted towards the validation of a specific hypothesis. In this section we explainsome of the techniques used to initialize atomic positions, as velocities are usuallyrandomly sampled from a Maxwell-Boltzmann distribution.2.4.1 Liquid SystemsLiquid systems can be initialized by placing the molecules on an evenly spacedgrid. This method can be effective when the objects to be placed are of similar sizeso that they are able to fit onto a grid without overlapping with each other.Another option is to insert particles at random and to perform a neighbordistance test to avoid particle collision. However, as the box becomes filled, randominsertions will likely cause collisions and, in high density systems, the many attemptsrequired can make this method impractically slow. Despite this drawback, thismethod is very general as it doesn’t require specification of a particular grid size,and for this reason it was used to initialize the liquid phase simulations in this thesis.232.4.2 CrystalsCrystals can be generated given information about the unit cell. Typically,the parameters required to specify the initial structures are the unit cell lengths a, band c, the unit cell angles, α, β and γ, the fractional coordinates of the atoms, andthe space group.This information can be usually found electronically as Crystallographic In-formation Files,76 that contain the data required to generate the unit cell. The unitcell can then be repeated and translated along the three crystal axes to yield arbi-trarily large crystals. Spherical crystals can be obtained by creating a large crystaland discarding all molecules or ions beyond a certain cutoff radius. Crystals canalso be extracted from molecular simulation after nucleating and growing a crystal,however, as such crystals are often imperfect, care must be taken to include thesame number of positive and negative ions necessary to maintain electroneutrality.2.4.3 Crystal-liquid InterfacesCrystal-liquid interfaces can be generated by placing a crystal in a simulationbox and then positioning liquid particles in the empty spaces using the techniquesdescribed above. Notice that, when using the crystals generated using the procedurementioned above, the results may vary greatly depending on the crystal face exposedand/or the presence of defects.242.5 Time EvolutionThe equations of motion described by Equation (2.1) can be integrated nu-merically to yield a trajectory. The method used in the simulations performed inthis thesis is the leap-frog method, termed such because the velocities are updatedat a different time step than the positionsri(t+ ∆t) = ri(t) + ∆tvi(t+ ∆t/2), (2.17)vi(t+ ∆t/2) = vi(t−∆t/2) + ∆tai(t), (2.18)ai(t) =fi(t)mi, (2.19)fi(t) = −∇ui(ri). (2.20)This method is particularly effective in molecular dynamics because it is time re-versible (integrating forward n steps and backward n steps brings the system to thesame initial position) and is able to keep the total energy constant to a good approx-imation. Care should be taken in choosing a time step ∆t that is small comparedto the motion in the system (a typical value is ∼ 1−2 femtoseconds).One drawback of the leapfrog integrator is that it doesn’t calculate velocitiesat integer time steps, and these need to be approximated from the velocities at halftime steps as followsvi(t) =vi(t+ ∆t/2) + vi(t−∆t/2)2. (2.21)252.5.1 ConstraintsThe vibration of chemical bonds and bond angles can be approximated inmolecular simulation using a harmonic potential. However, as the vibrational fre-quency is often high compared to the intermolecular motion, in certain models suchas the SPC/E water model, bonds and angles are kept at a fixed value (the moleculeis modeled as a rigid body). In simulation programs, rigid bonds and angles are usu-ally kept fixed using a constraining algorithm.A constraint is a function that introduces a dependency between coordinates,for example, a constraint that imposes a distance d between two sites l and m canbe expressed as followsg(r, t) = |rl − rm| − d = 0 (2.22)r = {r1 . . . rN}. (2.23)If the constraints are holonomic (i.e., they depend only on the coordinates r andtime t), the equations of motion for the constrained system can be obtained fromthe LagrangianL = K− V+ λ · g, (2.24)K =12Mv · v, (2.25)V =∑i∑j>iu(rij), (2.26)where λ = {λ0, ..., λK} are the K Lagrange multipliers associated with the K con-straints g = {g0, ..., gK}, K and V are, respectively, the kinetic and potential energyof the system. M is a diagonal (3N, 3N) matrix containing the mass of each particle26repeated three times (one for each of the x, y, z coordinates)M = diag(m1,m1,m1,m2,m2,m2, . . . ,mN ,mN ,mN). (2.27)The resulting equations of motion can be expressed in matrix notation77r¨ = (I−TB)M−1f −TB˙r˙, (2.28)T = M−1BT (BM−1BT )−1, (2.29)where B is the (K, 3N) matrix containing the partial derivatives with respect to theparticle coordinates of the constraintsBki =∂gk(r, t)∂xi. (2.30)The matrix I − TB projects the acceleration vector M−1f , and the last term ofEquation (2.28) represents centripetal forces due to rotation.In our simulations we adopted the LINCS algorithm77, which is a methodto efficiently integrate the constrained equations of motion discussed above. Theintegration of Equation (2.28) using the leap-frog algorithm can be performed byfirst conducting an unconstrained step to obtain the new coordinates and velocities(termed ru and vu), and then by projecting them using the matrix I − TB. Incomputer implementations, an additional term T(Br−d)/∆t is necessary to preventaccumulation of numerical errors77v(t+ ∆t/2) = (I−TB)vu(t+ ∆t/2)− T(Br(t)− d)∆t, (2.31)r(t+ ∆t) = (I−TB)ru(t+ ∆t) + Td. (2.32)27The slow step in this algorithm is the calculation of the matrix T, as itcontains the inversion of (BM−1BT )−1. In the LINCS method, the inversion isperformed efficiently by rearranging the matrix BM−1BT into a I − A form andapproximating its inverse using the series expansion(I−A)−1 = I + A + A2 + ... (2.33)The matrix BM−1BT can be transformed into the I−A form by using the diagonalmatrix S, defined as the inverse square root of the diagonal elements of BM−1BTas follows(BM−1BT )−1 = S(SBM−1BTS)−1S = S(I−A)−1S. (2.34)The advantages of the series expansion in Equation (2.33) are a high storage ef-ficiency (since A is usually sparse) and a faster computation time, compared tostandard matrix inversion. Nevertheless, Equation (2.33) is only valid when the ab-solute value of all the eigenvalues is less than 1, a condition that is not verified whenthere is a high degree of bond connectivity. The structure of matrix A also allowsfor easy parallelization, making the algorithm scalable to very large simulations.2.6 Temperature and Pressure ControlA system described by the equations of motion presented in Equation (2.1),will conserve energy and, if we enclose the system in a container or apply periodicboundary conditions, will conserve its volume (or density). Many systems of interest,like a laboratory environment, are best described by the conditions of constanttemperature and pressure. In the subsections below we will describe some of the28algorithms used to keep temperature and pressure to a fixed value in moleculardynamics simulations.2.6.1 Velocity Rescaling ThermostatThe instantaneous temperature in a molecular simulation can be calculatedfrom the kinetic energy of the system KT =2KkBNf, (2.35)where Nf is the number of degrees of freedom which, for a system subject to con-straints, is equal to the number of coordinates minus the number of constraints.A very simple approach, called the Berendsen thermostat,78 consists of fixingthe temperature by rescaling the velocities so that the kinetic energy assumes atarget value K′ consistent with the target temperature T ′. The scaling factor α,that is applied to the velocities is defined asα =√K′K, (2.36)K′ =12kBNfT′. (2.37)The factor α is usually applied with a predetermined frequency. The procedureoutlined above will keep the kinetic energy constant but will not reproduce fluctua-tions characteristic of the canonical ensemble, it was, however, used in some of thecalculations in this thesis to relax the systems to our target temperatures beforerunning the simulations for data collection.In order to retrieve the correct sampling for the canonincal ensemble, a dif-29ferent method, named the velocity rescaling thermostat,79 was used. In this method,the rescaling factor can be chosen to obtain a target kinetic energy that is sampledfrom the canonical ensemble distributionP (Kt) ∝ KNf2−1t e− KtkBT . (2.38)While, in principle, the kinetic energy could be rescaled to match a value sampledfrom this distribution at each time step, in practice this causes abrupt changes invelocity. In the velocity rescaling thermostat,79 this value is computed by evolvingthe kinetic energy according to an auxiliary stochastic process. The kinetic energyis rescaled to the target value gradually by updating the previous kinetic energyvalue according to the equationK(t+ ∆t) = K(t) + ∆K (2.39)∆K = (K′ −K(t))∆tτ+ 2√K′K(t)NfdW√τ, (2.40)where dW is a Wiener noise (a random displacement) and τ is a parameter, withthe dimensions of time, that can regulate how fast the system equilibrates to thetarget kinetic energy.The velocity rescaling thermostat was used to control the temperature in allof the simulations performed in this thesis.302.6.2 Berendsen BarostatFor a system consisting of N particles, volume V , and at an instantaneoustemperature T , the microscopic pressure can be calculated asp = kBTNV+13V∑ifi · ri, (2.41)where fi is the force acting on the ith particle at position ri. Control of the pressurecan be applied, in a manner similar to the velocity-rescaling technique, by rescalingthe particle positions and volume of the system by a factor αV˙ = 3αV, (2.42)x˙ = v + αx. (2.43)In the case of the Berendsen barostat,78 in order to bring the system pressure to avalue p′, the coefficient can be calculated asα = cβ(p′ − p)/(3τp), (2.44)where cβ is the isothermal compressibility of the system and the τp is a parameterwith the units of time. As cβ and τp appear as a ratio, only their relative magnitudeis important, and they affect how fast the system will equilibrate to the targetpressure.The Berendsen barostat does not produce the correct energy fluctuationsfor the isobaric-isoenthalpic ensemble, but is often used for its simplicity and toquickly bring the system to a target pressure. In this thesis, this barostat was used31when fluctuations generated by other methods would make the system unstable(more details are provided in Chapter 4), and was also used in other instances, aftertesting that it did not produce different results than more advanced methods suchas the Parrinello-Rahman barostat described below.2.6.3 Parrinello-Rahman BarostatThe Parrinello-Rahman barostat80 is another method to control pressure ina molecular dynamics simulation. The equations of motion are modified such asto allow the simulation cell vectors to evolve in time, in a way that produces thecorrect statistical properties of the isobaric-isoenthalpic ensemble.If we identify the cell vectors with a, b and c, we can construct a matrix hby stacking the vectors as its rowsh = {a,b, c}. (2.45)The equations of motion for the system subject to an isotropic pressure p can beobtained from the LagrangianL = K− V+ 12MTr(h˙T h˙)− pV, (2.46)where the additional terms can be interpreted as the kinetic and potential energyrelated to a piston. The term 12MTr(h˙T h˙) represents the kinetic energy of a pistonwith mass M that is compressing the cell, while pV is the potential energy drivingthe compression.The resulting equations of motion can be expressed as a function of the scaled32coordinates of the particles, si = h−1ri,s¨i = −∑i 6=jdudrij1mirij(si − sj)−G−1G˙s˙i, (2.47)h¨ =1M(pi − p)σ. (2.48)The other quantities defined in the above equations are σ, a matrix contain-ing the direction of the reciprocal vectors, G also termed the metric tensor, and pi,the stress tensor defined as followsσ = {b× c, c× a, a× b}, (2.49)V pi =∑imivivTi −∑i∑j>idudrij1rijrijrTji, (2.50)G = hTh. (2.51)If the system is far from equilibrium the Parrinello-Rahman barostat cangenerate large volume fluctuations and, in this case, the Berendsen barostat can beused to quickly bring the system to equilibrium before switching to the Parrinello-Rahman barostat. The Parrinello-Rahman barostat was used to control pressurein some of the simulations in this thesis. As will be discussed in Chapter 4, theBerendsen barostat was preferred for simulations where the fluctuations producedby the Parrinello-Rahman barostat caused the system to become unstable.33Chapter 3The Birth of NaCl Crystals:Insights from MolecularSimulations3.1 IntroductionAs stated in Chapter 1, the initial process by which one phase transitionsto another is termed nucleation. Nucleation has a central role in a wide range ofphysically interesting processes such as crystal and rock formation,81,82 drug de-velopment,83,84 and the formation of atmospheric aerosols.85,86 There have been anumber of interesting experimental studies of crystal nucleation,28,29,87–89 but thereremain significant challenges for current experimental methods, particularly for theanalysis of early prenucleation events. Therefore, to augment experimental stud-ies, molecular simulations are being increasingly used to gain physical insight intocrystal nucleation and growth.43–45,47,49,50,66,9034The usual approach to nucleation, classical nucleation theory (CNT), is asubject of current debate.26,91,92 In particular, the assumption that early stage po-tential nuclei can be characterized by properties of the nascent bulk phase (plus asurface term) is being challenged by both experiments and molecular simulations.One such example is provided by CaCO3, where experiments suggest that prenu-cleation clusters serve as precursors to nucleation.25 There is also evidence thatcrystallization can occur via a two-step mechanism consisting of an initial density“transition”, followed by a slower ordering transition whereby the dense liquid regionbecomes geometrically ordered.49,93,94In the present Chapter we use direct molecular dynamics simulations to in-vestigate the nucleation of NaCl nanocrystals. Supersaturated aqueous NaCl solu-tions are a good choice for detailed study because relatively accurate models exist,1and previous work has shown that spontaneous nucleation occurs on simulationtimescales.47,49,50,66 Earlier direct simulation studies of NaCl crystallization fromsupersaturated solution suggest a process by which less ordered and more hydratedNaCl clusters evolve with time into a largely anhydrous crystalline arrangement.49,59There is also some evidence that the two-step mechanism applies with a large density(concentration) fluctuation preceding any spatial ordering.49Nucleation and crystal growth have also been considered employing indirectsimulation techniques. In an early investigation Zahn47 approached the problemusing trajectory sampling methods. He noticed that small NaCl clusters tended toprefer a Na+ ion at the center, and suggested that such clusters might be importantin nucleation. A recent study reported by Zimmerman et al.66 used a seeded trajec-tory approach to determine the sizes of critical nuclei, ion attachment frequencies,and nucleation rates for NaCl in supersaturated solutions.35However, despite earlier efforts, the factors that influence the formation ofcritical nuclei have not been determined, and this is particularly true of the veryearly stages or the birth of the crystal. For example, it is interesting to ask if sizealone determines the relative stability of small NaCl clusters, or if other factors suchas the cluster geometry are important. This is one of the questions addressed in thepresent Chapter.Spontaneous nucleation events for NaCl are typically rare (on simulationtimescales) except at very high solute concentration, and multiple nucleation eventsmust be observed in order to draw any meaningful conclusions. Additionally, it isnecessary to define, detect, and follow in time a great many NaCl clusters that donot achieve nucleation. Here we carry out multiple direct simulations sufficient togenerate a number of nucleation events. We develop methodology to define anddetect crystal-like NaCl clusters, and to follow them in time from very early stages(∼ 6 ions) until nucleation is achieved or, much more frequently, the cluster dissolvesback into solution. An important conclusion of our analysis is that cluster size isnot the only factor influencing cluster lifetime and the probability of nucleation. Byintroducing a new parameter called the cluster “crystallinity” we show that clustergeometry is also a very significant factor influencing cluster lifetime and nucleationprobability; moreover, this is true for clusters as small as six ions.The remainder of this Chapter is divided into three parts. The models andsimulation methods are described in Section 3.2, results are presented in Section3.3, and our conclusions are summarized in Section 3.4.36σ (nm) (kJ mol−1)Na+ 0.2160 1.4754533Cl− 0.4830 0.0534924Table 3.1: Joung-Cheatham Lennard-Jones parameters for NaCl.13.2 Models and Methods3.2.1 Simulation DetailsIn our simulations we adopt the Joung-Cheatham parameter set1 for theNa+ and Cl− ions paired with the SPC/E water model.2 For this parameter set,the saturation mole fraction of NaCl, xNaCl, at 298 K was calculated60,64,95 to be∼ 0.06 using the chemical potential route, while results from direct coexistencemethods1,96–98 report a value of ∼ 0.09−0.11 (experimentally,99 the value is ∼ 0.10).With this force field all nonbonded interactions consist of Lennard-Jones (LJ) pluselectrostatic terms such that the site-site pair potentials have the formu(rij) = 4ij[(σijrij)12−(σrij)6]+14piε0qiqjrij, (3.1)where ε0 is the permittivity of free space, qi and qj are the partial charges on sites iand j, and σij and ij are the usual LJ length and energy parameters, respectively.The Joung-Cheatham parameters for Na+ and Cl− are given in Table 3.1, and σij andij are obtained from these and the SPC/E parameters2 using the Lorenz-Berthelotcombining rules, σij = (σi + σj)/2, and ij =√ij.Simulations were carried out under NPT conditions employing the GRO-MACS100 molecular dynamics package, version 4.5.4. The temperature was con-trolled using the velocity-rescale thermostat79 with a relaxation time of 0.1 ps, and37the pressure was kept constant at 1 bar using a Berendsen barostat78 with a com-pressibility of 4.5 × 10−5 bar−1 and a relaxation time of 1.0 ps. A timestep of 2 fswas used in all simulations. A spherical cutoff of 0.9 nm was applied to the pairpotentials, and the long-range electrostatic interactions were calculated using theparticle mesh Ewald (PME) method.74Details of all simulations performed are summarized in Table 3.2. Initialsimulations were carried out for xNaCl ranging from 0.20 to 0.30, and the lowestconcentration where a nucleation event was observed within 200 ns at 300 K wasxNaCl = 0.22. This time frame is feasible for an investigation employing directmolecular dynamics simulations, so we focus on this concentration. More nucle-ation events occur at higher concentrations, but the occurrence of many simultane-ous nucleations can complicate the analysis. Note that xNaCl = 0.22 is, based onthe estimates provided above, between two and four times the saturation value at298 K. The atomic coordinates for the system were collected at intervals of 0.04nanoseconds, and the coordinate sets collected with this time interval are referredto as time frames elsewhere in the paper.The analysis given in this paper is based on five “replica” simulations atxNaCl = 0.22 carried out as follows. Initially, 7040 NaCl pairs and 24960 watermolecules were distributed randomly on a lattice within the simulation cell. Thesystem was then “equilibrated” at 400 K and five configurations were extracted atintervals of ∼ 10 ns. Five simulations, initiated with these configurations, werecooled to 300 K, and allowed to evolve for ∼ 500 ns, as summarized in Table 3.2.The number of nucleation events observed in each simulation is also given in Table3.2.38Water Ions xNaCl Length (ns) Nuclei25600 12800 0.20 157 024960 14080 0.22 187 324000 16000 0.25 515 522400 19200 0.30 202 n.a.24960 14080 0.22 523 1547 2532 3530 4550 3Table 3.2: Summary of the simulations performed. The number of water molecules(Water), the number of ions (Ions), the mole fraction xNaCl, the run length (Length),and the number of clusters that achieved nucleation (Nuclei) are indicated. AtxNaCl = 0.3, multiple nucleations occurred resulting in a polycrystal. The fivesimulations at xNaCl = 0.22 listed at the bottom of the table were used for thestatistical analysis given in the text.3.2.2 Cluster DetectionThe aim of the present work is to detect NaCl clusters and follow their evo-lution in time in order to understand the factors influencing (or not) the probabilityof achieving crystal nucleation. From the perspective of classical nucleation theory(CNT), these clusters can be viewed as subcritical nuclei unless and until nucle-ation is achieved. Since we are trying to observe the very early stages of a phasetransition, it is necessary to find some measure to distinguish local liquid-like andsolid-like structures, and to determine which ions are part of the same solid-likecluster. Generally, such measures are based on local environments, and several ap-proaches to this problem have been employed. These include methods based onion connectivity,58,59 local ion density,101 local solvent density,47 and local bond or-der parameters.102,103 The last approach has been used extensively in studies of icenucleation.104–10639In the present paper, the computational procedures used to detect (define) ionclusters and follow their time evolution consists of three steps: filtering, clustering,and entity resolution, which are described in detail below.FilteringThe purpose of this step is to identify for each configuration along a trajectorythe ions that are most likely to be part of a local solid-like structure. To do thiswe follow the general bond-orientational order parameter approach of Steinhardt et.al.107 This method defines potential order parameters of the formql =[4pi2l + 1l∑m=−l|qlm|2]1/2,whereqlm =1NrN∑ri=r1Y ml (θ(ri), φ(ri)).In the above equation, Y ml is a spherical harmonic, ri is the position vector ofthe ith neighbor with respect to a central ion, and θ(ri), φ(ri) are, respectively,the polar and azimuthal angles with respect to an arbitrary frame of reference. Theactual choice of reference frame is irrelevant as the order parameter ql is rotationallyinvariant.The objective is to find the ql that best distinguishes between ions withliquid-like and solid-like local environments. To do this we examined ql distributionsfor l = 2, 4, 6, 8, and 10, using a spherical crystal of ∼ 2000 ions to represent thesolid phase, and a supersaturated solution at xNaCl = 0.20 (that did not nucleate)as representative of the liquid phase. Both reference systems were held at 300K. Based on this investigation, the order parameter q8 provides the best separation40Figure 3.1: Density of q8 values in a sample crystal (∼ 2000 ions)(blue) and a su-persaturated solution (xNaCl = 0.20) (green). The two distributions show excellentseparation. The vertical orange line indicates the selected threshold (q8 > 0.35).between liquid and solid phase distributions (Figure 3.1). Also, by testing a referencesolution of lower concentration (xNaCl ≈ 0.025), we found that the q8 distributionis largely independent of concentration (unlike order parameters based on local iondensities101). In the remainder of this chapter we refer to q8 as the bond orderparameter.The order parameter q8 is used to filter out ions that are part of liquid-likedisordered structures. From Figure 3.1, we see that the best separation point be-tween liquid-like and solid-like structures occurs at q8 ≈ 0.4. However, we selectq8 = 0.35 as the threshold for our analysis to capture structures from the upper endof the solution distribution. As clusters that achieve nucleation form spontaneouslyfrom solution, selecting a lower threshold value allows us to trace nucleation pro-cesses from the very initial stages. Selecting a very low threshold would include too41many transient, short-lived structures, while a very high threshold would filter outinteresting prenucleation stages. The value of 0.35 selected proved to be a goodcompromise between these two factors. We verified that our results and conclusionsare not affected by small variations in the selected threshold value. Note that withthe 0.35 threshold only ∼ 10% of the ions are classified as solid-like for each con-figuration. It is also important to emphasize that in our simulations all clustersthat eventually achieved nucleation originated within the set of ions identified assolid-like; thus our filtering process does not eliminate any interesting events.Cluster IdentificationWhile order parameter filtering efficiently detects ions in solid-like local en-vironments, another algorithm is necessary to identify which of the selected (solid-like) ions are part of the same cluster or aggregate. For this purpose we adopted thedensity-based cluster algorithm DBSCAN108 as implemented in the Python libraryscikit-learn.109 In the DBSCAN algorithm, clusters are identified by means of twoparameters, that in the present context are a distance () and the number of neigh-boring ions (η) within that distance. To find ions connected to each other, an ion isselected at random, and the number of ions within are counted. If the number ofions within is equal to or greater than η, then the central ion is labelled as core andits neighbors are labelled as border, otherwise, the ion is labelled as noise (Figure3.2). The procedure is repeated by picking an ion from the border. If that ion is core,it will generate a new border and the procedure is repeated recursively. Eventually,no more points are reachable and we label all of the collected core and border ionsas part of the same cluster. The procedure is then repeated by choosing another ionat random (that might belong to a new cluster) until all ions are processed.42Figure 3.2: Schematic diagram illustrating the DBSCAN algorithm for η = 4. Topleft panel: an ion is chosen at random and ions within are counted; the value is 4so the ion is marked as core (red). Its neighbours are labeled as border (pink). Topright panel: the procedure is repeated for one of the border ions which becomes acore ion. Bottom panel: all ions in the cluster have been processed, and no moreions are reachable from the core ions. The isolated ion is labeled as noise (blue).43The value of is the minimum distance for which two ions, are consideredpart of the same aggregate. A high value of causes some clusters that are closebut not connected to each other to be aggregated into a single cluster. On theother hand, a small value of would fail to connect ions that are effectively bonded(every ion would be detected as a single cluster). The value of η acts as a filterfor very small clusters, as every cluster that has fewer than η ions is discarded. Avery low value of η detects too many short-lived fluctuations that are difficult tofollow and analyse, while a high value would detect only large clusters providing noopportunity to investigate early stages of cluster formation. In the present analysisthe values = 0.3 nm and η = 6 ions are used. The value 0.3 nm is approximatelythe NaCl bond length, and η = 6 ions was chosen to be the minimum size cluster.These values proved to work well in practice, avoiding the problems noted above.Time Resolution of ClustersWe wish to follow the time evolution of all clusters identified as describedabove. Therefore, it is necessary to determine how clusters are connected in time.This would be a simple question if we had only a single cluster that grew or dissolvedas time progressed. However, in the present case we have many clusters that contin-ually change (divide, grow, dissolve) as the simulation advances in time. Therefore,we need to clearly define how we identify and follow clusters in time. Our generalapproach is to find “similar” clusters in successive configurations (using an appro-priate similarity index), and build a time connectivity graph using the followingprocedure.First we initialize an empty graph where each cluster detected appears as anode labelled ci(s), where i indicates the cluster and s the time frame (configuration).44Figure 3.3: Schematic diagram illustrating the cluster time resolution algorithm.At t = 0 the clusters are labeled employing the DBSCAN algorithm. A clusteris compared with clusters in the following frame, and a connection is made if theclusters are sufficiently similar, judged by the Jaccard index as described in thetext. The procedure is repeated recursively for the following frames. By performingthis procedure for all clusters, it is possible to obtain a connectivity graph (bottom)where each connected component represents the time evolution of a single cluster.Then for all consecutive frames pairwise similarities, dij(s) = sim(ci(s), cj(s + 1)),are calculated between clusters using the similarity measure defined below. If twoclusters in different frames are detected to be sufficiently similar (satisfy an appro-priate threshold) they are deemed to be the same cluster at different points in time,and the nodes representing these clusters are connected through an edge. The out-come of this procedure is a graph that connects clusters in different time frames,such that every connected component of this graph is the trajectory of a particularcluster as illustrated in Figure 3.3.45Several similarity measures could be used to look for matching clusters, andhere we employ the Jaccard index,110 defined as the size of the intersection dividedby the size of the union of two sets,J(A,B) =|A ∩B||A ∪B| .Here A and B represent two clusters in consecutive time frames, and clearly J(A,B)is one, if both clusters contain exactly the same ions, and zero if they do not share anyions. The threshold value of J(A,B) must be selected with care to avoid ambiguitiesin cluster identity. For example, if two clusters, a and b, detected in frame one mergeinto a single cluster in frame two, it is not clear if the merged cluster is a, b, or anentirely new cluster c. Similar ambiguities can arise from other possible clusterevolution scenarios. In order to assign a unique well-defined time line to a cluster,it is necessary to remove all ambiguities. A natural way to do this is to connectclusters in consecutive frames only if they satisfy the threshold J(A,B) > 0.5. Bydefinition of the Jaccard index, this threshold ensures that a cluster will match atmost one cluster in the following frame, giving a set of unique cluster trajectories.One potential issue with this method of obtaining cluster trajectories is howto take account of splits and merges that happen over multiple time frames. Forexample, if we have a large cluster (e.g., 500 ions) and a piece (e.g., 100 ions) becomestemporarily detached, the connectivity with the original cluster would be lost andthe detached piece would appear as a newborn cluster of 100 ions. This problemsis handled by detecting abnormally large newborn clusters and merging them withthe original cluster. Fortunately, such events are rare (it happened once in our fivesimulations) and the post processing step was sufficient to solve the problem.46Another possible issue is that small clusters can partially redissolve (in thesense that some of their ions become less ordered), fall below the 6 ion threshold,and go undetected for a a few time frames. To solve this problem we modified thealgorithm to look for a match in 10 consecutive frames following the initial clusterdetection. Newborn clusters that had no match in the following 10 frames wereassumed to be short-term fluctuations, and discarded from the analysis.3.3 Results and DiscussionIonic clusters form and disintegrate continuously in the supersaturated solu-tion. These clusters are detected and followed in time as described above. We areinterested in understanding what features (if any) of ionic clusters influence theirability to survive and eventually achieve nucleation. We considered several clusterattributes or properties, but only two had a significant easily observable influenceon the nucleation probability of a cluster. One of those, cluster size (calculated asthe number of ions in the cluster), is of rather obvious importance, and is in fact thecrucial parameter in the usual application of CNT, which assumes spherical nucleithat become critical at a certain radius. The other cluster property that proved veryinfluential, especially so in smaller clusters, we call q¯8, or the cluster crystallinity.This property is simply the average value of the q8 bond order parameter definedabove, taken over all ions in a cluster.Other cluster properties considered include volume, surface area, sphericity,average neighbor count, hydration, and radius of gyration. These properties areprecisely defined and discussed in Appendix A. For smaller clusters (e.g., 10 ions),none of these properties showed a significant correlation with the probability of nu-47cleation. For larger clusters (e.g., 30 ions), we do find noticeable correlations forthree properties, with smaller surface areas, higher sphericities, and smaller radii ofgyration all favoring nucleation. These attributes are all measures of cluster com-pactness, and it is not surprising that they are correlated with nucleation probabilityfor larger clusters.We note that all results reported below are based on the five simulationsof 523 − 550 ns at xNaCl = 0.22, as listed in Table 3.2. In our five simulations atotal of 13 nucleation events were observed. In these simulations we noticed thatclusters with a lifetime of at least 30 ns that reached a size greater than ∼ 50 ionsnever redissolved, and continued to grow during the simulation. These clusters weredeemed to have achieved nucleation. Note that we would expect this observational“definition” of nucleation to be concentration dependent.At inception (defined as the earliest time a cluster is detected) the ionicclusters span a rather wide range of size and crystallinity. The size and crystallinitydistributions of all clusters detected are shown in Figure 3.4, the top left and topright panels, respectively. Note that on average our algorithm detects ∼ 2 newclusters per time frame, giving a total of ∼ 80000 clusters detected over the courseof the five simulations. It is also worth noting that the filtering step of the algorithmremoves ∼ 90% of the ions, indicating that only ∼ 10% of the ions present in solutionbelong to a cluster of any sort.From the histogram of initial sizes (Figure 3.4, top left panel) we see thatthe peak occurs at 6 ions (the smallest cluster detected by our algorithm), and thatthe frequency of larger sizes follows an exponential-like decay. Although very rare,clusters of more than 30 ions are found but, as noted below, these correspond toelongated, amorphous-looking structures. The distribution of crystallinities (Figure48Figure 3.4: Size (top left) and crystallinity (top right) histograms for clusters at firstdetection. The bottom panel shows the joint distribution of size and crystallinity.Note that very few clusters begin with both high size and high crystallinity, suggest-ing that this feature develops as clusters evolve in time. The red points indicate thethirteen clusters that achieve nucleation, and we note that these show no obviouspreference for any region of the joint distribution.493.4, top right panel) is skewed to the right, but one does not observe the high valuescharacteristic of a bulk crystal (Figure 3.1). A scatter plot of cluster crystallinityversus size is shown in the bottom panel of Figure 3.4. It is apparent that at clusterinception the correlation between these variables is modest at most. For smallerclusters, there are large fluctuations of crystallinity for a given cluster size, withlarge and small values occurring with high frequency. For larger cluster sizes, thecrystallinity fluctuations decrease markedly, which is not surprising, since statisti-cally variance is expected to decrease with increasing sample size, which here is thenumber of ions in the cluster.Ionic clusters exhibit a range of crystallinities and the examples given inFigure 3.5 provide an idea of how cluster structure and crystallinity are related.Newly detected small clusters can have low (bottom left panel) or high (bottom rightpanel) crystallinity, with high values being associated with more regular structures.For larger clusters (more than 30 ions) we do not find any very high crystallinities atthe point of initial detection, as such clusters tend to have elongated shapes (Figure3.5, top left panel.)3.3.1 Cluster LifetimesAs discussed above, clusters are born with certain characteristics, and, aftersome time, they will either achieve nucleation (very rare) and continue to grow asa crystal or, much more commonly, dissolve back into solution. Given that onewould expect a connection between cluster lifetime and nucleation, it is of interestto investigate which cluster characteristics influence their lifetime, defined here asthe total time a cluster is detected (by the algorithm described above) as a distinctentity in solution. As, by definition, clusters that nucleate have infinite lifetimes,50Figure 3.5: Examples of clusters of different sizes and crystallinities (given in paren-thesis below each cluster). Na+ and Cl− ions are colored purple and green, respec-tively. At inception larger clusters tend to be elongated and of lower crystallinity(top left), whereas larger high crystallinity clusters (not observed at early stages)tend to be more compact (top right). For small clusters structural features as-sociated with low (bottom left) and high crystallinity (bottom right) are not soapparent.51Figure 3.6: Histograms of “failed” cluster lifetimes (left) and the overall survivalfunction (right) for all clusters detected in the simulations.we excluded these from our statistical analysis. Note that in the five simulations(Table 3.2) included in the analysis only 13 clusters achieved nucleation, whereas∼ 80000 failed to nucleate. Focusing on the large number of failures provides uswith excellent statistics, and clearly, the influence of the very few that did nucleatewould have a negligible effect on the lifetime and survival distributions discussedbelow.The overall lifetime distribution of failed clusters (uncontrolled for size orcrystallinity) is shown in the left panel of Figure 3.6. This plot shows that the greatmajority of clusters are very short lived (lifetimes less than a nanosecond), but thereare a significant number of clusters that live for more than 10 ns.A more instructive way to represent lifetime information is to construct sur-vival functions. The survival function, S(t), is defined as the probability that acluster has a lifetime, τ , larger than t, orS(t) = Pr(τ > t).The survival function starts at 1 (all clusters have a non-zero lifetime) and52decreases with t, depending on the surviving fraction of the population. The survivalfunction can be written as the cumulative productS(t) = R(1)R(2) . . . R(t) ,withR(t) =Population(t)Population(t− 1) ,where Population(t− 1) represents the population at frame t− 1, and Population(t)is the surviving population at the next frame t. R(t) is the probability of survivingfor one time interval at t. We note that this method of calculating S(t) is commonlyused, and is known as the Kaplan-Meier estimator.111 The overall survival functionincluding all clusters is shown in the right panel of Figure 3.6, and we note that over90% of all clusters survive less than 2 ns.Survival functions can be used to investigate if and how particular clustercharacteristics influence their probability of survival. To isolate the influence ofcrystallinity, we divide clusters of fixed size into groups of high (> 0.40) and low(≤ 0.40) crystallinity, and compare the survival functions. Note that the survivalfunctions are constructed using the size and crystallinity of a cluster at first detec-tion. Obviously, neither of these cluster characteristics remains fixed as the clusterevolves in time. Results for clusters of 10 ions are shown in Figure 3.7, and it isobvious that crystallinity has a large effect on the survival probability. Clusters inthe higher crystallinity group have a substantially higher survival probability thanclusters in the lower crystallinity group. This is true for all cluster sizes, remarkably,even in the 6 ion case, which is the smallest cluster size we consider. This is illus-trated in Figure 3.8, where we plot the ratio (high/low) of the median lifetimes for53Figure 3.7: Survival curves for clusters of size 10 grouped by high (> 0.40) and low(≤ 0.40) crystallinity.clusters of high and low crystallinity as a function of cluster size. We see that theratio is significant (∼ 1.5) for clusters of six ions, and rises rapidly with increasingcluster size, reaches a plateau at clusters of ∼ 10 ions, then rises again for clusterscontaining more than ∼ 18 ions. Results for clusters larger than 30 ions are notincluded because for the larger clusters the statistical analysis is less reliable due tosmall sample sizes. In an analogous manner, one can fix crystallinity and calculatehow survival functions vary with size. This is less interesting, simply showing that,as expected, larger clusters have a larger survival probability.While cluster survival is not a direct measure of nucleation probability, onewould expect both phenomena to be related to cluster stability. Therefore, based onthe above analysis, we would expect cluster crystallinity to be an important factorinfluencing the probability that a cluster will achieve nucleation. Specifically, for54Figure 3.8: Ratio of median lifetimes for high versus low crystallinity clusters ofdifferent size. The high crystallinity clusters always have the higher median lifetime.fixed cluster size higher crystallinities should favor higher nucleation probabilities.The relationship between cluster crystallinity and nucleation is treated directly inthe following section, and we show that the behavior is indeed as we would anticipatebased on the survival analysis. We remark that the survival analysis has the greatadvantage that the large number of clusters (events) included provides convincingstatistics. On the other hand, we observe only a few nucleation events, so a statisticalanalysis based on nucleation alone would be less convincing.3.3.2 Crystallinity and NucleationTo get a qualitative idea of the influence of cluster size and crystallinity onnucleation, it is useful to plot cluster trajectories in (size, crystallinity) space. Thisis done in Figure 3.9 for a cluster that nucleates (blue) and one that does not (red).55Figure 3.9: Trajectory shown in (size, crystallinity) space for a nucleated (blue) anda failed (red) cluster as they evolve in time. In early stages both clusters oscillatein size up to ∼ 30 ions, but the crystallinity of the cluster that eventually achievesnucleation reaches higher values than the failed case. This suggests that it is acombined effect of both size and crystallinity that promotes nucleation.It is easy to see that after a brief residency at low sizes and low crystallinities, thecluster that nucleates appears to surpass a “critical region” in the (size, crystallinity)space and never falls back. This suggests that both size and crystallinity influencethe probability of nucleation.This can also be seen in the early stage time evolution profiles shown in Fig.3.10 for the nucleating and failing clusters. We see that while relatively long lived(∼ 12 ns), and sometimes exceeding 30 ions in size, for most of the trajectory thecrystallinity of the failing cluster lies below that of the nucleating case. Snapshots(not shown) suggest that the failing cluster grows in an elongated manner with aresulting lower crystallinity. The growth profile for a nucleating cluster is shown56in Fig. 3.11. It is interesting to note that the growth rate increases with size,attributable to the growing surface area of the nucleus. The crystallinity profileinitially grows but plateaus as the cluster becomes larger in size, and the crystallinityapproaches that of a perfect crystal.The influence of cluster size and crystallinity on nucleation probability can beexplored more quantitatively by comparing statistically these properties for clustersthat nucleate with those that fail. If cluster size were the only factor influencingnucleation, then two clusters of the same size would have the same probability of nu-cleation. Therefore, if we fix the cluster size and observe a difference in crystallinitybetween nucleations and failures, we can isolate the effect of crystallinity.One way to fix cluster size is to follow the trajectory of each cluster in timeuntil it reaches a particular size, at which point the cluster crystallinity is recorded.Values obtained in this way can be reasonably assumed to be independent becausethe clusters exist at different points in time or in different simulations, and weregister only a single crystallinity value for each cluster. Crystallinity distributionsobtained as described are shown in Figure 3.12 for clusters of different size. Thecrystallinities of the clusters that achieved nucleation are also indicated in the plots,and we see that for all cluster sizes, the crystallinities of the nucleated clusters fallmainly in the upper half of the crystallinity distribution. In other words, the clustersthat eventually nucleate have on average higher crystallinities than clusters that fail,and it is remarkable that this distinction exists even for clusters as small as six ions.For comparison, similar plots for other properties are reported in Appendix A.The influence of cluster crystallinity on nucleation is further demonstratedin Figure 3.13, where the probability of nucleation, P (N), for clusters of N ions isshown for clusters of “high” and “low” crystallinities. If the cluster crystallinity is57Figure 3.10: Comparison of growth and crystallinity profiles for nucleated (colored)and failed (gray) clusters. After a short period, the nucleated cluster manages toreach quite a high crystallinity and maintain its size. In contrast, the failed cluster,while maintaining its size, experiences a steady decrease in crystallinity.58Figure 3.11: Size (blue) and crystallinity (orange) profiles for a nucleated cluster.The fluctuation in crystallinity is fairly high, and, as the crystal increases in size,the crystallinity reaches a plateau at ∼ 0.47. The snapshots represent the cluster at0 ns (left), 40 ns (center) and 80 ns (right).59Figure 3.12: Crystallinity distributions (blue histogram) of failed clusters of differentsize (number of ions). Clusters that achieved nucleation are indicated by singleorange lines. Note that clusters that achieve nucleation come preferentially fromthe upper part of the crystallinity distribution.60> 0.40, it is labelled high, otherwise it is labelled low. The nucleation probabilitiesP (N |high) and P (N |low) are shown in Figure 3.13 for N ranging from 6 to 35 ions.Results for larger clusters are not shown because the sample size is too small to givegood estimates of the probabilities. We note that P (N |high) is much higher thanP (N |low) for all cluster sizes. The effect is most pronounced for small clusters (whichhave a very low total probability of nucleation) where P (N |high) is up to eighttimes larger than P (N |low). The ratio P (N |high)/P (N |low) generally decreaseswith increasing N , but remains at ∼ 3 for N = 35.This analysis highlights the fact that nucleation is substantially influencedby factors other than cluster size. For a given cluster size, geometrical order, ascaptured by the crystallinity parameter, increases cluster lifetime and the probabilityof achieving nucleation.Cluster Binding EnergyIt is interesting to ask why high crystallinity leads to increased cluster life-times and influences the probability that a cluster will achieve nucleation. Onepossibility that immediately comes to mind is that in high crystallinity clusters thedirect ion-ion interactions lead to lower cluster energies and hence increased stability.To investigate this possibility we calculated the energiesUk =12Nk∑iNk∑ju(rij) , (3.2)where Nk is the total number of ions in cluster k, and u(rij) is the ion-ion interactionas defined in Section 3.2. One potential problem with this analysis is that largefluctuations can occur in the apparent cluster energy depending on whether or not61Figure 3.13: Probabilities of high (q¯8 > 0.40) and low (q¯8 < 0.40) crystallinitiesachieving nucleation. The error bars (vertical black lines) represent one standarddeviation. The lower panel shows the results for small clusters on an expanded scale.Note that for all cluster sizes P (N |high) is always larger than P (N |low), and thehigh/low ratio is shown as a gray line in the upper panel.62a cluster has a net charge. These fluctuations can be created by a single ion movingin or out across the defined cluster boundary, and hence can be rather arbitraryand artificial. Therefore, to avoid this problem we include only electrically neutralclusters in this analysis.Joint distributions of crystallinity and cluster energy for clusters of six andten ions are shown in Figure 3.14. For clusters of six ions the correlation coefficientr = −0.074 indicating only a very weak correlation between energy and crystallinity.For clusters of ten ions, r = −0.22 showing that the correlation increases with clustersize, but still remains rather weak. Moreover, we see from Figure 3.14, that theclusters that eventually nucleate come mainly from the high crystallinity tail of thecrystallinity distribution, but are nearly evenly spread over the energy distribution.These observations show that crystallinity is not a measure of increased clusterstability due to lower cluster energies coming through the direct ion-ion interactions.Clearly, other factors must be involved.The Role of Na+ and Cl− Ions in Small CrystalsIt is also of interest to more closely examine the structural nature of smallNaCl clusters, since these clusters represent a fundamental step in the nucleationprocess. To do this we divide all ions in the system into two classes, high q8 (≥ 0.4)and low q8 (< 0.4), and calculate running counterion coordination numbers for Na+and Cl−.Results for single frames at 40 and 240 ns are shown in Figure 3.15. Forthe high q8 class we see that both Na+ and Cl− have more counterions (3-4) in thefirst coordination shell (∼ 0.35 nm) than the solution average (∼ 1). This is notsurprising because the high q8 ions generally belong to ionic clusters. However, from63Figure 3.14: Joint distribution of crystallinity and cluster energy (kJ mol−1) forclusters of size 6 (left) and size 10 (right). Values corresponding to clusters thateventually nucleated are displayed in orange. The correlation coefficient is negativeand very low at size 6, and moderately low at size 10, indicating that crystallinitycarries different information than energetic stability.Figure 3.15 we also notice that high q8 Na+ ions have more counterion neighborsthan high q8 Cl− ions, suggesting that Na+ ions tend to lie more within the “interior”of small clusters than Cl− ions, consistent with an earlier observation of Zahn.47 Apossible reason for this is that the smaller Na+ ions interact more strongly with watermolecules, and hence more first shell counterions are needed to compensate for watermolecules lost from the first hydration shell when ionic clusters are formed. Notethat the gap between the coordination numbers of high q8 Na+ and Cl− decreasesin later frames as larger clusters develop in the simulation.To obtain additional evidence that Na+ and Cl− do not contribute symmet-rically to the structure of small clusters, for each cluster we calculate the parameter∆q8 = 〈q+8 〉 − 〈q−8 〉,64Figure 3.15: Running coordination numbers for: high q8 Na+ with all Cl− (blue);high q8 Cl− with all Na+ (green); and all Na+ and Cl− (red). In the first coordinationshell (∼ 0.35 nm) both high q8 ions are surrounded by more counterions (3-4) thanaverage (1). Also, high q8 Na+ are surrounded by more counterions than high q8Cl−, suggesting that on average they occupy positions deeper within the clusters.The difference in the first shell coordination number decreases at longer times, aslarger crystals develop in the simulation, and the internal preference of Na+ becomesless noticeable.where 〈q+8 〉 and 〈q−8 〉 are the average values q8 for the Na+ and Cl− ions in thecluster.Figure 3.16 shows ∆q8 as a function of cluster size. These results are from thesecond replica simulation (Table 3.2). We note that for small clusters ∆q8 is positiveon average indicating that Na+ tends to lie in more ordered environments than Cl−.As we would expect, the distinction between Na+ and Cl− grows smaller as clustersgrow in size, and ∆q8 approaches zero for large clusters. Survival functions andnucleation probabilities were analysed separately for positive and negative valuesof ∆q8, but no significant differences were found. Thus, while there is a structuraldistinction for small clusters, this effect does not have a significant influence oncluster survival or the probability of nucleation.65Figure 3.16: ∆q8 as a function of cluster size. The shaded area indicates one standarddeviation. The average value is slightly positive indicating that Na+ ions tend tooccupy ordered environments than Cl− ions.663.4 Summary and ConclusionsDirect molecular dynamics simulations have been employed to identify andinvestigate factors that influence the nucleation of NaCl crystals in model super-saturated aqueous solutions. We develop methods that allow potential nuclei tobe detected as small clusters (∼ 6 ions), and followed in time until nucleation isachieved, which occurs very rarely, or the cluster dissolves back into solution.Our analysis clearly demonstrates that cluster size is not the only propertythat has an important influence on the expected lifetime and nucleation probabilityof a particular cluster. We show that the geometric arrangement of the ions in thecluster, as measured by a single parameter which we call the cluster “crystallinity”is also very influential. For example, for clusters of ten ions the median lifetimefor clusters of high crystallinity is double that of those of low crystallinity, andtheir probability of achieving nucleation is ∼ 8 times greater. Similarly, smallerand larger clusters also have significantly longer lifetimes and greater probabilitiesof nucleation.Physically, it is not entirely clear why crystallinity (as measured by our pa-rameter) has such a large influence on cluster lifetime and nucleation probability,and this is especially true for small clusters of six or ten ions. One possibility is thatthere is a connection with the binding energy of a cluster. However, we did not finda strong correlation between binding energy and our crstallinity parameter, indi-cating that the crystallinity parameter is not merely a proxy for energetic stability.For small clusters, we did find that Na+ ions had some preference for a “central”position, in accordance with Zahn’s observation,47 but such cluster arrangementshad no influence on the probability of nucleation.67Given that there is no obvious energetic explanation for the increased stabilityof small clusters of higher crystallinity, it appears that the advantage must lie inthe microscopic dynamics of cluster growth and/or disintegration. One possibilityis that ions in more ordered environments are less exposed to water molecules, andhence high crystallinity clusters are less susceptible to disintegration.68Chapter 4Crystallization of Lithium Halidesfrom Molecular Simulation.4.1 IntroductionAlkali halides are inorganic compounds with high melting points that areusually found in nature as crystalline solids. Due to their importance in industrialand scientific applications such as corrosion,112 desiccants, nanotube preparation,113but also because of their simple structure, they have been widely studied for manyyears both in experiments and theory.114,115 More recently, sodium chloride has beenused as a model compound in nucleation from solution,47,49 providing importantinsights in the nucleation process. However, not much has been done regardingother compounds, partly because their phase behavior in simulation has not beenthoroughly examined.In nature, the stable crystal structures for alkali halides are rock salt (FCC)and CsCl (BCC), however both experimental116–118 and theoretical studies show69the presence of wurtzite-type crystal structures,119–123 especially in lithium halides.Using ab-initio methods, the wurtzite-type structure was shown to be metastable atstandard pressure for most of the alkali halides,122 and its existence was confirmedthrough low temperature deposition experiments (aimed at capturing metastablestructures) for LiCl,116 LiBr117,118 and LiI.118 Four coordinated structures were alsohypothesized by Pauling,124 who derived the radius ratio rules using a simple hardsphere model.Classical simulations of LiCl crystals also showed that, under the Tosi-Fumi125potential, the preferred structures achieved by melting and freezing small clusterswere described as having hexagonal motifs.121 These structures were found to beenergetically favorable (with respect to the rock salt structure) for sizes up to 64ions.121 The competition between the hexagonal and cubic structure of LiCl waslater found to be also present in larger systems.119Despite the evidence of the existence of a metastable wurtzite structure in al-kali halides, a systematic review of the most commonly used classical potentials withrespect to this structure has not been performed. Such a study is important for thedetermination of solid and solution properties (e.g. solubilities) which could greatlyaffect the nucleation mechanism from the melt and/or from aqueous solutions.Through molecular simulation we study the role of the wurtzite structure forlithium halides in both infinite crystals and finite size clusters. By comparing twoestabilished pair potentials used62–65,126–130 in classical simulations, the Tosi-Fumi125(TF) and Joung-Cheatham1 (JC) parameter sets, we find that the wurtzite structureis the lowest potential energy minimum for some of the lithium compounds. Thissuggests that current models need to be adjusted to reflect the higher stability ofthe rock salt structure, in order to agree with experiments.704.2 Models and MethodsAs noted above, in this Chapter we consider two commonly used classicalmodels for alkali halides. The TF potential125,131 is of the Born-Mayer-Huggins,132–134consisting of a short-range repulsion combined with electrostatic interactions, andattractive dispersion terms of dipole-dipole and dipole-quadrupole order. The inter-action uij(r) between ions i and j can be expressed asuij(r) =14piε0qiqjr+Bije−αijr − Cijr6− Dijr8. (4.1)The parameters required are the charges qi and qj, Bij, Cij, and Dij, which dependon the particular pair interaction (with different sets for different salts), and αijwhich is a constant depending only on the salt. The TF potential was obtainedby fitting experimental data (e.g. equation of state and volume derivatives) forpure salts. Because of its focus on pure salts, this model has been mainly used forcrystalline and molten salts,102,126–130 but it has also been adapted for use in aqueoussolution.98The JC model combines electrostatic and Lennard-Jones (LJ) interactionsand has the formuij(r) =14piε0qiqjr+ 4ij[(σijr)12−(σijr)6], (4.2)where σij and ij are the usual LJ length and energy parameters. Unlike the TFpotentials, the JC interactions were tuned to reproduce aqueous solution properties(e.g. solvation free energies), as well as solid state properties (e.g. lattice constantsand energies). However, we note that all solid state calculations assumed the rock71salt structure, and no attempt was made to determine if a more energetically sta-ble crystal structure existed. The parameter sets given by JC for a particular saltvary a little depending on the water model employed in calculations of the solutionproperties. Here we use the set given for the SPC/E water model, which is themodel we employ throughout this thesis. We note that the rock salt lattice energiesobtained for these parameters are in good agreement with experiments.1 JC param-eter sets have been used extensively to study solubility and solid-liquid transitionsin solution.60–65Lattice energies for LiX (X = F−, Cl−, I−, Br− and I−) salts are obtainedfor both rock salt and wurtzite crystal structures. The unit cells (see Figure 4.1) forthe two crystal types were generated by calculating the unit cell parameters (a, b, c),angles (α, β, γ), and atomic fractional positions given the parameter a (the unit cellconfiguration for both the wurtzite and rock salt structures is completely determinedby a, as the other parameters are constrained by symmetry.) The crystal structures,as well as the final crystals, were generated using the chemlab python library.135 Forthe models considered, the lattice energy can be divided into short-range (van derWaals) and Coulombic contributions, denoted Ew and Ec, respectively, such thatthe total energy E = Ew + Ec.The short-range interaction felt by ion i is given byEwi =nmax(rc)∑n=−nmax(rc)N∑jusrij (|sij + n|) with i 6= j, when n = {0} , (4.3)where usrij is the short-range part of the pair interaction, sij is the difference betweenthe position vectors of of ions i and j in the unit cell, and n is the displacement ofthe periodic cell. The first sum in Eq. (4.3) is over periodic repetitions of the unit72Figure 4.1: Unit cells used for the rock salt (right) and wurtzite (left) crystal struc-tures. The configurations are dependent on a single cell parameter a, all otherparameters being constrained. In the rock salt structure, a is the length of anyedge, while for the wurtzite structure a is the length of any of the two short edges.cell, and nmax is the minimum cell displacement needed to include all ions within acutoff radius, rc, of the central ion. In the present calculations, rc is taken to be 2.4nm.The electrostatic interactions were calculated using the Ewald summationmethod73 with conducting boundary conditions. Note that conducting and vacuumboundary conditions will give the same lattice energy for rock salt because, for thatcrystal structure, the unit cell does not have a dipole moment, but this is not truefor wurtzite where a dipole moment is present. The Ewald parameter κ (see Ref.73) was set at ∼ 5/a (a is the cell parameter), and 1000 wave vectors were used inthe Fourier space sum. Short-range interactions in the real space part of the Ewaldsum were truncated at a. The lattice energies obtained were in excellent agreement(within ∼ 0.001 kJ mol−1) with values reported in the literature.136Since in both rock salt and wurtzite crystals there are only two types of lattice73site (positive and negative ion sites), the total lattice energy can be written as thesum of two terms E = E+ +E−, where E+ and E− are contributions from the ionson different sites interacting with the remainder of the lattice. This decompositionof E gives some useful insight into the relative stabilities of rock salt and wurtzitestructures for various salts (see Section 4.3).Molecular dynamics simulations were carried out for both finite clusters andinfinitely periodic systems. Starting configurations were generated by repeatingthe rock salt unit cell along the 3 crystal axes, 3, 4 and 5 times for finite clusters(corresponding to 256, 512, and 1000 ions, respectively), and 5 times for infinitelyperiodic systems.Simulations were conducted using the GROMACS molecular dynamics pack-age version 4.6.100 Finite clusters were simulated using a cubic simulation cell oflength L that was ∼4-5 times larger than the cluster, such that interactions with itsperiodic images were negligibly small. All interactions were spherically truncated atL/2, and the cell volume was kept constant during the simulation. The temperaturewas regulated using the velocity-rescale thermostat.79 These simulation conditionseffectively correspond to a finite size crystal at zero external pressure.Periodically infinite systems (both liquid and solid) were simulated by apply-ing periodic boundary conditions in the usual manner.73 All short-range interactionswere spherically truncated at 0.9 nm, and the long-range electrostatic interactionswere calculated using the particle mesh Ewald (PME) method.74 The pressure wasfixed at 1 bar using the Berendsen barostat,78 and the temperature was regulatedusing the velocity rescale thermostat.79 The Berendsen barostat was chosen becauseit produced less drastic fluctuations during convergence, making the simulation lesslikely to become unstable (especially with the TF potential). Some tests were per-74formed with the Parrinello-Rahman barostat137 to verify that the choice of barostatdid not influence the results.We are mainly interested in determining crystal structures at finite tem-perature by spontaneous nucleation and growth. Therefore, the simulations wereperformed in the following manner. The simulations were started at 300 K in a rocksalt configuration and equilibrated for ∼ 0.1 − 0.5 ns. The temperature was thenraised to ∼ 1200 − 2000 K (depending on the melting point of the crystal) over atime period of ∼ 1 ns, causing the salt to melt. The molten salt was then frozen bycooling the system to ∼ 500− 600 K over a time interval of ∼ 5 ns. The solid wasthen maintained at the freezing temperature for ∼ 5− 20 ns to allow the system toanneal (the difference in annealing times is because small crystals require less timesto anneal). Finally, the temperature was lowered to 300 K to collect properties forcomparison with the initial rock salt crystal equilibrated at the same temperature.After this procedure, the final structure was determined by visual examina-tion, and by comparing radial distribution functions of the initial and final frames ofthe trajectory, as illustrated in Figure 4.2. A summary of the simulations performed,along with the resulting crystal structures, is reported in Table 4.4 and 4.3.4.3 Results and Discussion4.3.1 Lattice EnergiesTo illustrate how the TF and JC potentials differ, we compare the pair inter-actions for LiCl in Figure 4.3. One notices (see the top and bottom panels of Figure4.3) that the Cl− ion is slightly “larger” for the JC model. Also, the repulsive part75Figure 4.2: Radial distribution functions for an infinite LiCl crystal for the TF model(left panel) and JC model (right panel). The initial rock salt configuration and theliquid phase are colored in blue and red, respectively, and do not show noticeabledifferences between models. The rdf corresponding to the crystallized structures aredisplayed in green and correspond to wurtzite on the left and rock salt on the right.The inset plots show the rdf for ideal wurtzite on the left and ideal rock salt on theright (black lines) superimposed onto the simulation result (green).76of the potential, which dominates at short range, grows less sharply for the TF thanfor the JC model (this is due to the different functional form, ∼ e−r vs ∼ r−12,employed in the potentials). For the attractive cation-anion interaction (middlepanel of Figure 4.3), there is a difference in the well depth and in the position ofthe minimum, with the TF interaction more attractive than JC. At longer range thepotentials are equally dominated by the electrostatic interaction (r−1).Of course simply comparing pair potentials in this way does not provide muchinformation on interactions in particular crystal structures, therefore we now turnto the lattice energies. For lithium halides occupying rock salt or wurtzite lattices,the electrostatic (Coulombic) interaction is the same for both positive and negativeion sites such that 2Ec+ = 2Ec− = Ec, where Ec is the total electrostatic contributionto the lattice energy. For the crystal structures we considerEc =14pi0M±r0, (4.4)where r0 is the distance of an ion from its nearest neighbor, and M± is the Madelungconstant, which depends only on the crystal structure. Since M rocksalt± > Mwurtzite± ,the electrostatic contribution to the lattice energy, for fixed r0, is always lower forthe rock salt structure (a consequence of the higher coordination number). Thewurtzite structure can only become electrostatically more stable than rock salt byallowing the ions to pack more tightly, thus increasing the electrostatic interactionby reducing r0.Equilibrium ion-ion distances and the minimum lattice energies are signifi-cantly influenced by the short-range repulsive part of the pair potential. In the caseof oppositely charged ions, the short-range repulsion counteracts the Coulombic77Figure 4.3: TF and JC potentials for the pairs Li+-Cl− (top panel), Li+-Li+ (middlepanel) and Cl−-Cl−(bottom panel). The TF model has a shorter contact radius forthe Cl− ion, and is more slowly varying function at short range. Notice also how atvery short range (∼ 0.1 nm), the TF model becomes attractive (but this does notinfluence the simulations).78attractive forces creating a well defined “contact” distance (Figure 4.3, top panel).The interaction between ions of like charge is always repulsive, with a sharp increaseat short range that strongly penalizes close contact (Figure 4.3, bottom panel). Thebalance between the two competing interactions determines the nearneighbor equi-librium distance and the coordination number.In general the lattice-positive-ion interaction, E+, will be different from thelattice-negative-ion interaction, E−, because the short-range, like-charge interactionsare different. To obtain the lattice energy it is necessary to find the structure thatminimizes the sum of those two contributions.Lattice energies for perfect lithium halide crystals at 0 K are given in Tables4.1 and 4.2 for the JC and TF models, respectively. We see from that for the JCmodel (Table 4.1) the rock salt structure has the lower lattice energy for LiF andLiCl, but that the wurtzite energy is lower for LiBr and LiI. For the TF model,(Table 4.2) the wurtzite lattice energy is lower than that of rock salt for all fourlithium halide salts. Thus, only the JC results for LiF and LiCl are in accord withthe rock salt structure experimentally observed for all lithium halides.138 Latticeenergies for NaCl are included in Tables 4.1 and 4.2 for comparison purposes, andwe note that in this case the rock salt structure has the lower lattice energy for boththe JC and TF models.One can gain some insight into this behavior from Figure 4.4, where thelattice energies for lithium halides are plotted as functions of the cell parameter a.The lattice-positive-ion and lattice-negative-ion contributions, respectively, are alsoshown. Note that in the figure the actual plots represent 2E+ and 2E−; this preservesthe same vertical scale such that features of all three curves (e.g. the positions of theminima) can be easily compared. For the JC model (left column, Figure 4.4) we see79Compound Erock(kJ mol−1) Ewurtz(kJ mol−1) ∆E (kJ mol−1)LiF -1052.56 -1030.42 -22.14LiCl -870.31 -863.49 -6.82LiBr -826.04 -826.67 0.63LiI -761.83 -770.08 8.25NaCl -793.15 -769.67 -23.48Table 4.1: Lattice energies for the JC model for the wurtzite (Ewurtz) and rock salt(Erock) structures. The difference Erock- Ewurtzis reported as ∆E. The NaCl latticeenergies are added for comparison.Compound Erock (kJ mol−1) Ewurtz (kJ mol−1) ∆E (kJ mol−1)LiF -1043.04 -1051.83 8.78LiCl -844.65 -846.66 2.02LiBr -795.60 -796.48 0.87LiI -720.01 -729.20 9.18NaCl -777.73 -762.51 -15.22Table 4.2: Lattice energies (TF model) for the wurtzite (Ewurtz) and rock salt (Erock)structures. The difference Erock- Ewurtzis reported as ∆E. The NaCl lattice energiesare added for comparison.that as the anion becomes larger the repulsive interactions become more important,and the minima become less negative for both rock salt and wurtzite structures. Theeffect is larger for rock salt, hence the crossover to the wurtzite structure at LiBr.The apparent reason for this is the fact that for rock salt the minima for E+ (bluelines) and E− (green lines) are found at very different a values, or, in other words,the most favorable cell length for the lattice-positive-ion interaction is different fromthat of the lattice-negative-ion interaction. A similar effect occurs, but is much lesspronounced for the wurtzite structure, hence wurtzite becomes more stable thanrock salt as the ion size asymmetry becomes larger. For the TF model (right panel,Figure 4.4), the wurtzite is the more stable structure for all four salts.80Figure 4.4: Lattice energies as a function of the cell length a for the JC (left column)and TF (right column) potentials. The solid red lines indicate the energy of thewurtzite structure, while the dashed red lines indicate the rock salt structure. Theenergy minima are annotated below each curve. Blue and green lines represent thecontributions of the positive and negative ion sites, respectively.81Compound JC TFLiF rock salt wurtziteLiCl rock salt wurtziteLiBr wurtzite wurtziteLiI wurtzite wurtziteTable 4.3: Final structures obtained by melting and freezing periodic crystals. Thefinal temperature is 300 K, and the periodic box contained 1000 ions. The resultsare consistent with the lattice energy calculations (Tables 4.1 and 4.2).4.3.2 Infinitely Periodic System SimulationsThe lattice energy calculations given above determine the stable crystal struc-tures at T = 0 K. It remains interesting to ask if these are in fact the structuresthat crystallize spontaneously from the melt at finite temperature. To answer thisquestion we carried out MD simulations, first heating a rock salt crystal to a tem-perature above the melting point, then cooling the molten salt until spontaneouscrystallization occurs, following the protocol described in Section 4.2. The resultsobtained for both JC and TF models are summarized in Table 4.3, and we see thatin all cases the salt crystallized into the structure predicted by the lattice energycalculations.The lattice structures obtained upon spontaneous crystallization can be iden-tified by inspecting configurational snapshots, or, more quantitatively, by comparingradial distribution functions with expected results for rock salt and wurtzite crystalstructures. As an example, various rdfs for LiCl/JC and LiCl/TF are shown inFigure 4.2. Comparing the rdfs, it is apparent that LiCl/JC crystallizes as rock saltand LiCl/TF as wurtzite, consistent with the lattice energy predictions.In addition to identifying crystal structures, one can make some other ob-servations based on the rdfs in Figure 4.2. For the initial rock salt structure both82models give qualitatively similar rdfs (blue curves), as expected, but the peaks tendto be lower and broader for the TF model. We attribute this to the softer short-range repulsion of the TF potential, which allows the ions to have larger oscillationsaround their equilibrium positions. We see also that the molten salts give verysimilar rdfs for both models.4.3.3 Finite Size Cluster SimulationsWe have also carried out melting/freezing simulations on finite size clustersof 216, 512, and 1000 ions. This is of interest in the context of crystal nucleation,for example, from aqueous solution. It is important to know if one can expect smallclusters to nucleate in the stable structure of the bulk phase, or in some other size-dependent structure. The final structures obtained for finite clusters are summarizedin Table 4.4.We see from Table 4.4 that for salts where the rock salt and wurtzite latticeenergies differ by ∼ 2 kJ mol−1 or more, the finite cluster structures are consistentwith the infinite system results discussed above, and with predictions based on thelattice energies. One important difference is that for cases where rock salt is thestable structure, such as LiCl/JC and LiF/JC, even small clusters achieve well-defined rock salt structures, whereas this is not true of salts that prefer wurtzite.Rather, in the wurtzite case, one obtains structures that have a hexagonal motifcharacteristic of wurtzite, but are by no means perfect wurtzite structures, and inthe small cluster case can be mostly hollow inside. We refer to these structures aswurtzite-like.An example structure obtained for LiCl/TF with 216 ions is shown in Fig-83Compound Ions JC TFLiF 216 rock salt wurtzite-like512 rock salt wurtzite-like1000 rock salt wurtzite-likeLiCl 216 rock salt wurtzite-like512 rock salt wurtzite-like1000 rock salt wurtzite-likeLiBr 216 rock salt wurtzite-like512 mixed wurtzite-like1000 wurtzite-like wurtzite-likeLiI 216 wurtzite-like wurtzite-like512 wurtzite-like wurtzite-like1000 wurtzite-like wurtzite-likeTable 4.4: Final structures obtained by melting and freezing finite size clusters usingthe JC and TF models, the final temperature is 300 K. A wurtzite-like structurecorresponds to a structure that is not perfectly crystalline but shows hexagonalmotifs similar to that shown in the left panel of Figure 4.5.ure 4.5 (left panel). We note that these hexagonal structures are similar to thosereported in previous cluster simulations of LiCl employing the TF potential.121 Infact, based on considering finite and infinite hexagonal structures, Rodrigues andSilva Fernandes119 suggested that the TF potential does not produce the experi-mental rock salt structure for LiCl, and our lattice energy and simulations confirmthis conclusion.It is also interesting to compare the rdf of a perfect wurtzite and a finite sizewurtzite-like crystal. In Figure 4.6, the Li-I rdf for the JC model is shown for botha perfect, infinite, wurtzite crystal, and a finite size crystal of 1000 ions, both at 300K. To avoid the influence of surface ions, the rdf is computed by using only the Li+ion closest to the cluster’s geometric centre. One notices that the rdf peaks for thefinite size structure are broad and, with the exception of the contact peak, by nomeans a perfect match for the perfect wurtzite structure.84Figure 4.5: Spontaneously crystallized LiCl/TF (left panel) and LiCl/JC (rightpanel) clusters of 216 ions. For the cluster on the left, the structure has hexagonalmotifs characteristic of the wurtzite structure, but is mostly hollow inside. On theright, the cluster crystallizes into a rock salt structure, without substantial surfacedeformations. In this figure, Li+ is represented in purple and Cl− is in green.The reason why small wurtzite crystals are not stable is because the wurtziteunit cell has a non-zero dipole moment directed along the c axis. For a unit cell witha = 0.45 the dipole moment is ∼ 26 D, therefore a perfect wurtzite crystal will havea surface dipole with charges of the same sign occurring at the surface. This meansthat the surface configuration of a finite wurtzite crystal is energetically unfavorable,and will undergo surface reconstruction. For small clusters, one observes ratherhollow structures (most of the cluster is surface) with hexagonal patterns on thesurface. This effect is confirmed by noticing that the dipole moments per ion ofthe finite 216 ion crystals obtained in our simulations are . 0.1 D, which is muchsmaller than the value ∼ 6.5 D expected for a perfect wurtzite crystal.Surface reconstruction of a larger (2048 ion) wurtzite crystal of LiI/JC isillustrated in Figure 4.7. Simulated in vacuum at 300 K, the reconstruction of the(0001) face is apparent in the before and after snapshots shown in Figure 4.7. We85Figure 4.6: LiI/JC rdfs at 300 K for the perfect wurtzite crystal (red), and for theinterior of a finite size crystal of 1000 ions (blue). The peaks are well defined forthe perfect wurtzite crystal, while for the finite size crystal there is wide broadeningand merging of adjacent peaks.86Figure 4.7: Example of surface reconstruction for a LiI crystal (Li+ and I− ions arecolored in purple and grey, respectively). On the left, surface charges are presenton the wurtzite crystal face (0001). After simulating the crystal in vacuum at 300K, the resulting surface (on the right) rearranges to reduce the dipole moment. Thedipole moment per ion is reduced from ∼ 6.3 D for the structure on the left, to∼ 0.7 D for the structure on the right.87Figure 4.8: Spontaneously crystallized LiBr/JC cluster (1000 ions) using theLiBr/JC model (Li+ in purple, and Br− in brown). One notices both hexagonalmotifs, characteristic of the wurtzite structure, as well as square motifs characteris-tic of the rock salt structure.note that the reconstruction reduces the dipole moment of the crystal from ∼ 6.3D per ion initially to ∼ 0.7 D after simulated reconstruction at 300 K.An interesting situation occurs when the rock salt and wurtzite lattice en-ergies are very close together, such as LiBr/JC. The rock salt structure crystallizesspontaneously at small sizes but, as size increases, it crosses over to wurtzite, givingrise to mixed crystal structures, as shown in Figure 4.8. This confirms that therock salt structure has a lower surface tension (due to the lack of surface charges)than the wurtzite structure and crystallizes more easily. Also, since the standardprocedure we used involved ∼20 ns annealing time at high temperature, (see Section4.2 for details) we can exclude that this is due to dynamical effects.88Temp (K) Er Ew Ew† ∆E (kJ mol−1)0 -894.53 -889.66 n.a. -4.86300 -886.34 -881.70 -878.71 -4.65Table 4.5: LiY lattice energies at 0 K and 300 K. The value Ew was obtained bysimulating a box to initialized from the perfect wurtzite structure, while Ew† is theenergy of the crystal obtained from freezing of the molten salt.4.3.4 Temperature EffectsBy constructing a fictitious system, we studied how temperature favors oneor the other crystal structure. In our lithium halide set, we obtained cases where thewurtzite and rock salt lattice energies are spaced well apart and, for LiBr/JC, wherethe wurtzite structure is marginally more stable than the rock salt structure, butno system was found to have a marginally stable rock salt structure (the compoundclosest to this condition is LiCl/JC, where rock salt is ∼ 6.82 kJ mol−1 more stablethan wurtzite).A fictitious LiY/JC system, based on the LiBr/JC parameter set with asmaller σBr (in our case, reduced by 0.04 nm) was constructed to have a rock saltcrystal structure that is slightly more stable than the wurtzite structure. Accordingto lattice energy calculations, the minimum energies for the rock salt and wurtziteconfigurations for the LiY/JC were found to be, respectively, −894.53 and −889.66kJ mol−1, with the rock salt structure being more stable by ∼ 4.87 kJ mol−1.However, the simulations of crystallization from the melt at 300 K show that,despite its higher energy (Table 4.5), the preferred structure is wurtzite, suggestingthat the wurtzite structure is favored by entropy.894.4 Summary and ConclusionsIn this Chapter, we investigated the relative stability of the wurtzite androck salt crystal structures for a series of lithium halides using the popular Joung-Cheatham1 and Tosi-Fumi125 potentials.While wurtzite crystals are metastable in real systems, the TF model wasfound to crystallize preferentially into the wurtzite structure. In contrast, the JCmodel correctly predicts the rock salt crystal structure for LiCl and LiF but not forLiBr and LiI. The origin of this effect can be attributed to the short-range part ofthe pair potential and the excessive destabilization of the rock salt structure as ionsize asymmetry increases.Our simulations of finite size clusters in vacuum suggest that the wurtzitestructure is hard to form at small sizes because of the surface rearragements neededto avoid surface charge formation. We also found that if the lattice energies ofthe rock salt and wurtzite structures are close to each other, like in the case ofthe LiBr/JC model, the system may crystallize as rock salt at small sizes, and aswurtzite at larger sizes.The behavior of molten salts is interesting in its own right, and an under-standing of finite size cluster structures can be related to homogeneous nucleationfrom aqueous solution. Further work is required to show if the finite size crys-tal structure could prevent the stable bulk structure from forming or if there arestructure rearrangements during the nucleation process.90Chapter 5Crystallization of LithiumFluoride in Aqueous solution fromMolecular Simulation.5.1 OverviewMolecular simulation is employed to study nucleation and growth of lithiumfluoride crystals in aqueous solution. It was possible to observe nucleation on sim-ulation timescales only at a temperature of 500 K, and a fixed density of ∼ 1.24g/cm3 to prevent water evaporation. We find that the growth rate is temperaturedependent, and follows Arrhenius’ law, with an activation energy of ∼ 50 kJ mol−1.Since the activation energy for the diffusion of Li+ and F− ions in solution was foundto be much smaller (∼ 20 kJ mol−1) the growth barrier is likely related to the highenergy required to remove water from the first solvation shell of the ions. The solu-bility for the model salt was estimated to be ∼ 0.03 mole fraction at 300 K, which91is much higher than the experimental value (∼ 0.001 mole fraction). Additionally,a decrease of solubility with temperature was observed, likely due to a decrease inthe dielectric constant of the solvent.5.2 IntroductionThe study of alkali halide phase transitions with molecular simulation hasbeen an intense subject of study because of the simplicity of the compounds, andthe availablity of models designed to reproduce both properties of their crystals andaqueous solutions.1 Particular attention has been devoted to NaCl models,1,125,139for which exist nucleation studies,47,49,66,140 as well as estimates of their solubility.60Lithium halides crystallization from the melt has also been examined,121,127 but thestudy of their nucleation from solution is very limited.141In this chapter we investigate LiF nucleation from aqueous solution usingmolecular simulation. LiF, with an experimental solubility142 of ∼ 0.135 g/100 mL(∼0.001 mole fraction) at 298 K and 1 bar, is the least soluble of the lithium halides,and molecular dynamics simulations of its nucleation process provides insights intohow the nucleation mechanism may differ from that of a more soluble salt such asNaCl. Also, LiF does not introduce the complication of stable hydrate structuresthat are common for compounds such as LiCl.143We find that, within simulation timescales, nucleation from solution can onlybe obtained by heating the solution to ∼ 500 K, and that, for LiF, crystal growth isan activated process with an activation energy of ∼ 50 kJ mol−1. In the remainderof this Chapter we explore the origin of this temperature dependence by examiningproperties of the LiF solution such as diffusion, solubility, as well as the distribution92σ (nm) (kJ mol−1)Li+ 0.1409 1.4088967F− 0.4022 0.0309637Table 5.1: Joung-Cheatham Lennard-Jones parameters for LiF1 adapted for theSPC/E water model.2and lifetime of LiF clusters in solution.5.3 Models and MethodsThe model used to simulate LiF aqueous solutions are the Joung-Cheatham1parameters for the ions, combined with the SPC/E2 parameters for water. In thismodel the pairwise interactions are represented by a Lennard-Jones term plus anelectrostatic contribution (cf. Equation (3.1)). The Joung-Cheatham parametersfor Li+ and F− are given in Table 5.1.All simulations employed the GROMACS version 4.5.5 software package.100To keep the temperature fixed in our simulation we used the velocity rescalingthermostat,79 and to keep pressure fixed we used the Parrinello-Rahman barostat.137The time step used in our simulations was 2 fs. The short-range interactions werecalculate up to a cutoff of 0.9 nm while the electrostatic interactions were calculatedusing the particle mesh Ewald method.74A series of LiF solutions at different concentrations were initialized by ran-domly placing water molecules, together with Li+ and F− ions, in a cubic simulationcell of length ∼ 4 nm. The system was relaxed at a fixed temperature of 300 K anda pressure of 1 bar until its density, the Li+-F− radial distribution function, and theion cluster distribution (see below for explanation) were stabilized, which required∼ 20 ns.93NucleationT (K) Length (ns) ttn (ns) LiF H2O m.f. Pressure (bar)300 60.60 n/a 419 7977 0.05 1300 134.12 n/a 878 7903 0.10 1300 147.20 n/a 1677 7640 0.18 1300 200.00 n/a 2095 7430 0.22 1300 158.56 n/a 1369 7761 0.15 1350 199.74 n/a 1500400 102.12 n/a 3100450 136.82 n/a 4800500 132.84 15 6500500 80.00 5500 80.00 3GrowthT (K) Length (ns) LiF H2O m.f. Pressure (bar)300 80.00 1369 7761 0.15 1350 80.00 1500400 80.00 3100450 249.30 4800-6000500 24.64 6500-7200Table 5.2: Summary of the simulations performed to study LiF nucleation andgrowth. The temperature at which a simulation was conducted is indicated byT, the length of the simulation is indicated by Length, the concentration in molefraction is abbreviated as m.f., and the time to nucleation (ttn) is the time at whichthe first nucleus was observed. LiF and H2O indicate the number of moleculesof each type used in the simulation. The pressure values for the crystal growthsimulations are reported as a range for some simulations because, when substantialcrystal growth occurs, the pressure drifts to higher values.94Starting from the ∼ 0.15 mole fraction LiF equilibrated solution (which cor-responds to ∼ 5 to 15 times the saturation concentration for our model, dependingon the temperature), we simulated the system over a range of temperatures (be-tween 300 and 500 K) and observed crystal nucleation at a temperature of 500 K.In all simulations the volume was kept constant (with density ∼ 1.24 g/cm3) tomaintain the solution in the liquid state. Three replica simulations at 500 K werealso performed to have an estimate of the time necessary to achieve the first nucleus,the replicas showed a substantial variability resulting in a time to first nucleationranging from 3 to 15 ns. Lower concentrations were also tested but we were notable to obtain nucleation at any temperatures (between 300 and 500 K), concentra-tions higher than ∼ 0.15 mole fraction immediately yielded multiple nucleations ata temperature of 500 K. A summary of the simulations performed is given in theleft panel of Table 5.2. It is useful to notice that the experimental solubility for LiFat 300 K is ∼ 0.001 mole fraction,142 which is an order of magnitude lower than ourestimate for the model employed (see next section).To study crystal growth of a LiF crystal in solution, a configuration con-taining a partially grown crystal made of ∼ 300 ions was extracted from one of thesimulations at 500 K and 0.15 mole fraction. Beginning with this configuration,NV T simulations were performed at the temperatures 300, 350, 400, 450 and 500K (cf. Table 5.2).The amount of crystal phase present in the growth simulations was estimatedby calculating the Steinardth-Nelson bond order parameter107 q8 for each ion andby counting the number of ions with q8 > 0.40 (cf. Chapter 3 for further details).95The order parameter ql is defined asql =[4pi2l + 1l∑m=−l|qlm|2]1/2,whereqlm =1NrN∑ri=r1Y ml (θ(ri), φ(ri)),where {r1 . . . rN} are the positions of the first N neighboring ions relative to thecentral ion, Y ml are spherical harmonics, and θ and φ are the polar and azimuthalangles corresponding to the vector r. The method and parameters used are thesame as those described in Chapter 3 (which correspond to parameters N = 12 andl = 8).Another clustering method, equivalent to that employed by Hassan,58 wasused to study dissolution and the properties of the solution, including informationabout pairs and other ion aggregates that do not arrange with an octahedral geom-etry. We will refer to this method as connectivity based clustering.In this method, given a system configuration, all the pairwise distances be-tween positive and negative ions are first calculated. Every pair of oppositely chargedions that are closer than a certain threshold are deemed to be in contact and placedin the same cluster. The natural threshold for this procedure is a number slightlylarger than the contact radius (which can be found from the first peak of the Li-Fradial distribution function), and it was chosen to be 0.3 nm for LiF. The final resultof the above procedure is a set of clusters where every ion can be reached from everyother ion through a path of consecutive connections. Notice also that free ions areconsidered to be in a cluster of their own. See Figure 5.1 for a diagram.96Figure 5.1: Clustering based on neighbor distance. The dashed lines representdistances lower than a fixed threshold, and the colors represent different clusters.Free particles form a cluster of their own.5.4 Results and DiscussionA series of LiF solutions at different concentrations and temperatures weresimulated at constant volume and it was found that spontaneous, homogeneousnucleation within ∼ 3 to 15 ns can be obtained at a high temperature (500 K) andat a concentration ∼ 0.15 mole fraction and higher. Simulations of crystal growthrate also show an increase with temperature and, as we will show below, the growthrate dependency on temperature follows an Arrhenius law.The Arrhenian temperature dependency of the crystal growth simulationssuggests that there is an activation barrier to crystal growth. The origin of thebarrier could be due to the very slow diffusion of ions towards the surface of thecrystal. Another hypothesis is that solvation shells are strongly attached to the ions,and there is a high barrier to displace water molecules for ion attachment to occur.Lastly, a temperature-dependent crystal growth rate could also be explained by an97Figure 5.2: Schematic diagram of the free energy landscape at different temperaturesfor the crystal growth process. If the driving force (which is the difference in freeenergy of initial and final stages) increases and the barrier is weakly affected bythe increase in temperature, the activation barrier (which is the height of the curvebetween inital and final stages) of the process decreases.increase in the driving force at high temperature. If the solution is less stable athigh temperature, while the energy of the “transition” state may be weakly affected,the activation barrier might be reduced (see Figure 5.2 for a diagram).Since we did not observe crystal nucleation at different temperatures, simi-lar statements are not directly applicable to the nucleation process. However, it isreasonable to assume that nucleation is affected by growth in the sense that nucle-ation cannot be observed without crystal growth. In this sense, we would expect theactivation energy for growth to constitute a lower bound for the crystal nucleationactivation energy.5.4.1 DiffusionOne might expect diffusion to influence the speed at which crystal growthoccurs by decreasing the overall mobility of the ions in solution. The diffusionconstants for the Li+ and F− ions were estimated by calculating the mean square98T (K) D+ (10−5cm2/s) D− (10−5cm2/s)300 0.0190 (± 0.0007) 0.0205 (± 0.0002)350 0.1348 (± 0.0050) 0.1350 (± 0.0106)400 0.3842 (± 0.0174) 0.4216 (± 0.0133)450 0.7918 (± 0.0313) 0.8329 (± 0.0115)500 1.3514 (± 0.0239) 1.4928 (± 0.0089)Table 5.3: Diffusion constants for Li+ (D+) and F− (D−) at different temperatures.displacement between 1 and 9 ns of our simulations, while making sure that nonucleation is taking place, and the results are reported in Table 5.3.We calculated the diffusion coefficient at different temperatures to estimatethe activation energy of the process using the Arrhenius equation:ln(D) = − EaRT+ ln(A), (5.1)where D is the diffusion coefficient, Ea is the activation energy, R is the gas constant,and A is the preexponential factor. The activation energy for diffusion can beestimated by performing a linear regression of ln(D) versus 1T(see Figure 5.3 foran example), the resulting Ea for both Li+ and F− diffusion was found to be ∼ 26kJ mol−1. Interestingly, this value is less than half the activation energy for crystalgrowth (see below), suggesting that diffusion is not a controlling factor in the growthprocess.5.4.2 Growth BarrierBy measuring the growth rate at different temperatures it is possible toestimate an activation energy for the growth process, similarly to the estimation ofthe activation energy for diffusion.99Figure 5.3: Linear regression to determine the activation energy for the diffusionprocess for Li+ ions.The growth simulations summarized in Table 5.2, were initialized by extract-ing a configuration containing a small cluster of ∼ 300 ions from the simulation at500 K and ∼ 0.15 mole fraction. This seed was surrounded by the supersaturatedsolution, simulations were carried out at the temperatures 300, 350, 400, 450 and500 K, and the growth profiles were monitored using a threshold on the q8 orderparameter as described in Chapter 5.3. The initial growth rate, rg (over ∼ 1 ns)was obtained for each temperature (Table 5.4) and fitted to the Arrhenius equation(Figure 5.4). The activation energy was estimated to be ∼ 50 kJ mol−1.If we compare the activation energy for growth with the activation energyfor diffusion, it is easy to see that the activation energy for growth is much slowerthan that of diffusion, suggesting that the growth rate is not diffusion limited. Theorigin of such high activation barrier probably lies in the solvent-ion interaction, for100T(K) rg (ion ns−1)300 -0.05350 0.23400 1.66450 11.41500 32.40Table 5.4: Growth rates at different temperaturesFigure 5.4: Arrhenius fit to the growth ratesthis reason we further explored the structure and dynamics of the supersaturatedsolution.5.4.3 Solution Structure and DynamicsOne important property of the solution is the degree of association of the ions.To study ion association, we performed connectivity based clustering, as describedin Chapter 5.3. Using this method, every ion is assigned to a cluster, and it is101Figure 5.5: Fraction of ions in connectivity based clusters of a certain size at differenttemperatures.possible to study how the distribution of free ions, ion pairs, ion triples, and higherorder aggregates varies at different temperatures.The fraction of ions that are part of a cluster of a certain size at differenttemperatures is shown in Figure 5.5. At 300 K, about 40 % of the ions are notassociated (in clusters of size 1), about 35% of the ions are part of clusters of size2, and the remaining 15% are part of larger clusters up to ∼ 16 ions. As thetemperature increases, the distribution gradually shifts towards clusters of highersize and, at 500 K only ∼ 15% of the ions are not associated. This shows that astemperature increases the ions have a tendency to associate into larger aggregates,thus favoring nucleation.Further evidence can be found from the Li+F− radial distribution functions ofthe solution at different temperatures shown in Figure 5.6. At higher temperatures,102Figure 5.6: Li+F− radial distribution functions at different temperatures. The con-tact peak, characteristic of ion pairs increases with temperature, while the secondpeak (solvent-separated pairs) decreases.the LiF contact peak increases substantially, while the second peak, correspondingto solvent-separated pairs, decreases in intensity. This result is in agreement withthe ion association plots (Figure 5.5) as more ion association increases the numberof contact pairs.An important point is that the ion association distributions shown in Figure5.5 describe the metastable “equilibrium” of the solution at different temperatures.This means that, before nucleation happens, the system spontaneously relaxes tohave a certain ion association pattern. This effect, as well as the change in solu-bility, is likely to be due to the changing dielectric constant of water at differenttemperatures and pressures.Ion association distributions, however, do not describe the dynamics of thesolution. To assess how clusters exchange ions, it is useful to introduce a newmeasure that we will refer to as surviving pairs.At a given time frame, the ions are associated in a number of clusters. Twodifferent frames may (and usually) have a different number of clusters and, in orderto define a lifetime, it is necessary to match the clusters in two different frames103Figure 5.7: Diagram of surviving pairs. Initially, the {1, 2, 3, 4} cluster contains 6unique pairs, after the cluster splits, 2 of the pairs are still associated.(this process is usually called entity resolution and is described in detail in Chapter3). The matching can be quite problematic when there is substantial rearrangingbetween the frames, as there is a large redistribution of ions between the clusters.Another approach to the problem is to count the number of pairs that areassociated in a certain time frame, and count the number of pairs that are stillassociated at a different time frame. As shown in Figure 5.7, the surviving pairs isperfectly defined even when clusters split and, while it doesn’t provide informationabout the lifetime of individual clusters, can be used as a global measure of howlong lived are ion clusters in solution.By plotting the surviving pairs as a function of time delay between frames itis possible to obtain a plot that indicates how many pairs are still associated up totime t (Figure 5.8). As temperature increases, the initial number of pairs increases(corresponding to time delay t = 0 ) as there is a higher degree of association,however, their decay with t is much more pronounced. The 300 K surviving pairscurve shows that, at this temperature, the associated structures are very stable, anddecay slowly over the course of ∼ 20− 30 nanoseconds (the range is too wide to be104Figure 5.8: Surviving pairs at a time delay of t ns for a range of different tempera-tures. As the temperature increases, there are more associated pairs that decay ata faster rate.shown on the plot). In contrast, temperatures higher than 400 K show that after∼ 1 ns most of the pairs have decayed, suggesting fast rearrangements of the ionaggregates present in solution.A cluster analysis similar to the one presented in Chapter 3 was performed forthe LiF supersaturated solution. The number of clusters detected per frame, usingthe methodology described in Chapter 3, is reported in Table 5.5. An importantfinding is that, at low temperatures, no clusters are detected using this method. Themain reason for this effect is that the aggregates that form in the LiF supersaturatedsolution do not possess high enough crystallinity to pass through the filtering step(cf. Chapter 3). At higher temperatures, there is a sharp increase in number ofclusters that, presumably, are a prerequisite for crystal nucleation. In contrast, the105Compound T (K) nc × 105 m.f.LiF 300 0 0.15350 0400 0.15450 1.33500 27.65NaCl 300 142.33 0.22Table 5.5: Average number of clusters (nc) per time frame at different temperatures,for LiF and NaCl supersaturated solutions (see Chapter 3 for details). nc is nor-malized by the number of ion pairs in the simulations to allow comparisons betweensimulations with different numbers of ions. The concentration in mole fraction isreported in the m.f. column.NaCl solution, as discussed in Chapter 3, is composed of a large quantity of clusterswith longer lifetimes.Lifetime analysis was not performed on the cluster data produced from theLiF simulations because at high temperatures the time interval at which configura-tions were sampled (0.02 ns) is large compared to the time it takes a typical clusterto redissolve, hence it was not possible to follow the time evolution of most clusters.5.4.4 SolubilityIf the solubility decreases substantially at high temperature, one would ex-pect nucleation and growth rates to increase due to a higher driving force (since oursimulations are all at the same concentration, the solubility determines the degree ofsupersaturation). Experimentally, the solubility of LiF increases with temperature(at least in the range ∼ 300 to 350 K at 1 atm).142 However, it is still possible forour LiF model to show a decrease in solubility with temperature due to changes inthe solvent, for example the dielectric constant.106We estimated the solubility of LiF at 300 K by dissolving a crystal untildissolution stops, and at 500 K by both dissolving a crystal and, by growing a nucleusin a supersaturated solution until growth stops. More specifically, the dissolutionand growth were assumed to be completed when there was no change in the solutionconcentration which occurs after ∼ 100 to 200 ns.The solution concentration was calculated by detecting connectivity basedclusters (see Chapter 5.3) and counting the number of ions that are part of clusterswith size less than 100 ions. As the simulation contains a single cluster of size∼ 2000 ions, and transient clusters are of size definitely smaller than 100 ions (seeFigure 5.5), these numbers represent a good estimate of the ions present in solution.The dissolution and growth profiles are shown in Figure 5.9.It is useful to notice that the starting point of the dissolution at 300 K andat 500 K is larger than 0 ns. This is due to the fact that, since the initial crystalwas extracted from a particular configuration, it possessed a small net charge dueto ions fluctuating in proximity of the interface. To avoid net charges, extra ionswere added to the solution to retain the electroneutrality. Geometry optimizationwas also required and produced slightly different starting condition for the twotemperatures. The fluctuations observed for the dissolution simulation at 500 Kshow that the initial conditions are quite close to equilibrium; in fact, after aninitial rapid dissolution, the cluster grows again until the solution stabilizes arounda concentration of ∼ 0.01 mole fraction.The solubilities estimated using these methods yielded values in the range of∼ 0.01 to 0.03 mole fraction, which is much larger than the experimental solubilityof ∼ 0.001 mole fraction. Interestingly, at high temperature there is a decrease insolubility and this could be attributed to a change in solvent nature. To test this107Figure 5.9: Solution concentration profiles obtained by dissolving a crystal at 300K (blue), by growing a crystal at 500 K (red, top line), and dissolving a crystal at500 K (red, bottom line). The estimated solubilities are ∼ 0.025 mole fraction at300 K and ∼ 0.01 at 500 K.hypothesis, we measured the dielectric constant of SPC/E water to be ∼ 71 at 300K and ∼ 40 at 500 K. The decrease in dielectric constant of the solvent is in linewith the decrease in solubility at high temperature.The degree of supersaturation can certainly influence crystal growth rate asit is the driving force of the process. However, as it was shown above, the crystalgrowth rate shows an Arrhenian dependency on temperature, and is not likely due toa change in solubility, because we would expect the associated change in activationenergy to yield a nonlinear Arrhenius plot.Similar effects can also be expected for the nucleation rate; in fact, increasingthe degree of supersaturation by pushing the concentration to ∼ 22 mole fraction108didn’t result in nucleation at 300 K, showing that an increase in driving force is notenough to trigger nucleation at that temperature.5.5 Summary and ConclusionsIn this Chapter we focused on the nucleation and growth of LiF crystalsfrom solution. Our results suggest that ions in the supersaturated solution arestrongly solvated and increasing the temperature helps ions to overcome the barrierto dehydration, favoring growth.While we could not directly measure nucleation rate at different tempera-tures, we were able to measure growth rate of an already nucleated crystal, and wefound that the temperature dependence of the growth rate follows the Arrheniuslaw, indicating that the process is activated and has an activation energy of ∼ 50 kJmol−1. Such a high activation energy can not be explained by ion diffusion, whichhas a much smaller activation energy (∼ 20 kJ mol−1).At high temperature, the LiF solution contains a large quantity of ion clustersthat form and dissolve very quickly but, despite the fast decay of these structures,the solution is able to achieve nucleation. This situation is in contrast with whatwas observed in Chapter 3, where the NaCl solution was found to have many longlived clusters that failed to achieve the size and crystallinity for nucleation to occur.A mechanism consistent with our results for LiF is a two-step nucleationprocess.144 The first step involves cluster generation from solution and, in the caseon LiF, clusters appear at high temperature and pressure, where they are easier toform because dehydration is fast enough to achieve kinetically, and ion clusters arealso favored because the dielectric constant of the water decreases. Once a sufficient109Figure 5.10: Possible free energy landscape for LiF nucleation from solution. A firstbarrier is required to go from the free ions (solution) to the cluster aggregates. Thecluster aggregates are metastable, continuously form and disappear from solution,and are more easily found at high temperature. The second barrier reflects the factthat most clusters do not nucleate but dissolve back into solution, depending oncluster related properties such as size, crystallinity, and surface tension. In absenceof nucleation rate measurements, it is not possible to estimate the relative magnitudeof the two barriers, and two possible options are displayed as solid and dashed lines.number of clusters are generated from the saturated solution, there is another bar-rier, similar to that postulated by CNT, that is likely related to cluster propertiessuch as size and crystallinity as we found for NaCl (Chapter 3). Unfortunately, thedata we collected is not sufficient to determine the magnitude of the second barrier,as it requires precise measurement of nucleation rates sampled at different temper-atures. A diagram of possible free energy landscapes of the LiF nucleation processis shown in Figure 5.10.The solubility of LiF was estimated to be ∼ 0.03 mole fraction at 300 K, andwas found to decrease at high temperature to reach a value of ∼ 0.01. The variationis likely due to the decreasing dielectric constant of the solvent. The calculatedvalue of the solubility is substantially different from the experimental solubility of∼ 0.001 mole fraction, which likely indicates that, in this respect, the JC model isnot accurate for LiF, and further model improvements are needed.110Chapter 6Molecular Dynamics Simulation ofNaCl Dissolution.6.1 OverviewMolecular dynamics simulations are used to investigate the dissolution ofNaCl nanocrystals (containing ∼ 2400 ions) in water. We focus on systems undersink conditions at 300 K, but the influences of concentration and temperature arealso investigated. Cubical, spherical, tablet-shaped, and rod-shaped nanocrystalsare considered, and it is shown that the initial shape can influence the dissolutionprocess. Dissolution is observed to occur in three stages: an initial period where themost exposed ions are removed from the crystal surface, and the crystal takes on asolution-annealed shape which persists throughout the second stage of dissolution;a second long intermediate stage where dissolution roughly follows a fixed rate law;and a final stage where the small residual crystal (. 200 ions) dissolves at an everincreasing rate until it disappears. The second stage of dissolution which applies for111most of the dissolution process is well described by classical rate equations whichsimply assume that the dissolution rate is proportional to an active surface areafrom which ions are most easily detached from the crystal. The active area dependson the initial crystal shape. We show that for our model NaCl nanocrystals the ratedetermining step for dissolution under sink conditions is ion detachment from thecrystal, and that diffusion layers do not exist for these systems.6.2 IntroductionSolid dissolution is an important process in many physical systems and sit-uations, with one obvious example being in the area of drug development.145,146Many drugs are administered in solid form and dissolve in the gastrointestinal tractbefore being absorbed by the body. As the dissolution properties are related tothe drug bioavailablilty,147,148 there is considerable current interest in dissolutionprocesses.149During the past century, several dissolution models based on different as-sumptions have been suggested147,149,150 to describe and interpret dissolution rates,and identify the factors that affect them. The original model was put forward byNoyes and Whitney,151 and further generalized by Brunner and Tolloczko152 to ac-count for varying solid surface area during particle dissolution. For spherical solidsamples, the integrated form of the Brunner-Tolloczko equation is the well-knowncube root law, whereby the cube root of the particle weight decreases linearly in timeduring dissolution, first derived by Hixson and Crowell.153 Another early model,based on the assumption that the dissolution rate is controlled by solute diffusionfrom a concentrated layer of solution adjacent to the solid surface, was suggested by112Nernst154 and Brunner.155 Under some conditions the Nernst-Brunner equation alsoleads to the cube root law. Additionally, square root156 and two-thirds root156,157dissolution laws have been obtained and used in dissolution research. These modelsare all based on theoretical assumptions or empirical observation, depending on thesubstance being dissolved and its environment; however, little is known from a mi-croscopic, mechanistic point of view. Furthermore, the classical models do not takeany account of shape and finite size effects (other than surface area) that one mightexpect to be of some importance in nanoscale crystals.Thanks to large increases in computational power, molecular simulation hasstarted to gain traction as a tool in dissolution research. Studies have been con-ducted to uncover the mechanics of the first steps of dissolution and growth of smallcrystalline compounds.158–160 Gao and Olsen159 have reported simulations of thedrug acetaminophen in water, giving insights into the initial stages of the dissolu-tion process. Several studies have also been done to analyze the interactions thattake place at the interface between water and NaCl.161–165In the present paper, we investigate the complete dissolution of NaCl nanopar-ticles in water employing molecular dynamics simulations. NaCl was chosen becauseof its simple structure and because of the availability of force fields developed to re-produce both solid and solution phase properties.1 Our results provide an atomisticview of the dissolution process at every stage, from initiation, through a fixed ratelaw regime, to the finial disintegration of the residual crystal. Additionally, we ex-amine the influence of particle shape on dissolution, and compare our results withsuitable, shape-adapted, classical dissolution models. The influences of solutionconcentration and temperature on the dissolution rate are considered. The temper-ature dependence allows us to estimate the activation energy associated with ion113σ (nm) (kJ mol−1)Na+ 0.2160 1.4754533Cl− 0.4830 0.0534924Table 6.1: Joung-Cheatnam NaCl parameters1 for the Lennard-Jones potential.detachment, which for NaCl crystals controls the dissolution rate.The remainder of this Chapter is divided into three parts. The models andmethods are given in Section 6.3, our results are described and discussed in Section6.4, and our conclusions are summarized in Section 6.5.6.3 Models and MethodsSeveral simulations of NaCl nanocrystals dissolving in water are reported. Inall simulations, the nonbonded, site-site interactions u(rij) consist of Lennard-Jones(LJ) plus electrostatic terms such thatu(rij) = 4ij[(σijrij)12−(σijrij)6]+14piε0qiqjrij, (6.1)where qj and qi are the charges on sites i and j, rij is the site-site separation, σijand ij are the LJ length and energy parameters, and ε0 is the permittivity of freespace.To model Na+ and Cl− ions we adopt the parameters of Joung and Cheat-nam1 (Table 6.1), used together with the SPC/E water model.2 The Joung andCheatnam parameters were developed to reproduce both solid and solution prop-erties, such as the solubility.1 The usual Lorenz-Berthelot combining rules, σij =(σi + σj)/2, ij =√ij, were used to calculate the LJ interactions.114Figure 6.1: A representation of the crystals simulated. The number of NaCl ionpairs is given below each illustration.In order to investigate the possible influence of particle shape (number ofcrystal layers, surface area, etc.) on dissolution rate, NaCl crystals of differentshapes were generated by repeating the NaCl unit cell along the three crystal axesdirections, to produce the following four structures (Fig. 6.1):Cubic Crystal (1372 ion pairs) generated by repeating the unit cell 7 times in allcrystal axes directions.Spherical Crystal (1256 ion pairs) generated by taking a cubic crystal of ∼ 2000ion pairs and removing all ions pairwise (such as to maintain a charge neutralcrystal) whose distance from the center is greater than 2.28 nm.Rod-Shaped Crystal (1152 ion pairs) generated by repeating the unit cell 18, 4,and 4 times along the crystal axes.Tablet-Shaped crystal (1200 ion pairs) generated by repeating the unit cell 10,10 and 3 times along the crystal axes.All dissolution simulations were carried out under NPT conditions with the pressurefixed at 1 bar. The numbers of NaCl and water molecules in each simulation,together with the temperatures considered are summarized in Table 6.2.For the present model the saturation concentration estimates reported byAragones et al.98 correspond to NaCl mole fractions in the range 0.08− 0.09. Mostof our simulations were carried out at NaCl mole fractions. 0.015, or at∼ 1/6 of the115NaCl H2O T (K)Sphere1256 97956 3001256 97956 3201256 97956 3401256 75917 3001256 48238 300Cube 1372 92838 300Rod 1152 77952 300Tablet 1200 81200 300Table 6.2: Summary of crystal shape, the number of ion pairs, the number of watermolecules, and the temperatures used in the simulations.saturation concentration. This meets the so-called “sink condition” which refers tosolutions sufficiently dilute that ion reattachments do not significantly influence thedissolution rate.166 Note that sink conditions are usually attained when the volumeof solvent is 3 − 10 times greater than the saturation volume.166 We tested thatour systems did obey sink conditions by directly monitoring ion attachment anddetachment events, as well as by running simulations at different concentrations, asdescribed below (Section 6.4).All simulations were performed employing the GROMACS100 molecular dy-namics package version 4.5.4. The temperature was controlled using the velocity-rescale thermostat79 with a relaxation time of 0.1 ps. The pressure was kept constantby applying the Berendsen barostat78 with a compressibility of 4.5×10−5 bar−1 anda relaxation time of 1.0 ps. The time step chosen for all the simulations was 2 fs. Pe-riodic boundary conditions with the usual minimum image convention were applied,and all short-range interactions were spherically truncated at 0.9 nm. The long-range electrostatic interactions were accounted for using the particle mesh Ewald(PME) method74 with a Fourier spacing of 0.20 nm. Constraints on all bonds were116maintained using the LINCS algorithm.77Initial configurations were generated by placing the NaCl crystal at the centerof the simulation cell surrounded by water molecules randomly placed on a gridspaced by 0.3 nm. The system was equilibrated for 1 to 2 ns to obtain a stablepressure, temperature, and potential energy during which time no substantial part ofthe crystal dissolved. The system was then evolved in time until complete dissolutionof the crystal was achieved. For the systems considered here, this required timesranging from ∼ 100 to ∼ 700 ns, depending on the nature of the crystal, the amountof water present, and the temperature. Note that all simulations were at 300 K,except for the spherical crystal, where two additional simulations were performedat 320 K and 340 K, and used to estimate an activation energy for the dissolutionprocess.An order parameter based on the number of neighboring ions is used to clas-sify the ions as being part of the crystal or part of the solution. Ions belonging tothe crystal have a substantially higher number of neighbors than those in solution,and this difference is sufficient to define a relatively simple order parameter. Theorder parameter for a particular ion is calculated by counting the number of neigh-boring ions (both positive and negative) within a specified radius. The actual valueof the neighbor search radius did not greatly affect the result as long as there is asubstantial difference between the number of neighbors around ions in the crystaland in the solution. In the present calculations we selected a search radius of 0.6 nm.For this radius, ions associated with the crystal have on average 20 neighbors, whilethose in solution solution have 6. Therefore, an intermediate value of 11 was chosenas the classification threshold; ions having 11 or more neighbors were classified asbeing part of the crystal, otherwise they were classified as being part of the solution.117Figure 6.2: An example of the order parameter applied to a surface ion. The centralatom (yellow) is surrounded by 14 neighbors (red). The 0.6 nm cutoff applied tothe central atom is highlighted with a gray transparent sphere.An illustration of how this order parameter detects an ion at the crystal surface fora particular configuration is given in Fig. 6.2. Using this order parameter we couldclosely follow the dissolution of each crystal independent of its shape.6.4 Results and Discussion6.4.1 Stages of DissolutionWe first examine in detail the stages of the dissolution process. Simulationswere used to obtain dissolution profiles for several model NaCl crystals. A dissolu-118tion profile is a plot of the crystal weight, or, equivalently, the number of ions in thecrystal, as a function of time. The time derivative of this profile is the dissolutionrate, or the change in crystal weight (or ion number) per unit time. The dissolutionprofile can be influenced by the solution concentration and structural properties ofthe crystal, such as size and shape. Qualitatively, we observed that the dissolu-tion process occurs in three stages that can be roughly described as follows: InitialStage: detachment of the most exposed ions located on any sharp edges. Fixed RateRegime: crystal dissolution after the edges are removed appears to closely follow afixed rate law. Rapid Dissolution: the crystal reaches a “stability limit” and quicklydisappears.Initial StageDuring the first few nanoseconds, water interacts with the most exposed ionsand these move readily from crystal to solution. The detachment of these ions leadsto defects on the surface of a geometrically shaped crystal. In the initial stage, thedissolution rate depends on the crystal shape, as this influences the arrangement ofthe most exposed ions.Dissolution profiles for initially cubic and spherical crystals are plotted in Fig.6.3, and snapshots of the cubic crystal at various points on the dissolution curveare shown in Fig. 6.4. For the cubic crystal, one observes (see the magnificationof the short-time region in Fig. 6.3) an initial steep slope between 0 and ∼ 20 ns;dissolution then quickly slows down, and enters a nearly linear regime (the actualtime dependence is discussed below). During the rapid initial stage ∼ 200 ions leavethe crystal, and as can be seen from Fig. 6.4 (snapshot at 20 ns), these come mainlyfrom the corners and edges of the crystal. From Fig. 6.4 we see that after the119Figure 6.3: Dissolution profiles for the cubical and spherical crystals displayed asnumber of ions in the crystal vs time. The starting and final sections of the profilesare detailed in the zoomed-in plots.corners and edges of the initial crystal are gone, more ions are detached from thesites generated by their removal, and water continues to gradually consume crystaledges until they are completely rounded off. At 185 ns (Fig. 6.4), the dissolvingcrystal is almost spherical in shape, and its structure and dissolution rate are verysimilar to those of the initially spherical crystal. Note from Fig. 6.3 that, after theinitial stage, the dissolution profiles of the cubic and spherical crystals are essentiallyparallel to each other. Note also, that the initial rate is much faster for the cubethan for the sphere, which has no corners and edges. Nevertheless, for the sphereas well, dissolution is faster at the beginning, as the more weakly attached ions areswept from the surface.120Figure 6.4: Snapshots corresponding to different points in the dissolution profile ofthe cubic crystal. After 20 ns, only ions located at edges and corners are removed.The edges and corners gradually get consumed (100 ns) and the crystal becomesroughly spherical after about 185 ns. After that point the shape doesn’t changeuntil the final stage of dissolution.The dissolution profiles of rod-shaped and tablet-shaped crystals also showfast initial stages followed by slower nearly linear behavior (see Fig. 6.5, and thediscussion below). The overall qualitative picture can be seen from snapshots takenfrom rod and tablet dissolution trajectories shown in Figs. 6.6 and 6.7, respectively.For the rod (Fig. 6.6), ions are preferentially removed from the ends, and for thetablet (Fig. 6.7), water peels layers from the thin sides, while none are removedfrom the interior of the larger flat surfaces. In general, as the dissolution continues,the detachments progressively round the corners of the crystals, giving both thetablet-shaped and rod-shaped crystals a more cylindrical appearance. After theinitial stage, for the rod and tablet crystals we are left with structures that resemblelong and flat cylinders, respectively. From that point on, the dissolution takes placefrom the side in the case of the flat cylinder (originally a tablet) and from the endsin the case of the long cylinder (originally a rod). These dissolution induced shapechanges are included in the dissolution rate laws discussed below.Physically, it is not difficult to understand the detachment patterns discussedabove. We would expect ion detachment to be an activated process (see below forconfirmation) and the activation energy of a particular ion will depend on its location121Figure 6.5: Comparison of the dissolution profiles for different shapes expressed interms of fraction of crystal dissolved fcry(t), the derivatives dfcry/dt (on the order of−10−3 ns−1) are depicted in the inset. The color coding is: dark blue for the sphere;green for the cube; red for the tablet; light blue for the rod. The rod-shaped andtablet-shaped crystals show a higher overall dissolution rate.Figure 6.6: Snapshots along the rod-shaped crystal trajectory. Ions are detachedmainly from the ends of the rod.122Figure 6.7: Snapshots along the tablet-shaped crystal trajectory. The top andbottom surfaces are never attacked by water molecules.in the crystal. We would expect ions at corners and edges to be more weakly boundto the crystal and hence to have a lower activation energy. Moreover, the edgesand corners are more exposed to water molecules, which should aid in getting overactivation barriers inherent in the crystal.To further examine this reasoning, we simulated a system with a perfect NaCl(001) crystal face (no defects, corners, or edges) in contact with water. The crystalface was constructed by repeating the unit cell 7, 7, and 4 times (for a total of 1578ions) along the x, y, and z axes, and placing the resulting crystal wall in a simulationcell periodic in the x and y directions. The crystal wall was in contact with 2855water molecules, with all conditions, parameters, and potentials identical to ourother dissolution simulations at 300 K. During an ∼ 200 ns simulation, not a singleion left the crystal wall, demonstrating a much lower probability of detachment forions located in large flat crystal surfaces, and highlighting again the importance ofcorners, edges etc., in the dissolution process.The picture that emerges is that the detachment of ions depends on the local123environment on the surface. To quantify this effect we performed a calculation ofthe probability of detachment (within a certain time frame) given the number ofions present (0 to 6) in the first coordination shell.The ions detached from the surface of the spherical crystal within 1.6 ns wasrecorded and averaged in a time frame between 20 and 40 ns. The distribution ofthe order parameter showed that the detached ions assume a value no larger than22 neighbors. Given that the detached ions come from the surface, this thresholdvalue was used to classify ions as part of the surface.The ions on the surface were further classified based on the number of imme-diate neighbors (within a radius of 0.37 nm), corresponding to the first coordinationshell. This step was necessary to analyse the relationship between the local envi-ronment and the probability of detachment.The total number of ions on the surface, given the local environment is re-ported in Fig. 6.8, top right corner. The number of detached ions is reported inthe bottom right corner, by calculating the fraction of detached ions for each localenvironment, it is possible to obtain the probability of detachment (reported on theleft).The probabilities of detachment show that, as the coordination number in-creases (the ion is more internal), the probability of detachment greatly decreases.However, the total number of ions detached shows that most of the ions dissolvedcome from an environment containing 3 neighbors. This is due to the fact that, evenif the probability of detachment is very low (∼ 0.08), the larger supply of those ionson the surface makes them, on the 1.6 ns time frame, the most common kind of ionto detach.124Figure 6.8: Three quantities obtained for the spherically shaped crystal as functionsof the number of neighbor in the first coordination shell (within 0.37 nm), Leftpanel: The relative probability of detachment (within 1.6 ns) of a suface ion, giventhe number of neighbors. Upper right panel: The total number of ions on the surfaceat t = 0. Lower right panel: The total number of ions detached over 1.6 ns. Allquantities are averaged over 50 time slices betwen 20 and 40 ns of the dissolutionsimulation.125Fixed Rate RegimeAfter the initial stage, the dissolution of all crystals follows a similar route.The dissolution profile is almost linear (the actual behavior is discussed in moredetail below), until the remaining crystal (∼ 200 ions) becomes very unstable andthere is a sudden and continuing increase in dissolution rate as shown in Figure 6.3.During this stage of dissolution, the crystal surfaces are well annealed bythe water, and dissolution takes place evenly over reasonably well defined “active”areas of the crystal surface. Roughly, the active areas are spherical surfaces for thecube and sphere, the cylindrical bases for the rod, and the cylindrical side for thetablet. For these active surfaces that develop in the partially dissolved crystals, theprobability of detachment appears essentially uniform such that dissolution followsa fixed rate law until the final rapid dissolution stage.We next consider more closely the influence of crystal shape, and compareour dissolution curves with classical models of the type briefly discussed in Section6.2.Influence of the crystal shape. The crystals we consider differ not only inshape, but also a little in size. Therefore, in order to compare dissolution profilesin Fig. 6.5 we plot the fraction of the crystal dissolved, fcry(t) = (N0 − N(t))/N0,where N(t) is the number of ions in the crystal at time t, and N0 is the numberpresent initially. It is apparent that for the crystals we consider the tablet and roddissolve substantially faster (∼ 6 ions ns−1) than the cube and sphere (∼ 4 ionsns−1). For the present examples, this is likely explained by the fact that ions onthe active surfaces of the tablet and rod interact more weakly with the bulk crystalthan those on the surface of the water-annealed cube or sphere. Note that in their126smallest dimension the tablet and rod have only 8 and 6 crystal layers, respectively,whereas the cube has 14 layers initially. Thus the faster rates we observe for thetablet-shaped and rod-shaped crystals, while true for the present examples, are notonly functions of shape, and likely will not apply as general rules.Comparison with dissolution models. As discussed in Section 6.2, dissolutionprocesses have been studied for well over a century, and a number of models havebeen put forward in the form of differential and integrated rate laws. Recent reviewsof these rate laws and their history are given in refs. 147 and 149. Here we findthat simple rate laws give a reasonably good description of our results in the fixedrate law regime. If we assume that sink conditions apply, and that detachment fromthe surface is the rate-determining step, we would expect the dissolution rate to begiven bydN(t)dt= −kSactive(t) , (6.2)where N(t) is the number of ions remaining in the crystal at time t, Sactive(t), is the“active surface area” where ion detachments occur, and k is a constant. This ratelaw that takes account of the changing surface area was first suggested by Brunnerand Tolloczko,152 as an extension of the original Noyes-Whitney equation151 (i.e.,dN(t)/dt = −k) which applies to experimental situations where the surface areais kept fixed during dissolution. The full Brunner-Tolloczko and Noyes-Whitneyequations also apply if sink conditions are not obeyed, but here we consider onlythe sink condition limit. The applicability of the sink condition to our simulationsis confirmed below (Section 6.4).We also remark that another model of dissolution proposed by Nernst154 and127Brunner155 assumes that solute diffusion from a surface layer where the concen-tration approaches the value at saturation, Cs, is the rate determining step in thedissolution process.149 If we consider a spherical diffusion layer of surface area S(t),set S(t) = Sactive(t), and k = DCs/δ, where D is a solute diffusion coefficient, and δis the thickness of the layer, then Eq. (6.2) becomes the Nernst-Brunner equation.In principle, both D and δ could be time dependent, but if these quantities areregarded as constant, then at least for spherical crystals, Eq. (6.2) and the Nernst-Brunner equation are formally identical, even if the underlying physical mechanismsand rate constants are different. Below, we show explicitly that the diffusion layermodel does not apply to the dissolution of NaCl.Integrated rate laws can be obtained by specifying particular forms for Sactive(t).Above, we argued that in the fixed rate law regime, the active surface areas are ap-proximately spherical for the sphere and the cube, the area of a nearly cylindricalwall for the tablet, and mainly the area of the circular base for the rod which an-neals into a roughly cylindrical shape. Assuming these active surface areas, andusing appropriate area to volume relationships, it is easy to deduce that for thesphere and the cube Sactive(t) ∝ N2/3(t), for the tablet Sactive(t) ∝ N1/2(t), and forthe rod Sactive is a constant independent of t. From Eq. (6.2), we then obtain theintegrated rate laws3√N(t) = 3√N0 − kt (Sphere,Cube) , (6.3)√N(t) =√N0 − kt (Tablet) , (6.4)N(t) = N0 − kt (Rod) , (6.5)128where of course k represents a different constant in each rate law. The cube rootlaw applicable to spherical crystals is just the well-known Hixson-Crowell result153originally obtained by integrating the Brunner-Tolloczko equation (or, equivalently,the Nernst-Brunner equation with DCs/δ held fixed) assuming a uniform sphericalsurface that decreases with time in proportion to w2/3, where w is the weight of thecrystal. We note that a square root law has been obtained empirically by Niebergalland Goyan,156 but this is unlikely related to the tablet expression obtained here.Diffusion controlled assumptions have also been used to obtain two-thirds157,167 andsquare root167 laws.We tested all suggested rate laws by attempting linear least squares fits tothe dissolution profiles. The initial and final portions of the dissolution profiles wereremoved from the fit as they constitute a deviation from the laws. The goodnessof fit parameters R2 are given in Table 6.3. Results for a two-thirds root law arealso included since this model is sometimes used in analysis of experimental data.From Table 6.3 we see that even though for some crystal shapes (sphere and tablet)the best fits correspond to the rate laws derived above, in general the differencesamongst the goodness of fit parameters are not large enough to draw any firmconclusions. Essentially, different power laws can fit a particular dissolution profilenearly equally well, and profiles on much longer time scales would be necessary tomake meaningful distinctions.We can also easily check to see if a diffusion layer can reasonably explain theobserved dissolution rates. By calculating the ∆N/∆t in for the spherical crystalin the fixed rate regime (at ≈ 200 ns), we estimate the dissolution rate dN/dt ≈−2.36 × 10−3 ions ps−1. If a diffusion layer is rate determining, then dN/dt =−D ∗ S ∗ Cs/δ, and for our model Cs ≈ 3.29 ions nm−3, S at ≈ 200 ns (assuming129R2 values for the rate lawsSphere Cubelinear 0.98873 0.99846cube root 0.99738 0.99291square root 0.99598 0.99635two thirds root 0.99407 0.99838Rod Tabletlinear 0.98962 0.98798cube root 0.99775 0.99805square root 0.99808 0.99885two thirds root 0.99679 0.99733Table 6.3: Goodness of fit parameters R2 obtained for different rate laws. The bestfits are indicated in bold.spherical geometry) is 58.08 nm2, and at saturation we obtain D+ ≈ 0.36×10−3 andD− ≈ 0.42× 10−3 nm2 ps−1. This means that the thickness of the diffusion layer δwould need to be ∼ 29 nm in order to explain the observed dissolution rate. This isobviously much too thick to be physically relevant, since it greatly exceeds the sizeof our crystal (radius ∼ 2 nm) and system. The ion density as a function of distancefrom the crystal center at ∼ 300 ns is shown in Fig. 6.9, and it is obvious that theion density in solution is nearly uniform and considerably less than that implied byCs, except possibly for a very narrow region near the crystal surface. Therefore, thediffusion layer model clearly does not apply to NaCl dissolution.Rapid DissolutionDissolution profiles (Fig. 6.3) show that in the final stage of crystal disso-lution there is a sharp increase in the rate at which ions are lost to solution. Thecrystal becomes very unstable and disappears rapidly once it has been sufficientlyreduced in size, in all likelihood due to the decreased lattice energy of the small crys-130Figure 6.9: Radial ion density profile starting from the center of a spherical crystalat 300 ns. The dashed horizontal line in the inset indicates the density at saturation.The calculated radius of the crystal is ∼ 2 nm, and a concentration gradient (thatcould be related to a diffusion layer) is not observed.tal. The dissolution profiles indicate that the onset of increased instability occurswhen the crystal contains ∼ 200 ions.A closer view of the final stage of dissolution of the spherical crystal is pro-vided by the snapshots shown in Fig. 6.10. We see that the crystal holds its structuredown to ∼ 64 ions, then deforms and disintegrates completely when only ∼ 12 ionsremain. As one would expect, the onset of the final stage marked by increasinginstability is independent of the initial crystal size. For example, spherical crystalsinitially containing 532 and 240 ion pairs also enter an instability region at ∼ 200ion. A similar picture pertains for crystals of other shapes.It is worth mentioning a possible connection to the opposite process of crys-tal nucleation from supersaturated solution. For model NaCl solution, the criticalnucleus just above saturation is estimated to contain ∼ 75 ions,50 which is reason-ably close to the size where the residual NaCl crystal becomes very unstable in thedissolution process.131Figure 6.10: Snapshots of the last stage of dissolution of the spherical crystal. Whenthe crystal reaches ∼64 ions, the structure begins to disintegrate. In the finalsnapshot, water quickly penetrates the structure and breaks apart the remainingcrystal nucleus.6.4.2 Concentration EffectsIn order to confirm that the dissolution profiles discussed above correspondto sink conditions and are not strongly influenced by increasing salt concentration,dissolution simulations of the spherical crystal were carried out at three salt molefractions (calculated at complete dissolution), 0.0254, 0.0163, and 0.0127, keepingall other conditions fixed. The dissolution profiles are shown in Fig. 6.11, and wenote that the results are very similar for 0.0163, and 0.0127, indicating that any saltconcentration effects are small, and we are indeed in the sink regime in these systems.However, at the highest concentration considered, 0.0254, the dissolution rate clearlyslows down after the initial ∼ 120 ns. There are several possible explanations for theobserved concentration effect at the highest salt mole fraction, and it is interestingto examine some possibilities in more detail.One possible reason is ion reattachments. When the solution surroundingthe crystal reaches a certain concentration, one might expect some ions to reattachto the crystal, hence slowing down the net dissolution rate. This possibility canbe eliminated by counting the number of detachments and attachments over time.132Figure 6.11: Concentration effects on rate. Curves (a), (b), and (c) are for NaClmole fractions (calculated assuming complete dissolution of the entire crystal) of0.0254, 0.0163, and 0.0127. Note that curve (a) shows a sharp decrease in rate at∼120 ns, while there is no substantial difference between the profiles (b) and (c).Figure 6.12 shows that the reattachments are approximately constant at a rate of2 events/ns, while detachments show a sharp rate change at about 120 ns, corre-sponding to the change in the dissolution profile noted above (Fig. 6.11). Sincethere is no substantial variation in the attachment rate, reattachments are not re-sponsible for the change in the dissolution rate. We further remark that since thereattachment rate oscillates but does not increase with time (concentration), mostof the reattachments counted are likely due to ions at the crystal-solution interfacemoving back and forth across the artificially sharp boundary introduced by our or-der parameter. Another possibility is that the ions in solution change the propertiesof the solvent, for example, by effectively “sequestering” a significant fraction of thewater molecules, hence reducing the dissolution rate. To isolate the effect of theions in solution on the dissolution rate, we took the system at 120 ns (before thecrossover) and simply removed the dissolved ions from the solution, while the par-tially dissolved crystal was left in the system. This procedure completely removed133Figure 6.12: Detachment and reattachment events for the spherical crystal. Dur-ing the course of dissolution, the detachment rate systematically slows while theattachment rate remains substantially constant.the crossover in the dissolution rate, confirming that the drop in rate is caused bythe presence of ions in solution, and not by any change in the crystal.6.4.3 Temperature Dependence and Activation EnergyWe would expect ion detachment to be an activated process, and in orderto estimate the activation energy dissolution simulations of spherical crystals (∼1256 ion pairs) were carried out at three temperatures, 300, 320, and 340 K. Rateconstants are found by fitting the cube root law [Eq. (6.2)] in the fixed rate lawregime (i.e., eliminating the initial and final regions of the dissolution profiles).Excellent linear fits are obtained as shown in Fig. 6.13. The rate constants obtainedfollow the Arrhenius equation (Fig. 6.14)ln(k) = ln(A)− Ea/RT , (6.6)134Figure 6.13: Fits to the cube root law for the spherical crystal at temperatures of300, 320 and 340 K.Figure 6.14: Fits to the Arrhenius equation of the rate constants obtained for thespherical crystal at 300, 320, and 340 K. Error bars represent one standard deviationof uncertainty.135giving an activation energy of 33.46± 1.46 kJ mol−1 (SD).The standard deviation in the activation energy quoted above and thoseincluded as error bars in Figure 6.14) were estimated as follows. The trajectoriesobtained for the three temperatures considered were divided into four equal parts(time slices), and rate constants for each part were obtained by fitting to the cuberoot law (eq. 6.2). The four different estimates were then used to obtain theerror bars shown in Figure 6.14. Similarly, by fitting the four different sets of rateconstants to the Arrhenius equation, four estimates of the activation energy wereobtained and used to estimate the standard deviation in the activation energy.6.5 Summary and ConclusionsMolecular dynamics simulations were employed to examine in detail thedissolution of NaCl nanocrystals of different shape. Specifically, cubic, spherical,tablet-like and rod-like crystals were considered with the dissolution carried outunder so-called sink conditions, where the solution is always very dilute comparedto saturation. In all cases, dissolution was found to occur in three distinct stages.Initially, the more exposed and/or most weakly bound ions are quickly detachedfrom the crystal. The initial crystal shape greatly affects the first stage as it mainlyinvolves ions at edges and corners. After the rapid initial period, the crystals takeon solvent-annealed shapes (cubic crystals become nearly spherical and rods andtablets roughly cylindrical) that persist until the final stage of dissolution. Duringthe long intermediate stage the dissolution appears to closely follow a fixed rate law.In the final stage (. 200 ions) the crystal becomes very unstable and dissolves atan ever increasing rate until it disappears.136In the intermediate fixed rate law regime, the dissolution process is well de-scribed by assuming that the rate is proportional to an active surface area from whichthe ions are preferentially detached. For the solvent-annealed cubic and sphericalcrystals the active area is the entire spherical surface, and we obtain the classicalcube root law often employed in dissolution studies of macroscopic crystals, assuminga spherical shape. We show that the cube root law also provides a good descriptionof the nanoscopic crystals considered here in the intermediate stage. We also showthat ion detachment from the surface, and not the existence of a Nernst-Brunnerdiffusion layer,149,154,155 is the rate determining step in the dissolution process forNaCl. For the spherical crystal, simulations were done at three temperatures andthe rate constants, determined by fitting to the cube root law, closely follow theArrhenius equation, giving an activation energy of ∼ 33.5 kJ mol−1.After the initial stage, both the tablet-shaped and rod-shaped crystals haveacquired roughly cylindrical shapes, and we observed that ions did not leave uni-formly over the entire surface area for these geometries. Rather, for the annealedtablet and rod we identified the active surface areas to be the cylinder walls, and thecylinder base, respectively, giving a square root law for the tablet and a linear lawfor the rod. These rate laws give good fits to the dissolution profiles, but the totaldissolution times for the nanocrystals considered here are not long enough to makeclear, unambiguous distinctions amongst the different rate laws. Nevertheless, oursimulations show that, apart from the initial and final stages of dissolution, modelsof the classical type give a good description of the dissolution of NaCl nanocrystals.137Chapter 7Summary and Conclusions7.1 SummaryThe microscopic mechanism of nucleation is still largely unknown becausethe high temporal and spatial resolution needed to observe the process is difficult toreach using current experimental techniques.32 A considerable body of research hasbeen conducted using molecular simulation that, thanks to the availability of parallelcomputational power, can be used to model nucleation for a range of compounds ofscientific and industrial interest.32 In this thesis, we investigated different aspects ofnucleation using alkali halide salts as a model, and designed data analysis techniquesable to better exploit the results produced by molecular dynamics simulations.In Chapter 2, the key algorithms and methods of molecular dynamic sim-ulations are introduced, as well as the models used to represent alkali halides inaqueous solutions.Chapter 3 describes an investigation of nucleation of NaCl from aqueoussolution. We found that NaCl nucleation is observable within simulation timescales138(10 to 200 ns), at 300 K and 1 bar, starting from a solution at about twice thesaturation concentration for the model. A major challenge to the investigation ofthe early stages of nucleation was the lack of a reliable method to detect smallclusters in solution. To address this issue, we developed a novel three-step methodto detect cluster formation, and follow cluster evolution over time. Thanks to thenew methodology, it was possible to collect and monitor several properties of thelarge number of clusters that continuously form and disappear from solution.By analyzing a variety of cluster properties we showed, for the first timeby direct measurement, that the lifetime and probability of nucleation of a clusterdepend not only on size, but also on the specific geometrical arrangement of thecluster. The importance of this result lies in its contrast with the common CNTassumption that the probability of nucleation depends only on cluster size. A wayto extend CNT to include structural effects could be the addition of a structuredependent free energy in the exponential term of the nucleation rate. It is also usefulto notice that factors such as geometrical arrangement are not merely surrogates ofsurface tension, as they do not scale with surface area.In Chapter 4 we investigated crystallization from the melt of the TF and JCmodels for lithium halides, partially to assess their probable behavior in simulationsof crystallization from aqueous solutions. We found that the TF model incorrectlypredicts the wurtzite crystal structure as the stable structure for lithium salts, whilethe more recent JC model reproduces the correct (rock salt) stable crystal structurefor LiF and LiCl, but not for LiBr and LiI. Furthermore, we found that the wurtzitecrystal structure is highly irregular in finite size crystals, and that this can be likelyattributed to surface rearrangements that are necessary to avoid the formation of asurface dipole moment.139In Chapter 5 we further explore crystal nucleation from solution focusing onlithium fluoride. For the JC model, which has the correct rock salt crystal structure,nucleation and growth can be observed only at high temperature and a high degreeof supersaturation. By measuring the growth rate at different temperatures wefound that growth requires a high activation energy, which is likely associated withthe high barrier required to remove water molecules from the first solvation shellof an ion. While we couldn’t measure an activation energy for nucleation, it isreasonable to assume that the activation energy for nucleation is equal to, or greater,than the activation energy for growth. We also found that the structure of themetastable solution is quite different from that observed for NaCl. In LiF solutions,the metastable ionic clusters tend to be much shorter lived and less regular in shape.Based on our observations, in supersaturated LiF aqueous solutions, a crucialstep to nucleation is the ability of an ion to lose at least part of its solvation shelland join a cluster. In this system, the applicability of CNT is in question becausethe mechanism seems to invalidate some of the foundational assumptions of CNT.For example, the LiF supersaturated solution is characterized by the presence ofmetastable ionic clusters, that, compared to those found in NaCl solutions, aresmall and do not resemble the bulk crystal. Interestingly, these metastable speciesincrease in number at higher temperatures and, their formation likely constitutes anecessary step to nucleation.In Chapter 6, the dissolution of NaCl was investigated to gain insight intothe stability of finite size ionic crystals in water, and the dynamics of ion attachmentand detachment. It was found that NaCl dissolution can be modeled as a three-stage process where the crystal initially loses ions around the edges, then followsan essentially fixed dissolution rate law until it reaches a certain size (. 200 ions),140where it becomes very unstable and quickly dissolves into solution. By measuringthe dissolution rate of NaCl at different temperatures, it was found that the processis activated, with a barrier of ∼ 33 kJ mol−1. Additionally, we found that it ispossible to obtain rate laws from a simple expression involving the surface area ofthe crystal.7.2 Future DirectionsOur results suggest that, in general, the mechanism of homogeneous nucle-ation commonly assumed in CNT does not correctly model the species that precedethe formation of a critical nucleus. In the case of NaCl nucleation, we found thatsize is not the only factor affecting nucleation, and geometric arrangement also influ-ences the probability of nucleation, as well as the survival of prenucleation clusters.Despite the influence of geometric regularity, the process observed resembles theCNT process in that there is a distribution of clusters that stochasticly form or dis-appear in solution, even though the clusters do not have the same properties as thebulk phase, since they are mostly made of interfacial ions. The presence of small,metastable clusters can be also seen as evidence of a multi-step nucleation process,where the prenucleation clusters are intermediate species that convert to the bulkphase as the cluster grows.A mechanism that explains the data obtained from LiF simulations of Chap-ter 5 is a process with two kinetic barriers. The first, is a barrier to cluster formation.In order to form metastable clusters, it is necessary to surpass a barrier that appearsto be related to the solvent-solute interactions. The second barrier, is a nucleationbarrier that likely depends on cluster related properties (e.g. size, surface tension141and crystallinity). Characterization of nucleation rates at different temperatures(through simulation) could help probe the nature of the second barrier, and lead tocomparison, validation, and extension of multi-step nucleation theories.144,168Further speculating, a compound where the barrier to generating prenu-cleation clusters is very weak could have a behavior similar to that exhibited byCaCO3,25 where there is aggregation into an amorphous phase, followed by rear-rangement into a crystal structure. Investigation on more limiting cases would beuseful in formalizing a theory that encompasses multiple nucleation mechanisms.In this respect, a kinetic formulation involving multiple steps such as diffusion, de-hydration, and aggregation would possibly constitute a more realistic model fornucleation.Model development for solute and solvent interactions is also of crucial impor-tance. Models, in order to give realistic results in simulations of phase transitions,should reproduce both solution and solid properties such as solubilities and latticeenergies. This is especially important for compounds that form hydrates (e.g. LiCl)where a model that at least supports the correct stable structure at the pressureand temperature of interest is necessary. Additionally, increasing attention shouldbe devoted to the study of finite-size structures since, as we have shown throughoutthis thesis, crystallization nuclei can be very small and have different propertiesthan the bulk crystal phase.As we have shown, by using medium- and large-scale molecular dynamics sim-ulations combined with advanced data analysis techniques, it is possible to directlyobserve individual nucleation events, degree of ion association, cluster lifetimes andother properties of interest. More experiments are definitely needed to validate themodels and, while some work has been done for NaCl,29 a better coverage of simple142compounds would be effective in evaluating the results obtained through molecularsimulation.A combination of realistic models, experiments, and data analysis techniqueswill ultimately lead to a better understanding of crystal nucleation and a refinementof classical theories.143Bibliography[1] I. S. Joung and T. E. Cheatham, J. Phys. Chem. B 112, 9020 (2008).[2] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91,6269 (1987).[3] J. Angus, Structure and Thermochemistry of Diamond, The Physics of Dia-mond, OS Press, Amsterdam, 1997.[4] A. L. DeVries and D. E. Wohlschlag, Science 163, 1073 (1969).[5] F. Franks, Biophysics and biochemistry at low temperatures, Cambridge Uni-versity Press, 1985.[6] M. K. Yau and R. Rogers, A short course in cloud physics, Elsevier, 1996.[7] M. Kumai, J. Meteor. 18, 139 (1961).[8] M. Kumai and K. E. Francis, J. Atmos. Sci. 19, 474 (1962).[9] P. G. Debenedetti, Metastable Liquids: Concepts and Principles, PrincetonUniversity Press, 1996.[10] M. Volmer and A. Weber, Z. Phys. Chem. 119, 277 (1926).[11] L. Farkas, Z. Phys. Chem. 125 (1927).[12] R. Becker and W. Do¨ring, Ann. Phys. (Berlin) 416, 719 (1935).144[13] J. B. Zeldovich, Acta physicochim. URSS 18, 1 (1943).[14] J. Frenkel, Kinetic Theory of Liquids, Dover Publications, Dover, 1955.[15] J. L. Schmitt, G. W. Adams, and R. A. Zalabsky, J. Chem. Phys. 77, 2089(1982).[16] A. Kacker and R. H. Heist, J. Chem. Phys. 82, 2734 (1985).[17] R. Strey, T. Schmeling, and P. E. Wagner, J. Chem. Phys. 85, 6192 (1986).[18] C. Hung, M. J. Krasnopoler, and J. L. Katz, J. Chem. Phys. 90, 1856 (1989).[19] J. Wlk, , and R. Strey*, J. Phys. Chem. B 105, 11683 (2001).[20] S. Auer and D. Frenkel, Nature 409, 1020 (2001).[21] J. L. Harland and W. van Megen, Phys. Rev. E 55, 3054 (1997).[22] A. Heymann, A. Stipp, C. Sinn, and T. Palberg, J. Colloid Interface Sci. 207,119 (1998).[23] P. Wette and H. J. Scho¨pe, Phys. Rev. E 75, 051405 (2007).[24] J. D. Rodriguez-Blanco, S. Shaw, and L. G. Benning, Nanoscale 3, 265 (2011).[25] D. Gebauer, A. Vo¨lkel, and H. Co¨lfen, Science 322, 1819 (2008).[26] D. Gebauer, M. Kellermeier, J. D. Gale, L. Bergstro¨m, and H. Co¨lfen, Chem.Soc. Rev. 43, 2348 (2014).[27] W. Ostwald, Z. Phys. Chem. 22, 289 (1897).[28] M. Sleutel, J. Lutsko, A. E. Van Driessche, M. A. Dura´n-Olivencia, andD. Maes, Nat. Commun. 5 (2014).145[29] J. Desarnaud, H. Derluyn, J. Carmeliet, D. Bonn, and N. Shahidzadeh, J.Phys. Chem. Lett. 5, 890 (2014).[30] H.-S. Na, S. Arnold, and A. S. Myerson, J. Cryst. Growth 139, 104 (1994).[31] R. J. Davey, S. L. M. Schroeder, and J. H. terHorst, Angew. Chem. Int. Ed.52, 2166 (2013).[32] G. C. Sosso, J. Chen, S. J. Cox, M. Fitzner, P. Pedevilla, A. Zen, andA. Michaelides, Chem. Rev. (2016).[33] J.-P. Hansen and L. Verlet, Phys. Rev. 184, 151 (1969).[34] J. Q. Broughton, G. H. Gilmer, and K. A. Jackson, Phys. Rev. Lett. 49, 1496(1982).[35] R. O. Rosenberg and D. Thirumalai, Phys. Rev. A 36, 5690 (1987).[36] D. S. Corti and P. G. Debenedetti, Chem. Eng. Sci. 49, 2717 (1994).[37] F. de Wette, R. Allen, D. Hughes, and A. Rahman, Phys. Lett. A 29, 548(1969).[38] M. J. Mandell, J. P. McTague, and A. Rahman, J. Chem. Phys. 64, 3699(1976).[39] P. Rein ten Wolde, M. J. RuizMontero, and D. Frenkel, J. Chem. Phys. 104,9932 (1996).[40] P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, J. Chem. Phys. 110,1591 (1999).[41] B. O’Malley and I. Snook, Phys. Rev. Lett. 90, 085702 (2003).146[42] E. B. Moore and V. Molinero, Phys. Chem. Chem. Phys. 13, 20008 (2011).[43] J. Anwar and D. Zahn, Angew. Chem. Int. Ed. 50, 1996 (2011).[44] M. Salvalaglio, C. Perego, F. Giberti, M. Mazzotti, and M. Parrinello, Proc.Natl. Acad. Sci. 112, E6 (2015).[45] F. Jimenez-Angeles and A. Firoozabadi, J. Phys. Chem. C 118, 11310 (2014).[46] H. Ohtaki and N. Fukushima, Pure Appl. Chem. 63 (1991).[47] D. Zahn, Phys. Rev. Lett. 92, 040801 (2004).[48] F. Giberti, G. A. Tribello, and M. Parrinello, J. Chem. Theory Comput. 9,2526 (2013).[49] D. Chakraborty and G. N. Patey, Chem. Phys. Lett. 587, 25 (2013).[50] D. Chakraborty and G. N. Patey, J. Phys. Chem. Lett. 4, 573 (2013).[51] J. Yan and G. N. Patey, J. Phys. Chem. Lett. 2, 2555 (2011).[52] J. Yan and G. N. Patey, J. Phys. Chem. A 116, 7057 (2012).[53] J. Yan, S. Overduin, and G. N. Patey, J. Chem. Phys. 141, 074501 (2014).[54] S. A. Zielke, A. K. Bertram, and G. N. Patey, J. Phys. Chem. B 119, 9049(2015).[55] S. A. Zielke, A. K. Bertram, and G. N. Patey, J. Phys. Chem. B 120, 1726(2015).[56] S. A. Zielke, A. K. Bertram, and G. N. Patey, J. Phys. Chem. B 120, 2291(2016).147[57] D. Murphy, Geophys. Res. Lett. 30 (2003).[58] S. A. Hassan, Phys. Rev. E 77, 031501 (2008).[59] S. A. Hassan, J. Chem. Phys. 134, 114508 (2011).[60] A. Benavides, J. Aragones, and C. Vega, J. Chem. Phys. 144, 124504 (2016).[61] I. S. Joung and T. E. Cheatham, J. Phys. Chem. B 113, 13279 (2009).[62] G. A. Orozco, O. A. Moultos, H. Jiang, I. G. Economou, and A. Z. Pana-giotopoulos, J. Chem. Phys. 141, 234507 (2014).[63] J. L. Aragones, M. Rovere, C. Vega, and P. Gallo, J. Phys. Chem. B 118,7680 (2014).[64] F. Moucˇka, I. Nezbeda, and W. R. Smith, J. Chem. Phys. 138, 154102 (2013).[65] J. Aragones, E. Sanz, C. Valeriani, and C. Vega, J. Chem. Phys. 137, 104507(2012).[66] N. E. Zimmermann, B. Vorselaars, D. Quigley, and B. Peters, J. Am. Chem.Soc. 137, 13352 (2015).[67] T. A. Halgren, J. Comput. Chem. 17, 490 (1996).[68] P.-O. A˚strand, P. Linse, and G. Karlstro¨m, Chem. Phys. 191, 195 (1995).[69] H. Sun, S. J. Mumby, J. R. Maple, and A. T. Hagler, J. Am. Chem. Soc. 116,2978 (1994).[70] W. Eckhardt, A. Heinecke, R. Bader, M. Brehm, N. Hammer, H. Huber, H.-G. Kleinhenz, J. Vrabec, H. Hasse, M. Horsch, M. Bernreuther, C. W. Glass,148C. Niethammer, A. Bode, and H.-J. Bungartz, 591 TFLOPS Multi-trillionParticles Simulation on SuperMUC, pages 1–12, Springer Berlin Heidelberg,Berlin, Heidelberg, 2013.[71] S. W. de Leeuw, J. W. Perram, and E. R. Smith, Proceedings of the RoyalSociety of London A: Mathematical, Physical and Engineering Sciences 373,27 (1980).[72] P. P. Ewald, Ann. Phys. (Berlin) 369, 253 (1921).[73] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, ClarendonPress: Oxford, United Kingdom, New York, NY, USA, 1989.[74] T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993).[75] J. W. Cooley and J. W. Tukey, Math. Comp. 19, 297 (1965).[76] S. R. Hall, F. H. Allen, and I. D. Brown, Acta Crystallogr. Sect. A 47, 655(1991).[77] B. Hess et al., J. Comput. Chem. 18, 1463 (1997).[78] H. J. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. Di Nola, andJ. Haak, J. Chem. Phys. 81, 3684 (1984).[79] G. Bussi, D. Donadio, and M. Parrinello, J. Chem. Phys. 126, 014101 (2007).[80] M. Parrinello and A. Rahman, J Appl Phys 52, 7182 (1981).[81] H. Blatt, R. Tracy, and B. Owens, Petrology: Igneous, Sedimentary, andMetamorphic, Macmillan, 2006.149[82] K. L. Webber, A. U. Falster, W. B. Simmons, and E. E. Foord, J. Petrol. 38,1777 (1997).[83] P. Arosio, M. Vendruscolo, C. M. Dobson, and T. P. Knowles, Trends Phar-macol. Sci. 35, 127 (2014).[84] B. E. Rabinow, Nat. Rev. Drug Discovery 3, 785 (2004).[85] A. Laaksonen, V. Talanquer, and D. W. Oxtoby, Annu. Rev. Phys. Chem.46, 489 (1995).[86] M. Kulmala, T. Peta¨ja¨, M. Ehn, J. Thornton, M. Sipila¨, D. Worsnop, andV.-M. Kerminen, Annu. Rev. Phys. Chem. 65, 21 (2014).[87] D. Aastuen, N. Clark, L. Cotter, and B. J. Ackerson, Phys. Rev. Lett. 57,1733 (1986).[88] Z. Wang, F. Wang, Y. Peng, and Y. Han, Nat. Commun. 6 (2015).[89] N. Pienack and W. Bensch, Angew. Chem. Int. Ed. 50, 2014 (2011).[90] M. Matsumoto, S. Saito, and I. Ohmine, Nature 416, 409 (2002).[91] M. Kellermeier, A. Picker, A. Kempter, H. Co¨lfen, and D. Gebauer, Adv.Mater. 26, 752 (2014).[92] J. Merikanto, E. Zapadinsky, A. Lauri, and H. Vehkama¨ki, Phys. Rev. Lett.98, 145702 (2007).[93] P. G. Vekilov, Cryst. Growth Des. 4, 671 (2004).[94] V. J. Anderson and H. N. Lekkerkerker, Nature 416, 811 (2002).150[95] Z. Mester and A. Z. Panagiotopoulos, J. Chem. Phys. 143, 044505 (2015).[96] H. M. Manzanilla-Granados, H. Saint-Mart´ın, R. Fuentes-Azcatl, and J. Ale-jandre, J. Phys. Chem. B 119, 8389 (2015).[97] K. Kobayashi, Y. Liang, T. Sakka, and T. Matsuoka, J. Chem. Phys. 140,144705 (2014).[98] J. Aragones, E. Sanz, and C. Vega, J. Chem. Phys. 136, 244508 (2012).[99] Merck index 14th ed. whitehouse station, nj: Merck and co, 2006.[100] B. Hess, C. Kutzner, D. Van der Spoel, and E. Lindahl, J. Chem. TheoryComp. 4, 435 (2008).[101] G. Lanaro and G. N. Patey, J. Phys. Chem. B 119, 4275 (2015).[102] J. R. Espinosa, C. Vega, C. Valeriani, and E. Sanz, J. Chem. Phys. 142,194709 (2015).[103] W. Lechner and C. Dellago, J. Chem. Phys. 129, 114707 (2008).[104] R. Radhakrishnan and B. L. Trout, J. Am. Chem. Soc. 125, 7743 (2003).[105] A. Reinhardt, J. P. Doye, E. G. Noya, and C. Vega, J. Chem. Phys. 137,194504 (2012).[106] T. Li, D. Donadio, G. Russo, and G. Galli, Phys. Chem. Chem. Phys. 13,19807 (2011).[107] P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784(1983).151[108] M. Ester, H.-P. Kriegel, J. Sander, and X. Xu, A density-based algorithm fordiscovering clusters in large spatial databases with noise., in Kdd, volume 96,pages 226–231, 1996.[109] F. Pedregosa et al., J. Mach. Learn. Res. 12, 2825 (2011).[110] P. Jaccard, New Phytol 11, 37 (1912).[111] E. L. Kaplan and P. Meier, J. Am. Stat. Assoc. 53, 457 (1958).[112] C. P. O’Hagan, B. J. O’Brien, F. Griffin, B. Hooper, S. B. Leen, and R. F. D.Monaghan, Energy & Fuels 29, 3082 (2015).[113] A. R. Kamali and D. J. Fray, Carbon 77, 835 (2014).[114] J. R. Hardy and A. M. Karo, The Lattice Dynamics and Static of Alkali HalideCrystals, Springer Science + Business Media, 1979.[115] T. Martin, Phys. Rep 95, 167 (1983).[116] A. Bach, D. Fischer, and M. Jansen, Z. Anorg. Allg. Chem. (2009).[117] Y. Liebold-Ribeiro, D. Fischer, and M. Jansen, Angew. Chem. 120, 4500(2008).[118] D. Fischer, A. Mueller, and M. Jansen, ChemInform 36 (2005).[119] P. C. Rodrigues and F. M. S. Fernandes, EPJ D 44, 109 (2007).[120] P. C. Rodrigues and F. M. S. Fernandes, J Mol Struc-THEOCHEM 946, 94(2010).[121] T. Croteau and G. N. Patey, J. Chem. Phys. 124, 244506 (2006).152[122] Zˇ. P. Cˇancˇarevic´, J. C. Schoen, and M. Jansen, Chem Asian J 3, 561 (2008).[123] A. Aguado, A. Ayuela, J. M. Lo´pez, and J. A. Alonso, Phys. Rev. B 56, 15353(1997).[124] L. Pauling, The Nature of the Chemical Bond and the Structure of Moleculesand Crystals: an Introduction to Modern Structural Chemistry, volume 18,Cornell university press, 1960.[125] M. Tosi and F. Fumi, J. Phys. Chem. Solids 25, 45 (1964).[126] J. R. Espinosa, C. Vega, C. Valeriani, and E. Sanz, J. Chem. Phys. 144,034501 (2016).[127] J. Wang, J. Wu, Z. Sun, G. Lu, and J. Yu, J. Mol. Liq. 209, 498 (2015).[128] C. Valeriani, E. Sanz, and D. Frenkel, J. Chem. Phys. 122, 194501 (2005).[129] T. Zykova-Timan, D. Ceresoli, U. Tartaglino, and E. Tosatti, J. Chem. Phys.123, 164701 (2005).[130] A. B. Belonoshko, R. Ahuja, and B. Johansson, Phys. Rev. B 61, 11928(2000).[131] F. Fumi and M. Tosi, J. Phys. Chem. Solids 25, 31 (1964).[132] M. Born and J. E. Mayer, Z. Phys. 75, 1 (1932).[133] M. L. Huggins and J. E. Mayer, J. Chem. Phys. 1, 643 (1933).[134] J. E. Mayer, J. Chem. Phys. 1, 270 (1933).[135] G. Lanaro, chemlab 0.4, https://chemlab.readthedocs.io, 2015.153[136] Y. Sakamoto, J. Chem. Phys. 28, 164 (1958).[137] M. Parrinello and A. Rahman, Phys. Rev. Lett. 45, 1196 (1980).[138] D. B. Sirdeshmukh, L. Sirdeshmukh, and K. Subhadra, Alkali Halides: AHandbook of Physical Properties, volume 49, Springer Science & BusinessMedia, 2013.[139] D. E. Smith and L. X. Dang, J. Chem. Phys. 100, 3757 (1994).[140] G. Lanaro and G. N. Patey, J. Phys. Chem. B 120, 9076 (2016).[141] L. Le and V. Molinero, J. Phys. Chem. A 115, 5900 (2010).[142] C. B. Stubblefield and R. O. Bach, J. Chem. Eng. Data 17, 491 (1972).[143] W. A. Hart, O. Beumel, and T. P. Whaley, The Chemistry of Lithium, Sodium,Potassium, Rubidium, Cesium and Francium: Pergamon Texts in InorganicChemistry, volume 7, Elsevier, 2013.[144] D. Erdemir, A. Y. Lee, and A. S. Myerson, Acc. Chem. Res. 42, 621 (2009),PMID: 19402623.[145] J. B. Dressman, G. L. Amidon, C. Reppas, and V. P. Shah, Pharm. Res. 15,11 (1998).[146] O. Anand, X. Y. Lawrence, D. P. Conner, and B. M. Davit, AAPS J. 13, 328(2011).[147] A. Dokoumetzidis and P. Macheras, Int. J. Pharm. 321, 1 (2006).[148] S. S. Jambhekar and P. J. Breen, Drug Disc. Today 18, 1173 (2013).154[149] J. Siepmann and F. Siepmann, Int. J. Pharm. (2013).[150] P. Costa and J. M. Sousa Lobo, Eur. J. Pharm. Sci. 13, 123 (2001).[151] A. A. Noyes and W. R. Whitney, J. Am. Chem. Soc. 19, 930 (1897).[152] L. Brunner and S. Tolloczko, Z. Phys. Chem , 283 (1900).[153] A. Hixson and J. Crowell, Industrial & Engineering Chemistry 23, 923 (1931).[154] W. Nernst, Z. Phys. Chem 47, 52 (1904).[155] E. Brunner, Z. Phys. Chem 43, 56 (1904).[156] P. Niebergall, G. Milosovich, and J. Goyan, J. Pharm. Sci. 52, 236 (1963).[157] W. Higuchi and E. Hiestand, J. Pharm. Sci. 52, 67 (1963).[158] S. Piana and J. D. Gale, J. Am. Chem. Soc. 127, 1975 (2005).[159] Y. Gao and K. W. Olsen, Mol. Pharm. 10, 905 (2013).[160] N. H. De Leeuw, S. C. Parker, and J. H. Harding, Phys. Rev. B 60, 13792(1999).[161] H. Shinto, T. Sakakibara, and K. Higashitani, J. Phys. Chem. B 102, 1974(1998).[162] R. Bahadur, L. M. Russell, S. Alavi, S. T. Martin, and P. R. Buseck, J. Chem.Phys. 124, 154713 (2006).[163] O. Engkvist and A. J. Stone, J. Chem. Phys. 112, 6827 (2000).[164] J. Klime, D. R. Bowler, and A. Michaelides, J. Chem. Phys. 139, (2013).155[165] N. Holmberg, J. C. Chen, A. S. Foster, and K. Laasonen, Phys. Chem. Chem.Phys. (2014).[166] C. M. Riley and T. W. Rosanske, Development and Validation of AnalyticalMethods, volume 3, Elsevier, 1996.[167] J. Wang and D. R. Flanagan, J. Pharm. Sci. 88, 731 (1999).[168] W. Pan, A. B. Kolomeisky, and P. G. Vekilov, J. Chem. Phys. 122, 174905(2005).[169] T. J. Richmond, J. Mol. Biol. 178, 63 (1984).[170] W. E. Lorensen and H. E. Cline, Marching cubes: A high resolution 3d surfaceconstruction algorithm, in ACM siggraph computer graphics, volume 21, pages163–169, ACM, 1987.[171] D. Xu and Y. Zhang, PloS one 4, e8140 (2009).156Appendix ADistribution of Various Propertiesfor Nucleated and Failed Clusters.In Figures A.1 and A.2, distributions of various properties for nucleated andfailed clusters of size 10 and 30 are reported. A description of such properties isgiven below:• Average neighbor count is calculated by counting for each ion of the clusterthe number of neighbors within 0.6 nm. The counts are then averaged for everyion in the cluster.• Volume is approximated by encasing a cluster in a cuboid box, subdividedinto 32× 32× 32 grid points, and counting the number of grid points that areinside of the cluster, the volume is obtained by multiplying this number bythe volume of the cell (Figure A.3). The radius that each ion covers is chosento be the ionic radius plus the solvent radius, according to the definition ofexcluded volume.169• Surface area is approximated by extracting the surface from the volumet-ric representation described above using the marching cube algorithm,170 theprocedure is inspired by the work of Xu et al.171157• Sphericity is the ratio of the surface area of a a sphere with equivalent volumeand the surface area of the cluster, and is defined aspi1/3(6V )2/3A,where V is volume and A is surface area. Its value is 1 when the object isperfectly spherical and assumes lower values for less spherical objects. It isoften used to represent the compactness of a three dimensional object.• Hydration is a measure of water content in the cluster and is calculated bycounting for each ion the number of water molecules within 0.3 nm, the countsare then averaged over all the ions in the clusters.• Radius of gyration is defined asR2g =1N∑a(xa − x¯a)2,where the sum is over the N ions in the cluster, xa is the position of the ion,and x¯a is the geometric center of the cluster. It is commonly used as a measureof compactness.At size 30, the results are consistent with the intuition that more compactclusters are favored for nucleation. Surface area, sphericity, and the radius of gy-ration of nucleated cluster tend to assume values compatible with compact shapes.The effect, however, is not noticeable at size 10, suggesting that those measures arenot sensitive enough to describe the fine structure of smaller clusters, and that theseaspects do not influence the early stages of nucleation.158Figure A.1: Distributions of various properties for failed clusters (blue histogram) ofsize 10. Clusters that achieved nucleation are indicated by single orange lines. Notethat the property values of the clusters that achieve nucleation are found aroundthe mean of the distribution.159Figure A.2: Distributions of various properties for failed clusters (blue histogram)of size 30. Clusters that achieved nucleation are indicated by single orange lines.At this size, some preference for low radius of gyration, low surface area and highersphericity can be observed.160Figure A.3: Two dimensional representation of the algorithm used to estimate thevolume. The cluster is encased in a grid, and the grid points that lie within thecluster are represented in red.161
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Molecular simulation of nucleation and dissolution...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Molecular simulation of nucleation and dissolution of alkali halides Lanaro, Gabriele 2017
pdf
Page Metadata
Item Metadata
Title | Molecular simulation of nucleation and dissolution of alkali halides |
Creator |
Lanaro, Gabriele |
Publisher | University of British Columbia |
Date Issued | 2017 |
Description | The process of crystal nucleation, despite being so fundamental and ubiquitous in industrial and natural processes, is still not fully understood because of its stochastic nature, and the high spatial and temporal resolution needed to observe it through experiments. This thesis investigates several aspects of nucleation through the use of molecular dynamics, a computational technique that is able to simulate systems up to ~10¹² atoms (as of today's computational power). The projects in this thesis focus on the nucleation from aqueous solution of alkali halide salts, with supplementary studies on the related processes of dissolution in water, and crystallization from the melt. The mechanism of NaCl nucleation from solution is examined in Chapter 3 by direct simulation. The NaCl supersaturated solution was found to contain many small ionic clusters that continuously form and disappear from solution until one (or more) of them nucleates and grows irreversibly. An original method was developed to detect and follow clusters in time, producing results useful in the study of their characteristics and lifetimes. Most importantly, it was found that the lifetime of transient clusters is about ~1 ns, and that both the cluster lifetime and nucleation probability are significantly higher if the cluster is more geometrically ordered. The dissolution of NaCl crystals was also investigated. The process was found to happen in stages, is characterized by an activation barrier, and can be described by a simple rate law. The crystal nucleation of LiF from supersaturated solution was observed, in our simulations, only at high pressure and temperature. The growth rate for an already nucleated crystal was found to have a temperature dependence that follows the Arrhenius law, and further evidence suggests that the reason for such behavior is the high activation energy required to dehydrate the ions. The crystallization from the melt of the Joung-Cheatham and Tosi-Fumi models for lithium halides was also investigated. We found that, for the Tosi-Fumi model, all lithium halides crystallize as wurtzite. For the Joung-Cheatham model, LiF and LiCl crystallize as rock salt, while LiBr and LiI crystallize as wurtzite. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2017-03-15 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-ShareAlike 4.0 International |
DOI | 10.14288/1.0343232 |
URI | http://hdl.handle.net/2429/60909 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2017-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-sa/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2017_may_lanaro_gabriele.pdf [ 10.67MB ]
- Metadata
- JSON: 24-1.0343232.json
- JSON-LD: 24-1.0343232-ld.json
- RDF/XML (Pretty): 24-1.0343232-rdf.xml
- RDF/JSON: 24-1.0343232-rdf.json
- Turtle: 24-1.0343232-turtle.txt
- N-Triples: 24-1.0343232-rdf-ntriples.txt
- Original Record: 24-1.0343232-source.json
- Full Text
- 24-1.0343232-fulltext.txt
- Citation
- 24-1.0343232.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0343232/manifest