Coarse-Grained Simulation of CaCO3 Aggregation and Crystallization

Subscriber access provided by TULANE UNIVERSITY

C: Physical Processes in Nanomaterials and Nanostructures

Coarse-Grained Simulation of CaCO3 Aggregation and Crystallization Made Possible by Non-Bonded Three-Body Interactions Michael King, Simon Pasler, and Christine Peter J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b09604 • Publication Date (Web): 16 Jan 2019 Downloaded from http://pubs.acs.org on January 17, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Coarse-Grained Simulation of CaCO3 Aggregation and Crystallization Made Possible by Non-Bonded Three-Body Interactions Michael King, Simon Pasler, and Christine Peter∗ Department of Chemistry, University of Konstanz, 78464 Konstanz, Germany E-mail: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Calcium-containing minerals are key model systems to investigate fundamental principles of nucleation and mineral formation both experimentally and by simulation. Due to the rare event character of nucleation, the different dimensions of pre- and postnucleation stages and the possible relevance of non-classical nucleation pathways, such investigations require advanced sampling techniques and simulation models on a range of resolution levels. To this end we have developed coarse-grained (CG) models for calcium carbonate. We present a strategy to devise CG parameters – including nonbonded angular-dependent terms – such that the model correctly represents the calcite phase along with properties of the constituents in solution. We show how the CG interactions affect the crystallization pathways by stabilizing different intermediates spanning a wide range of degrees of crystallinity and water content. This will allow us to investigate contributions to crystallization transitions and link them to experimentally observed non-classical crystallization intermediates.

Introduction Calcium-containing minerals, such as calcite or hydroxyapatite, are not only of major importance as structural materials both in living and inanimate objects 1,2 , they are also key model systems to investigate how nanostructured, multi-component mineral-based materials form. In recent years experimental and computational studies have shed new light on the processes and structures that occur during the early stages of mineralization, and the traditional view of classical nucleation theory and classical growth models was significantly expanded. New pathways and intermediates have been investigated, involving prenucleation clusters, liquid-liquid phase separation, or metastable amorphous precursor phases 3–6 . In this context, classical molecular dynamics simulation can provide a valuable microscopic view on the processes, transition states, mechanistic principles and driving forces. Unfortunately, due to the rare event character of the nucleation and phase transition steps 2


Page 2 of 25

Page 3 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


involved, simulations at full atomistic resolution are usually limited to separate stages or intermediates of the total mineralization process, such as the formation of prenucleation clusters 4,7 , the stability of nanoparticles 8–10 or the growth/dissolution of surfaces 8,11 . In spite of the use of enhanced sampling techniques 9,12 , the time- and – most importantly – length-scale limitations are so severe, that appropriate scale-bridging techniques are merited if one wants to cover the entire mineralization process. Here, reduced-resolution models, in particular particle-based coarse-grained (CG) models are of central importance. An unfortunate consequence of reducing the resolution is, however, that CG models inevitably do not represent all structural and thermodynamic properties of a given system equally well. Moreover, they are very state-point dependent, i.e. are usually not simultaneously applicable for very different chemical environments or phases 13,14 . This is the particular challenge if one wants to design a CG model suitable to simulate mineralization: the model needs to represent very different state points, including ions in solution at varying concentrations, ionic solutions in contact with solid mineral phases, and the bulk crystalline phases themselves. In the latter case, the need to account for anisotropic entities and interactions poses yet another difficulty. Here, we present a suite of CG models for calcium carbonate that address these challenges, i.e. that are designed to represent both the pre-nucleation solution phase as well as the early stages and transitions during nanoparticle formation and growth. These models are based on reference data from atomistic simulations of solutions and solid phases and are tightly connected to the atomistic scale, which will allow us in the future to reinsert atomistic details and back-calculate experimentally accessible properties. Due to the use of multiple reference state points and parametrization targets, and due to the large number of different CG interactions which are highly interdependent, the design of a CG model is not unambiguous. We present a sequential strategy to determine the different types of CG interaction potentials and ultimately three CG models that reproduce both the thermodynamically stable crystalline phase (calcite) as well as properties of the aqueous ionic solution 3



Page 4 of 25

such as ion solvation structure and ion pairing. We then simulate aggregation from supersaturated solution and monitor the formation of metastable amorphous clusters or crystalline nanoparticles. Interestingly, the models exhibit different aggregation pathways and different types of nanoparticles, representing various mineralization stages that have been found in experimental studies. This demonstrates how CG models not only give access to larger system sizes and processes on longer timescales but also – by systematic tuning of interactions – open up the possibility to tackle fundamental questions regarding classical and non-classical nucleation mechanisms.

Computational Methods All MD simulations were performed using the LAMMPS code 15 . Atomistic reference simulation were carried out with the CaCO3 force field of Raiteri et al. 4 and the SPC/Fw water model. NPT simulations were run with a 1 fs time step, a Nosé-Hoover chain thermostat with a chain length of 5 and a 100 fs relaxation time, and a barostat of chain length 5 and a 1000 fs relaxation time 16,17 . Electrostatic interactions were computed with the PPPM algorithm 18 with an accuracy of 10−5 . CG simulations were carried out with a 5 fs time step and the same thermostat and barostat settings as the atomistic simulations. In the following the general procedure of parameterizing the CG model is described; more details and the full set of CG interaction parameters are presented the SI. For the analysis of the crystalline order in nanoparticles, a local q6 Steinhardt parameter was calculated by using PLUMED 19 with the switching function used by Quigley et al. 9 , as shown in fig. S5. The crystallinity of a cluster or nanoparticle was calculated by eq. (1),

χ=

Ncrystalline Ncluster

(1)

with Ncrystalline being the number of crystalline-like particles in a cluster, i.e. atoms with a local q6 larger than 0.3, and Ncluster being the number of particles in the cluster. Further 4


Page 5 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


technical and implementation details are given in the SI.

Parameterization of the CG model For the present CG model we chose a level of resolution where water molecules are represented as single CG beads and used the monoatomic water (mW) model of Molinero and coworkers 20 . This model, which relies on short range three-body interactions, reproduces structural and thermodynamic properties of water in an excellent manner and should be an ideal basis to represent not only aqueous ionic solutions but also water rich ion clusters that may play a role during mineralization. Ca2+ and CO32– ions are then represented as non-charged singlebead particles, i.e. electrostatic interactions and anisotropy of ions and ion-ion interactions need to be represented by appropriately tuned two- and three-body non-bonded potentials. In the following, these CG-beads will be referred to as particles. In contrast, aggregated, supramolecular species that are observed in the simulations will be referred to as clusters (in analogy to the experimentally frequently used term prenucleation clusters). Clusters of somewhat larger size and temporal stability will be referred to as nanoparticles. The large number of two- and three-body interactions between ions and water were parameterized in a sequential procedure illustrated in fig. 1.

5



Figure 1: Illustration of sequential CG parameterization strategy, with parameterized interactions shown next to the respective parameterization target. Black arrows indicate twobody interactions, gray three-body. A set of ion-ion interactions is determined to reproduce bulk phase properties of calcite, with pair interactions mostly determining the unit cell dimensions (stage 1) and three-body interactions the full crystal structure (stage 2). In the next stage, ion-water interactions and some of the direct ion-ion interactions are refined in iterative cycles against structural and thermodynamic properties of the aqueous solution (local solvation structure, ion-ion association, and the mineral-water interface). First, most of the ion-ion interactions were tuned against bulk phase properties of calcite, most importantly density and the calcite crystal structure. In a second step, ion-water 6


Page 6 of 25

Page 7 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


interactions and some of the direct ion-ion interactions were refined against structural and thermodynamic properties of the aqueous solution (obtained from atomistic reference simulations and experimental data), most importantly local solvation structure, ion-ion association, and finally the behavior of the mineral-water interface. Due to the combination of such diverse parameterization targets and the interdependencies of the parameters one needs to iteratively refine parameters by going through several optimization cycles. In consequence, one does not obtain a single CG model but rather a range of models with slightly different parameter combinations, illustrating once again the difficulty of developing a CG model for such a multi-component, multi-phase problem. We have selected three of these CG models which vary in the strength of water-calcium interactions and in the stability of the bulk phase, and which give us the possibility to study different mineralization phenomena. For a short summary which specifically highlights the differences between the models see table 1, the full parameter tables can be found in the SI. In the following, the parameterization will be discussed in more details. Table 1: Brief overview of the CG models used property Ca−mW 2-bodyCaCO3 3-bodyCaCO3 +3-bodyCa−mW −Ca

model1 0.55 LJ strong no

model2 3.0 buck weak yes

model3 5.0 buck weak yes

In a first step, ion-ion interactions were tuned to reproduce the bulk phase properties of CaCO3 (calcite). The Ca-CO3 interactions were modeled with a Lennard-Jones potential in model 1 and a Buckingham potential in model 2 and 3. For the Ca-Ca and CO3-CO3 interactions short-ranged repulsive Weeks-Chandler-Andersen pair potentials were used in model 1, in models 2 and 3 Buckingham potentials with a very strong repulsive term were employed. In order to reproduce the full crystal structure, these interactions were complemented by non-bonded three-body terms as energy penalty functions to introduce anisotropy in the system. Here, different types of three-body potentials were tested, starting from the 7



Page 8 of 25

Stillinger-Weber potential as used by Molinero 21 . However, the distorted octrahedral coordination of the C positions of carbonates around Ca2+ and vice versa (for an illustration see fig. S1b) can not be represented by a harmonic term as in the Stillinger-Weber potential. Therefore, a new type of three-body term needed to be introduced. Starting from the ansatz of a periodic cosine potential, as used in the DREIDING force field 22 , we complemented a Stillinger-Weber like three-body term with an angle-dependent term that was specifically modified to be capable of describing the angles between ion sites in calcite, in particular allowing for multiple minima, see eq. (3). γik σik γij σij ) exp( ) φ3 (rij , rik θijk ) = λij ij G(θijk ) exp( rij − aij σij rik − aik σik # " (f cos(θ))2 + b cos(2θ) + d G(θ) = g 1 − cos(4θ) + a exp − 2c2

(2) (3)

The function in eq. (3) was implemented in the form of eq. (4) as a tabulated three-body potential using the implementation of the polymorphic extension of the Stillinger-Weber potential by Zhou et al. 23 in LAMMPS. "

# 2 (f t) + b(2t2 − 1) + d G0 (t) = g 1 − (8t4 − 8t2 + 1) + a exp − 2c2

(4)

with: t = cos(θ) These three-body interactions were used as short range energy penalty functions in conjunction with the above mentioned Lennard-Jones or Buckingham two-body interaction potentials. The targets of the this first optimization step of these ion-ion pair and three-body potentials were the cell parameters, density, radial distribution function and the distance dependent angle-distribution in calcite (for an illustration of the latter target, see for example fig. S1b). Due to the large number of adjustable parameters, symmetries that are present in the calcite structure were enforced during parameter search and optimization. Neverthe8


Page 9 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


less, due to the vast number of possible parameter combinations, a conjugate-gradient based parameter optimization did not converge to a model that reproduced the calcite structure. Thus, the parameter search was carried out by a numerical sensitivity-analysis type algorithm. Here, one to four parameters were varied simultaneously in forward and backward steps while the others were kept unaltered. For each set of parameters, a simulation of a calcite bulk system was performed for 20000 steps and the standard errors of all target observables with respect to the calcite reference data was computed and the best set chosen. The parameter search was iteratively performed with decreasing step size until a minimum step size of 1 % per parameter. The bulk phase properties of the final models (on the level of the solid phase model2 and model3 are identical) after a simulation of 1.0 · 107 steps at 300 K and a subsequent energy minimization are presented in table 2. Figure 2 shows the distance-dependent CO3 Ca CO3 and Ca Ca Ca angle distributions found in the simulations in comparison to the calcite reference, showing that the models reproduce the full calcite structure. (Note that the shown data refer to the final models after the full cycle of also parameterizing in the solution phase and iterating between solid and solution phase as described below.)

9



Page 10 of 25

Figure 2: Angle vs. distance distributions for the symmetric CO3 Ca CO3 and Ca Ca Ca triangles (with the distance being the distance between the central and the two outer atoms. The reference data obtained from the perfect crystal positions in calcite are shown as red crosses, while the data from the CG simulations are shown as density plots (red indicating high, blue indicating low probability density).

Table 2: Selected properties of bulk & solution phase properties a [Å] b [Å] c [Å] Volume [Å3 ] Density [g/cm3 ]

experiment atomistica model1 Calcite bulk phase 4.990 24 4.934 5.070 24 4.990 4.938 5.070 24 17.062 17.082 16.562 367.961 24 360.356 368.713 2.710 24 2.767 2.705 Solution properties

model2

model3

5.008 5.008 16.463 357.615 2.788

5.008 5.008 16.463 357.615 2.788

Ca rCa-OW b [Å] 2.33-2.44 25 2.35 2.35 2.35 2.35 c 25 CN 6-10 7.33 7.15 7.32 7.52 CO3 rCCO3 -OW b [Å] 3.35 26 3.24 3.31 3.1 3.08 c 26 CN 9.1 12.0 ∼ 15 ∼ 15 ∼ 15 a 4 b atomistic reference simulations are performed with Raiteri model ; based on first peak in g(r); c coordination number (CN) of first solvation shell. For full technical details see SI. In the second stage of the parametrization, we focused on the interactions of ions in solution. These interactions are highly interdependent, therefore, an iterative cycle between 10


Page 11 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


different parameterization targets was needed. The Lennard-Jones parameters for the Ca mW, respectively CO3 mW, interactions were optimized against the position of the first peak in the ion water radial distribution function and the coordination number of the first hydration shell (CNion−mW ). A detailed overview of the solution properties is shown in table 2. In the three models, three different Ca mW interaction strengths were used with Ca−mW = 0.55, 3.0, 5.0

Kcal mol

resulting in CNCa−mW

values of 7.2, 7.3, 7.5 (experimental reference values are in the range of 6-10 25 while atomistic simulations yielded 7.2 4 ). Here, one should mention that the values obtained for the coordination numbers from integrating the radial distribution function in the simulation are not completely unambiguous. Here, an integration cutoff in the minimum after the first peak of the radial distribution function was used. Application of a switching function, as described by Raiteri et al. 4 , results in too small values for the coordination numbers in the case of the CG models due to the broadening of the peaks. The radial distribution functions obtained for the three CG models are shown in fig. 3. The discrepancies in the position of the second peak are due to the (mW) water model used and could not be shifted by adjusting the two-body interactions between the ions and the water beads.

(a)

(b)

Figure 3: Radial distribution function of a single ion in water: a) Ca - mW/OW; b) CO3 mW/OW Since the interaction between the Ca2+ and CO32– ions is a major driving force of the ion 11



association and subsequent aggregate formation, the CG Ca-CO3 interaction was tuned after the ion-water interactions had been fixed. To this end the potential of mean force (PMF) between a Ca2+ and CO32– ion was calculated for the atomistic model as a reference. The Lennard-Jones potential in model1 and the Buckingham potential in model2 and model3 (which had before been determined based on the solid calcite phase) were readjusted so that the overall association strength (approximate width and depth) of the atomistic PMF was reproduced (see fig. 4). In this level of CG representation the clear solvent structuring of the ion pairing with defined solvent-shared and solvent-separated states is smeared out, which is in part due to the solvent structure around the ions (see fig. 3). Note that this modification of ion-ion interaction potentials required readjusting (re-scaling) all other ion-ion interactions to preserve the correct solid-phase properties of calcite. This turned out to be possible due to the fact that the models for the solid phase mostly rely on structural properties and are therefore scalable to some extent.

Figure 4: Potential of mean force for the distance between Ca and CO3 as a single ion pair in solution. In addition to the ion-water structure and the ion-ion PMF in solution, the water structure on top of the thermodynamically stable (104) surface of calcite was evaluated. Here it turned out that in several CG simulations with Ca-mW interaction strengths in the range of Ca−mW = 3.0 to 5.0

Kcal , mol

a penetration of mW particles into the mineral phase was observed.

This indicates that the barriers for water beads passing the first crystal layer is compara12


Page 12 of 25

Page 13 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tively low due to the spherical shape of the CO3 beads in contrast to the disc shaped ions in the atomistic model. To prevent this coarse-graining artifact, an additional Stillinger-Weber three-body term for the Ca-mW-Ca interactions, with an equilibrium angle of 0°, was added to the models. This term leads to an additional repulsive contribution at the surfaces, in particular upon penetrating the surfaces, while does not affect the solution properties of the individual ions and the oppositely charged ion pairs. It only minimally interferes with higher ion-ion species (in dilute solution) by effectively slightly increasing the Ca-Ca repulsion. This Ca-mW-Ca interaction term was parameterized in such a way that the lateral distribution of the water molecules on the calcite surface better reproduces the atomistic reference data (and no interpenetration artifact is observed any more). The resulting water structure both perpendicular and in the first layer parallel to the calcite surface is shown in fig. S2 and fig. S3, respectively. Note that while no higher order solution species or higher-dimensional potentials of mean force between ions and water beads were used during the parametrization of the models, data regarding such properties have been evaluated for validation purposes and will be shown in the results section where the formation of aggregates from solution is discussed. The aim to derive a CG model that is transferable from a very dilute ion solution to the thermodynamically stable mineral phase has motivated the presented parameterization process. Given the known transferability and representability problems of CG models and the challenges involved in simulating a multiphase system, compromises in reproducing the – potentially conflicting – reference properties need to be made and an easy, unique solution for a model may remain elusive. It was found before that the use of bottom-up coarse graining approaches – which are designed to generate CG potentials that very accurately reproduce the given targets – may result in particularly severe limitations with regards to predicting thermodynamic properties and the behavior at thermodynamic state points that are not part of the reference data. In particular it is very challenging to devise bottom-up models for applications at a wide range of state points and concentrations or across phase 13



transitions or a phase separation 13,14,27–29 . Thus, in the present work we have not aimed for a single set of CG parameters that very accurately reproduces a specific set of target properties, but we have rather used those target properties to identify parameter ranges and those types of interactions that are important to simultaneously reasonably well reproduce of the two reference state points. We present now results for three different resulting CG models which were chosen such that the effect of certain interaction types on the aggregation and mineralization process and the intermediates that form can be investigated. Model 1 relies on CaCO3 potentials with hard repulsive interaction terms (Lennard-Jones and comparatively strong three-body terms) which lead to a stronger ordering tendency compared to model2 and 3. The latter two models use less repulsive CaCO3 terms (Buckingham potentials that can be combined with weaker three-body interactions) allowing for a softer behavior of the mineral. The three models further differ in their Ca-mW interaction strength, in increasing order from model 1 to model 3, allowing to study the effect of the solvation strength of ions on the formation of ionic clusters from solution.

Results and Discussion Aggregation of CaCO3. With these three models, we simulated an ion solution at a concentration of ∼ 0.35

mol L

(400 CaCO3 formal units and 63.200 water beads in a cubic box

of 124 Å side length) for 100 or 250 million time steps (formally corresponding to 500 and 1250 ns, not accounting for the typical speedup gained from the smoothness of CG potentials 30 ). Details regarding the simulation protocol are given below and in the SI. Figure 5 shows snaphots that illustrate the state of the system at different times during the simulations and the sizes, crystallinity and wetness of the clusters. Initially, in all three models small clusters of a few ions form which are similar to the ones found in atomistic simulations by Demichelis et al. 31 . A more detailed comparison of the local structure of these solution species is shown in the SI, where the geometry (angles versus distances) of different three-ion

14


Page 14 of 25

Page 15 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


species is compared across the CG and atomistic simulation models. Next, these solution species merge into small amorphous clusters with differing water content depending on the CG model. These subsequently merge into bigger amorphous clusters. Here, the behavior found in the three models starts to diverge. In the case of model1, which relies on LJ-based pair potentials, strong intra-mineral three-body terms, and comparatively weak water-Ca interactions, the amorphous clusters, which have a very low water content at all times, ultimately crystallize into a perfect calcite structure. In the case of model2 and model3, which both use softer Buckingham potentials and stronger Ca-water interactions, no transformation into calcite is observed. Due to the differences between these two models in the strength of the Ca-water interaction, the water content of the clusters is drastically different. For model2, the amorphous clusters initially contain water molecules that are gradually expelled with increasing nanoparticle size until finally a large, amorphous nanoparticle is obtained. In contrast, the water content of the clusters remains high at all times for model3, which has the strongest water-Ca interactions. The cluster has a water content of 30.5 wt% in the final structure, which is slightly above experimentally reported values between 15 wt% 32 and 24 wt% 33,34 . This qualitative description of the different pathways is corroborated by the data presented at the bottom of fig. 5, where time series of the size and crystallinity of the respective largest clusters and the water content of all clusters are presented.

15



Figure 5: Illustration of CaCO3 aggregation observed with three different coarse-grained models. Snapshots of simulations are shown for every model, with bulk water molecules omitted for clarity. Insets show the following properties of the respective snapshots: N is the number of Ca & CO3 beads in the respective largest cluster, χ is the crystallinity of the largest cluster (using local Steinhardt parameters 35 , see also computational details), H2O is the water content of the clusters in weight percent. In the graphs at the bottom, timeseries of the following properties are shown: the number of atoms, N, in the largest cluster (green lines), the number of crystalline atoms, #lq6, in the largest cluster (red lines), and the number of water molecules, #H2O, in the clusters (blue line). 16


Page 16 of 25

Page 17 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In addition, further simulations at different system sizes and ion concentrations were carried out to investigate whether the different pathways depend on the degree of supersaturation. To this end, in a second set of simulations the concentration was decreased by a factor of 8 while keeping the total number of ions constant (resulting in a total system size of 512000 beads) and in a third set the total number of beads was kept constant (at 64000) while the number of ions was decreased. A full description of all simulated systems and resulting structures is presented in the SI. In short, they consistently confirm the reported mechanisms and structures (crystallinity, water content) for the three models as described above, irrespective of the degree of supersaturation or system size. This in particular shows that the amorphous states (metastable intermediates) which are found for model2 and 3 not only occur due to very high supersaturation but also occur at rather low ion concentrations. Furthermore, model1 apparently somewhat overestimates the crystallinity of nanoparticles since both amorphous and crystalline clusters are found at small sizes of around 2.0-2.4 nm (i.e. smaller than the previously reported stability crossover at a size of ∼4 nm for crystalline nanoparticles 8,10 . Figure 6 shows the coexistence of such small crystalline structures and amorphous structures and a merging event which at first results in a particle with an amorphous and a crystalline domain. The crystalline domain seems to induce crystallization in the amorphous one resulting in a single crystalline nanoparticle, see cyrstallinity graph (graph B in fig. 6). In other simulations also the merging of two crystalline nanoparticles was observed, in this case twinned crystals were found (data not shown).

17



Figure 6: Merging event for model1 for 400 CaCO3 formal units in a 64000 bead system. Left panel: timeseries of the number of atoms, N, in the largest cluster (green lines), the number of crystalline atoms, #lq6, in the largest cluster (red lines), and the number of water molecules, #H2O, in the clusters (blue line). Right panel: Average crystallinity (per ion) in the largest cluster as function of time. Note: A different simulation is shown compared to the one in fig. 5. Note that still all concentrations investigated here are higher than the experimental concentrations where nucleation occurs. All simulations that are shown are in the spinodal regime. In the future, the CG models will allow us to study systems of larger sizes and to probe lower saturation states. We will perform extensive comparisons between atomistic and CG simulations at different concentrations in order to better understand where the models differ regarding the aggregation/phase separation thermodynamics. This will aid the future development of the CG models towards a better transferability across concentrations and phases. Nanoparticles. Given the very different types of nanoparticles that formed in the different models, the question arises whether for model2 and 3 the simulations are merely arrested in amorphous metastable intermediates towards calcite or whether these nanoparticles are stable. Therefore, MD simulations were started from solvated calcite nanoparticles of two sizes 18


Page 18 of 25

Page 19 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(22 Å and 48 Å) with (104) surfaces (i.e. the thermodynamically stable surface of calcite). The sizes were chosen since atomistic simulation studies had predicted that nanoparticles smaller than ∼4 nm should become amorphous while larger particles should be crystalline 8,10 . The crystallinity was monitored via the local q6 Steinhardt parameter for each atom (see fig. 7 and fig. S6). In all models, the larger nanoparticle remained crystalline. Remarkably, in the case of model3, the smaller nanoparticle underwent a spontaneous transition to an amorphous structure, in good accordance with the atomistic data (see fig. 7).

Figure 7: Histogram of the local Steinhardt parameter q6 for calcite nanoparticles of 22 Å (panel a) and 48 Å (panel b) using model3 as a function of simulation time (108 steps, formally corresponding to 500 ns).

Conclusions In summary, we can show that the use of short-range three-body interactions allow us to generate CG models that are able to represent calcite in water. Our models can reproduce structural and thermodynamic properties of the calcite bulk phase and the ions in solution fairly well, while gaining both computational efficiency (between a factor of 20 to 40 depending on the system size) and speeding up the dynamics of the different processes involved. (Benchmarking and a detailed discussion of all aspects of computational efficiency and speedup of the CG dynamics compared to atomistic simulations can be found in the SI). 19



Thus, the models allow to cheaply access bigger time scales and larger system sizes compared to what is easily accessible to atomistic simulations. This approach allows us now to simulate and study the aggregation and crystallization of CaCO3 in solution. The usage of these CG models provides easy test systems for complex mechanisms, which can be useful in method developing for more detailed or complex systems. By varying CG interactions we could observe the formation of qualitatively different ion clusters and nanoparticles which reflect the spectrum of experimentally proposed intermediates and mineralization mechanisms very well. Thus the suite of CG models will allow us to investigate different pathways towards mineral formation and address general questions about nucleation mechanisms in the future that can be transferred to other minerals. Furthermore, we have demonstrated that it is possible to generate CG models for minerals by describing the spatial environment of the constituents in the solid phase with specifically tuned non-bonded three-body interactions that are able to represent local anisotropy in a coordination shell and angle potentials with multiple minima. This can be extended to other minerals and other geometric motives and future investigations of the formation of complex solid materials.

Acknowledgement This work was supported by grants from the German Research Foundation (SFB1214). The computations were performed on computational resources funded within the bwHPC program by the state of Baden-Württemberg.

Supporting Information Available • supplemental.pdf Containing the full force field parameter set, description of the parameterization procedure and targets, monitored properties and additional information about the stability of nanoparticles and concentration dependency of the nucleation process. 20


Page 20 of 25

Page 21 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References (1) Weiner, S.; Sagi, I.; Addadi, L. Choosing the Crystallization Path Less Traveled. Science 2005, 309, 1027–1028. (2) Ma, Y.; Weiner, S.; Addadi, L. Mineral Deposition and Crystal Growth in the Continuously Forming Teeth of Sea Urchins. Adv. Funct. Mater. 2007, 17, 2693–2700. (3) Gebauer, D.; Kellermeier, M.; Gale, J. D.; Bergström, L.; Cölfen, H. Pre-Nucleation Clusters as Solute Precursors in Crystallisation. Chem. Soc. Rev. 2014, 43, 2348–2371. (4) Raiteri, P.; Demichelis, R.; Gale, J. D. Thermodynamically Consistent Force Field for Molecular Dynamics Simulations of Alkaline-Earth Carbonates and Their Aqueous Speciation. J. Phys. Chem. C 2015, 119, 24447–24458. (5) Sebastiani, F.; Wolf, S. L. P.; Born, B.; Luong, T. Q.; Cölfen, H.; Gebauer, D.; Havenith, M. Water Dynamics from THz Spectroscopy Reveal the Locus of a LiquidLiquid Binodal Limit in Aqueous CaCO3 Solutions. Angew. Chem. Int. Ed. 2017, 56, 490–495. (6) Henzler, K.; Fetisov, E. O.; Galib, M.; Baer, M. D.; Legg, B. A.; Borca, C.; Xto, J. M.; Pin, S.; Fulton, J. L.; Schenter, G. K. et al. Supersaturated Calcium Carbonate Solutions are Classical. Sci. Adv. 2018, 4 . (7) Smeets, P. J. M.; Finney, A. R.; Habraken, W. J. E. M.; Nudelman, F.; Friedrich, H.; Laven, J.; De Yoreo, J. J.; Rodger, P. M.; Sommerdijk, N. A. J. M. A Classical View on Nonclassical Nucleation. PNAS 2017, (8) Raiteri, P.; Gale, J. D. Water Is the Key to Nonclassical Nucleation of Amorphous Calcium Carbonate. J. Am. Chem. Soc. 2010, 132, 17623–17634. (9) Quigley, D.; Freeman, C. L.; Harding, J. H.; Rodger, P. M. Sampling the Structure

21



of Calcium Carbonate Nanoparticles with Metadynamics. J. Chem. Phys. 2011, 134, 044703. (10) Bano, A. M.; Rodger, P. M.; Quigley, D. New Insight into the Stability of CaCO3 Surfaces and Nanoparticles via Molecular Simulation. Langmuir 2014, 30, 7513–7521. (11) De La Pierre, M.; Raiteri, P.; Stack, A. G.; Gale, J. D. Uncovering the Atomistic Mechanism for Calcite Step Growth. Angew. Chem. Int. Ed. 2017, 56, 8464–8467. (12) Quigley, D.; Rodger, P. M.; Freeman, C. L.; Harding, J. H.; Duffy, D. M. Metadynamics Simulations of Calcite Crystallization on Self-Assembled Monolayers. J. Chem. Phys. 2009, 131 . (13) Mukherjee, B.; Delle Site, L.; Kremer, K.; Peter, C. Derivation of Coarse Grained Models for Multiscale Simulation of Liquid Crystalline Phase Transitions. J. Phys. Chem. B 2012, 116, 8474–8484. (14) Dunn, N. J. H.; Foley, T. T.; Noid, W. G. Van der Waals Perspective on CoarseGraining: Progress toward Solving Representability and Transferability Problems. Acc. Chem. Res. 2016, 49, 2832–2840. (15) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1 – 19. (16) Shinoda, W.; Shiga, M.; Mikami, M. Rapid Estimation of Elastic Constants by Molecular Dynamics Simulation Under Constant Stress. Phys. Rev. B 2004, 69, 134103. (17) Tuckerman, M. E.; Alejandre, J.; López-Rendón, R.; Jochim, A. L.; Martyna, G. J. A Liouville-Operator Derived Measure-Preserving Integrator for Molecular Dynamics Simulations in the Isothermal-Isobaric Ensemble. J. Phys. A: Math. Gen. 2006, 39, 5629.

22


Page 22 of 25

Page 23 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(18) Hockney, R. W.; Eastwood, J. W. Computer Simulation Using Particles, 1st ed.; Bristol ; Philadelphia : A. Hilger, 1988. (19) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New Feathers for an old Bird. Comput. Phys. Commun. 2014, 185, 604 – 613. (20) Molinero, V.; Moore, E. B. Water Modeled As an Intermediate Element between Carbon and Silicon. J. Phys. Chem. B 2009, 113, 4008–4016. (21) Stillinger, F. H.; Weber, T. A. Computer Simulation of Local Order in Condensed Phases of Silicon. Phys. Rev. B 1985, 31, 5262–5271. (22) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. DREIDING: a Generic Force Field for Molecular Simulations. J. Phys. Chem. 1990, 94, 8897–8909. (23) Zhou, X. W.; Foster, R., M. E. Jones; Yang, P.; Fan, H.; Doty, F. P. A Modified Stillinger-Weber Potential for TlBr, and Its Polymorphic Extension. J. Mater. Sci. Res. 2015, 4, 15–32. (24) Graf, D. L. Crystallographic Tables for the Rhombohedral Carbonates. Am. Mineral. 1961, 46, 1283–1316. (25) Ohtaki, H.; Radnai, T. Structure and Dynamics of Hydrated Ions. Chem. Rev. 1993, 93, 1157–1204. (26) Kameda, Y.; Sasaki, M.; Hino, S.; Amo, Y.; Usuki, T. Neutron Diffraction Study on the Hydration Structure of Carbonate Ion by Means of 12C/13C Isotopic Substitution Method. Physica B 2006, 385-386, Part 1, 279 – 281. (27) Villa, A.; Peter, C.; van der Vegt, N. F. A. Transferability of Nonbonded Interaction Potentials for Coarse-Grained Simulations: Benzene in Water. J. Chem. Theory Comput. 2010, 6, 2434–2444.

23



(28) Potestio, R.; Peter, C.; Kremer, K. Computer Simulations of Soft Matter: Linking the Scales. Entropy 2014, 16, 4199–4245. (29) Lu, J.; Qiu, Y.; Baron, R.; Molinero, V. Coarse-Graining of TIP4P/2005, TIP4P-Ew, SPC/E, and TIP3P to Monatomic Anisotropic Water Models Using Relative Entropy Minimization. J. Chem. Theory Comput. 2014, 10, 4104–4120. (30) Fritz, D.; Koschke, K.; Harmandaris, V. A.; van der Vegt, N. F. A.; Kremer, K. Multiscale Modeling of Soft Matter: Scaling of Dynamics. Phys. Chem. Chem. Phys. 2011, 13, 10412–10420. (31) Demichelis, R.; Raiteri, P.; Gale, J. D.; Quigley, D.; Gebauer, D. Stable Prenucleation Mineral Clusters are Liquid-Like Ionic Polymers. Nat. Commun. 2011, 2, 590. (32) Ihli, J.; Wong, W. C.; Noel, E. H.; Kim, Y.-Y.; Kulak, A. N.; Christenson, H. K.; Duer, M. J.; Meldrum, F. C. Dehydration and Crystallization of Amorphous Calcium Carbonate in Solution and in Air. Nat. Commun. 2014, 5, 3169. (33) Michel, F. M.; MacDonald, J.; Feng, J.; Phillips, B. L.; Ehm, L.; Tarabrella, C.; Parise, J. B.; Reeder, R. J. Structural Characteristics of Synthetic Amorphous Calcium Carbonate. Chem. Mater. 2008, 20, 4720–4728. (34) Radha, A. V.; Forbes, T. Z.; Killian, C. E.; Gilbert, P. U. P. A.; Navrotsky, A. Transformation and Crystallization Energetics of Synthetic and Biogenic Amorphous Calcium Carbonate. PNAS 2010, 107, 16438–16443. (35) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Icosahedral Bond Orientational Order in Supercooled Liquids. Phys. Rev. Lett. 1981, 47, 1297–1300.

24


Page 24 of 25

Page 25 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry

25


Coarse-Grained Simulation of CaCO3 Aggregation and Crystallization

Recommend Documents