Computational Insight into CalciumâSulfate Ion Pair Formation - The

Article Cite This: J. Phys. Chem. C XXXX, XXX, XXX-XXX

pubs.acs.org/JPCC

Computational Insight into Calcium−Sulfate Ion Pair Formation Emily H. Byrne, Paolo Raiteri,* and Julian D. Gale Curtin Institute for Computation, The Institute for Geoscience Research (TIGeR), and Department of Chemistry, Curtin University, PO Box U1987, Perth, Western Australia 6845, Australia S Supporting Information *

ABSTRACT: The thermodynamics of ion pair formation between Ca2+ and SO42− has been studied using a rigid ion force field, the polarizable AMOEBA force field, and ab initio molecular dynamics simulation. The results obtained from the three methods are remarkably similar and consistent with the available experimental data and show that the ion association is driven by an increase in entropy, which can be related to the release of water molecules as previously found for Ca2+ and CO32−. Two new rigid ion force fields targeting different solvation free energies for sulfate have been developed. The comparison between static and dynamic properties of the solvated anion, as well as the pairing free energy with Ca2+, suggest that the model with the strongest solvation is more realistic, which may help to resolve the inconsistency in the current literature.

1. INTRODUCTION Calcium sulfate (CaSO4·xH2O) and its’ hydrates are a family of minerals commonly found in nature, which can take the form of one of three main phases; gypsum (x = 2), bassanite (x = 0.5), and anhydrite (x = 0).1,2 Each of these three phases also has commercial uses, such as construction materials,3 desiccants,4 and medical applications.5 Their precipitation can also be less beneficial since they can be produced as byproducts in industrial processes, such as the formation of gypsum scale in pipes used for desalination and oil/gas extraction.1,2 As shown above, the primary difference between each of the phases is the water content per formula unit of calcium sulfate,1,5,6 with gypsum having the highest water content and anhydrite the lowest. Due to this distinction, the three calcium sulfate phases can be formed through dehydration or hydration of one of the other phases.2,3 Despite the importance of these materials, very little is known about the nucleation and growth of calcium sulfate in water at the molecular level. Water clearly plays a key role in the nucleation process of calcium sulfate as it is often incorporated into the structure of the growing crystal. Although there have been some previous studies relating to the nucleation and growth of these minerals, the majority of them are experimental, and so the information provided on the molecular mechanisms underpinning the nucleation processes is often indirect.2 Knowledge of calcium sulfate nucleation in water would be beneficial not only for a better understanding of the nucleation of minerals in water in general but also from an industrial perspective, where it could be used to develop scale inhibition techniques. Other minerals, such as calcium carbonate,4−6 calcium phosphate,6,7 and silica,6,7 have been studied considerably more and have been found to follow complex indirect pathways involving precursors2,4−10 that then evolve to form the final crystal. © XXXX American Chemical Society

The aim of this work is to take an initial step toward understanding calcium−sulfate association in aqueous solution by providing a thorough description of ion pair formation using both classical and ab initio methods. The first step was to develop a force field capable of capturing the thermodynamics of the calcium and sulfate ions in water, as well as the solubilities of the hydrated crystalline phases. In order to achieve this, we followed a similar approach to our previous studies of calcium carbonate.11,12 This involves the use of lattice dynamics (LD) and molecular dynamics (MD) to calculate the free energies of the crystalline phases and the solvation free energies of the ions in water, respectively. After optimization of the force field parameters to best reproduce the known experimental data, we have studied the thermodynamics of Ca− SO4 ion pair formation using classical MD, including with the polarizable AMOEBA force field. Ab initio molecular dynamics simulations were then also used to calculate the ion pairing free energy using metadynamics13 and to validate the simulations obtained with the classical force fields.

2. METHODS 2.1. Classical Molecular Dynamics. All molecular dynamics simulations with the rigid ion force field developed in this study were performed using the LAMMPS code,14 with a time step of 1 fs in the NVT ensemble at 300 K, except for the solvation free energy and ion pairing free energy simulations, which were performed in the NPT ensemble over the temperature range 290−350 K. A Nosé−Hoover chain thermostat of length 5 was used in all simulations, with the addition of the MTK barostat15 for the NPT simulations. All simulations were performed using a cubic box of side length Received: October 4, 2017 Revised: November 1, 2017 Published: November 1, 2017 A

DOI: 10.1021/acs.jpcc.7b09820 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C ∼50 Å containing 4127 SPC/Fw16 water molecules plus the ions. The free energy perturbation method17,18 was used to calculate the ion solvation free energies, as implemented in LAMMPS. The solute−solvent interactions were gradually switched off in 40 stages of 5 ns each; 20 stages for the electrostatic interactions followed by a further 20 stages for the van der Waals interactions; finite size19 and standard state corrections20 have been included, while the volume correction required for the NPT runs turned out to be negligible for the chosen system size. The ion pairing free energy was calculated from well-tempered21 and multiple walker22 metadynamics simulations using the PLUMED 2.2 plug-in.23 The collective variables used were the calcium−sulfur distance and the calcium−water oxygen coordination number, which was calculated using a switching function of the form n

( ) f (d ) = 1−( ) 1−

d − d0 r0

d − d0 r0

solubility of the sulfate phases was as well reproduced as possible. All intramolecular parameters for the sulfate anion were kept fixed at the values proposed by Allan et al.29 2.3. Ab Initio Molecular Dynamics. In order to validate the classical simulations, ab initio molecular dynamics (AIMD) has also been performed. Here we use Gaussian plane-wave approach as embodied in the QuickStep module of the CP2K code.30 All atoms were represented using the pseudopotentials of Goedecker−Teter−Hutter form31 (with Ca having 10 valence electrons) and a triple-ζ double polarized Gaussian basis set (TZ2P), while the electron density was expanded in an auxiliary basis set with a plane-wave cutoff of 400 Ry. Calculations were performed using the BLYP-D3 functional32,33 (i.e., GGA with Grimme’s dispersion corrections). To partially correct for the systematic overstructuring of water, an elevated temperature of 330 K was employed with deuterium replacing hydrogen. This combination of functional, temperature, and atomic mass has been found to give a more accurate representation of the experimental water radial density distribution function at ambient conditions.34 Self-consistent field was achieved using the orbital transformation algorithm.35 Molecular dynamics was performed in the NVT ensemble using a stochastic thermostat and a time step of 0.5 fs. A cell dimension of 16 Å was chosen containing the CaSO4 ion pair and 134 water molecules for evaluation of the free energies, while a smaller system with a cell of 14.611 Å with SO4 and 100 waters was used for the hydration behavior of the sulfate ion alone. Initial configurations were taken from classical simulations and then re-equilibrated with the ab initio molecular dynamics. Task farming was used to run metadynamics for the ion pair with 10 walkers, each of which was run for at least 110 ps, resulting in a total simulation time of 1.14 ns. Two collective variables were sampled, namely, the Ca−S distance and the Ca−O(water) coordination number using parameters of n = 16 and m = 32 for the powers of distances and an r0 of 3 Å (note that compared to eq 1 there is no d0 in the CP2K implementation of the coordination number). It is important to note that the use of the water coordination number for Ca2+ is essential for determining the ion pairing free energy landscape in AIMD because the rate of water exchange at this cation is still slow compared with the time scale that is accessible to each walker, even if it is relatively fast on an experimental time scale. For the unbiased dynamics of the sulfate anion in water, a simulation of 54 ps was performed. 2.4. Ab Initio Solvation Free Energy. To provide a further benchmark for the solvation free energy of the sulfate anion, this quantity has also been determined using molecular quantum mechanical methods. Here, calculations were performed using QChem 4.136 at the M06/6-31+G** level of theory within the SM8 continuum solvation model.37 The solvation free energy taken is that for the geometry of the sulfate anion after optimization in the presence of the solvation model.

m

(1)

where d0 = 2.1 Å, r0 = 1.0 Å, n = 4, and m = 10. In order to reconstruct the free energy surface, Gaussians of widths 0.2 Å and 0.1 along the distance and coordination number, respectively, were deposited every 1 ps with an initial height of 0.026 eV. A total of 25 walkers were used, resulting in an aggregate run time of 250 ns. A bias factor of 5 was used to gradually reduce the Gaussian heights and ensure convergence of the free energy. All the simulations with the AMOEBA force field24 have been performed with GULP25 using an in-house interface to routines from the Tinker program26 in order to provide the energy and forces for the AMOEBA model. In order to improve the efficiency, the calculation of the AMOEBA energy and forces was parallelized using MPI. Parameters for the sulfate− water system have been taken from the work of Lambrecht et al.27 where the results of the AMOEBA model were benchmarked against those from both density functional and ab initio methods. For the Ca2+ ion, the parameters were taken from the AMOEBA09 library, though we note that other parameter sets for Ca2+ in water are also available.28 Full selfconsistency was used for the multipole moments at each step of the trajectory. All simulations were initially equilibrated in the NPT ensemble for an isotropic cubic cell of approximately 25 Å in size containing 520 water molecules plus the ions using a stochastic barostat and thermostat with a time constant of 0.1 ps. The PLUMED 2.2 plug-in23 was also used to calculate the ion pairing free energy as a function of the same two collective variables using task farming of multiple walkers within GULP. 2.2. Lattice Dynamics. Calculations of the bulk properties of the solid phases were performed using lattice dynamics within the GULP code.25 The free energies of the crystalline phases were computed using the quasi-harmonic approximation at 298.15 K with phonons sampled using a shrinking factor of 6 to define the Monkhorst−Pack mesh. For consistency with the target application in molecular dynamics simulations, the zeropoint energy contribution was excluded from the free energy. Single-point calculation of the free energies was used at the fully optimized bulk structure with respect to the internal energy. Relaxed fitting was used to determine the calcium−sulfate interaction parameters using the anhydrite structures and bulk moduli, where available, as observables. The free energy was then iteratively refined as a further constraint to ensure that the

3. FORCE FIELD VALIDATION As a starting point for the development of a thermodynamically accurate force field for CaSO4 in water, we took our previous force field for Ca2+ with SPC/Fw water12,38 and the parameters for sulfate ions developed by Allan et al.29 The calcium−water interaction parameters were fitted to reproduce the solvation free energy of the ion and consisted of a simple Lennard-Jones potential: B


Article

The Journal of Physical Chemistry C ⎡⎛ σ ⎞12 ⎛ σ ⎞6 ⎤ Uvdw = 4ε⎢⎜ ⎟ − ⎜ ⎟ ⎥ ⎝r⎠ ⎦ ⎣⎝ r ⎠

41

However, in a later publication Marcus suggested that the enthalpic (−1035 kJ/mol) and entropic (−200 J/mol/K) terms are more reliable than the free energy itself, and the reported value, ΔG = −975 kJ/mol at standard conditions, was indeed simply obtained using the definition of the Gibbs free energy and not actually measured. Because of this ambiguity in the literature, we decided to develop two sets of force field parameters (hereafter labeled as FF1 and FF2) that target the two extreme sets of values of the hydration free energy (−1080/−1090 and −975 kJ/mol) and use both to study the Ca−SO4 ion pairing process. The comparison between the classical force fields, the polarizable force field (AMOEBA) and the ab initio calculations reported below, together with experimental data, does however suggest that the more negative solvation free energy is a more reasonable value based on the properties of the sulfate ion in water that result from this choice. The calculated solvation free energies for the sulfate dianion are shown in Figure 1 as a function of temperature and

(2)

The sulfate intramolecular interactions are described by a Morse bond stretching term and a harmonic angle bending term; UMorse bond = D[1 − e−α(rij − r0)]2

(3)

1 kθ(θijk − θ0)2 2

(4)

Uharmonic angle =

where D is the potential well depth, α is the stiffness parameter, rij is the distance between atoms i and j, r0 is the equilibrium bond length, kθ is the angle-bending force constant, θijk is the angle between the atoms at the pivot atom j, and θ0 is the equilibrium angle. The remaining interactions between sulfate−water and Ca− sulfate were described using a Buckingham potential: Ubuck = A e−rij / ρ −

(6)

ΔG = ΔH − T ΔS

C rij 6

(5)

The parameters for these interactions were fitted to reproduce the experimental solvation free energy and the dissolution free energy of calcium sulfate crystalline phases, respectively. The atomic charges (in atomic units) were set to qCa = +2.0, qS = +1.36, qO = −0.84 and qHw = +0.41, qOw = −0.82. For completeness, in Table 1 we report all the parameters used for the classical MD simulations carried out for this work including those taken from previous work. Table 1. Final Force Field Parameters Used in the Present Worka kb (eV/Å2)

harmonic bonds Ow Morse bonds S

Hw D (eV)

O harmonic angles

Ow Hw O S Lennard-Jones 12−6 Ow Ow Buckingham FF1 FF1 FF2 FF2

O Ow O Ow O

Ow Ca O O Ca O Ca

5.0 Hw O

r0 (Å)

45.93 α (Å−1)

1.012 r0 (Å)

1.2 kθ (eV/rad2)

1.505 θ0 (deg)

3.29136 15.0 ε (eV)

0.0067400 0.0009500 A (eV) 103585.02 12534.455133 1815.6986 12634.455133 2345.6986

ρ (Å) 0.200 0.246 0.283373 0.2649 0.283474

Figure 1. Plot showing the temperature dependence of the solvation free energy for FF1 (bottom line, open squares) and FF2 (top line, open circles). Symbols indicate the individual data points, while the lines represent the least-squares fits to the data over the range.

113.24 109.47 σ (Å) 3.165492 3.350000 C (eV Å6)

Table 2. Solvation Free Energy of SO42− in Water and Corresponding Enthalpy and Entropy Contributions Calculated Using the Two Force Fields Developed in This Work (FF1 and FF2) and from Molecular Quantum Mechanical Calculations in the Presence of Continuum Solvent and Previous Literaturea

0.0 0.0 0.0 0.0 0.0

a

FF1 and FF2 are the two sulfate−water Buckingham potentials developed in this work. Here O and Ow denote the oxygens of sulfate and water, respectively. All intermolecular potentials are truncated at 9 Å using a taper function over a range of 3 Å.

force field FF1 FF2 DFT (SM8) Marcus39 Marcus40 Marcus41

3.1. Hydration Free Energy. The works by Marcus are generally regarded as being among the most reliable sources for the solvation free energy of ions in water. At the time this research started, the values for ΔG (−1080/−1090 kJ/mol), ΔH (−1035 kJ/mol), and ΔS (−219/−200 J/mol/K) reported by Marcus39,40 were inconsistent with the definition of the Gibbs free energy:

ΔG (kJ/mol) @ 298.15 K

ΔH (kJ/mol)

ΔS (J/mol/K)

−1081 −967 −1118

−1145 ± 11 −1016 ± 16

−213 ± 34 −164 ± 49

−1090 −1080 −975

−1035 −1035 −1035

−200 −219 −200

a Literature values that are inconsistent with the definition of the Gibbs free energy are shown in italics. The reported errors in the calculated values are the uncertainties of the fitting parameters obtained during a least squares fit.

C


Article

The Journal of Physical Chemistry C summarized in Table 2 for both force fields along with the enthalpies and entropies of solvation obtained from a linear fit of the Gibbs free energy. Naturally the agreement for the solvation free energy at ambient conditions is good since this was the objective of the parameter fitting process. However, the enthalpy and entropy terms, obtained from the fit to the results in Figure 1, are significantly different between the two force fields and show more variable agreement. In particular, FF1 reproduces extremely well the entropic term, while FF2 better reproduces the experimental hydration enthalpy. It is interesting to note that the difference between FF1 and FF2 is not just a rigid shift in the enthalpic term, but there is also a significant change in the temperature dependence of the solvation free energy. The solvation entropy of FF1 (−213 J/ mol/K) is perfectly consistent with the values reported by Marcus for sulfate and for anions of similar charge and shape (SeO42−, CrO42−, MoO42−, and WO42−). Although the values for some of those compounds vary between different papers, they are all in the range −186/−230 J/mol/K.39,40 Therefore, we would argue that FF1 appears to be a better representation of the aqueous sulfate ion and that the solvation enthalpy is probably slightly more negative than the experimentally reported value (−1035 kJ/mol). In order to further discriminate between the two choices of solvation thermodynamics we turned to ab initio calculations and determined the solvation free energy using the parametrized continuum solvent model SM8, which gave ΔG = −1118 kJ/mol. While SM8 is generally more accurate for neutral species, and one or more explicit shells of water may be required to obtain a more reliable free energy of solvation, the difference between the value obtained and the less exothermic hydration of free energy of Marcus is more than twice the typical mean unsigned error for the method.42 This is therefore suggestive that the more exothermic experimental free energy is more consistent with molecular quantum mechanical calculations with continuum solvation. 3.2. Water Coordination Shell. In order to further differentiate between the two force fields developed, we turned our attention to the water structure around SO42− and calculated the radial pair distribution functions (Figure 2) and the 3D density map of water around the anion (Figure 3). In all cases (rigid ion force fields, AMOEBA, and BLYP-D3) the first peak of the S−Ow pair distribution function is between 3.7 and 3.8 Å, with the only exception being for FF2 where it is shifted to 4 Å. The sulfate by water coordination number is also consistent across all models (12.2−13.4 molecules within the first minimum of the g(r)) with FF2 again being the slight exception, having a larger solvation shell that can accommodate one extra water molecule; these numbers compare well with the X-ray diffraction studies reviewed by Ohtaki and Radnai,43 which estimated the S−Ow distance to be in the range 3.7−3.9 Å and the water coordination number in the range 7−12. Our results are also consistent with previous simulations using classical force fields,44−46 polarizable force fields,45 and ab initio47 methods, which report the first peak of the S−Ow pair distribution function to be located at 3.7−3.8 Å and the coordination number to be in the range 12−14. It is particularly instructive to compare the results reported by Williams et al.46 with our FF2 simulations as they have both been developed to give a hydration free energy of approximately −970 kJ/mol. Although the two force fields have similar O−Ow distances (Figure SI1), the S−Ow distance is much shorter for the Williams et al.46 force field, which

Figure 2. Pair distribution functions for (a) S−Ow, (b) S−Hw, (c) O− Ow, and (d) O−Hw obtained using FF1 and FF2 compared to results obtained with the AMOEBA force field and quantum mechanical (BLYP-D3) simulations.

results in a significantly different 3D arrangement of the water molecules around the ion (Figure SI2). Indeed, by using the Williams et al.46 force field we obtained a 3D water structure where there is a high density of water oxygen (red) and hydrogen (white) atoms around the planes that bisect the O− S−O angle. Therefore, the water Ow−Hw bond points preferentially toward the S atom, suggesting that only a very weak hydrogen bond is formed. In the case of FF2 the regions of high probability of finding a water oxygen (red) and hydrogen (white) atoms are spread out over a large area, and there is an almost homogeneous water distribution around the sulfate anion, which again is a reflection of a weak hydrogen bond. However, the BLYP-D3, AMOEBA and FF1 simulations clearly show that the water forms well-defined, strong hydrogen bonds with sulfate and that the water oxygen atoms are localized in a toroid around each of the oxygens of sulfate. This picture is more consistent with the generally accepted structure of water around sulfate and is also consistent with previously published AIMD simulations by Pegado et al.45 where the water oxygen is preferentially located around the sulfate oxygen atoms. 3.3. Self-Diffusion Coefficient. Besides the structural properties, the dynamic behavior of the ions in water is also a key element in the crystallization process. We therefore examined the self-diffusion coefficient of SO42− in water as calculated using both force fields and AIMD simulations and compared it to the literature values (Table 3 and Figure SI3). For the classical and polarizable AMOEBA force fields, the selfdiffusion coefficient at infinite dilution, D0, was calculated by extrapolating the results obtained from three simulations using 3D periodic water boxes of 25, 50, and 75 Å to the limit of infinite box size. The diffusion coefficient from the AIMD was instead calculated only with a box size of 14.611 Å due to the higher computational cost of this method. A correction was then applied to the AIMD results based on the size dependence of the force field results. It should also be noted that the use of D


Article

The Journal of Physical Chemistry C

Table 3. Self-Diffusion Coefficients of SO42− and H2O in Water at Infinite Dilution Computed with the Force Fields and ab Initio MD in the Present Study, and the Comparison to Literature Valuesa self-diffusion coefficient (10−5 cm2/s)

SO42−

H2O

FF1 FF2 AMOEBA BLYP-D3 at 330 K Krynicki et al.64 (exp.) Yuan-Hui and Gregory48 (exp.) Williams and Carbone46 (sim.)

1.191 1.263 1.1b 1.17c

2.86 2.86 2.75 2.7 2.3

1.070 1.489d

2.80

a

The self-diffusion coefficient has been calculated by extrapolating the values obtained using three box sizes: 25, 50, and 75 Å to infinite size. b Obtained using a 25 Å box and corrected by adding the difference between the water self-diffusion coefficient in a 25 Å box and an infinite system (0.610−5 cm2/s). cObtained using a 14.6 Å box. d Williams and Carbone cite a value of 1.20 × 10−5 cm2/s, which has been scaled by a factor of 0.82 that corresponds to the ratio of the experimental value of the self-diffusion coefficient of water and the one obtained from the SPC/E water model.

where the sulfate model that forms the strongest hydrogen bonds has the slower diffusion coefficient. 3.4. Structure and Solubility of the CaSO4 Phases. Having fitted the sulfate−water model, the Ca2+−SO42− interaction parameters for the rigid ion force fields were then optimized to best reproduce the experimental structure and solubilities of the calcium sulfate solid phases (Table 4, Table SI1). Because there are only four partially correlated parameters still to be fitted, namely, A and ρ of the Buckingham potential for the sulfate−sulfate and calcium−sulfate interactions, a degree of compromise is necessary. As per our previous work on the alkaline-earth carbonate force field derivation, the solvation free energy of the solid phases had the largest weight during the fitting process. Hence, some discrepancy in the lattice parameters and bulk moduli are to be expected. The lattice parameters of the five CaSO4 solid phases analyzed here predicted by FF1 are typically within ±2.7% of the experimental values, although for both hydrates the computed crystalline structure appears to be slightly elongated in c and compressed in the other directions. These discrepancies are likely to be related to the fact that the SPC/Fw water potential has been parametrized to reproduce bulk properties of water, and it does not necessarily perform well when the molecules are in a crystalline or cluster environment. We note that the dissolution free energies for gypsum and bassanite calculated using FF1 are in reasonable agreement with the experimental estimates. However, there is a noticeable discrepancy for the anhydrous phases, and unlike in the experiments, these are predicted to be less soluble than gypsum. For the case of FF2, which has a quite different hydration free energy for sulfate, the prediction of the solubilities of the calcium sulfate phases becomes even worse. Indeed, we could only reproduce one phase at a time, and therefore, we report here only the set of parameters that best reproduced gypsum, which we then used for the ion pairing free energy. The reason for the inability of the present classical force field models to reproduce the solubility of all the solid phases simultaneously almost certainly lies in the choice of a rigid ion model and the lack of any polarization effects. The absence of this contribution does not allow the force fields to capture the rearrangement of the charge distribution in the molecules when they go from

Figure 3. Representations of the 3D water structure around the sulfate ion as calculated with FF1 (a), FF2 (b), AMOEBA (c), and BLYP-D3 (d). Iso-surfaces for O and H atoms of water are represented in red and white, respectively. The oxygen and sulfur atoms of sulfate are shown in ball and stick representations, colored in red and yellow, respectively. For all structures, the isovalues of the Ow and Hw atoms are 0.108 and 0.216 atoms/Å3, respectively.

a higher temperature (330 K) and deuterium instead of hydrogen in the BLYP-D3 simulations means that care is needed when making comparisons to the other values. A direct comparison of the self-diffusion coefficient with experiments is somewhat difficult because it is often measured for electrolyte solutions with a mixture of ions. Yuan-Hui and Gregory48 experimentally measured the self-diffusion coefficient to be 1.07 × 10−5 cm2/s for SO42− in seawater at 25 °C, which is very close to our estimated value. However, it is worth noting that the water mobility is a major contributing factor to the diffusion of the solute atoms, and because the water models used here overestimate the self-diffusion by 25%, the sulfate dynamics relative to the water would be slightly slower than in the experiments. The values obtained here are generally better than the uncorrected value reported by Williams and Carbone,46 1.488 × 10−5 cm2/s, and follow an expected trend E


Article


Table 4. Comparison with Experiment of the Structural Properties of the Crystalline Calcium Sulfate Phases Calculated from Lattice Dynamicsa gypsum I12/C1

bassanite I121

insoluble anhydrite Cmcm

soluble anhydrite 1 C222

soluble anhydrite 2 P6222

3

volume (Å ) a (Å) b (Å) c (Å) β (deg) K (GPa) ΔGdiss (kJ/mol) volume (Å3) a (Å) b (Å) c (Å) β (deg) K (GPa) ΔGdiss (kJ/mol) volume (Å3) a (Å) b (Å) c (Å) K (GPa) ΔGdiss (kJ/mol) volume (Å3) a (Å) b (Å) c (Å) K (GPa) ΔGdiss (kJ/mol) volume (Å3) a(Å) c (Å) K (GPa) ΔGdiss (kJ/mol)

experiment

FF1

FF2

AMOEBA

488.829 5.674 15.105 6.491 118.513 43.5 24.91/27.4257 1056.074 12.032 6.927 12.671 90.270

498.087 5.517 15.168 6.656 116.577 42.3 25.38 1087.560 12.169 6.955 12.854 91.435 41.4 19.25 306.810 6.908 6.367 6.976 63.8 39.82 535.398 11.992 6.924 6.448 59.1 31.33 267.699 6.924 6.448 59.1 31.19

565.825 5.702 16.174 6.884 116.980 30.2 25.00 1182.490 12.600 7.217 13.008 91.206 37.3 53.05 327.522 7.056 6.460 7.186 57.8 89.09 570.146 12.331 7.119 6.495 54.3 80.45 285.073 7.119 6.495 54.3 80.45

492.137 5.886 14.596 6.494 118.091 53.9

20.07/19.11 305.481 6.993 6.245b 6.995b 54.9 23.68/26.157 530.856 12.078 6.972 6.304 15.31 265.149 6.969 6.300 10.87

1060.329 11.964 6.800 13.034 90.269 59.1 257.603 6.470 6.115 6.511 162.4 549.439 12.125 7.023 6.453 71.9 225.506 6.470 6.222 150.0

a

The experimental values are taken from Schofield et al.65 (gypsum), Bezou et al.66 (bassanite, C222), Hawthorne et al.67 (Cmcm), and Lager et al.68 (P6222). The experimental values for the bulk moduli (K) of gypsum and Cmcm are taken from Comodi et al.69 and Schwerdtner et al.,70 respectively. The dissolution free energies are from Wagman et al.,43 unless otherwise specified. bThe cell parameters b and c have been switched relative to the original experimental data, and the space group changed from Amma to the alternative setting of Cmcm.

field, FF1, overall, and slightly superior in terms of the cell volume and angles (where not right angles). However, for the anhydrous phases of calcium sulfate, there are substantial errors in the cell parameters, in some cases of up to 7.5%. Interestingly, the results for anhydrous calcium sulfate in the soluble anhydrite 1 structure are quite accurate, while both of the other polymorphs show large errors. This suggests that further refinement of the AMOEBA parameters would be necessary for use in the solid-state under dry conditions. Specifically, it may be necessary to derive individual pairwise terms for the Lennard-Jones interactions rather than rely on combination rules. Beyond structure, the AMOEBA model also seems to systematically overestimate the mechanical hardness of phases based on the current parameters, as judged using the bulk moduli. At present, we cannot comment on the bulk solubilities for AMOEBA for reasons discussed above. However, it should be noted that even first-principles methods have been shown to struggle to yield accurate thermodynamics for hydrated mineral phases.51

solution to a crystalline environment. This problem is particularly acute for hydrate salts where the water model is optimized primarily for the liquid phase and often fails to capture the different properties of water when embedded in an ionic solid. This also has an effect on the structural properties of the CaSO4 phases predicted by FF2 that show errors of up to 7.5% in the lattice parameters and up to 15.75% in the predicted volume of the solid phases. In this regard, both the AMOEBA force field and the AIMD simulations should be superior and better able to reproduce the solubilities of all the structures simultaneously. Because of the much larger computational cost of using these methods, we are unable to obtain solvation free energies at present, thus preventing the calculation of the solubilities of mineral phases. However, recent work has demonstrated that the determination of free energies of solvation for ions using AMOEBA49 and AIMD is becoming possible.50 When analyzing the performance of the AMOEBA model for the structure of calcium sulfate bulk phases, the performance appears to be somewhat variable. Given that the parameters were developed with a focus on species in water, it is perhaps unsurprising that the reproduction of the structure of the hydrate phases is best, being comparable to the rigid ion force

4. ION PAIRING RESULTS Having developed two new candidate classical force fields for CaSO4 and thoroughly benchmarked them against the F


Article


Figure 4. Ca−SO4 ion pairing free energy calculated using FF1 (a), FF2 (b), AMOEBA (c), and BLYP-D3 (d) MD as a function of the Ca−S distance and of the Ca by water coordination number. The free energy maps have been shifted along the y-axis to correct any artifacts in the calculation of the Ca coordination number due to use of two slightly different switching functions, which can either partially include neighbors in the second shell or consider only fractional amounts for the water in the first solvation shell.

AMOEBA force field and ab initio MD, where possible, we now turn our attention to the ion pairing free energy as a thermodynamic property that is not explicitly fitted during the force field parametrization. Here the free energy landscape has been calculated as a function of the calcium−sulfur distance and the calcium−water coordination number. As in our previous work on the ion pairing of alkaline earth metals with carbonate,12 the analysis of the 2D free energy maps (Figure 4) provides a clear picture of the ion pair formation. Although the finer details of the free energy profiles are different, all the methods used here clearly show that the Ca2+ ion in solution has a preferred water coordination number of 7 and that the formation of a contact ion pair (CIP) involves a progressive desolvation of the cation to be surrounded by six and five water molecules to form a monodentate and bidentate ion pair, respectively. The results reported here are qualitatively consistent with previous DFT calculations of the Ca−Cl ion pair formation by Baer and Mundy, which also considered the hydration of the Ca2+ ion alone.52 In both the BLYP-D3 results of this work and the force fields of the present study, the preferred Ca by water coordination number is also 7. However, Baer and Mundy52 found a Ca2+ coordination number of 6 to be more stable than 8 as the secondary favored hydration state, while our force field calculations suggest the reverse (Figure SI4). Due to the limited range of Ca−S distance sampled in our

AIMD, it is not possible to unambiguously determine the relative free energy of the Ca water coordination states for a truly isolated Ca2+ ion. In order to better appreciate the differences between the classical force field and the AIMD, it is instructive to look at the 1D pairing free energy obtained by integrating out the water coordination number (Figure 5). Because of the lower computational cost of the classical and polarizable MD simulations, it is feasible to use large enough cells to converge the long-range tail of the free energy and align them onto the analytic solution for the pairing free energy of two point charges interacting via a screened electrostatic potential (dashed line). However, for the AIMD simulations, we could only explore the contact ion pair (CIP) and solvent-shared ion pair (SSHIP) states, and therefore, we aligned the free energy based on the position of the SSHIP minimum. All methods used in this work clearly indicate that the SSHIP state has a lower free energy than the CIP and has similar barriers for the formation of the CIP and for the dissociation of the complex. However, there are significant differences between the predicted stability of the CIP, with FF1 and the AMOEBA calculation indicating that the CIP is only marginally higher in energy than the SSHIP (∼2 kJ/mol), while BLYP-D3 shows a larger difference (∼7 kJ/mol), and FF2 only has a shallow minimum for the CIP, which is approximately 15 kJ/mol less G


Article


Figure 6. Temperature dependence of the ion pairing free energy compared to the experimental results of Ainsworth,54 Nordstrom et al.,57 Bešter-Rogač,56 and Wachter and Riederer.55

Table 5. Ion Pairing Free Energy of Ca2+−SO42− at 298.15 K (ΔG0), along with the Corresponding Enthalpic and Entropic Contributionsa

Figure 5. One-dimensional Ca−SO4 ion pairing free energy as a function of the Ca−S distance calculated using the classical force fields developed in this work (FF1 and FF2), the AMOEBA force field, and ab initio MD simulations.

FF1 FF2 AMOEBA Ainsworth54 Nordstrom et al.57 Wachter and Riederer55 Bešter-Rogač56

favorable than the SSHIP. Another apparent qualitative difference between the classical and AMOEBA force fields (with the exception of FF1) compared to the AIMD results is the seeming absence of a stable bidentate ion pair minimum in the 1D free energy profiles. However, careful examination of the 2D free energy landscape indicates that sometimes this can be deceptive. For example, for AMOEBA there is a minimum corresponding to 5-fold water coordination by calcium as expected for the bidentate case, but because the monodentate state is lower in free energy over the full distance range, there is no indication of the existence of the higher energy state in the 1D profile, as would be readily apparent if a crossover in stability as a function of distance was to occur. The bottom line remains that the bidentate CIP is unlikely to be of particular importance. The pairing free energy calculations (except BLYP-D3) were carried out at different temperatures in order to extract the enthalpic and entropic contributions. The classical force field simulations have been run at intervals of 5 K between 290 and 350 K (Figure SI5), and the AMOEBA calculations at four different temperatures between 300 and 345 K. The ion-pairing dissociation/association constant and association free energy have then been calculated by integrating the potential of mean force;12,53 K a −1 = K d = C 0

∫R

R1 0

exp(−βϕ(r ))4πr 2 dr

ΔG (kJ/mol)

ΔH (kJ/mol)

−12.55 −12.78 −12.21 −13.18 −13.2 −13.0 −12.9

11.3 ± 0.7 5.7 ± 0.8 7.1 ± 1.8 9.5 ± 0.9 6.9 6.7 5.4

ΔS (J/mol/K) 80 ± 62 ± 64 ± 76 ± 67.0 66.0 61.5

2 2 5 3

The results calculated using the force fields are compared with experimental results by Ainsworth.54 The two ion pairing free energies of FF1 and FF2 correspond to the values obtained from the line of best fit. The entropic term for the CaSO4 association by Nordstrom et al.57 was obtained from the free energy and enthalpy of formation. a

therefore, we cannot directly calculate the association constant for this method. However, by assuming that it connects with the force field free energy curves at the SSHIP, it is possible to estimate that the association constant predicted by AIMD is likely to fall between the values obtained from FF1 and FF2. Both the classical and AMOEBA force fields used in this work predict a positive association enthalpy (5.7/11.3 kJ/mol) and a pairing free energy at standard conditions in the range −12.2/−12.8 kJ/mol, which is perfectly consistent with the available experimental data (Table 5).54−57 These results are also consistent with the similar case of Ca2+ and CO32− association, which we analyzed in detail in ref 58, where the ion pair formation was shown to be driven by the entropy of the solvent molecules that are released during ion pair formation, rather than by the electrostatic interaction between the ions. From the analysis of the metadynamics trajectories obtained using FF1, it is possible to extract the average water coordination number for the ion pair (17 water molecules), and by comparing it with that of the isolated species (7 for Ca and approximately 13 for SO4), we find that on average three water molecules are released into solution (Figure SI6). Analogously, for the AMOEBA and FF2 force fields, four water molecules are released back into solution during the ion pair formation. Hence, the entropic contribution to ion pair formation per water molecule released is in the range 16−26 J/mol/K, which is also consistent with the experimental and computational

(7)

where C0 is a constant used to convert to standard concentration and ϕ(r) = ΔG + kBT ln(4πr2) is the potential of mean force. The resulting temperature-dependent free energies are plotted in Figure 6, while the enthalpy and entropy contributions obtained by fitting these curves are given in Table 5. As noted above, because of the small simulation cell used in the AIMD simulations, it is impossible to align the pairing free energy to the analytic solution in the asymptotic region, and H


Article

The Journal of Physical Chemistry C values for the CaCO3 system (29.3 J/mol/K)58 and other experimental estimates of the water entropy change due its localization in crystalline hydrated phases (32/33.5 J/mol/ K).59,60 This phenomenon of entropy-driven association, which has been experimentally observed in a variety of systems61 and extensively studied with computational methods,62,63 may have important consequences for the thermodynamic stability and abundance of prenucleation clusters, beyond that expected from classical nucleation theory. Among the many computational works in the literature that deal with the formation of ion pairs between ionic species, a few recent ones have started looking at the contribution of water and entropy to the pairing free energy.62,63 In particular, the works by Ballard and Dellago62 and by Shi and Beck,63 which both focused on alkali halide ion pairs, conclude that the ion association is driven by an increase in entropy due to the release of water molecules from the ions’ coordination shell into the solution. This picture is perfectly consistent with the results presented here for CaSO4 and with our previous work on the formation of CaCO3 ion pairs, where ion pair formation is driven by an increase in entropy proportional to the number of water molecules released into solution and opposed by a positive enthalpic term. The unfavorable enthalpy change can be ascribed to the loss of electrostatic interactions with the released water that is not compensated by the cation−anion electrostatic interaction.

formation of the contact ion pair might be a prerequisite to subsequent formation of any pre- or postnucleation species. Future work is required to probe the further ion association processes of calcium sulfate, and this is now possible based on the force field (FF1) developed within this study.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b09820. Comparison between radial pair distribution functions computed with the potential described in this work and literature ones; 3D water density map computed using a force field from the literature; self-diffusion coefficient as a function of the box size; molecular arrangement in the various states of the ion pair formation; all the pairing free energy curves for the two force fields discussed in this work; table with the dissolution free energy and solubilities of the calcium sulfate phases (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: + 61 8 9266 2687. ORCID

5. CONCLUSIONS In conclusion, we have developed and tested two rigid ion force fields by targeting the two available experimental estimates of the hydration free energy of the sulfate anion and the solubility of the most common hydrated solid phases of CaSO4. Because the two available hydration free energies differ by about 100 kJ/ mol, the two force fields gave markedly different predictions for most of the structural and dynamic properties computed in this work. To discriminate between these force field models, we have additionally performed calculations using both a more sophisticated, polarizable force field using the AMOEBA model and ab initio molecular dynamics with a dispersion-corrected GGA functional. Through comparison of the two rigid ion force fields against the results of the other methods, as well as other literature data, it is possible to identify FF1, which gives the more exothermic hydration free energy for sulfate, as being the superior representation of aqueous calcium sulfate systems. In particular, the distribution of water around the sulfate ion is markedly different between FF2 and all of the other results. Having identified the best rigid ion force field model for the system, it was then possible to examine the ion pairing of calcium sulfate in aqueous solution. Based on this force field, as well as AMOEBA and BLYP-D3, there is a consistent picture of the ion pairing process. Consistent with previous estimates of the ion pairing free energy, we find exothermic binding under standard conditions of approximately −13 kJ/mol. By decomposing the free energy into the enthalpy and entropy, based on the temperature dependence, it is found that ion pairing is entropically driven, which is consistent with the release of an average of three water molecules observed during binding. The contact ion pair is found to be between 2 and 8 kJ/mol less stable than the solvent-shared state, with a barrier of at least 14 kJ/mol to reach a distance at which direct coordination between the ions occurs. Therefore, it appears unlikely that the contact state plays a substantial role in the initial ion pairing. However, the balance between these two states may change as further ion association occurs, and so

Paolo Raiteri: 0000-0003-0692-0505 Author Contributions

All authors contributed to performing the calculations, analyzing the data, and writing the manuscript. All authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS E.H.B. would like to acknowledge the contribution of an Australian Government Research Training Program Scholarship in supporting this research, and P.R. and J.D.G. thank the Australian Research Council for financial support through the Discovery Programme (FT130100463, DP160100677). This work was also supported by resources provided by the Pawsey Supercomputing Centre and National Computational Infrastructure with funding from the Australian Government and the Government of Western Australia.

■

REFERENCES

(1) Azimi, G.; Papangelakis, V. G.; Dutrizac, J. E. Modelling of Calcium Sulphate Solubility in Concentrated Multi-Component Sulphate Solutions. Fluid Phase Equilib. 2007, 260 (2), 300−315. (2) Van Driessche, A. E. S.; Benning, L. G.; Rodriguez-Blanco, J. D.; Ossorio, M.; Bots, P.; García-Ruiz, J. M. The Role and Implications of Bassanite as a Stable Precursor Phase to Gypsum Precipitation. Science 2012, 336 (6077), 69−72. (3) Pina, C. M. Nanoscale Dissolution and Growth on Anhydrite Cleavage Faces. Geochim. Cosmochim. Acta 2009, 73 (23), 7034−7044. (4) 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. (5) Gebauer, D.; Voelkel, A.; Coelfen, H. Stable Prenucleation Calcium Carbonate Clusters. Science 2008, 322 (5909), 1819−1822. (6) Gebauer, D.; Kellermeier, M.; Gale, J. D.; Bergström, L.; Cölfen, H. Pre-Nucleation Clusters as Solute Precursors in Crystallisation. Chem. Soc. Rev. 2014, 43 (7), 2348.

I


Article


(29) Allan, N. L.; Rohl, A.; Gay, D.; Catlow, C. R. A.; Davey, R.; Mackrodt, W. Calculater Bulk and Surface Properties of Sulfates. Faraday Discuss. 1993, 95, 273−280. (30) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. Comput. Phys. Commun. 2005, 167 (2), 103−128. (31) Krack, M. Pseudopotentials for H to Kr Optimized for GradientCorrected Exchange-Correlation Functionals. Theor. Chem. Acc. 2005, 114 (1−3), 145−152. (32) Becke, A. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38 (6), 3098−3100. (33) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132 (15), 154104−154120. (34) Bankura, A.; Karmakar, A.; Carnevale, V.; Chandra, A.; Klein, M. L. Structure, Dynamics, and Spectral Diffusion of Water From FirstPrinciples Molecular Dynamics. J. Phys. Chem. C 2014, 118 (50), 29401−29411. (35) VandeVondele, J.; Hutter, J. An Efficient Orbital Transformation Method for Electronic Structure Calculations. J. Chem. Phys. 2003, 118 (10), 4365−4369. (36) Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown, S. T.; Gilbert, A. T. B.; Slipchenko, L. V.; Levchenko, S. V.; O’Neill, D. P.; et al. Advances in Methods and Algorithms in a Modern Quantum Chemistry Program Package. Phys. Chem. Chem. Phys. 2006, 8 (27), 3172−3191. (37) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Perspective on Foundations of Solvation Modeling: the Electrostatic Contribution to the Free Energy of Solvation. J. Chem. Theory Comput. 2008, 4 (6), 877−887. (38) Wu, Y.; Chen, H.; Wang, F.; Paesani, F.; Voth, G. A. An Improved Multistate Empirical Valence Bond Model for Aqueous Proton Solvation and Transport. J. Phys. Chem. B 2008, 112 (2), 467− 482. (39) Marcus, Y. Ion Solvation; Wiley, 1985. (40) Marcus, Y. A Simple Empirical Model Describing the Thermodynamics of Hydration of Ions of Widely Varying Charges, Sizes, and Shapes. Biophys. Chem. 1994, 51 (2), 111−127. (41) Marcus, Y. Ions in Solution and Their Solvation; John Wiley & Sons, Inc., 2016. (42) Cramer, C. J.; Truhlar, D. G. A Universal Approach to Solvation Modeling. Acc. Chem. Res. 2008, 41 (6), 760−768. (43) Ohtaki, H.; Radnai, T. Structure and Dynamics of Hydrated Ions. Chem. Rev. 1993, 93 (3), 1157−1204. (44) Cannon, W.; Pettitt, B.; McCammon, J. A. Sulfate Anion in Water - Model Structural, Thermodynamic, and Dynamic Properties. J. Phys. Chem. 1994, 98 (24), 6225−6230. (45) Pegado, L.; Marsalek, O.; Jungwirth, P.; Wernersson, E. Solvation and Ion-Pairing Properties of the Aqueous Sulfate Anion: Explicit Versus Effective Electronic Polarization. Phys. Chem. Chem. Phys. 2012, 14 (29), 10248. (46) Williams, C. D.; Carbone, P. A Classical Force Field for Tetrahedral Oxyanions Developed Using Hydration Properties: the Examples of Pertechnetate (TcO4−) And Sulfate (SO42−). J. Chem. Phys. 2015, 143 (17), 174502−174509. (47) Vchirawongkwin, V.; Rode, B. M.; Persson, I. Structure and Dynamics of Sulfate Ion in Aqueous Solution-an Ab Initio QMCF MD Simulation and Large Angle X-Ray Scattering Study. J. Phys. Chem. B 2007, 111 (16), 4150−4155. (48) Yuan-Hui, L.; Gregory, S. Diffusion of Ions in Sea Water and in Deep-Sea Sediments. Geochim. Cosmochim. Acta 1974, 38 (5), 703− 714. (49) Mohamed, N. A.; Bradshaw, R. T.; Essex, J. W. Evaluation of Solvation Free Energies for Small Molecules with the AMOEBA Polarizable Force Field. J. Comput. Chem. 2016, 37 (32), 2749−2758.

(7) Navrotsky, A. Energetic Clues to Pathways to Biomineralization: Precursors, Clusters, and Nanoparticles. Proc. Natl. Acad. Sci. U. S. A. 2004, 101 (33), 12096−12101. (8) Raiteri, P.; Gale, J. D. Water Is the Key to Nonclassical Nucleation of Amorphous Calcium Carbonate. J. Am. Chem. Soc. 2010, 132 (49), 17623−17634. (9) De Yoreo, J. J.; Vekilov, P. G. Principles of Crystal Nucleation and Growth. Rev. Mineral. Geochem. 2003, 54, 57−93. (10) Erdemir, D.; Lee, A. Y.; Myerson, A. S. Nucleation of Crystals From Solution: Classical and Two-Step Models. Acc. Chem. Res. 2009, 42 (5), 621−629. (11) Raiteri, P.; Gale, J. D.; Quigley, D.; Rodger, P. M. Derivation of an Accurate Force-Field for Simulating the Growth of Calcium Carbonate From Aqueous Solution: a New Model for the Calcite− Water Interface. J. Phys. Chem. C 2010, 114 (13), 5997−6010. (12) 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 (43), 24447−24458. (13) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99 (20), 12562−12566. (14) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117 (1), 1−19. (15) Martyna, G.; Tobias, D.; Klein, M. L. Constant Pressure Molecular Dynamics Algorithms. J. Chem. Phys. 1994, 101, 4177. (16) Wu, Y.; Tepper, H.; Voth, G. A. Flexible Simple Point-Charge Water Model with Improved Liquid-State Properties. J. Chem. Phys. 2006, 124 (2), 024503. (17) Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. J. Chem. Phys. 1935, 3 (5), 300. (18) Zwanzig, R. High-Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases. J. Chem. Phys. 1954, 22 (8), 1420. (19) Hummer, G.; Pratt, L. R.; Garcia, A. Ion Sizes and Finite-Size Corrections for Ionic-Solvation Free Energies. J. Chem. Phys. 1997, 107 (21), 9275−9277. (20) Kastenholz, M. A.; Hünenberger, P. H. Computation of Methodology-Independent Ionic Solvation Free Energies From Molecular Simulations. II. the Hydration Free Energy of the Sodium Cation. J. Chem. Phys. 2006, 124 (22), 224501. (21) Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: a Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100 (2), 20603. (22) Raiteri, P.; Laio, A.; Gervasio, F. L.; Micheletti, C.; Parrinello, M. Efficient Reconstruction of Complex Free Energy Landscapes by Multiple Walkers Metadynamics. J. Phys. Chem. B 2006, 110 (8), 3533−3539. (23) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. Computer Physics Communications. Comput. Phys. Commun. 2014, 185 (2), 604−613. (24) Ponder, J. W.; Wu, C.; Ren, P.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A., Jr.; et al. Current Status of the AMOEBA Polarizable Force Field. J. Phys. Chem. B 2010, 114 (8), 2549−2564. (25) Gale, J. D. GULP: Capabilities and Prospects. Z. Kristallogr. Cryst. Mater. 2005, 220 (5−6), 552−554. (26) Ponder, J. W.; Richards, F. M. An Efficient Newton-Like Method for Molecular Mechanics Energy Minimization of Large Molecules. J. Comput. Chem. 1987, 8, 1016. (27) Piquemal, J.-P.; Perera, L.; Cisneros, G. A. S.; Ren, P.; Pedersen, L. G.; Darden, T. A. Towards Accurate Solvation Dynamics of Divalent Cations in Water Using the Polarizable Amoeba Force Field: From Energetics to Structure. J. Chem. Phys. 2006, 125 (5), 054511− 054517. (28) Lambrecht, D. S.; Clark, G. N. I.; Head-Gordon, T.; HeadGordon, M. Exploring the Rich Energy Landscape of Sulfate. J. Phys. Chem. A 2011, 115, 11438−11454. J


Article

The Journal of Physical Chemistry C (50) Bell, D. R.; Qi, R.; Jing, Z.; Xiang, J. Y.; Mejias, C.; Schnieders, M. J.; Ponder, J. W.; Ren, P. Calculating Binding Free Energies of Host-Guest Systems Using the AMOEBA Polarizable Force Field. Phys. Chem. Chem. Phys. 2016, 18 (44), 30261−30269. (51) Demichelis, R.; Raiteri, P.; Gale, J. D.; Dovesi, R. Examining the Accuracy of Density Functional Theory for Predicting the Thermodynamics of Water Incorporation Into Minerals: the Hydrates of Calcium Carbonate. J. Phys. Chem. C 2013, 117 (34), 17814−17823. (52) Baer, M. D.; Mundy, C. J. Local Aqueous Solvation Structure Around Ca 2+During Ca 2+···Cl − Pair Formation. J. Phys. Chem. B 2016, 120 (8), 1885−1893. (53) Fuoss, R. M.; Kraus, C. A. Properties of Electrolytic Solutions. III. the Dissociation Constant. J. Am. Chem. Soc. 1933, 55 (3), 1019− 1028. (54) Ainsworth, R. G. Dissociation Constant of Calcium Sulphate From 25 to 50°C. J. Chem. Soc., Faraday Trans. 1 1973, 69 (0), 1028. (55) Wachter, R.; Riederer, K. Properties of Dilute ElectrolyteSolutions From Calorimetric Measurements. Pure Appl. Chem. 1981, 53 (7), 1301−1312. (56) Bester-Rogac, M. Determination of the Limiting Conductances and the Ion-Association Constants of Calcium and Manganese Sulfates in Water From Electrical Conductivity Measurements. Acta Chim. Slov. 2008, 55 (1), 201−208. (57) Nordstrom, D. K.; Plummer, L. N.; Langmuir, D.; Busenberg, E.; May, H. M.; Jones, B. F.; Parkhurst, D. L. Revised Chemical Equilibrium Data for Major Water−Mineral Reactions and Their Limitations. In Chemical Modeling of Aqueous Systems II; ACS Symposium Series; American Chemical Society: Washington, DC, 2009; Vol. 416, pp 398−413. (58) Kellermeier, M.; Raiteri, P.; Berg, J. K.; Kempter, A.; Gale, J. D.; Gebauer, D. Entropy Drives Calcium Carbonate Ion Association. ChemPhysChem 2016, 17 (21), 3535−3541. (59) Dunitz, J. D. The Entropic Cost of Bound Water in Crystals and Biomolecules. Science 1994, 264 (5159), 670−670. (60) Konigsberger, E.; Konigsberger, L.; Gamsjager, H. LowTemperature Thermodynamic Model for the System Na2Co3MgCO3-CaCO3-H2O. Geochim. Cosmochim. Acta 1999, 63 (19−20), 3105−3119. (61) De Visscher, A.; Vanderdeelen, J.; Königsberger, E.; Churagulov, B. R.; Ichikuni, M.; Tsurumi, M. IUPAC-NIST Solubility Data Series. 95. Alkaline Earth Carbonates in Aqueous Systems. Part 1. Introduction, Be and Mg. J. Phys. Chem. Ref. Data 2012, 41 (1), 013105. (62) Ballard, A. J.; Dellago, C. Toward the Mechanism of Ionic Dissociation in Water. J. Phys. Chem. B 2012, 116 (45), 13490−13497. (63) Shi, Y.; Beck, T. Deconstructing Free Energies in the Law of Matching Water Affinities. J. Phys. Chem. B 2017, 121 (9), 2189−2201. (64) Krynicki, K.; Green, C. D.; Sawyer, D. W. Pressure and Temperature Dependence of Self-Diffusion in Water. Faraday Discuss. Chem. Soc. 1978, 66, 199−208. (65) Schofield, P. F.; Knight, K. S.; Stretton, I. C. Thermal Expansion of Gypsum Investigated by Neutron Powder Diffraction. Am. Mineral. 1996, 81 (7−8), 847−851. (66) Bezou, C.; Nonat, A.; Mutin, J. C.; Christensen, A. N.; Lehmann, M. S. Investigation of the Crystal Structure of Γ-CaSO4, CaSO4· 0.5 H2O, and CaSO4· 0.6 H2O by Powder Diffraction Methods. J. Solid State Chem. 1995, 117, 165−176. (67) Hawthorne, F. C.; Ferguson, R. B. Anhydrous Sulphates; II, Refinement of the Crystal Structure of Anhydrite. Can. Mineral. 1975, 13, 289−292. (68) Lager, G. A.; Armbruster, T.; Rotella, F. J.; Jorgensen, J. D.; Hinks, D. G. A Crystallographic Study of the Low-Temperature Dehydration Products of Gypsum, CaSO4· 2h2O - Hemihydrate CaSO4· 0.50h2O, and Γ-CaSO4. Am. Mineral. 1984, 69 (9−10), 910− 919. (69) Comodi, P.; Nazzareni, S.; Zanazzi, P. F.; Speziale, S. HighPressure Behavior of Gypsum: a Single-Crystal X-Ray Study. Am. Mineral. 2008, 93 (10), 1530−1537.

(70) Schwerdtner, W. M.; Tou, J.; Hertz, P. B. Elastic Properties of Single Crystals of Anhydrite. Can. J. Earth Sci. 1965, 2 (6), 673−683.

K


Computational Insight into CalciumâSulfate Ion Pair Formation - The

Recommend Documents

Computational Insight into CalciumâSulfate Ion Pair Formation - The