Derivation of an Accurate Force-Field for Simulating the Growth of

J. Phys. Chem. C 2010, 114, 5997–6010

5997

Derivation of an Accurate Force-Field for Simulating the Growth of Calcium Carbonate from Aqueous Solution: A New Model for the Calcite-Water Interface Paolo Raiteri,*,† Julian D. Gale,*,† David Quigley,‡ and P. Mark Rodger§ Nanochemistry Research Institute, Department of Chemistry, Curtin UniVersity of Technology, GPO Box 1987, Perth, WA 6845, Australia, Department of Physics and Centre for Scientific Computing, UniVersity of Warwick, Gibbet Hill Road, CoVentry CV4 7AL, United Kingdom, and Department of Chemistry and Centre for Scientific Computing, UniVersity of Warwick, Gibbet Hill Road, CoVentry CV4 7AL, United Kingdom ReceiVed: NoVember 19, 2009; ReVised Manuscript ReceiVed: January 14, 2010

The performance of existing force-field models for the calcium carbonate - water system has been critically assessed with particular reference to the thermodynamic consequences. It is demonstrated that all currently available parametrizations fail to describe the calcite-aragonite phase transition, and the free energies of solvation for the calcium cation are also considerably in error leading to a poor description of the dissolution enthalpy for calcite. A new force-field, based on rigid carbonate ions, has been developed that corrects these deficiencies and accurately describes the thermodynamics of the aqueous calcium carbonate system within molecular dynamics simulations. Not only does this new model lead to quantitative changes in the properties of the calcite (101j4) surface in contact with water, but also significant qualitative differences. With this more accurate model it is found that calcium ions do not adsorb at the pristine basal plane of calcite, while carbonate ions only weakly bind. Carbonate diffusion across the surface is found to occur only when the anion is solvent separated from the underlying surface, with there being an equal tendency to readsorb or migrate into the bulk liquid. Introduction Carbonates are important minerals in the field of geochemistry from a number of perspectives. In the area of geosequestration, the reaction of dissolved carbon dioxide with metal-containing compounds forms the basis of mineral trapping via carbonate formation.1 In the presence of living organisms, the crystallization of carbonates can lead to complex structures, via the process of biomineralization,2 with a variety of important functions in the natural world. Conversely, the crystallization of carbonates can be problematic when induced by polymer membranes, leading to fouling during desalination.3 Given the significance of the above processes, there is a desire to fully understand the chemistry that occurs during the crystallization of carbonates from aqueous solution or conversely their dissolution. Considerable information concerning the structure of the carbonate-water interface4 and growth mechanisms has been obtained by experiment.5,6 In particular, the advent of atomic force microscopy has allowed dynamical investigations of in situ growth to be performed under a variety of conditions.7 While such studies can provide a wealth of valuable details, there are limits to the resolution achievable and the certainty of the interpretation, and therefore computer simulations of the aqueous-mineral interface offer the potential to be a useful, complementary technique to achieve further insights. Because of the complexity of the interfacial region and the need to sample many possible configurations, such computa* To whom correspondence should be addressed. E-mail: [email protected]; [email protected]. † Curtin University of Technology. ‡ Department of Physics and Centre for Scientific Computing, University of Warwick. § Department of Chemistry and Centre for Scientific Computing, University of Warwick.

tional studies will require the use of molecular dynamics. Despite recent advances that make large-scale quantum mechnical simulations feasible,8 the challenging nature of studying growth processes at interfaces precludes the realistic use of these methods for the foreseeable future. Consequently, the key to being able to perform reliable simulations of the aqueouscarbonate interface is the existence of an accurate force-field model to describe this system. It is the contention of this article that all currently available models have deficiences and therefore a new force-field is essential. It will be demonstrated that this is not just a matter of small quantitative discrepancies but that inadequacies in a model can lead to a qualitatively incorrect picture of this important class of interface. Although there are many interesting families of carbonate minerals, we choose to focus on arguably the most studied exemplar; calcium carbonate (CaCO3). This material is widely studied due to its relative abundance and significance in terms of occurrence and usage. Calcium carbonate commonly occurs as three distinct anhydrous minerals, namely calcite, aragonite, and vaterite, as well as further polymorphs at high pressure.9 Calcite and aragonite have well-characterized structures, with the former transforming to the latter at moderate pressure. Vaterite is metastable with respect to calcite and aragonite and will readily transform in contact with solution in many environments. However, it can be important during nucleation and growth as an intermediate phase between amorphous calcium carbonate and the formation of the other two more common polymorphs. In the last few decades several force-fields to study CaCO3, both in anhydrous and aqueous conditions, have been developed. Arguably one of the first attempts to create a force-field for calcium carbonate was that of Yuen et al. in 1978,10 who analyzed the charge distribution with the aim of extracting lattice energies for calcite and aragonite. However, it has not been

10.1021/jp910977a  2010 American Chemical Society Published on Web 02/18/2010

5998

J. Phys. Chem. C, Vol. 114, No. 13, 2010

widely applied to more general properties. The year 1992 saw considerable activity within this area, with three separate forcefields being published due to Dove et al.,11 Pavese et al.,12 and Jackson and Price,13 in chronological order of manuscript acceptance. All three studies took similar approaches, in that the carbonate group was treated as a molecular entity with harmonic angle-bending and improper torsions to supplement the ionic model. The work of Pavese also considered the use of a shell model (SM) description for oxygen, in addition to a rigid ion model (RIM). Pavese et al.14 also published a revised shell model four years later, in which a Morse potential replaces the Buckingham form for the intramolecular C-O interaction, to enable them to describe the thermal expansion behavior of calcite within the quasi-harmonic approximation. Subsequent to this period, there have been several other attempts to refine the potential model. Fisler et al.15 retained the shell model approach of Pavese et al. but changed the description of the carbonate group to be in the spirit of molecular mechanics (i.e., distinct interatomic potentials for intra- versus intermolecular interactions with cut-offs determined by bonding connectivity). Further refinements were later made to this model by Rohl et al.16 to yield one of the most refined force-fields for structural and vibrational properties. This model was able to verify the (2 × 1) reconstruction of the (101j4) surface of calcite under dry conditions.5 Similarly, Hwang et al.17 have also explored a conventional molecular mechanics representation based on the Dreiding force-field. A more sophisticated approach was taken by Thackeray and Siders,18 who used a four-state valence bond approach to capture the distortions and polarizability of the carbonate anion. However, this has also limited its wider use. Recently, the attention has shifted from deriving improved models for calcium carbonate alone to creating force-fields that are capable of describing the interface of this material with water, and even organics, with a view to addressing biomineralization. de Leeuw and Parker19 developed a shell model description of water and combined this with the shell model of Pavese et al. to consider the hydrated surfaces of calcite, aragonite and vaterite.20 Kerisit and Parker21 later introduced a revised water model in an attempt to overcome the tendency of the water to freeze at room temperature.22 Arguably the best force-field to date for the simulation of calcium carbonate in water is that of Bruneval et al.23,24 Here they blended parameters from both of the models of Rohl et al.16 and Archer et al.25 (both of which are a revision of the Fisler et al. model15), while adjusting a few further quantities and combining it with a welltested description of water in the form of the polarizable SWM4NDP force-field.26 For the calcium and carbonate to water interactions they maintained the parameters of de Leeuw and Parker. One further force-field is worthy of note, which is that of Freeman et al.,27 where parameters were developed with the specific aim of being able to study biomineralization. Again several existing models were combined, with TIP3P being employed for water and de Leeuw and Parker for the crossspecies interactions. For calcium carbonate itself, the revised model of Pavese et al. from 1996 was chosen, though for reasons of efficiency the shells were omitted. Furthermore, the oxygen-oxygen interaction was modified by removing the dispersive attraction and adjusting the repulsive parameters to compensate. As described above, there has been considerable effort directed toward the development of the “ultimate” force-field for calcium carbonate. However, the focus of the majority of these previous studies has been to obtain an accurate description

Raiteri et al. of the solid-state properties of this substance, with particular emphasis on the structure, mechanical hardness and phonon spectra. In the present work, we stress the importance of thermodynamics while developing a new force-field designed to specifically model the crystallization and growth of calcium carbonate, in contrast to previous studies. Based on this, we explore the consequences for the free energy of ion adsorption and water structure at the most widely studied (101j4) surface of calcite. A particular goal of future simulations will be to reproduce and understand the selection of the metastable aragonite polymorph during nucleation and growth. We therefore examine the location of the calcite-aragonite phase boundary to assess the extent to which the relative energetics of these two phases have been captured. Methodology Evaluation of Existing Force-Fields for Calcium Carbonate. The objective of the present work is to obtain a force-field that accurately describes the thermodynamics of the calcium carbonate system. Hence, we need to examine how the existing force-fields fare against this criterion. Experimentally, the thermodynamics of the phase transitions of CaCO3 have been widely studied, though with considerable variation in the results obtained. Wolf et al.28,29 have performed a series of evaluations to produce recommend values that are arguably the most precise available at present. The calcite-aragonite phase transition is a particularly interesting one, since aragonite has a lower enthalpy than calcite, but the entropic contribution makes the latter phase more stable at ambient conditions. The phase boundary for the calcite-aragonite transition as a function of pressure has received particular attention, having proved difficult to define precisely due to its highly nonlinear nature. Estimates of the transition pressure at ambient temperature vary between 0.3 and 0.4 GPa.30-32 In order to examine how the existing force-fields perform for the thermodynamics of the calcite-aragonite phase transition at 298 K, we have optimized the structures of both phases with respect to the internal energy using a range of models. Given that the majority of manuscripts do not specify the cutoff parameters used, there will be slight discrepancies from the numbers specified in the original papers. Furthermore, in the case of the model of Archer et al.,25 the cutoff had to be increased to 25 Å from the value in the manuscript in order to achieve smooth convergence. For the Jackson and Price model, we further assume that the parameters S and n for the torsion term are 1 and 2, respectively, as per the model of Dove et al., since the parameters are omitted from the original manuscript. Based on the optimized structures, the free energy is evaluated within the quasi-harmonic approximation for a temperature of 298 K using a Monkhorst-Pack mesh to numerically integrate the phonons across the Brillouin zone.33 Shrinking factors of 6 × 6 × 4 and 6 × 6 × 2 were employed for aragonite and calcite, respectively. For selected models the relaxation of the structures with respect to the free energy was examined, though this only led to changes in the free energy differences of less than 0.1 kJ/mol. It should be noted that the models of Dove et al.11 and Freeman et al.27 both exhibit imaginary phonon frequencies at the gamma point for aragonite, corresponding to rotation of the carbonate group away from the experimental symmetry. If we examine how the existing force-fields perform, as shown in Table 1, the accurate description of thermodynamic quantities can be seen to be a major weakness, despite the generally good description of structural and mechanical properties shown in the original papers. If we consider the enthalpy differences, then

Simulating the Growth of Calcium Carbonate

J. Phys. Chem. C, Vol. 114, No. 13, 2010 5999

TABLE 1: Calculated Enthalpy and Free Energy Differences (at 298 K) between Calcite and Aragonite for Literature Force-Fields for CaCO3a force-field 11

Dove et al. Jackson and Price13 Pavese et al. I RIM12 Pavese et al. I SM12 Freeman et al.27 Rohl et al.16 Archer et al.25 Bruneval et al.23 experiment28

∆H [kJ mol-1]

∆G [kJ mol-1]

-12.67 4.32 -5.36 2.74 -7.33 -3.75 1.76 4.82 0.44

-12.54 1.01 -6.34 0.56 -6.91 -3.98 -0.47 3.24 -0.84

a Here a positive value indicates that aragonite is more stable than calcite. Free energy differences are evaluated based on the structure optimized with respect to the enthalpy.

there is no consistent prediction that aragonite has a lower value; indeed nearly all force-fields are in error for this quantity by an order of magnitude with respect to experiment. Consequently, when the free energy is evaluated there is no possibility that a reversal of stability can occur. Clearly a new parametrization is required to rectify this situation if it is to be feasible to describe the competitive nucleation of different polymorphs. As an aside, the question may be posed as to whether quantum mechanical calculations could improve upon the force-field description of the calcite-aragonite energy difference. In order to test this, the structures of calcite and aragonite have been optimized at constant pressure using density functional theory within the linear-scaling SIESTA methodology.8 Here the Kohn-Sham equations are solved using a strictly localized pseudoatomic orbital basis set in combination with a normconserving pseudopotential representation of the potential due to the core electrons and nuclei. In the present calculations, a small core pseudopotential was used for Ca, such that the 3s and 3p electrons are explicitly included in the valence. Softconfined basis functions were used with a potential of 100 Ry and an inner radius of 0.95 times the radial cutoff, which was set to 6 au for Ca and C, while a larger value of 7 au was required for O. Both double-ζ polarized (DZP) and triple-ζ, double polarized (TZ2P) basis sets were used with a split norm parameter of 0.15. An auxiliary Cartesian mesh was used to expand the density with a mesh cutoff of 600 Ry. The Brillouin zone was integrated using a Monkhorst-Pack mesh with a K-grid cutoff of 16 Å. Calculations were performed for both the PBE34 and PBEsol35 functionals. Our recent work36 on the phases of iron sulphide has demonstrated that the new generation of GGA functionals generated specifically for the solid-state, such as PBEsol,35 Wu and Cohen37 and AM05,38 are much more accurate for relative phase stability, and this appears also to be the case for calcium carbonate. Based on the widely used PBE functional, the transition from aragonite to calcite is found to have an exothermic enthalpy of -9.3 kJ mol-1 at the DZP level of basis set, while the PBEsol yields a smaller value of -2.8 kJ mol-1. However, when using the PBEsol functional with a TZ2P basis set aragonite becomes more stable than calcite by 4.9 kJ mol-1 in terms of the enthalpy, in accord with experiment. Although this value requires further correction for the zero point energy, it indicates that density functional theory is likely to yield a calcite-aragonite phase boundary that is in error by as much as current force-field models, even for the best available GGA functionals. Therefore, we can conclude that empirical repa-

rameterization of force-field models will be more appropriate that the use of density functional theory results for the present system. Derivation of an Improved Force-Field for the AqueousCarbonate System. Fitting Procedure for Mineral Phases. Before attempting any refitting of a force-field against the energy difference between two polymorphs, it is important to consider what the correct quantity to fit against is. For two solid polymorphs at standard conditions, an accurate quasi-harmonic lattice dynamical description should be able to capture both the enthalpy and entropy differences, thereby leading to a correct prediction of the free energy. While this is possible, in principle, it is a very demanding requirement. Of the force-fields examined, only the Pavese rigid ion model comes close to capturing the entropic contribution correctly (despite not offering the best description of the carbonate anion vibrational modes), while the model of Archer et al. exhibits the correct change of sign, but with the wrong magnitude. Although lattice dynamics is the correct method to use for computing the free energy differences of solids at moderate temperature, our ultimate objective here is study the interface between carbonate and aqueous systems, which requires the use of molecular dynamics. It is therefore essential to consider the effects of this scenario on the prediction of the thermodynamics. Because molecular dynamics is a classical technique, there is no quantization of the nuclear motions and thus no zero point vibrational contribution either. This significantly changes the nature of the free energy contributions in the solid state. To analyze the consequences of the use of molecular dynamics for the solid-state, let us assume that the system will approximate one where the equipartition theorem is applicable (which is effectively the same as allowing the energy level spacing in the harmonic oscillator tend to zero and neglecting the zero point energy). In this case the entropy at a temperature, T, for a vibrational frequency ω, is given by

(

( ))

S(T) ) kB 1 - ln

hω kBT

(1)

which is different to the equivalent expression for the quantized lattice dynamical case. The above expression can be evaluated for the phonons of both calcite and aragonite, allowing for integration across the Brillouin zone. Based on the model of Rohl et al.,16 this leads to an entropy of transition of only 0.164 J mol-1 K-1, as compared to the recommended experimental value of 4.3 J mol-1 K-1.28 If the intramolecular vibrations of the carbonate anion are excluded, as would be the case for a rigid molecule description, then this is further reduced to 0.0026 J mol-1 K-1. The reason for this large reduction is that the low frequency phonon modes are dispersed around the value of ambient thermal energy at room temperature, leading to strong cancelation of positive and negative logarithmic contributions. Based on this analysis, it is estimated that the use of molecular dynamics will lead to an underestimation of the entropy of transition, potentially by orders of magnitude for the model we will introduce later. Starting from the perspective that the most important objective is that the free energy difference between calcite and aragonite be correct at 298.15 K, such that calcite is indeed the most stable polymorph under ambient conditions, we have chosen to use the experimental free energy difference as an observable for fitting the enthalpy of transition within a static calculation. While determining the free energy difference between calcite and aragonite by molecular dynamics is possible, in principle, this is a demanding undertaking and not well suited to force-field

6000


Raiteri et al.

derivation. Given the likelihood that the entropic contribution would be negligible within this technique, we argue that the proposed approach is justified, as well as considerably more efficient. This is further validated by the Monte Carlo calculation of the calcite-aragonite phase equilibrium in our model, presented later. In reparameterizing a force-field for the accurate description of CaCO3 polymorphs, we closely follow the procedure previously described by Fisler et al. and Rohl et al. Here the relaxed fitting methodology39 is used within the program GULP40 to perform a least-squares refinement of the force-field against the structure and mechanical properties of calcite and aragonite, as well as the energy difference between these two phases (in contrast to previous work). Phonon stability was also enforced by including the frequency of the lowest mode at the gamma point as an observable. In addition, the lowest vibrational mode at the (1/2, 1/2, 0) point in the Brillouin zone for calcite was fitted to a value of 40 cm-1 with a low weight in order to avoid a potential imaginary mode in this region of the phonon dispersion. For large scale MD simulations polarizable force fields are extremely expensive and hinder the possibility of reaching adequate time- and length-scales to study problems such as the nucleation and growth of CaCO3 nanoparticles in an aqueous environment. With this in mind, it was decided to develop a simpler force field for calcium carbonate systems with efficient molecular dynamics as the goal. Consequently, it was necessary to make two simplifications with respect to some of the more sophisticated models previously discussed. First, all atoms are treated as rigid ions (i.e., no polarization model is included) since this accelerates simulations by up to an order of magnitude. Many previous studies have also worked within this rigid ion description and so this is a standard approximation.41 Second, instead of describing the intramolecular degrees of freedom of the carbonate ion, it is treated as a rigid molecule. When combined with a similar description of water, of which there are many to choose from, this allows a longer time step to be employed provided the integrator used is able to efficiently treat rigid bodies. Although carbonate flexibility is important for achieving an accurate description of the mechanical and vibrational properties of the solid-state phases of calcium carbonate, it is not critical for the successful reproduction of the thermodynamics. For the proposed rigid molecule force-field for calcium carbonate, the fitted parameters for the solid-state are reduced to the atomic charges and the coefficients of the interatomic two-body potentials, which are chosen to be of the Buckingham form:

( Fr ) - rC

U(r) ) A exp -

6

(2)

Here, the C coefficients, which represent the dispersive attraction between species are set to zero for all interactions involving the relatively unpolarizable calcium cation, and thus only fitted for the intermolecular oxygen-oxygen potential. For the interaction between the calcium ions and water, the LennardJones (12,6) potential was employed since this functional form is more robust for solution phase simulations:

[( σr ) - ( σr ) ]

U(r) ) 4ε

12

6

(3)

TABLE 2: Force Field Parameters for the Aqueous Calcium Carbonate Systema Buckingham Ca Ca O O O Ca

O C O OW HW OW

A [eV]

F [Å]

C [eV Å6]

3161.6335 120000000. 63840.199 12534.455133 396.320957 1186.492877

0.271511 0.120000 0.198913 0.215172 0.230006 0.297

0.0000000 0.0000000 27.899008 12.090225 0.0000000 0.0000000

Lennard-Jones

ε [eV]

σ [Å]

Ca Ca OW

0.0005 0.0010 0.0070576

3.33 3.25 3.16435

OW OW OW

Set 0

Set 1 Set 2

a

Here O represents the oxygen of a carbonate group, while OW is the oxygen of a water molecule. The two-body potentials have been multiplied by a tapering function (eq 4) to make them go exactly to zero at rc ) 10 Å rt was always set equal to 6 Å. Charges (in a.u.) for Ca, C, and O are +2.0, +1.123282 and -1.041094, respectively.

A tapering function was applied to all two-body interactions to ensure a smooth and continuous energy and force of interaction:

x)

r - rt ; rc - rt

{

f(x) )

1 (1 + 3x + 6x )(1 - x) 0 2

3

if r < rt if rt g r e rc if r > rc

(4)

Here the parameters rt and rc were set to be 6 and 10 Å, respectively. The final parameters for the present force field are reported in Table 2. Unless otherwise stated, the molecular dynamics simulations reported below have been performed with the DL_POLY 2.19 package42 which we have modified to include the above tapering function. The equations of motion were integrated by using the velocity Verlet algorithm with a 2 fs time step. The temperature was always controlled by the Nose´-Hoover thermostat with a 0.1 ps relaxation time and for the constant pressure simulations the time reversible barostat of Martyna et al.43 was used with a relaxation time of 2 ps. The long-range electrostatics have been calculated with the Smooth Particle Mesh Ewald algorithm as implemented in DL_POLY.42 Here the Ewald precision was specified to be 10-8, leading to a convergence parameter of 0.38377 Å-1. For simulations of ions in liquid water using a cubic periodic cell this corresponds to a reciprocal space grid dimension of 16, while for the studies of the calcite-water interface the corresponding values are 32 and 64 in the plane of the surface and normal to this, respectively. For calculations on charged systems, a neutralizing background correction was applied.44 We have also performed a series of Monte Carlo calculations to locate the calcite-aragonite phase boundary. This provides a valuable diagnostic for the extent to which our new potential has correctly captured the relative energetics of the two phases. We once again stress that we are interested in the phase behavior of our model when treated classically, as in a large-scale molecular dynamics simulation. A standard procedure to locate solid-solid phase equilibria is to compute the free energy change along a pathway connecting each phase to the Einstein crystal,45 for which the absolute



free energy can be evaluated analytically. Having performed this procedure for a reference temperature and pressure in each crystal phase, the free energy change along an isotherm or isochore can be computed via thermodynamic integration. Intersections in the two free energy surfaces can be used to locate phase equilibria. While perfectly reasonable, this approach suffers from accumulation of error when performing the integration. Any inaccuracy in the initial reference free energy is therefore amplified in the location of phase boundaries. A superior approach is the lattice-switching Monte Carlo methodology of Bruce et al.46,47 Our implementation of this method for mineral systems is described briefly below, and discussed in more detail in a forthcoming paper. The lattice-switching approach allows for accurate computation of the free-energy difference between two phases at a particular temperature and density or pressure by sampling from both pure phases in a single simulation without the need to traverse the complex mixed-phase region. Simulations follow a standard Metropolis Monte Carlo procedure, with the positions of each atom, and the components of the three cell-vectors stored relative to one of two reference configurations, A and B, which belong to calcite and aragonite, respectively. In our simulations ion displacement and rotation moves are attempted with equal probability. Cell vector displacement moves are attempted with a probability 2Nu times smaller than ion moves, where Nu is the number of CaCO3 units. The lattice-switching methodology introduces a fourth Monte Carlo move (the lattice-switch) in which the reference configuration is swapped from A to B or B to A. This is attempted with the same average frequency as cell-vector displacement moves. In this case the lattice-switch incorporates a volume-change which must be accounted for in the acceptance probability. In principle, the above recipe will lead to configuration samples with the correct population within each of the two phases. The probability of inhabiting each phase at a particular temperature and pressure can then be computed and Boltzmanninverted to obtain accurate free-energy differences. In practice, any simulation of a feasible length will fail to sample both phases adequately due to the vanishing probability of accepting the lattice-switch. To overcome this problem, the latticeswitching method requires an “overlap parameter”, which in our case we define as

µ ) (HA - HB)/kBT

(5)

where HA and HB are the enthalpies of the system when its coordinates are measured from reference configurations A and B, respectively. Configurations with negative µ favor phase A (calcite) and those with positive µ favor B (aragonite). We then seek a suitable weight function η(µ) for a multicanonical simulation48 which will lead to uniform sampling over the entire range of µ including those within kBT of µ ) 0 where the lattice switch is likely to be accepted. The histogram P(µ) from the multicanonical simulation can then be unbiased to give the desired population statistics. Furthermore, the samples can be reweighted49 to other temperature and pressures near that simulated to give population statistics over a range of thermodynamic conditions. All simulations reported here divide the appropriate range of µ into 800 bins between µ ) -2000 and +2000. An initial guess of η(µ) ) 0 is used and iteratively refined using a Wang-Landau inspired50 recursion procedure. We judge our weight function to be converged when the change at each

TABLE 3: Lattice Parameters of Anhydrous and Hydrated Phases of CaCO3 at 300 K and 1 atm as Calculated with NσT MD Simulationsa mineral calcite

aragonite

vaterite

monohydrocalcite

ikaite

3

volume (Å ) a (Å) b (Å) c (Å) γ (deg) volume (Å3) a (Å) b (Å) c (Å) volume (Å3) a (Å) b (Å) c (Å) γ (deg) volume (Å3) a (Å) b (Å) c (Å) γ (deg) volume (Å3) a (Å) b (Å) c (Å) R (deg) β (deg) γ (deg)

force field

expt.

363.301 ( 0.586 4.933 ( 0.007 4.933 ( 0.007 17.242 ( 0.036 119.999 ( 0.107 225.400 ( 0.417 5.003 ( 0.005 8.047 ( 0.014 5.659 ( 0.012 1126.287 ( 1.924 7.245 ( 0.012 7.245 ( 0.012 24.772 ( 0.028 119.999 ( 0.132 252.442 ( 0.584 6.263 ( 0.012 6.263 ( 0.012 7.430 ( 0.014 120.003 ( 0.136 720.998 ( 2.082 8.380 ( 0.027 8.359 ( 0.018 11.044 ( 0.030 90.108 ( 0.539 111.230 ( 0.337 90.083 ( 0.438

367.61 4.988 4.988 17.061 120.00 226.91 4.961 7.967 5.741 1126.5 7.16 7.16 25.47 120.02 242.57 6.09 6.09 7.54 120.00 754.07 8.792 8.310 11.021 90.00 110.53 90.00

a The experimental data are taken from the American Mineralogist database. For calcite the experimental structure at 297 K is quoted,74 while for aragonite the structure is that of de Villiers.75 The vaterite structure used in our simulations is the one reported by Wang and Becker51 which turned out to have a lower internal energy than the disordered76 and the other ordered ones77 and the vaterite structural parameters are those from Kamhi76 with a supercell of a ) b) 3×a’ and c ) 3 × c’. Standard deviations are also quoted for the simulation results.

iteration is less than 10-10 kBT. Typically this requires between 4 and 10 million Monte Carlo cycles of Nu trial moves per cycle. Validation of Force-Field for Mineral Phases. Structural Properties. In Table 3 we report the structural parameters of the known anhydrous polymorphs of CaCO3 averaged over 200 ps of molecular dynamics in the isothermal-isobaric ensemble with a fully flexible cell (NσT). The agreement with the experimental cell parameters is acceptable and the error is always less than 2.5%. While this is inferior to the performance of the model of Rohl et al., which achieves a maximum error of 0.7% for all cell parameters of calcite and aragonite under static conditions, it is a necessary compromise in order to achieve the improved description of the thermodynamics. The agreement in lattice parameters is best for calcite, with a deviation of 1%, since the rigid geometry of the carbonate ion was chosen to have the experimental C-O bond length (1.284 Å) found in this polymorph. There is some variation of C-O bond length between phases, with it exhibiting values of 1.277-1.285 Å in aragonite and 1.287-1.298 Å in vaterite.51 There is also some difference between the cell parameters from the fitting procedure, which was performed statically, and those obtained in the molecular dynamics. However, these differences are small compared with deviations from the experimental values. Indeed, significant thermal expansion is only observed perpendicular to the plane of the carbonate groups for both calcite and aragonite. Phase Diagram. We have located the calcite-aragonite phase equilibrium in our new model for temperatures between 260

6002


Figure 1. Multicanonical weight function η(µ) at 300 K and 0.51 GPa of pressure (dashed line). The unbiased histogram P(µ) (solid line) is obtain from the resulting multicanonical simulation. The reweighted histogram at the estimated coexistence pressure of 0.54 GPa (dotted lines) is shown for comparison.

Figure 2. Calcite-aragonite phase boundary located using latticeswitching Monte Carlo calculations. Each cross represents the temperature and pressure at which a simulation was performed. The correspondingly colored dashed line identifies phase equilibria located by reweighting of samples. Three calculations of the phase boundary based on experimental data are shown for comparison.30,31,84 In the case of Carlson the curve has been extrapolated from higher temperature (350 °C) data.

and 360 K, this covering the range relevant to biomineralisation. Our simulations use unit cells containing a number of formula units, Nu, equal to 72. Reference configurations are taken as the optimized zero temperature configuration of each phase within the average cell vectors determined at 300 K and 0.51 GPa. At this temperature and pressure our lattice switching calculation produced the weight function and histogram shown in Figure 1. This clearly shows that calcite is the thermodynamically stable phase under these conditions in our model. Using this set of multicanonical weights we have performed a simulation of 4 million Monte Carlo cycles. Reweighting of these samples to alternative temperatures and pressures locate the 300 K coexistence pressure at 0.54 GPa. Although exact in principle, the reweighting procedure assumes that we have sampled all configurations that are statistically important at 0.54 GPa during our original simulation at 0.51 GPa, which is not guaranteed. This value is therefore approximate. To estimate the range of uncertainty on this value we have performed additional lattice-switching calculations at other temperature and pressures identified by linearly extrapolating the estimated phase boundary at 300 K to higher and lower temperatures using the gradient computed from the Clausius-Clapeyron equation. The locations of the four lattice-switching calculations in the pressure-temperature plane, along with the portion of the phase boundary estimated in each case, are shown in Figure 2 along with various experimental estimates. The dislocation between the segments of the computed phase boundary indicate the error

Raiteri et al. in the histogram reweighting procedure. Although the predicted phase boundary is at higher pressure than the experimental values, the difference is of a similar magnitude to the experimental uncertainty. For static calculations, as used to derive the force-field, the pressure of transition is approximately 0.4 GPa, and so the higher value calculated using Monte Carlo lattice switching is as the result of thermal effects. It should be noted that the difference in transition pressure of this magnitude corresponds to an error of approximately 0.3 kJ/mol in the enthalpy difference and so this represents a demanding test. It should also be noted that the agreement in the slope of the phase boundary is particularly satisfying, given that this is not guaranteed to be the case. Indeed we developed the latticeswitching methodology for mineral systems using a rigidcarbonate version of the force-field of Freeman et al.27 for which the phase boundary was located at a 45° angle to the experimental data. DeriWation of Force-Field for Mineral-Water Interactions. The second important set of interactions that need to be fitted are those between the mineral and the water molecules. From the perspective of studying nucleation and growth for CaCO3 the most important properties that need to be accurately reproduced are the hydration free energy of the ions and the dissolution enthalpy of the mineral. Hence, these are the objective quantities used to refine the present force-field. In the present work the hydration free energy has been calculated by means of the free energy perturbation (FEP) technique (see, e.g., ref 52). Within this approach the hydration free energy of an ion can be calculated by making the ion appear (disappear) in the solvent, which is readily achieved by gradually switching on (off) the solute-solvent interaction terms and the ion charges

1 ∆G ) - ln〈exp[-β∆U]〉0 β

(6)

where ∆U is the difference in the potential energy between the target and the reference states and 〈...〉0 denotes an ensemble average over configurations sampled from the reference state. In some cases eq 6 cannot be successfully used because the ensemble average may suffer from convergence problems, especially when large perturbations are involved. This difficulty in applying FEP theory can be circumvented through a simple stratification strategy, also often called staging. It relies on constructing several intermediate states between the reference and the target state, such that the direct evaluation of the free energy difference between successive states, ∆Gi, i+1, does not suffer from the aforementioned convergence problems. Intermediate states do not have to be physically meaningful, i.e., they do not have to correspond to systems that actually exist. More generally, we can consider the Hamiltonian as a function of some parameter, λ. Without loss of generality, we can choose 0 e λ e 1, such that λ ) 0 and λ ) 1 for the reference and target states, respectively. If we create N - 2 intermediate states linking the reference and the target states such that λ1 ) 0 and λN ) 1 then eq 6 can be used sequentially to yield

∆G ) -

1 β

N-1

∑ i)1

ln〈exp[-β(U(λi+1) - U(λi))]〉λi

(7)

If we then swap the reference and target states, eq 6 can be rewritten as


∆G )

1 ln〈exp[β∆U]〉1 β


(8)

where the ensemble average is now taken over configurations sampled from the target state. This ultimately leads to the GibbsBogoliubov inequalities

∆G e 〈∆U〉0

(9)

∆G g 〈∆U〉1

(10)

which set the upper and lower bounds for the free energy as 〈∆U〉0 and 〈∆U〉1, respectively. Therefore performing two sets of calculations in which the target and reference states are swapped is extremely useful in providing an estimate of the uncertainty of the calculation. The FEP calculations have been performed using a locally modified version of DL_POLY 2.1942 adapted for this purpose. The calculations have been performed mainly in the direction of extracting the solute from the solvent. However, the reverse process has also been performed in a number of cases to estimate the uncertainty of the calculation. This has been found to be approximately 2.9 kJ mol-1, which is comparable to ambient thermal energy. In our calculations the electrostatic and dispersion interactions have been treated separately, with the Coulombic term perturbed to zero first (or switched on last), with ten linearly spaced values of λ ) 1.0, 0.9, ..., 0.1 for each type of interaction. For the reverse process, a much finer logarithmic spacing for λ has been used when the Lennard-Jones interaction was switched on. Once the cavity around the solute was properly generated (λLJ ) 0.1) we used the same δλ ) 0.1 spacing as before. For every value of λ, NVT simulation with 20 ps of equilibration and 1 ns of averaging was carried out. The system size was one ion immersed a 521 water molecule box (25.1 Å in side). A finite size correction for calculation of the hydration free energy of the ions was applied, as described in the work of Hummer et al.53 To check if there were any residual finite size effects, the hydration free energy of Ca2+ was also calculated for a 2449 water box (42 Å in side) and no difference was found. Aside from the size-dependent effects, there has been considerable debate regarding the asymptotic limit for the free energy of solvation of ions, relating to the difference between M-summation versus P-summation,54 also known as the C1-type correction.55,56 Here the difference relates to whether the realspace summation is truncated based on the position of a centroid of a molecule (M) or the position of each individual particle (P). Although some controversy remains, in the present work we adopt the P-summation limit as advocated by Hummer et al.54 For comparison with previous models, we also calculated the hydration free energy for the models of Freeman et al.27 and Bruneval et al.23 We should note that for the Freeman et al. potential we were forced to keep the carbonate molecule rigid since the calculations in the original article were performed with a locally modified version of DL_POLY in which the intramolecular electrostatic and van der Waals interactions were not excluded. We also note that for the model of Bruneval et al. the shell position was relaxed at every time step, which was reduced to 1 fs. The estimation of the computational error of the hydration free energies was done by splitting the trajectory into 5 intervals and by taking the average and standard deviation of the 5 measures.

TABLE 4: Free Energy of Hydration for Ca2+a ∆Ghyd [kJ mol-1] -1594. -1514. -1505. -1491. -1444. -1350. -1333. -1283. -1503. -1447.

( ( ( ( (

7. 1. 2. 2. 2.

year

ref

1965 1977 1991 2000 2001 this work this work this work this work this work

78 79 80 81 70 SWM4-NDP23 Freeman et al.27 Set 0 Set 1 Set 2

a The rows with a year specified are experimental values, while those labelled as being from this work are values computed by free energy perturbation as described in the text. For the simulation values, standard deviations are given to indicate the degree of statistical uncertainty.

As a starting point for fitting of the ion-water interactions we used the potentials that have been developed by de Leeuw and Parker20 and successfully used by Bruneval et al.23 This model proved accurate in reproducing the structural properties of the ions in solution and interesting results regarding the stability of the ion pair were previously obtained. Several water models have been used in prior works in combination with these ion-water potentials. Originally, de Leeuw and Parker employed a shell model description combined with a flexible water molecule, though this has subsequently been shown to have serious limitations. Bruneval et al. replaced this water force field with an alternative shell model description, based on rigid molecules, that has been more extensively tested and shown to give good results for the bulk liquid. In the present work, we have chosen to use an alternative water model in keeping with the objective of computational efficiency. Water-water interactions are therefore described using the four-site rigid model TIP4P-Ew.57 To evaluate the performance of the ion-water interaction parameters we have calculated the hydration free energy of the Ca2+ ion with the original polarizable potential of Bruneval et al. (SWM4-NDP), which turned out to significantly underestimate the hydration free energy (Table 4). As a first approximation for our rigid model we decided to use the same watermineral interaction parameters as for the polarizable water case; this set of parameters will be hereafter called Set 0. It is evident from Table 4 that Set 0 is significantly worse than the polarizable SWM4-NDP model. However, this is not surprising since the water-mineral interactions were fitted for a completely different water model. Despite the poor performance of this model, results are included since they illustrate the failings that will result from underestimating the hydration free energy. From Table 4 it is apparent that there is actually a wide variation in the experimental data for the hydration energy of a calcium cation, with three values around -1500 kJ mol-1 and one just above -1450 kJ mol-1. Given this uncertainty, we have produced two new sets of parameters (Table 2) that reproduce these two values (Table 4), hereafter called Set 1 and Set 2. Following the same procedure as described above, we calculated the hydration free energy of CO32- (Table 5). In this case the original parameters by de Leeuw and Parker20 gave a satisfactory agreement with the experimental value, and thus no further refinement was required. In the present work we shall calculate several properties of Ca2+, CO32-, and calcite to show what dramatic effects an under/ overestimation of the mineral-water interaction may potentially have on the CaCO3 nucleation and growth in aqueous solution.

6004


Raiteri et al.

TABLE 5: CO32- Hydration Free Energy As Determined Experimentally (First Line) and from Free Energy Perturbation Calculations Performed in the Present Worka ∆Ghyd [kJ mol-1] -1314. -1404. ( 7. -1175. ( 1. -1301. ( 2. a

exp.80 SWM4-NDP23 Freeman et al.27 this work

See text for details.

TABLE 6: Calcite Dissolution Enthalpy for the Three Sets of Force-Field Parameters Developed in the Present Work As Compared against Both the Experimental Value and Those Obtained Based on Two Literature Force-Fields ∆Hdiss [kJ mol-1] -12.54 102.28 139.91 164.03 -68.51 -12.54

ref 82

exp. SWM4-NDP23 Freeman et al.27 Set 0 Set 1 Set 2

Our results will also demonstrate that Set 2 is the one that best reproduces the energetics of the ions in water, and thereby leads to the most accurate description of the thermodynamics for the calcium carbonate-water system. Results Having derived a new force-field model for the aqueouscarbonate system that targets the correct thermodynamic properties, the application of this parametrization to a range of properties that were not explicitly fitted can now be performed, as described in the following sections. Calcite Dissolution Enthalpy. By combining the results of the simulations used during the parametrization, a quantity that can be readily extracted is the dissolution enthalpy of calcite CO3 calcite ∆Hdiss ) ∆HCa hyd + ∆Hhyd - ∆Hlat

(11)

Ca CO3 and ∆Hhyd are the hydration enthalpies of Ca2+ where ∆Hhyd 2and CO3 , respectively, and ∆Hcalcite is the lattice enthalpy of lat calcite. Calculated results for the enthalpy of dissolution are given in Table 6, as is the experimental value. Both of the literature force-fields shown, along with Set 0, which underestimates the solvation energy for calcium, predict that dissolution of calcite would exhibit a large positive enthalpy change, in contrast to experiment. On the other hand, Set 1 with its stronger affinity of calcium cations for water yields an enthalpy change that is too favorable. Somewhat fortuitously, the dissolution enthalpy for Set 2 exactly matches the experimental value to within the statistical fluctuations present. The significance of this good agreement, despite not having been explicitly fitted, is that it suggests that the choice of the most recent experimental value for the free energy of solvation of Ca2+ leads to greatest degree of internal consistency between the data for the individual ions and the overall dissolution enthalpy. Based on the fact that Set 2 reproduces both the solvation free energies of the individual ions and the enthalpy of dissolution for calcite, this is our recommend final parameter set. However, comparisons will also be presented for the other sets, in the case of some properties, since it will illustrate the consequences of either over- or under-estimating the affinity

Figure 3. Average Ca2+-OW distance (a) and total energy (b) for different numbers of water ligands around the Ca2+ ion. The DFT values are taken from ref 58, whereas the curves labeled SWM4-NDP and Freeman et al. are obtained using the parameters reported in refs 23 and 27, respectively.

of calcium cations for water on the description of the aqueous calcium carbonate system. Calcium-Water Binding in Vacuo. As a further examination of the solute-solvent interaction parameters, the binding of water to a Ca2+ ion in vacuo has been calculated as a function of cluster size. In Figure 3 we plot both the average Ca-O distance and the total binding energy as a function of the number of water molecules for the three sets of parameters derived in the present study. Included in Figure 3 as a reference are the values previous calculated using density functional theory at the BLYP/6-311+G** level by Bako et al.58 For further comparison, we have also included the data for the two most recent calcite-water models, namely those of Freeman et al.27 and Bruneval et al.23 All cluster geometries were optimized within DL_POLY 2.19 until the forces on the atoms were less than 1 × 10-6 kJ mol-1/Å. The force-field combination of Bruneval et al. produces Ca2+-OW distances in very good agreement with those obtained from DFT (Figure 3a). In contrast, all the rigid ion force-fields consistently predict shorter Ca2+-OW distances than DFT, with the water molecules becoming closer to the Ca2+ ion as the hydration free energy becomes more negative. This is to be expected due to the fact that the main contribution to the hydration free energy comes from the electrostatics; hence the closer the water molecules can approach, the larger the electrostatic interaction becomes. This effect is also evident in the cluster internal energy (3b), which shows the same trend as the hydration free energy. Although Set 2 consistently predicts shorter Ca-O distances than the published DFT values, it must be recognized that there is potential for error from both sides. Given that the DFT values



Figure 4. Ca2+-OW pair distribution function for the calcium cation in aqueous solution.

TABLE 7: Ca2+-OW Distance, rCa-O, Coordination Number, n, and Water Residence Time, τ, in the First Ca2+ Hydration Shella exp. avg theo. avg Set 0 Set 1 Set 2

rCa-O [Å]

n

2.43 2.47 2.42 (2.49 ( 0.057) 2.17 (2.21 ( 0.036) 2.31 (2.36 ( 0.045)

7.28 8.18 7.48 ( 0.53 6.07 ( 0.26 7.14 ( 0.44

τ [ps] ∼100 107 ( 2 650 ( 30 100 ( 11

a For the Ca2+-OW distance, two values are quoted; the first is the position of the first peak in the radial distribution function, whereas the second (in parentheses) is the average within the cut-off radius of 3.3 Å. The averages of experimental and theoretical data for the Ca2+-OW distance and the coordination number are taken from Table 1 in ref 59, whereas the water residence time is from ref 83.

were obtained with a GGA functional, it is likely that there is a systematic overestimation in the bond lengths caused by underbinding. Indeed, in the work of Bruneval et al. they contrast their force-field with a calculation at the MP2/6-311+G* level for a single water molecule and show that the Ca2+-OW distance is overestimated and that the binding energy is significantly underestimated, as would be anticipated. In 3b it can be readily seen that only Sets 1 and 2 give more exothermic binding energies than DFT in the limit of the coordination environment of Ca2+ cations found in bulk water, while the remaining models are the wrong side of this upper bound. Coordination Shell. Although cluster calculations represent a straightforward means of examining the quality of a forcefield, it is more important that the description of ions in bulk water is correct, rather than that of the gas phase. Consequently, extensive molecular dynamics simulations have been performed for the ions within a periodic box of explicit water molecules. Both experiment and computer simulations have reported a wide range of values for the coordination number (5.5-10), whereas there is marginally better agreement on the Ca2+-OW average distance (2.39-2.64 Å) (see ref 59 and references therein). The pair distribution function (Figure 4), the average coordination number (Table 7), and the water residence time in the first Ca2+ hydration shell (Table 7) have been calculated for our three sets of parameters during 100 ns of NPT simulation at 300 K and 1 atm. To determine whether a water molecule is within the first solvation shell of Ca2+ a cutoff distance of 3.3 Å was employed. This value lies within the broad minimum in the radial distribution function for Ca2+-OW where this function is close to zero. Varying the cutoff between 3.0 and 3.5 Å leads to a change in the coordination number of only 0.08. The water residence time, τ, has also been estimated by calculating the

Figure 5. Dissociation free energy profile for an ion pair in water. The dissociation free energies are 69.5 ( 2, 30.9 ( 1 and 39.6 ( 1 kJ mol-1 for Sets 0, 1, and 2, respectively.

distribution of the water-exchange time in the Ca2+ first solvation shell and then fitting it with an exponentially decaying function, C exp(-t/τ). As found for the cluster calculations, the first Ca2+ hydration shell for Set 1 is extremely compact and close to the ion, whereas Set 0 has a more open solvation shell, which implies that extra water molecules can be temporarily incorporated, increasing the average Ca2+-OW coordination number from the more intuitive 6 to almost 8. A close analysis of the MD trajectory for Set 1 revealed that the Ca2+ hydration shell is so compact and strongly bonded that it behaves almost like a supramolecular complex with the water molecules’ residence time being more than half a nanosecond (Table 7). In contrast, Set 2 shows good agreement with both the experimental averaged coordination number, and the residence time for water within the first solvation shell. Although the Ca2+-OW distance is underestimated by this parameter set, this is less significant than obtaining the correct thermodynamics and dynamics for the water environment. Furthermore, rigid water models, such as TIP4P-Ew, use an idealized geometry to fix the interaction centers, which inevitably leads to an offset in local geometric measurements. The same calculation has also been performed for CO32- in water, and we observe the presence of a hydrophobic region on top of the carbon atom, in agreement with previous first principles and classical MD simulations.23 Because the same carbonate-water interaction parameters are taken in the present study as used in previous works, it is unsurprising to find a very similar description of the solvation of CO32- with the change of water model only making a minor contribution. Ion Pairing in Water. The first step toward the understanding of CaCO3 nucleation and growth from aqueous solution is the aggregation of ions in water. We have therefore calculated the dissociation free energy of the charge neutral ion pair by performing a set of Umbrella Sampling simulations60 using the PLUMED utility.61 The free energy was then extracted via the weighted-histogram method, utilizing the program WHAM.62 As a reaction coordinate we used the Ca-C distance, with 25 windows equally spaced every 0.4 Å starting from 1 Å and with a spring constant of 2 eV/Å2. The three sets of parametrizations give similar free energy profiles, as shown in Figure 5. Here the deep minimum for the tight ion pair is split into two, with the two configurations corresponding to a bi- and monodentate binding of the CO32to the Ca2+ in order of increasing Ca-C distance. This feature is not present in the simulations with the polarizable force field of Bruneval et al.,23 though they do observe a point of inflection

6006


in this region of the free energy profile. In terms of the free energy of dissociation, Sets 1 and 2, with values of 31 and 40 kJ mol-1, respectively, are consistent with the previous calculations based on the polarizable force field, 35.7 kJ mol-1,23 whereas Set 0 predicts a value that is almost twice as large at 67.5 kJ mol-1. Despite the quantitative variations, all the models indicate that the tight ion pair is the most stable state and that if this associated configuration is formed in solution it is unlikely to dissociate. Of course, we note that the carbonate ion is in equilibrium with the bicarbonate species in solution, as a function of pH, and so dissociation can be assisted by proton transfer. Aqueous Interface with the (101j4) Calcite Surface. The morphology of calcite crystals is typically dominated by the (101j4) termination, and therefore this is most widely studied interface between calcium carbonate and water. Although the growth of calcite usually occurs at surface defects, such as screw dislocations, a proper description of this surface is still an important prerequiste for a model to be applicable to studies of the aqueous calcium carbonate system. Under vacuum conditions the (101j4) surface has been observed to undergo a (2 × 1) reconstruction; a phenomenon supported by the most accurate of the polarizable force-fields. However, we find that even with such force-field models the reconstruction is averaged out by heating to ambient conditions within molecular dynamics for the dry surface, and in the presence of water it is not observed at all (of course it is impossible to examine the low temperature behavior in the presence of water since the solvent will freeze). Although the present model, being rigid ion, fails to give a (2 × 1) reconstruction, this is of negligible consequence for the solvated (101j4) surface. We have studied the aqueous interface for the (101j4) calcite surface by performing several NVT simulations of a calcite slab in water at 300 K. Before starting these simulations the aqueous interface was equilibrated to 1 atm with the free energy minimization procedure described in ref 63. The calcite slab was composed of 864 CaCO3 formula units and approximately 3 nm thick. 3D periodic boundary conditions were applied with the calcite slab being separated from its images by 4.7 nm of water (2974 molecules) along the z direction. Validation studies were performed to ensure that the calcite slab and water layer are both sufficiently large that the system behavior is not influenced by these parameters. Simulations were performed for all the three sets of parameters described here. However, only minor differences were observed and only the results for Set 2, as the most accurate parametrization, are presented here. The similarity between different parametrizations of the Ca2+ cation-water interaction strength is not surprising because of the strong influence of the water hydrogen-bonding network on the structure of the interface, and the fact that the main surface interaction is with the carbonate groups, which remains unchanged between the three sets of parameters. Experimental data exists4 for the water density distribution as a function of the distance from the surface (Figure 6) against which comparison can be made in terms of their fitted model. Our results are found to be in broad agreement with the experimental model4 and with previous simulations of the (101j4) calcite surface,21,64,65 as well as data for similar systems.66 There is a slight difference from one previous simulation study by Perry et al.67 in that the innermost water peak exhibits a small splitting in their study, which is absent in the present work. According to our model the average distance of the water molecules of first water layer from the surface (2 Å) is approximately 0.3 Å shorter than the experimental model.4

Raiteri et al.

Figure 6. Density profiles for species normal to the aqueous interface with the (101j4) calcite surface. The distance is defined relative to the average z coordinate of the Ca2+ ions in the outermost layer, where the origin is indicated with a dashed line.

Figure 7. Water structure at the (101j4) surface. The top figure illustrates the complete aqueous interface with oxygen, carbon, calcium, and hydrogen colored red, green, blue, and white, respectively. The first solvation layer is highlighted in orange and shown in the middle figure, whereas the second layer is similarly highlighted in yellow in the bottom figure.

However, Geissbuhler et al.4 also noted that their model implies a Ca2+-OW distance of about 3 Å which is close to 0.5 Å longer than the value found for the ions in solution, whereas in our simulations, the distance is almost unchanged. A better insight into the water layering structure can be obtained by examining the atomic trajectories (see for example a snapshot in Figure 7) and from the water density map in planes parallel to the (101j4) surface (Figure 8). From the middle panel of Figure 7 it is evident how the water molecules belonging to



Figure 9. Dissolution free energy of a Ca2+ ion (a) and a CO32- ion (b) from within the (101j4) calcite surface. For ease of comparison the curves have been aligned to have a zero free energy for the nondefective surface.

Figure 8. Water density map in planes parallel to the (101j4) calcite surface at different distances along the surface normal. Panels a-c correspond to the first, second and third peaks, respectively, in the water density profile shown as the curve for OW in 6. Density is shaded from blue (low density) to red (high density). For the sake of clarity we show the underlying calcite surface only in half of the image.

the layer closest to the (101j4) surface have an OW-HW bond oriented parallel to the mineral surface, forming a weak hydrogen bond with the outermost CO32- molecules. On the other hand, the water molecules in the second water layer (lower panel of Figure 7) form strong directional hydrogen bonds with the CO32- of the calcite surface, consistent with previous MD simulations.21 The extreme ordering of the two water layers close to the mineral surface can be shown by looking at the averaged 2D water distribution in planes parallel to the (101j4) surface (Figure 8).

According to our simulations, the position of the water molecules in the second layer is strictly defined by the strong hydrogen bonds with the mineral surface, which is to an extent consistent with the model in ref 4 where their fitting accuracy was insensitive to the lateral position of the water molecules in that layer. Our simulations seem to indicate the presence of a third weakly ordered layer of water, as shown by the peaks in the density profile normal to the surface in Figure 6, which was also suggested by previous calculations.21 From the water density map for this layer (Figure 8c), the molecules display a slight preference for hexagonal order that resembles that of the first layer to which it exhibits some tendency to hydrogen bond. One qualitative difference with respect to experiment is the presence of two nonequivalent water positions, both in the first and second water layer, whereas the structural model by Geissbuhler et al. allowed for only one in each layer. However, because of the difficulty of achieving lateral resolution in the experiments, this does not imply an inconsistency. Dissolution of Ions from the (101j4) Calcite Surface. While recognizing that dissolution of calcite will most probably occur from defects, such as kink sites, the energetic profile for removal of ions from the flat (101j4) surface is a quantity that has been previously computed to examine the tendency for dissolution from this facet and therefore a quantity that can be contrasted between models. Hence, the free energy to remove an ion from a flat (101j4) calcite surface by dissolving it into water (Figure 9 and Table 8) has been determined for the present force-field that yields the proper balance between the thermodynamics of the solid and solution. To calculate this quantity we performed a set of umbrella sampling simulations using the z coordinate of the center of mass of the ions as a reaction coordinate. The calcite (101j4) surface was set to z ) 0, and its position was defined as the average z coordinate of the outermost Ca2+ ions.

6008


TABLE 8: Dissolution Free Energy [kJ mol-1] of Ca2+ and CO32- Ions Belonging to the Outermost (101j4) Calcite Layera Set 0 Set 1 Set 2 ref 65

Ca2+

CO32-

262.4 ( 7.7 118.7 ( 1.9 211.3 ( 12.5 ∼270

315.5 ( 1.0 251.8 ( 25.1 248.9 ( 1.9 ∼40

a Values from Spagnoli et al.65 are approximate as they are estimated from figures rather than quoted in the text.

Figure 10. Dissolution free energy of a Ca2+ ion (a) and a CO32- ion (b) adsorbed on the (101j4) calcite surface. For ease of comparison the curves have been aligned to have a zero free energy at 12 Å from the surface. The error-bars on the curves are between 0.5 and 7 kJ mol-1, depending on the position and simulations. For Set 2, the error-bar is always smaller than 2.9 kJ mol-1.

To prevent any surface displacements during the simulation, the central layer of the calcite slab was restrained to a fixed position. The starting configuration was the final configuration that we obtained after studying the aqueous interface for the (101j4) calcite surface. For the umbrella sampling, a set of 30 windows evenly spaced at 0.5 Å and a spring constant 1 eV/Å2 have been employed. In every window we performed a 2 ns NVT simulation with a Nose´-Hoover thermostat and a relaxation time of 0.1 ps; the first 50 ps were considered as equilibration. Consistent with previous simulations,65 the extraction free energy profile is structure-less and converges to a flat plateau between 5 and 10 Å from the (101j4) surface. However, the dissolution free energies estimated by our model greatly differ from those reported by Spagnoli et al..65 First, we observe the opposite behavior when comparing the two ions, with Ca2+ always being more easily extracted than CO32- which is consistent with the relative solvation free energies. Second, the free energy of dissolution of both ions is found to be considerably endothermic, consistent with the conclusion that direct dissolution of the flat surface is highly improbable. This contrasts

Raiteri et al. with the results of Spagnoli et al.,65 where the dissolution free energy for the CO32- anion is relatively modest and there is no significant additional kinetic barrier given the smoothly increasing nature of the free energy profile. Finally, it is evident that the free energy of dissolution for Ca2+ is sensitive to the parameter set, while that for the CO32- ion is relatively insensitive. This suggests that the solvation state of the ion in solution is the dominant factor, rather than the solvation of surface defect that is created. Dissolution of Ions Adsorbed on the (101j4) Calcite Surface. From the perspective of understanding the growth of CaCO3 in water, it is important to characterize the interaction between ions and the flat (101j4) surface. This will determine whether there is any tendency for nucleation on this facet and whether the surface plays a role in guiding ion diffusion to other growth sites, such as steps and kinks. Hence, the adsorption/ dissolution free energy of additional ions on the flat (101j4) surface (Figure 10) has been determined. The same umbrella sampling procedure was used here as described in the previous section. To maintain the chargeneutrality of the system, Ca2+ and CO32- ions were added to opposite sides of the calcite slab, while the free energy profile was performed separately for each ion. By changing the relative position of the ions for selected umbrellas we verified that the thickness of both the water layer and of the calcite slab were sufficient to screen any interaction between the ad-ions. In some cases extra windows were necessary to better resolve the free energy barriers; these were placed at 0.25 Å from the existing ones, with a spring constant of 3 eV/Å2. First of all we note that Set 0 predicts a deep absorption minimum both for Ca2+ and CO32- at about 3 Å from the surface, which progressively disappears in Set 2 and Set 1 while the rest of the free energy profile remains almost unchanged. This position corresponds to the ions sitting between the two water layers (see Figure 6) with the ions directly bonded to the surface with a water molecule having been displaced. Hence, this state can be considered as an inner sphere complex. It is interesting to note that Set 0 predicts for CO32- that the inner sphere complex has two distinct minima that correspond to a monodentate and bidentate binding of the carbonate ion to the surface. This feature is still partially visible for Set 2, while it is completely absent for Set 1. The case of CO32- is particularly interesting since the most stable configuration goes from the bidentated to the monodentated and eventually to the outersphere complex as the hydration free energy of Ca2+ becomes more exothermic. This appears to be a property where there is a strong interdependency of the Ca2+ and CO32- solvation parameters in the force-field. With the exception of Set 0, which is the least accurate model included to demonstrate the effect of overestimating the hydration free energy for calcium, there is only a shallow minimum for the Ca2+ ion adsorbed at the surface and this configuration is unfavorable with respect to the cation being in bulk water. This observation is consistent with the behavior previously seen for the mineral barite,68 where the molecular anion must already be present on the surface for the barium cation to be exothermically adsorbed. In contrast to the present data, the results of Kerisit and Parker21 suggest that the Ca2+ ion is almost isoenergetic when comparing the inner sphere complex to bulk solution, with a significant barrier to diffusion away from the surface. This qualitative difference can be traced back to the fact that their model underestimates the free energy of solvation for Ca2+ with respect to experiment and so their results most closely resemble those of Set 0 in the present study. From now on we will focus

Simulating the Growth of Calcium Carbonate only on the most reliable force-field parametrization, Set 2, and its consequences. In addition to the free energy profile normal to the surface, the diffusion pathways for the ions parallel to the surface plane are also important to consider because of the consequences of mass transport during crystal growth. For example, it might be imagined that growth at kink sites is accelerated by diffusion of ions that are adsorbed on the surface until they reach a step, which in turn will guide them to the active site for incorporation to the crystal, thereby reducing the need for ions to collide directly with the kink from solution. For such a mechanism to occur requires the ions to adsorb on the surface, and to be able to diffuse across it. Conventional molecular dynamics is unable to reach the time scales required to observe significant diffusion of ions at the calcite (101j4) surface. Consequently, metadynamics simulations69 have been performed using the x and y coordinates of the ad-ions as order parameters. While insufficient sampling has been achieved to obtain an accurate free energy surface, insights are available into the diffusion pathways across the approximate landscape. The key observation is that no diffusion occurs across the surface via an inner sphere complex. In the case of Ca2+, the ion readily leaves the local minimum near the surface and diffuses into the bulk solution. For CO32-, the ion jumps from the inner sphere complex to the broader minimum associated with the solvent-separated outer sphere complex. Once at this distance from the surface, the ion begins to diffuse readily parallel to the plane of the surface with only a low activation energy. From this outer sphere position, the barrier to return to the inner sphere complex is similar to, though slightly greater than, that required to migrate into the bulk solution. Therefore, little preference is exhibited between the two directions. The consequence of this is that the flat surfaces will not contribute to the reactive surface area of calcite, as there is no adsorption in this region, with implications for kinetic growth models.6 Conclusions In the present work we have examined the thermodynamics of the calcium carbonate-water system based on a force-field description of the interaction energy. While there have been numerous previous attempts to produce an accurate set of interatomic potentials to model the mineral phases of calcium carbonate and their aqueous interfaces, it has been demonstrated that none of these are able to reproduce the free energies of the component ions in solution, as well as the energy difference for the calcite-aragonite phase transition. In light of this, three new force-field parametrizations have been examined. Two of these models, Set 1 and Set 2, reproduce the experimental thermodynamics, but differ in the choice of the literature value for the free energy of solvation of the Ca2+ ion. By examining the dissolution enthalpy for calcite, which is exactly reproduced only by Set 2, it is possible to discriminate between the models and conclude that this parameter set is the most accurate. This implies that the choice of the most recent hydration free energy from David et al.70 for the calcium cation leads to be best internal consistency between the data. Given that Set 2 is the most accurate and is based on an intermediate free energy of solvation, relative to other experiment and theoretical models, the other two parameter sets act as an indication of the consequences of using a force-field that is in error for the thermodynamics toward either extreme. In the present work we have use a force-field that specifically aims to be efficient within the context of molecular dynamics

J. Phys. Chem. C, Vol. 114, No. 13, 2010 6009 simulations through the use of a rigid description of the carbonate anion and water molecule. In other contexts, the use of rigid species may prove less convenient and so the development of a flexible form of the present model is currently underway. A further limitation common to the present forcefield, as well as those that have gone before, is that the speciation of carbonate cannot be directly described. Hence, other methods are required to resolve the question of when carbonatebicarbonate interconversion occurs as calcium carbonate begins to crystallize from aqueous solution. The combination of a flexible form of the present force-field with multi-state empirical valence bond theory, as has been employed to simulate proton transfer in water,71,72 provides a means of overcoming this limitation in the future and is also being explored. Aside from quantitative differences between the present optimized model and previous work, several qualitative discrepancies become apparent that are demonstrated to have implications for the atomic-scale picture of the calcite-water interface. Unlike previous works, the inner sphere complex of Ca2+ at the (101j4) surface of calcite is found to be irrelevant to the structure of the interface, with the ion showing no tendency to associate with the clean surface in the absence of a previously deposited carbonate anion. Carbonate is also found to only weakly associated with the flat calcite basal plane. The destabilization of the inner sphere complex has significant implications for the calcite surface growth. Indeed, the traditional picture of surface growth, based on parallels with gas phase behavior, where the ions are adsorbed on the surface, followed by diffusion until they reach a step and then ultimately a kink, are not applicable to the (101j4) calcite surface. Both ions prefer to diffuse via the solution, though carbonate migration may be slightly enhanced when separated by a single solvation layer from the underlying surface. Although speculative at present, it is highly likely that the large changes observed here in the thermodynamic predictions for calcium carbonate in water will have important consequences for attempts to simulate the structure, nucleation and growth of all phases of this material from aqueous solution, including nanoparticles.73 Therefore, the newly develop model represents an important step toward reliably and efficiently simulating processes such as surface growth and dissolution, as well as ultimately biomineralization. Acknowledgment. P.R. and J.D.G. thank the Australian Research Council for Fellowships and support through a Discovery grant, as well as NCI and iVEC for the provision of computing resources. D.Q. and P.M.R. acknowledge the support of the EPSRC through Grant Nos. GR/S80127 and GR/S80103. The use of facilities at the University of Warwick Centre for Scientific Computing is also gratefully acknowledged. We also thank Dr. Colin Freeman for useful discussions and provision of some parameters. References and Notes (1) Lackner, K. S.; Wendt, C. H.; Butt, D. P.; Joyce, E. L.; Sharp, D. H. Energy 1995, 20, 1153–1170. (2) Biomineralization. Dove, P. M., de Yoreo, J. J., Weiner, S., Eds.; ReViews in Mineralogy and Geochemistry; Mineralogical Society of America: Washington, D.C., 2003; Vol. 54. (3) Al-Anezi, K.; Hilal, N. Sep. Purif. ReV. 2006, 35, 223–247. (4) Geissbuhler, P.; Fenter, P.; DiMasi, E.; Srajer, G.; Sorensen, L. B.; Sturchio, N. C. Surf. Sci. 2004, 573, 191–203. (5) Stipp, S. L. S. Mol. Simul. 2007, 107, 342–381. (6) Morse, J. W.; Arvidson, R. S.; Luttge, A. Chem. ReV. 2007, 107, 342–381. (7) Elhadj, S.; de Yoreo, J. J.; Hoyer, J. R.; Dove, P. M. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 19237–19242.

6010


(8) Soler, J. M.; Artacho, E.; Gale, J. D.; Garcia, A.; Junquera, J.; Ordejon, P.; Sanchez-Portal, D. J. Phys. Cond. Matter 2002, 14, 2745– 2780. (9) Carlson, W. D. ReV. Mineral. 1983, 11, 191–225. (10) Yuen, P. S.; Lister, M. W.; Nyburg, S. C. J. Chem. Phys. 1978, 68, 1936–1941. (11) Dove, M. T.; Winkler, B.; Leslie, M.; Harris, M. J.; Salje, E. K. H. Am. Mineral. 1992, 77, 244–250. (12) Pavese, A.; Catti, M.; Price, G. D.; Jackson, R. A. Phys. Chem. Miner. 1992, 19, 80–87. (13) Jackson, R. A.; Price, G. D. Mol. Simul. 1992, 9, 175–177. (14) Pavese, A.; Catti, M.; Parker, S. C.; Wall, A. Phys. Chem. Miner. 1996, 23, 89–93. (15) Fisler, D. K.; Gale, J. D.; Cygan, R. T. Am. Mineral. 2000, 85, 217–224. (16) Rohl, A. L.; Wright, K.; Gale, J. D. Am. Mineral. 2003, 88, 921– 925. (17) Hwang, S.; Blanco, M.; Goddard III, W. A. J. Phys. Chem. B 2001, 105, 10746–10752. (18) Thackeray, D. J.; Siders, P. D. J. Chem. Soc., Faraday Trans. 1998, 94, 2653–2661. (19) de Leeuw, N. H.; Parker, S. C. Phys. ReV. B 1998, 58, 13901– 13908. (20) de Leeuw, N. H.; Parker, S. C. Phys. Chem. Chem. Phys. 2001, 3, 3217–3221. (21) Kerisit, S.; Parker, S. C. J. Am. Chem. Soc. 2004, 126, 10152– 10161. (22) van Maaren, P. J.; van der Spoel, D. J. Phys. Chem. B 2001, 105, 2618–2626. (23) Bruneval, F.; Donadio, D.; Parrinello, M. J. Phys. Chem. B 2007, 111, 12219–12227. (24) Tribello, G. A.; Liew, C.; Parrinello, M. J. Phys. Chem. B 2009, 113, 7081–7085. (25) Archer, T. D.; Birse, S. E. A.; Dove, M. T.; Redfern, S. A. T.; Gale, J. D.; Cygan, R. T. Phys. Chem. Miner. 2003, 30, 416–424. (26) Lamoureux, G.; Harder, E.; Vorobyov, I. V.; Roux, B.; MacKerell, A. D., Jr. Chem. Phys. Lett. 2006, 418, 245–249. (27) Freeman, C. L.; Harding, J. H.; Cooke, D. J.; Elliot, J. A.; Lardge, J. S.; Duffy, D. M. J. Phys. Chem. C 2007, 111, 11943–11951. (28) Wolf, G.; Lerchner, J.; Schmidt, H. G.; Gamsjager, H.; Konigsberger, E.; Schmidt, P. J. Therm. Anal. 1996, 46, 353–359. (29) Wolf, G.; Konigsberger, E.; Schmidt, H. G.; Konigsberger, L. C.; Gamsjager, H. J. Thermal Anal. Calor. 2000, 60, 463–472. (30) MacDonald, G. J. F. Am. Mineral. 1956, 41, 744–756. (31) Carlson, W. D. Am. Mineral. 1980, 65, 1252–1262. (32) Redfern, S. A. T.; Salje, E. J. H.; Navrotsky, A. Contrib. Miner. Petrology 1989, 101, 479–484. (33) Gale, J. D. J. Phys. Chem. B 1998, 102, 5423–5431. (34) Perdew, J. P.; Burke, K.; Erznerhof, M. Phys. ReV. Lett. 1996, 77, 3865–3868. (35) Perdew, J. P. A. R.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constantin, L. A.; Zhou, X. L.; Burke, K. Phys. ReV. Lett. 2008, 100, 136406. (36) Spagnoli, D.; Refson, K.; Wright, K.; Gale, J. D. Phys. ReV. B. in press. (37) Wu, Z.; Cohen, R. E. Phys. ReV. B 2006, 73, 235116. (38) Armiento, R.; Mattsson, A. E. Phys. ReV. B 2005, 72, 085108. (39) Gale, J. D. Phil. Mag. B 1996, 73, 3–19. (40) Gale, J. D.; Rohl, A. L. Mol. Simul. 2003, 29, 291–341. (41) Quigley, D.; Rodger, P. M.; Freeman, C. L.; Harding, J. H.; Duffy, D. M. J. Chem. Phys. 2009, 131, 094703. (42) Smith, W. Mol. Simul. 2006, 32, 933–1121. (43) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117–1157. (44) Leslie, M.; Gillan, M. J. Phys. C, Solid State Phys. 1985, 18, 973– 982. (45) Frenkel, D.; Ladd, A. J. C. J. Chem. Phys. 1984, 81, 3188–3193.

Raiteri et al. (46) Bruce, A. D.; Wilding, N. B.; Ackland, G. J. Phys. ReV. Lett. 1997, 79, 3002–3005. (47) Bruce, A. D.; Jackson, A. N.; Ackland, G. J.; Wilding, N. B. Phys. ReV. E 2000, 61, 906–919. (48) Berg, B.; Neuhaus, T. Phys. Lett. B. 1991, 267, 249–253. (49) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1988, 61, 2635–2638. (50) Wang, F.; Landau, D. P. Phys. ReV. Lett. 2001, 86, 2050–2053. (51) Wang, J.; Becker, U. Am. Mineral. 2009, 94, 380–386. (52) Free Energy Calculations; Chipot, C., Pohorille, A., Eds.; Chem. Phys. Springer: New York, 2007; pp 1-43. (53) Hummer, G.; Pratt, L. R.; Garcia, A. E. J. Chem. Phys. 1997, 107, 9275–9277. (54) Hummer, G.; Pratt, L. R.; Garcia, A. E.; Berne, B. J.; Rick, S. W. J. Phys. Chem. B 1997, 101, 3017–3020. (55) Kastenholz, M. A.; Huenenberger, P. H. J. Chem. Phys. 2006, 124, 124106. (56) Kastenholz, M. A.; Huenenberger, P. H. J. Chem. Phys. 2006, 124, 224501. (57) Horn, H. W.; Scope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. J. Chem. Phys. 2004, 120, 9665–9678. (58) Bako, I.; Hutter, J.; Palinkas, G. J. Chem. Phys. 2002, 117, 9838– 9843. (59) Megyes, T.; Grosz, T.; Radnai, T.; Bako, I.; Palinkas, G. J. Phys. Chem. A 2004, 108, 7261–7271. (60) Torrie, G. M.; Valleau, J. P. J. Comp. Phys. 1977, 23, 187–199. (61) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A. Comput. Phys. Commun. 2009, 180, 1961–1972. (62) Grossfield, A. WHAM; World Wide Web electronic publication, 2008. (63) Karamertzanis, P. G.; Raiteri, P.; Parrinello, M.; Leslie, M.; Price, S. L. J. Phys. Chem. B 2008, 112, 4298–4308. (64) Kerisit, S.; Parker, S. C.; Harding, J. H. J. Phys. Chem. B 2003, 107, 7676–7682. (65) Spagnoli, D.; Kerisit, S.; Parker, S. C. J. Cryst. Growth 2006, 294, 103–110. (66) Spagnoli, D.; Gilbert, B.; Waychunas, G. A.; Banfield, J. F. Geochim. Cosmochim. Acta 2009, 73, 4023–4033. (67) Perry IV, T. D.; Cygan, R. T.; Mitchell, R. Geochim. Cosmochim. Acta 2007, 71, 5876–5887. (68) Piana, S.; Jones, F.; Gale, J. D. J. Am. Chem. Soc. 2006, 128, 13568–13574. (69) Laio, A.; Parrinello, M. Proc. Natl. Acad. Sci. 2002, 99, 12562– 12566. (70) David, F.; Vokhmin, V.; Ionova, G. J. Mol. Liq. 2001, 90, 45–62. (71) Wu, Y.; Chen, H.; Wang, F.; Paesani, F.; Voth, G. A. J. Phys. Chem. B 2008, 112, 467–482. (72) Kumar, R.; Christie, R. A.; Jordan, K. D. J. Phys. Chem. B 2009, 113, 4111–4118. (73) Quigley, D.; Rodger, P. M. J. Chem. Phys. 2008, 128, 221101. (74) Markgraf, S. A.; Reeder, R. Am. Mineral. 1985, 70, 590–600. (75) de Villiers, J. P. R. Am. Mineral. 1971, 56, 758–767. (76) Kamhi, S. R. Acta Crystallogr. 1963, 16, 770–&. (77) Meyer, H. J. Angew. Chem. 1959, 71, 678–679. (78) Rosseinsky, D. R. Chem. ReV. 1965, 65, 467–&. (79) Gomer, R.; Tryson, G. J. Chem. Phys. 1977, 66, 4413–4424. (80) Marcus, Y. J. Chem. Soc. Faraday Trans. 1991, 87, 2995–2999. (81) Schmid, R.; Miah, A. M.; Sapunov, V. N. Phys. Chem. Chem. Phys. 2000, 2, 97–102. (82) Handbook of Chemistry and Physics, 62nd ed.; Weast, R. C.; Astle, M. J., Eds.; CRC Press: Boca Raton, FL, 1981. (83) Ohtaki, H.; Radnai, T. Chem. ReV. 1993, 93, 1157–1204. (84) Jamieson, J. C. J. Chem. Phys. 1953, 21, 1385–1390.

JP910977A

Derivation of an Accurate Force-Field for Simulating the Growth of

Recommend Documents