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 The aim of this work is to take an initial step toward

understanding calcium−sulfate association in aqueous solution
Calcium sulfate (CaSO4·xH2O) and its’ hydrates are a family of
by providing a thorough description of ion pair formation using
minerals commonly found in nature, which can take the form of
both classical and ab initio methods. The first step was to
one of three main phases; gypsum (x = 2), bassanite (x = 0.5),
develop a force field capable of capturing the thermodynamics
and anhydrite (x = 0).1,2 Each of these three phases also has of the calcium and sulfate ions in water, as well as the
commercial uses, such as construction materials,3 desiccants,4 solubilities of the hydrated crystalline phases. In order to
and medical applications.5 Their precipitation can also be less achieve this, we followed a similar approach to our previous
beneficial since they can be produced as byproducts in studies of calcium carbonate.11,12 This involves the use of lattice
industrial processes, such as the formation of gypsum scale in dynamics (LD) and molecular dynamics (MD) to calculate the
pipes used for desalination and oil/gas extraction.1,2 As shown free energies of the crystalline phases and the solvation free
above, the primary difference between each of the phases is the energies of the ions in water, respectively. After optimization of
water content per formula unit of calcium sulfate,1,5,6 with the force field parameters to best reproduce the known
gypsum having the highest water content and anhydrite the experimental data, we have studied the thermodynamics of Ca−
lowest. Due to this distinction, the three calcium sulfate phases SO4 ion pair formation using classical MD, including with the
can be formed through dehydration or hydration of one of the polarizable AMOEBA force field. Ab initio molecular dynamics
other phases.2,3 simulations were then also used to calculate the ion pairing free
Despite the importance of these materials, very little is energy using metadynamics13 and to validate the simulations
known about the nucleation and growth of calcium sulfate in obtained with the classical force fields.
water at the molecular level. Water clearly plays a key role in
the nucleation process of calcium sulfate as it is often 2. METHODS
incorporated into the structure of the growing crystal. Although 2.1. Classical Molecular Dynamics. All molecular
there have been some previous studies relating to the dynamics simulations with the rigid ion force field developed
nucleation and growth of these minerals, the majority of in this study were performed using the LAMMPS code,14 with a
them are experimental, and so the information provided on the time step of 1 fs in the NVT ensemble at 300 K, except for the
molecular mechanisms underpinning the nucleation processes solvation free energy and ion pairing free energy simulations,
is often indirect.2 Knowledge of calcium sulfate nucleation in which were performed in the NPT ensemble over the
water would be beneficial not only for a better understanding of temperature range 290−350 K. A Nosé−Hoover chain
the nucleation of minerals in water in general but also from an thermostat of length 5 was used in all simulations, with the
industrial perspective, where it could be used to develop scale addition of the MTK barostat15 for the NPT simulations. All
inhibition techniques. Other minerals, such as calcium simulations were performed using a cubic box of side length
carbonate,4−6 calcium phosphate,6,7 and silica,6,7 have been
studied considerably more and have been found to follow Received: October 4, 2017
complex indirect pathways involving precursors2,4−10 that then Revised: November 1, 2017
evolve to form the final crystal. Published: November 1, 2017

© XXXX American Chemical Society A DOI: 10.1021/acs.jpcc.7b09820

J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C Article

∼50 Å containing 4127 SPC/Fw16 water molecules plus the solubility of the sulfate phases was as well reproduced as
ions. The free energy perturbation method17,18 was used to possible. All intramolecular parameters for the sulfate anion
calculate the ion solvation free energies, as implemented in were kept fixed at the values proposed by Allan et al.29
LAMMPS. The solute−solvent interactions were gradually 2.3. Ab Initio Molecular Dynamics. In order to validate
switched off in 40 stages of 5 ns each; 20 stages for the the classical simulations, ab initio molecular dynamics (AIMD)
electrostatic interactions followed by a further 20 stages for the has also been performed. Here we use Gaussian plane-wave
van der Waals interactions; finite size19 and standard state approach as embodied in the QuickStep module of the CP2K
corrections20 have been included, while the volume correction code.30 All atoms were represented using the pseudopotentials
required for the NPT runs turned out to be negligible for the of Goedecker−Teter−Hutter form31 (with Ca having 10
chosen system size. The ion pairing free energy was calculated valence electrons) and a triple-ζ double polarized Gaussian
from well-tempered21 and multiple walker22 metadynamics basis set (TZ2P), while the electron density was expanded in an
simulations using the PLUMED 2.2 plug-in.23 The collective auxiliary basis set with a plane-wave cutoff of 400 Ry.
variables used were the calcium−sulfur distance and the Calculations were performed using the BLYP-D3 functional32,33
calcium−water oxygen coordination number, which was (i.e., GGA with Grimme’s dispersion corrections). To partially
calculated using a switching function of the form correct for the systematic overstructuring of water, an elevated
n temperature of 330 K was employed with deuterium replacing

f (d ) =
1−( ) d − d0
r0
hydrogen. This combination of functional, temperature, and
atomic mass has been found to give a more accurate
m
1−( ) d − d0 representation of the experimental water radial density
r0 (1) distribution function at ambient conditions.34 Self-consistent
where d0 = 2.1 Å, r0 = 1.0 Å, n = 4, and m = 10. In order to field was achieved using the orbital transformation algorithm.35
reconstruct the free energy surface, Gaussians of widths 0.2 Å Molecular dynamics was performed in the NVT ensemble
and 0.1 along the distance and coordination number, using a stochastic thermostat and a time step of 0.5 fs. A cell
respectively, were deposited every 1 ps with an initial height dimension of 16 Å was chosen containing the CaSO4 ion pair
of 0.026 eV. A total of 25 walkers were used, resulting in an and 134 water molecules for evaluation of the free energies,
aggregate run time of 250 ns. A bias factor of 5 was used to while a smaller system with a cell of 14.611 Å with SO4 and 100
gradually reduce the Gaussian heights and ensure convergence waters was used for the hydration behavior of the sulfate ion
of the free energy. alone. Initial configurations were taken from classical
All the simulations with the AMOEBA force field24 have simulations and then re-equilibrated with the ab initio
been performed with GULP25 using an in-house interface to molecular dynamics.
routines from the Tinker program26 in order to provide the Task farming was used to run metadynamics for the ion pair
energy and forces for the AMOEBA model. In order to improve with 10 walkers, each of which was run for at least 110 ps,
the efficiency, the calculation of the AMOEBA energy and resulting in a total simulation time of 1.14 ns. Two collective
forces was parallelized using MPI. Parameters for the sulfate− variables were sampled, namely, the Ca−S distance and the
water system have been taken from the work of Lambrecht et Ca−O(water) coordination number using parameters of n = 16
al.27 where the results of the AMOEBA model were and m = 32 for the powers of distances and an r0 of 3 Å (note
benchmarked against those from both density functional and that compared to eq 1 there is no d0 in the CP2K
ab initio methods. For the Ca2+ ion, the parameters were taken implementation of the coordination number). It is important
from the AMOEBA09 library, though we note that other to note that the use of the water coordination number for Ca2+
parameter sets for Ca2+ in water are also available.28 Full self- is essential for determining the ion pairing free energy
consistency was used for the multipole moments at each step of landscape in AIMD because the rate of water exchange at
the trajectory. All simulations were initially equilibrated in the this cation is still slow compared with the time scale that is
NPT ensemble for an isotropic cubic cell of approximately 25 Å accessible to each walker, even if it is relatively fast on an
in size containing 520 water molecules plus the ions using a experimental time scale. For the unbiased dynamics of the
stochastic barostat and thermostat with a time constant of 0.1 sulfate anion in water, a simulation of 54 ps was performed.
ps. The PLUMED 2.2 plug-in23 was also used to calculate the 2.4. Ab Initio Solvation Free Energy. To provide a
ion pairing free energy as a function of the same two collective further benchmark for the solvation free energy of the sulfate
variables using task farming of multiple walkers within GULP. anion, this quantity has also been determined using molecular
2.2. Lattice Dynamics. Calculations of the bulk properties quantum mechanical methods. Here, calculations were
of the solid phases were performed using lattice dynamics performed using QChem 4.136 at the M06/6-31+G** level
within the GULP code.25 The free energies of the crystalline of theory within the SM8 continuum solvation model.37 The
phases were computed using the quasi-harmonic approximation solvation free energy taken is that for the geometry of the
at 298.15 K with phonons sampled using a shrinking factor of 6 sulfate anion after optimization in the presence of the solvation
to define the Monkhorst−Pack mesh. For consistency with the model.
target application in molecular dynamics simulations, the zero-
point energy contribution was excluded from the free energy. 3. FORCE FIELD VALIDATION
Single-point calculation of the free energies was used at the As a starting point for the development of a thermodynamically
fully optimized bulk structure with respect to the internal accurate force field for CaSO4 in water, we took our previous
energy. force field for Ca2+ with SPC/Fw water12,38 and the parameters
Relaxed fitting was used to determine the calcium−sulfate for sulfate ions developed by Allan et al.29 The calcium−water
interaction parameters using the anhydrite structures and bulk interaction parameters were fitted to reproduce the solvation
moduli, where available, as observables. The free energy was free energy of the ion and consisted of a simple Lennard-Jones
then iteratively refined as a further constraint to ensure that the potential:
B DOI: 10.1021/acs.jpcc.7b09820
J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C Article

⎡⎛ σ ⎞12 ⎛ σ ⎞6 ⎤ ΔG = ΔH − T ΔS (6)
Uvdw = 4ε⎢⎜ ⎟ − ⎜ ⎟ ⎥ 41
⎣⎝ r ⎠ ⎝r⎠ ⎦ (2)
However, in a later publication Marcus suggested that the
enthalpic (−1035 kJ/mol) and entropic (−200 J/mol/K) terms
The sulfate intramolecular interactions are described by a are more reliable than the free energy itself, and the reported
Morse bond stretching term and a harmonic angle bending value, ΔG = −975 kJ/mol at standard conditions, was indeed
term; simply obtained using the definition of the Gibbs free energy
and not actually measured. Because of this ambiguity in the
UMorse bond = D[1 − e−α(rij − r0)]2 (3) literature, we decided to develop two sets of force field
parameters (hereafter labeled as FF1 and FF2) that target the
1
Uharmonic angle = kθ(θijk − θ0)2 two extreme sets of values of the hydration free energy
2 (4) (−1080/−1090 and −975 kJ/mol) and use both to study the
where D is the potential well depth, α is the stiffness parameter, Ca−SO4 ion pairing process. The comparison between the
rij is the distance between atoms i and j, r0 is the equilibrium classical force fields, the polarizable force field (AMOEBA) and
bond length, kθ is the angle-bending force constant, θijk is the the ab initio calculations reported below, together with
angle between the atoms at the pivot atom j, and θ0 is the experimental data, does however suggest that the more negative
equilibrium angle. solvation free energy is a more reasonable value based on the
The remaining interactions between sulfate−water and Ca− properties of the sulfate ion in water that result from this
sulfate were described using a Buckingham potential: choice.
The calculated solvation free energies for the sulfate dianion
C
Ubuck = A e−rij / ρ − are shown in Figure 1 as a function of temperature and
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
harmonic bonds kb (eV/Å2) r0 (Å)
Ow Hw 45.93 1.012
Morse bonds D (eV) α (Å−1) r0 (Å)
S O 5.0 1.2 1.505
Figure 1. Plot showing the temperature dependence of the solvation
harmonic angles kθ (eV/rad2) θ0 (deg)
free energy for FF1 (bottom line, open squares) and FF2 (top line,
Hw Ow Hw 3.29136 113.24 open circles). Symbols indicate the individual data points, while the
O S O 15.0 109.47 lines represent the least-squares fits to the data over the range.
Lennard-Jones 12−6 ε (eV) σ (Å)
Ow Ow 0.0067400 3.165492
Ow Ca 0.0009500 3.350000
Buckingham A (eV) ρ (Å) C (eV Å6) Table 2. Solvation Free Energy of SO42− in Water and
O O 103585.02 0.200 0.0 Corresponding Enthalpy and Entropy Contributions
FF1 Ow O 12534.455133 0.246 0.0 Calculated Using the Two Force Fields Developed in This
FF1 O Ca 1815.6986 0.283373 0.0 Work (FF1 and FF2) and from Molecular Quantum
FF2 Ow O 12634.455133 0.2649 0.0 Mechanical Calculations in the Presence of Continuum
FF2 O Ca 2345.6986 0.283474 0.0 Solvent and Previous Literaturea
a
FF1 and FF2 are the two sulfate−water Buckingham potentials force field ΔG (kJ/mol) @ 298.15 K ΔH (kJ/mol) ΔS (J/mol/K)
developed in this work. Here O and Ow denote the oxygens of sulfate
FF1 −1081 −1145 ± 11 −213 ± 34
and water, respectively. All intermolecular potentials are truncated at 9
Å using a taper function over a range of 3 Å. FF2 −967 −1016 ± 16 −164 ± 49
DFT −1118
(SM8)
Marcus39 −1090 −1035 −200
3.1. Hydration Free Energy. The works by Marcus are Marcus40 −1080 −1035 −219
generally regarded as being among the most reliable sources for Marcus41 −975 −1035 −200
the solvation free energy of ions in water. At the time this
research started, the values for ΔG (−1080/−1090 kJ/mol), a
Literature values that are inconsistent with the definition of the Gibbs
ΔH (−1035 kJ/mol), and ΔS (−219/−200 J/mol/K) reported free energy are shown in italics. The reported errors in the calculated
by Marcus39,40 were inconsistent with the definition of the values are the uncertainties of the fitting parameters obtained during a
Gibbs free energy: least squares fit.

C DOI: 10.1021/acs.jpcc.7b09820
J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C Article

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). Figure 2. Pair distribution functions for (a) S−Ow, (b) S−Hw, (c) O−
In order to further discriminate between the two choices of Ow, and (d) O−Hw obtained using FF1 and FF2 compared to results
solvation thermodynamics we turned to ab initio calculations obtained with the AMOEBA force field and quantum mechanical
and determined the solvation free energy using the para- (BLYP-D3) simulations.
metrized continuum solvent model SM8, which gave ΔG =
−1118 kJ/mol. While SM8 is generally more accurate for results in a significantly different 3D arrangement of the water
neutral species, and one or more explicit shells of water may be molecules around the ion (Figure SI2). Indeed, by using the
required to obtain a more reliable free energy of solvation, the Williams et al.46 force field we obtained a 3D water structure
difference between the value obtained and the less exothermic where there is a high density of water oxygen (red) and
hydration of free energy of Marcus is more than twice the hydrogen (white) atoms around the planes that bisect the O−
typical mean unsigned error for the method.42 This is therefore S−O angle. Therefore, the water Ow−Hw bond points
suggestive that the more exothermic experimental free energy is preferentially toward the S atom, suggesting that only a very
more consistent with molecular quantum mechanical calcu- weak hydrogen bond is formed. In the case of FF2 the regions
lations with continuum solvation. of high probability of finding a water oxygen (red) and
3.2. Water Coordination Shell. In order to further hydrogen (white) atoms are spread out over a large area, and
differentiate between the two force fields developed, we turned there is an almost homogeneous water distribution around the
our attention to the water structure around SO42− and sulfate anion, which again is a reflection of a weak hydrogen
calculated the radial pair distribution functions (Figure 2) bond. However, the BLYP-D3, AMOEBA and FF1 simulations
and the 3D density map of water around the anion (Figure 3). clearly show that the water forms well-defined, strong hydrogen
In all cases (rigid ion force fields, AMOEBA, and BLYP-D3) bonds with sulfate and that the water oxygen atoms are
the first peak of the S−Ow pair distribution function is between localized in a toroid around each of the oxygens of sulfate. This
3.7 and 3.8 Å, with the only exception being for FF2 where it is picture is more consistent with the generally accepted structure
shifted to 4 Å. The sulfate by water coordination number is also of water around sulfate and is also consistent with previously
consistent across all models (12.2−13.4 molecules within the published AIMD simulations by Pegado et al.45 where the water
first minimum of the g(r)) with FF2 again being the slight oxygen is preferentially located around the sulfate oxygen
exception, having a larger solvation shell that can accommodate atoms.
one extra water molecule; these numbers compare well with the 3.3. Self-Diffusion Coefficient. Besides the structural
X-ray diffraction studies reviewed by Ohtaki and Radnai,43 properties, the dynamic behavior of the ions in water is also a
which estimated the S−Ow distance to be in the range 3.7−3.9 key element in the crystallization process. We therefore
Å and the water coordination number in the range 7−12. Our examined the self-diffusion coefficient of SO42− in water as
results are also consistent with previous simulations using calculated using both force fields and AIMD simulations and
classical force fields,44−46 polarizable force fields,45 and ab compared it to the literature values (Table 3 and Figure SI3).
initio47 methods, which report the first peak of the S−Ow pair For the classical and polarizable AMOEBA force fields, the self-
distribution function to be located at 3.7−3.8 Å and the diffusion coefficient at infinite dilution, D0, was calculated by
coordination number to be in the range 12−14. extrapolating the results obtained from three simulations using
It is particularly instructive to compare the results reported 3D periodic water boxes of 25, 50, and 75 Å to the limit of
by Williams et al.46 with our FF2 simulations as they have both infinite box size. The diffusion coefficient from the AIMD was
been developed to give a hydration free energy of instead calculated only with a box size of 14.611 Å due to the
approximately −970 kJ/mol. Although the two force fields higher computational cost of this method. A correction was
have similar O−Ow distances (Figure SI1), the S−Ow distance then applied to the AIMD results based on the size dependence
is much shorter for the Williams et al.46 force field, which of the force field results. It should also be noted that the use of
D DOI: 10.1021/acs.jpcc.7b09820
J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C Article

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 1.191 2.86
FF2 1.263 2.86
AMOEBA 1.1b 2.75
BLYP-D3 at 330 K 1.17c 2.7
Krynicki et al.64 (exp.) 2.3
Yuan-Hui and Gregory48 (exp.) 1.070
Williams and Carbone46 (sim.) 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
Figure 3. Representations of the 3D water structure around the sulfate
crystalline structure appears to be slightly elongated in c and
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 compressed in the other directions. These discrepancies are
and white, respectively. The oxygen and sulfur atoms of sulfate are likely to be related to the fact that the SPC/Fw water potential
shown in ball and stick representations, colored in red and yellow, has been parametrized to reproduce bulk properties of water,
respectively. For all structures, the isovalues of the Ow and Hw atoms and it does not necessarily perform well when the molecules are
are 0.108 and 0.216 atoms/Å3, respectively. in a crystalline or cluster environment.
We note that the dissolution free energies for gypsum and
a higher temperature (330 K) and deuterium instead of bassanite calculated using FF1 are in reasonable agreement with
hydrogen in the BLYP-D3 simulations means that care is the experimental estimates. However, there is a noticeable
needed when making comparisons to the other values. discrepancy for the anhydrous phases, and unlike in the
A direct comparison of the self-diffusion coefficient with experiments, these are predicted to be less soluble than gypsum.
experiments is somewhat difficult because it is often measured For the case of FF2, which has a quite different hydration free
for electrolyte solutions with a mixture of ions. Yuan-Hui and energy for sulfate, the prediction of the solubilities of the
Gregory48 experimentally measured the self-diffusion coefficient calcium sulfate phases becomes even worse. Indeed, we could
to be 1.07 × 10−5 cm2/s for SO42− in seawater at 25 °C, which only reproduce one phase at a time, and therefore, we report
is very close to our estimated value. However, it is worth noting here only the set of parameters that best reproduced gypsum,
that the water mobility is a major contributing factor to the which we then used for the ion pairing free energy. The reason
diffusion of the solute atoms, and because the water models for the inability of the present classical force field models to
used here overestimate the self-diffusion by 25%, the sulfate reproduce the solubility of all the solid phases simultaneously
dynamics relative to the water would be slightly slower than in almost certainly lies in the choice of a rigid ion model and the
the experiments. The values obtained here are generally better lack of any polarization effects. The absence of this contribution
than the uncorrected value reported by Williams and does not allow the force fields to capture the rearrangement of
Carbone,46 1.488 × 10−5 cm2/s, and follow an expected trend the charge distribution in the molecules when they go from
E DOI: 10.1021/acs.jpcc.7b09820
J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C Article

Table 4. Comparison with Experiment of the Structural Properties of the Crystalline Calcium Sulfate Phases Calculated from
Lattice Dynamicsa
experiment FF1 FF2 AMOEBA
3
gypsum volume (Å ) 488.829 498.087 565.825 492.137
I12/C1 a (Å) 5.674 5.517 5.702 5.886
b (Å) 15.105 15.168 16.174 14.596
c (Å) 6.491 6.656 6.884 6.494
β (deg) 118.513 116.577 116.980 118.091
K (GPa) 43.5 42.3 30.2 53.9
ΔGdiss (kJ/mol) 24.91/27.4257 25.38 25.00
bassanite volume (Å3) 1056.074 1087.560 1182.490 1060.329
I121 a (Å) 12.032 12.169 12.600 11.964
b (Å) 6.927 6.955 7.217 6.800
c (Å) 12.671 12.854 13.008 13.034
β (deg) 90.270 91.435 91.206 90.269
K (GPa) 41.4 37.3 59.1
ΔGdiss (kJ/mol) 20.07/19.11 19.25 53.05
insoluble anhydrite volume (Å3) 305.481 306.810 327.522 257.603
Cmcm a (Å) 6.993 6.908 7.056 6.470
b (Å) 6.245b 6.367 6.460 6.115
c (Å) 6.995b 6.976 7.186 6.511
K (GPa) 54.9 63.8 57.8 162.4
ΔGdiss (kJ/mol) 23.68/26.157 39.82 89.09
soluble anhydrite 1 volume (Å3) 530.856 535.398 570.146 549.439
C222 a (Å) 12.078 11.992 12.331 12.125
b (Å) 6.972 6.924 7.119 7.023
c (Å) 6.304 6.448 6.495 6.453
K (GPa) 59.1 54.3 71.9
ΔGdiss (kJ/mol) 15.31 31.33 80.45
soluble anhydrite 2 volume (Å3) 265.149 267.699 285.073 225.506
P6222 a(Å) 6.969 6.924 7.119 6.470
c (Å) 6.300 6.448 6.495 6.222
K (GPa) 59.1 54.3 150.0
ΔGdiss (kJ/mol) 10.87 31.19 80.45
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.

solution to a crystalline environment. This problem is field, FF1, overall, and slightly superior in terms of the cell
particularly acute for hydrate salts where the water model is volume and angles (where not right angles). However, for the
optimized primarily for the liquid phase and often fails to anhydrous phases of calcium sulfate, there are substantial errors
capture the different properties of water when embedded in an in the cell parameters, in some cases of up to 7.5%.
ionic solid. This also has an effect on the structural properties of Interestingly, the results for anhydrous calcium sulfate in the
the CaSO4 phases predicted by FF2 that show errors of up to soluble anhydrite 1 structure are quite accurate, while both of
7.5% in the lattice parameters and up to 15.75% in the the other polymorphs show large errors. This suggests that
predicted volume of the solid phases. In this regard, both the further refinement of the AMOEBA parameters would be
AMOEBA force field and the AIMD simulations should be necessary for use in the solid-state under dry conditions.
superior and better able to reproduce the solubilities of all the Specifically, it may be necessary to derive individual pairwise
structures simultaneously. Because of the much larger terms for the Lennard-Jones interactions rather than rely on
computational cost of using these methods, we are unable to combination rules. Beyond structure, the AMOEBA model also
obtain solvation free energies at present, thus preventing the seems to systematically overestimate the mechanical hardness
calculation of the solubilities of mineral phases. However, of phases based on the current parameters, as judged using the
recent work has demonstrated that the determination of free bulk moduli. At present, we cannot comment on the bulk
solubilities for AMOEBA for reasons discussed above.
energies of solvation for ions using AMOEBA49 and AIMD is
However, it should be noted that even first-principles methods
becoming possible.50
have been shown to struggle to yield accurate thermodynamics
When analyzing the performance of the AMOEBA model for
for hydrated mineral phases.51
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 4. ION PAIRING RESULTS
unsurprising that the reproduction of the structure of the Having developed two new candidate classical force fields for
hydrate phases is best, being comparable to the rigid ion force CaSO4 and thoroughly benchmarked them against the
F DOI: 10.1021/acs.jpcc.7b09820
J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C 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 AIMD, it is not possible to unambiguously determine the
turn our attention to the ion pairing free energy as a relative free energy of the Ca water coordination states for a
thermodynamic property that is not explicitly fitted during truly isolated Ca2+ ion.
the force field parametrization. Here the free energy landscape In order to better appreciate the differences between the
has been calculated as a function of the calcium−sulfur distance classical force field and the AIMD, it is instructive to look at the
and the calcium−water coordination number. As in our 1D pairing free energy obtained by integrating out the water
previous work on the ion pairing of alkaline earth metals coordination number (Figure 5). Because of the lower
with carbonate,12 the analysis of the 2D free energy maps computational cost of the classical and polarizable MD
(Figure 4) provides a clear picture of the ion pair formation. simulations, it is feasible to use large enough cells to converge
Although the finer details of the free energy profiles are the long-range tail of the free energy and align them onto the
different, all the methods used here clearly show that the Ca2+ analytic solution for the pairing free energy of two point
ion in solution has a preferred water coordination number of 7 charges interacting via a screened electrostatic potential
and that the formation of a contact ion pair (CIP) involves a (dashed line). However, for the AIMD simulations, we could
progressive desolvation of the cation to be surrounded by six only explore the contact ion pair (CIP) and solvent-shared ion
and five water molecules to form a monodentate and bidentate pair (SSHIP) states, and therefore, we aligned the free energy
ion pair, respectively. The results reported here are qualitatively based on the position of the SSHIP minimum. All methods
consistent with previous DFT calculations of the Ca−Cl ion used in this work clearly indicate that the SSHIP state has a
pair formation by Baer and Mundy, which also considered the lower free energy than the CIP and has similar barriers for the
hydration of the Ca2+ ion alone.52 In both the BLYP-D3 results formation of the CIP and for the dissociation of the complex.
of this work and the force fields of the present study, the However, there are significant differences between the
preferred Ca by water coordination number is also 7. However, predicted stability of the CIP, with FF1 and the AMOEBA
Baer and Mundy52 found a Ca2+ coordination number of 6 to calculation indicating that the CIP is only marginally higher in
be more stable than 8 as the secondary favored hydration state, energy than the SSHIP (∼2 kJ/mol), while BLYP-D3 shows a
while our force field calculations suggest the reverse (Figure larger difference (∼7 kJ/mol), and FF2 only has a shallow
SI4). Due to the limited range of Ca−S distance sampled in our minimum for the CIP, which is approximately 15 kJ/mol less
G DOI: 10.1021/acs.jpcc.7b09820
J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C Article

Figure 6. Temperature dependence of the ion pairing free energy

compared to the experimental results of Ainsworth,54 Nordstrom et
al.,57 Bešter-Rogač,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
Figure 5. One-dimensional Ca−SO4 ion pairing free energy as a
Entropic Contributionsa
function of the Ca−S distance calculated using the classical force fields ΔG (kJ/mol) ΔH (kJ/mol) ΔS (J/mol/K)
developed in this work (FF1 and FF2), the AMOEBA force field, and
FF1 −12.55 11.3 ± 0.7 80 ± 2
ab initio MD simulations.
FF2 −12.78 5.7 ± 0.8 62 ± 2
AMOEBA −12.21 7.1 ± 1.8 64 ± 5
favorable than the SSHIP. Another apparent qualitative Ainsworth54 −13.18 9.5 ± 0.9 76 ± 3
difference between the classical and AMOEBA force fields Nordstrom et al.57 −13.2 6.9 67.0
(with the exception of FF1) compared to the AIMD results is Wachter and Riederer55 −13.0 6.7 66.0
the seeming absence of a stable bidentate ion pair minimum in Bešter-Rogač56 −12.9 5.4 61.5
the 1D free energy profiles. However, careful examination of a
The results calculated using the force fields are compared with
the 2D free energy landscape indicates that sometimes this can experimental results by Ainsworth.54 The two ion pairing free energies
be deceptive. For example, for AMOEBA there is a minimum of FF1 and FF2 correspond to the values obtained from the line of
corresponding to 5-fold water coordination by calcium as best fit. The entropic term for the CaSO4 association by Nordstrom et
expected for the bidentate case, but because the monodentate al.57 was obtained from the free energy and enthalpy of formation.
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 therefore, we cannot directly calculate the association constant
1D profile, as would be readily apparent if a crossover in for this method. However, by assuming that it connects with
stability as a function of distance was to occur. The bottom line the force field free energy curves at the SSHIP, it is possible to
remains that the bidentate CIP is unlikely to be of particular estimate that the association constant predicted by AIMD is
importance. likely to fall between the values obtained from FF1 and FF2.
The pairing free energy calculations (except BLYP-D3) were Both the classical and AMOEBA force fields used in this
carried out at different temperatures in order to extract the work predict a positive association enthalpy (5.7/11.3 kJ/mol)
enthalpic and entropic contributions. The classical force field and a pairing free energy at standard conditions in the range
simulations have been run at intervals of 5 K between 290 and −12.2/−12.8 kJ/mol, which is perfectly consistent with the
350 K (Figure SI5), and the AMOEBA calculations at four available experimental data (Table 5).54−57 These results are
different temperatures between 300 and 345 K. The ion-pairing also consistent with the similar case of Ca2+ and CO32−
dissociation/association constant and association free energy association, which we analyzed in detail in ref 58, where the
have then been calculated by integrating the potential of mean ion pair formation was shown to be driven by the entropy of
force;12,53 the solvent molecules that are released during ion pair
R1 formation, rather than by the electrostatic interaction between
K a −1 = K d = C 0 ∫R 0
exp(−βϕ(r ))4πr 2 dr
(7)
the ions. From the analysis of the metadynamics trajectories
obtained using FF1, it is possible to extract the average water
where C0 is a constant used to convert to standard coordination number for the ion pair (17 water molecules), and
concentration and ϕ(r) = ΔG + kBT ln(4πr2) is the potential by comparing it with that of the isolated species (7 for Ca and
of mean force. The resulting temperature-dependent free approximately 13 for SO4), we find that on average three water
energies are plotted in Figure 6, while the enthalpy and molecules are released into solution (Figure SI6). Analogously,
entropy contributions obtained by fitting these curves are given for the AMOEBA and FF2 force fields, four water molecules are
in Table 5. released back into solution during the ion pair formation.
As noted above, because of the small simulation cell used in Hence, the entropic contribution to ion pair formation per
the AIMD simulations, it is impossible to align the pairing free water molecule released is in the range 16−26 J/mol/K, which
energy to the analytic solution in the asymptotic region, and is also consistent with the experimental and computational
H DOI: 10.1021/acs.jpcc.7b09820
J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C Article

values for the CaCO3 system (29.3 J/mol/K)58 and other formation of the contact ion pair might be a prerequisite to
experimental estimates of the water entropy change due its subsequent formation of any pre- or postnucleation species.
localization in crystalline hydrated phases (32/33.5 J/mol/ Future work is required to probe the further ion association
K).59,60 This phenomenon of entropy-driven association, which processes of calcium sulfate, and this is now possible based on
has been experimentally observed in a variety of systems61 and the force field (FF1) developed within this study.

■
extensively studied with computational methods,62,63 may have
important consequences for the thermodynamic stability and ASSOCIATED CONTENT
abundance of prenucleation clusters, beyond that expected * Supporting Information
S
from classical nucleation theory. Among the many computa-
The Supporting Information is available free of charge on the
tional works in the literature that deal with the formation of ion
ACS Publications website at DOI: 10.1021/acs.jpcc.7b09820.
pairs between ionic species, a few recent ones have started
looking at the contribution of water and entropy to the pairing Comparison between radial pair distribution functions
free energy.62,63 In particular, the works by Ballard and computed with the potential described in this work and
Dellago62 and by Shi and Beck,63 which both focused on alkali literature ones; 3D water density map computed using a
halide ion pairs, conclude that the ion association is driven by force field from the literature; self-diffusion coefficient as
an increase in entropy due to the release of water molecules a function of the box size; molecular arrangement in the
from the ions’ coordination shell into the solution. This picture various states of the ion pair formation; all the pairing
is perfectly consistent with the results presented here for CaSO4 free energy curves for the two force fields discussed in
and with our previous work on the formation of CaCO3 ion this work; table with the dissolution free energy and
pairs, where ion pair formation is driven by an increase in solubilities of the calcium sulfate phases (PDF)

■
entropy proportional to the number of water molecules
released into solution and opposed by a positive enthalpic AUTHOR INFORMATION
term. The unfavorable enthalpy change can be ascribed to the
loss of electrostatic interactions with the released water that is Corresponding Author
not compensated by the cation−anion electrostatic interaction. *E-mail: p.raiteri@curtin.edu.au. Phone: + 61 8 9266 2687.
ORCID
5. CONCLUSIONS Paolo Raiteri: 0000-0003-0692-0505
In conclusion, we have developed and tested two rigid ion force Author Contributions
fields by targeting the two available experimental estimates of All authors contributed to performing the calculations,
the hydration free energy of the sulfate anion and the solubility analyzing the data, and writing the manuscript. All authors
of the most common hydrated solid phases of CaSO4. Because have given approval to the final version of the manuscript.
the two available hydration free energies differ by about 100 kJ/
mol, the two force fields gave markedly different predictions for Notes
The authors declare no competing financial interest.

■
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 ACKNOWLEDGMENTS
sophisticated, polarizable force field using the AMOEBA model E.H.B. would like to acknowledge the contribution of an
and ab initio molecular dynamics with a dispersion-corrected Australian Government Research Training Program Scholar-
GGA functional. Through comparison of the two rigid ion force ship in supporting this research, and P.R. and J.D.G. thank the
fields against the results of the other methods, as well as other Australian Research Council for financial support through the
literature data, it is possible to identify FF1, which gives the Discovery Programme (FT130100463, DP160100677). This
more exothermic hydration free energy for sulfate, as being the work was also supported by resources provided by the Pawsey
superior representation of aqueous calcium sulfate systems. In Supercomputing Centre and National Computational Infra-
particular, the distribution of water around the sulfate ion is structure with funding from the Australian Government and the
markedly different between FF2 and all of the other results. Government of Western Australia.

■
Having identified the best rigid ion force field model for the
system, it was then possible to examine the ion pairing of REFERENCES
calcium sulfate in aqueous solution. Based on this force field, as (1) Azimi, G.; Papangelakis, V. G.; Dutrizac, J. E. Modelling of
well as AMOEBA and BLYP-D3, there is a consistent picture of Calcium Sulphate Solubility in Concentrated Multi-Component
the ion pairing process. Consistent with previous estimates of Sulphate Solutions. Fluid Phase Equilib. 2007, 260 (2), 300−315.
the ion pairing free energy, we find exothermic binding under (2) Van Driessche, A. E. S.; Benning, L. G.; Rodriguez-Blanco, J. D.;
standard conditions of approximately −13 kJ/mol. By Ossorio, M.; Bots, P.; García-Ruiz, J. M. The Role and Implications of
decomposing the free energy into the enthalpy and entropy, Bassanite as a Stable Precursor Phase to Gypsum Precipitation. Science
based on the temperature dependence, it is found that ion 2012, 336 (6077), 69−72.
pairing is entropically driven, which is consistent with the (3) Pina, C. M. Nanoscale Dissolution and Growth on Anhydrite
release of an average of three water molecules observed during Cleavage Faces. Geochim. Cosmochim. Acta 2009, 73 (23), 7034−7044.
binding. The contact ion pair is found to be between 2 and 8 (4) Demichelis, R.; Raiteri, P.; Gale, J. D.; Quigley, D.; Gebauer, D.
Stable Prenucleation Mineral Clusters Are Liquid-Like Ionic Polymers.
kJ/mol less stable than the solvent-shared state, with a barrier
Nat. Commun. 2011, 2, 590.
of at least 14 kJ/mol to reach a distance at which direct (5) Gebauer, D.; Voelkel, A.; Coelfen, H. Stable Prenucleation
coordination between the ions occurs. Therefore, it appears Calcium Carbonate Clusters. Science 2008, 322 (5909), 1819−1822.
unlikely that the contact state plays a substantial role in the (6) Gebauer, D.; Kellermeier, M.; Gale, J. D.; Bergström, L.; Cölfen,
initial ion pairing. However, the balance between these two H. Pre-Nucleation Clusters as Solute Precursors in Crystallisation.
states may change as further ion association occurs, and so Chem. Soc. Rev. 2014, 43 (7), 2348.

I DOI: 10.1021/acs.jpcc.7b09820
J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C Article

(7) Navrotsky, A. Energetic Clues to Pathways to Biomineralization: (29) Allan, N. L.; Rohl, A.; Gay, D.; Catlow, C. R. A.; Davey, R.;
Precursors, Clusters, and Nanoparticles. Proc. Natl. Acad. Sci. U. S. A. Mackrodt, W. Calculater Bulk and Surface Properties of Sulfates.
2004, 101 (33), 12096−12101. Faraday Discuss. 1993, 95, 273−280.
(8) Raiteri, P.; Gale, J. D. Water Is the Key to Nonclassical (30) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.;
Nucleation of Amorphous Calcium Carbonate. J. Am. Chem. Soc. 2010, Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density
132 (49), 17623−17634. Functional Calculations Using a Mixed Gaussian and Plane Waves
(9) De Yoreo, J. J.; Vekilov, P. G. Principles of Crystal Nucleation Approach. Comput. Phys. Commun. 2005, 167 (2), 103−128.
and Growth. Rev. Mineral. Geochem. 2003, 54, 57−93. (31) Krack, M. Pseudopotentials for H to Kr Optimized for Gradient-
(10) Erdemir, D.; Lee, A. Y.; Myerson, A. S. Nucleation of Crystals Corrected Exchange-Correlation Functionals. Theor. Chem. Acc. 2005,
From Solution: Classical and Two-Step Models. Acc. Chem. Res. 2009, 114 (1−3), 145−152.
42 (5), 621−629. (32) Becke, A. Density-Functional Exchange-Energy Approximation
(11) Raiteri, P.; Gale, J. D.; Quigley, D.; Rodger, P. M. Derivation of with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys.
an Accurate Force-Field for Simulating the Growth of Calcium 1988, 38 (6), 3098−3100.
Carbonate From Aqueous Solution: a New Model for the Calcite− (33) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and
Water Interface. J. Phys. Chem. C 2010, 114 (13), 5997−6010. Accurate Ab Initio Parametrization of Density Functional Dispersion
(12) Raiteri, P.; Demichelis, R.; Gale, J. D. Thermodynamically Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010,
Consistent Force Field for Molecular Dynamics Simulations of 132 (15), 154104−154120.
Alkaline-Earth Carbonates and Their Aqueous Speciation. J. Phys. (34) Bankura, A.; Karmakar, A.; Carnevale, V.; Chandra, A.; Klein, M.
Chem. C 2015, 119 (43), 24447−24458. L. Structure, Dynamics, and Spectral Diffusion of Water From First-
(13) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Principles Molecular Dynamics. J. Phys. Chem. C 2014, 118 (50),
Natl. Acad. Sci. U. S. A. 2002, 99 (20), 12562−12566. 29401−29411.
(14) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular (35) VandeVondele, J.; Hutter, J. An Efficient Orbital Trans-
Dynamics. J. Comput. Phys. 1995, 117 (1), 1−19. formation Method for Electronic Structure Calculations. J. Chem. Phys.
(15) Martyna, G.; Tobias, D.; Klein, M. L. Constant Pressure 2003, 118 (10), 4365−4369.
Molecular Dynamics Algorithms. J. Chem. Phys. 1994, 101, 4177. (36) Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.;
(16) Wu, Y.; Tepper, H.; Voth, G. A. Flexible Simple Point-Charge Brown, S. T.; Gilbert, A. T. B.; Slipchenko, L. V.; Levchenko, S. V.;
Water Model with Improved Liquid-State Properties. J. Chem. Phys. O’Neill, D. P.; et al. Advances in Methods and Algorithms in a Modern
2006, 124 (2), 024503. Quantum Chemistry Program Package. Phys. Chem. Chem. Phys. 2006,
(17) Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. J. Chem. 8 (27), 3172−3191.
Phys. 1935, 3 (5), 300. (37) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Perspective on
(18) Zwanzig, R. High-Temperature Equation of State by a Foundations of Solvation Modeling: the Electrostatic Contribution to
Perturbation Method. I. Nonpolar Gases. J. Chem. Phys. 1954, 22 the Free Energy of Solvation. J. Chem. Theory Comput. 2008, 4 (6),
(8), 1420. 877−887.
(19) Hummer, G.; Pratt, L. R.; Garcia, A. Ion Sizes and Finite-Size (38) Wu, Y.; Chen, H.; Wang, F.; Paesani, F.; Voth, G. A. An
Corrections for Ionic-Solvation Free Energies. J. Chem. Phys. 1997, 107 Improved Multistate Empirical Valence Bond Model for Aqueous
(21), 9275−9277. Proton Solvation and Transport. J. Phys. Chem. B 2008, 112 (2), 467−
(20) Kastenholz, M. A.; Hünenberger, P. H. Computation of 482.
(39) Marcus, Y. Ion Solvation; Wiley, 1985.
Methodology-Independent Ionic Solvation Free Energies From
(40) Marcus, Y. A Simple Empirical Model Describing the
Molecular Simulations. II. the Hydration Free Energy of the Sodium
Thermodynamics of Hydration of Ions of Widely Varying Charges,
Cation. J. Chem. Phys. 2006, 124 (22), 224501.
Sizes, and Shapes. Biophys. Chem. 1994, 51 (2), 111−127.
(21) Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered
(41) Marcus, Y. Ions in Solution and Their Solvation; John Wiley &
Metadynamics: a Smoothly Converging and Tunable Free-Energy
Sons, Inc., 2016.
Method. Phys. Rev. Lett. 2008, 100 (2), 20603. (42) Cramer, C. J.; Truhlar, D. G. A Universal Approach to Solvation
(22) Raiteri, P.; Laio, A.; Gervasio, F. L.; Micheletti, C.; Parrinello, M. Modeling. Acc. Chem. Res. 2008, 41 (6), 760−768.
Efficient Reconstruction of Complex Free Energy Landscapes by (43) Ohtaki, H.; Radnai, T. Structure and Dynamics of Hydrated
Multiple Walkers Metadynamics. J. Phys. Chem. B 2006, 110 (8), Ions. Chem. Rev. 1993, 93 (3), 1157−1204.
3533−3539. (44) Cannon, W.; Pettitt, B.; McCammon, J. A. Sulfate Anion in
(23) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Water - Model Structural, Thermodynamic, and Dynamic Properties. J.
Bussi, G. Computer Physics Communications. Comput. Phys. Commun. Phys. Chem. 1994, 98 (24), 6225−6230.
2014, 185 (2), 604−613. (45) Pegado, L.; Marsalek, O.; Jungwirth, P.; Wernersson, E.
(24) Ponder, J. W.; Wu, C.; Ren, P.; Pande, V. S.; Chodera, J. D.; Solvation and Ion-Pairing Properties of the Aqueous Sulfate Anion:
Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, Explicit Versus Effective Electronic Polarization. Phys. Chem. Chem.
R. A., Jr.; et al. Current Status of the AMOEBA Polarizable Force Phys. 2012, 14 (29), 10248.
Field. J. Phys. Chem. B 2010, 114 (8), 2549−2564. (46) Williams, C. D.; Carbone, P. A Classical Force Field for
(25) Gale, J. D. GULP: Capabilities and Prospects. Z. Kristallogr. - Tetrahedral Oxyanions Developed Using Hydration Properties: the
Cryst. Mater. 2005, 220 (5−6), 552−554. Examples of Pertechnetate (TcO4−) And Sulfate (SO42−). J. Chem.
(26) Ponder, J. W.; Richards, F. M. An Efficient Newton-Like Phys. 2015, 143 (17), 174502−174509.
Method for Molecular Mechanics Energy Minimization of Large (47) Vchirawongkwin, V.; Rode, B. M.; Persson, I. Structure and
Molecules. J. Comput. Chem. 1987, 8, 1016. Dynamics of Sulfate Ion in Aqueous Solution-an Ab Initio QMCF MD
(27) Piquemal, J.-P.; Perera, L.; Cisneros, G. A. S.; Ren, P.; Pedersen, Simulation and Large Angle X-Ray Scattering Study. J. Phys. Chem. B
L. G.; Darden, T. A. Towards Accurate Solvation Dynamics of 2007, 111 (16), 4150−4155.
Divalent Cations in Water Using the Polarizable Amoeba Force Field: (48) Yuan-Hui, L.; Gregory, S. Diffusion of Ions in Sea Water and in
From Energetics to Structure. J. Chem. Phys. 2006, 125 (5), 054511− Deep-Sea Sediments. Geochim. Cosmochim. Acta 1974, 38 (5), 703−
054517. 714.
(28) Lambrecht, D. S.; Clark, G. N. I.; Head-Gordon, T.; Head- (49) Mohamed, N. A.; Bradshaw, R. T.; Essex, J. W. Evaluation of
Gordon, M. Exploring the Rich Energy Landscape of Sulfate. J. Phys. Solvation Free Energies for Small Molecules with the AMOEBA
Chem. A 2011, 115, 11438−11454. Polarizable Force Field. J. Comput. Chem. 2016, 37 (32), 2749−2758.

J DOI: 10.1021/acs.jpcc.7b09820
J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C Article

(50) Bell, D. R.; Qi, R.; Jing, Z.; Xiang, J. Y.; Mejias, C.; Schnieders, (70) Schwerdtner, W. M.; Tou, J.; Hertz, P. B. Elastic Properties of
M. J.; Ponder, J. W.; Ren, P. Calculating Binding Free Energies of Single Crystals of Anhydrite. Can. J. Earth Sci. 1965, 2 (6), 673−683.
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 Thermody-
namics 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 Electrolyte-
Solutions 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. Low-
Temperature Thermodynamic Model for the System Na2Co3-
MgCO3-CaCO3-H2O. Geochim. Cosmochim. Acta 1999, 63 (19−20),
3105−3119.
(61) De Visscher, A.; Vanderdeelen, J.; Kö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. High-
Pressure Behavior of Gypsum: a Single-Crystal X-Ray Study. Am.
Mineral. 2008, 93 (10), 1530−1537.

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