Geometrically frustrated interactions drive structural complexity in amorphous calcium carbonate

Amorphous calcium carbonate is an important precursor for biomineralization in marine organisms. Key outstanding problems include understanding the structure of amorphous calcium carbonate and rationalizing its metastability as an amorphous phase. Here we report high-quality atomistic models of amorphous calcium carbonate generated using state-of-the-art interatomic potentials to help guide fits to X-ray total scattering data. Exploiting a recently developed inversion approach, we extract from these models the effective Ca⋯Ca interaction potential governing the structure. This potential contains minima at two competing distances, corresponding to the two different ways that carbonate ions bridge Ca2+-ion pairs. We reveal an unexpected mapping to the Lennard-Jones–Gauss model normally studied in the context of computational soft matter. The empirical model parameters for amorphous calcium carbonate take values known to promote structural complexity. We thus show that both the complex structure and its resilience to crystallization are actually encoded in the geometrically frustrated effective interactions between Ca2+ ions.


Introduction
Calcium carbonate is relatively unusual amongst simple inorganic salts in that it precipitates from aqueous solution in a metastable hydrated, amorphous form. 1 ACC-with a nominal composition of CaCO ternatively directed to crystallise into a number of different polymorphs by varying pH or temperature. 4ture exploits this complex phase behaviour in a variety of biomineralisation processes to control the development of shells and other skeletal structures. 5,6 ot only does the amorphous nature of biogenic ACC allow transformation to different crystalline CaCO 3 polymorphs, but it helps organisms fashion larger-scale hierarchical morphologies that are important in biomineral architectures. 7In seeking to develop bio-inspired crystal-engineering approaches for synthetic control over phase and morphology selection, there is an obvious need to understand why such a chemically simple system can exhibit such complex phase behaviour.
One domain in which a similar kind of phase complexity has been studied deeply from a theoretical perspective is that of soft-matter systems governed by multi-well pair potentials.Whereas isotropic particles that interact via a simple single-well potential (e.g.Lennard-Jones, LJ) self-assemble into structurally simple crystalline phases (e.g.][10][11] An elegant example in two dimensions is that of quasicrystal self-assembly for specific parameters of the double-well Lennard-Jones-Gauss (LJG) potential. 8,12 n three dimensions the same interaction model can be tuned to stabilise complex crystals with enormous unit cells, 11 and some combinations of the LJ and Gauss well positions even appear to frustrate crystallisation altogether. 9at competing length-scales might be relevant to ACC is a point raised by the observation of two preferred Ca• • • Ca distances dominating medium-range order in synthetic ACC. 13,14 hese two distances are directly evident in the experimental X-ray pair distribution function (PDF) of ACC and are attributed to different bridging modes of the carbonate ion, which can connect a pair of Ca 2+ ions either directly through one of the oxygen atoms (Ca-O-Ca pathway) or by inserting the carbonate ion fully between the cations (Ca-O-C-O-Ca pathway).
Here we explore the possibility that the structure of ACC is governed by effective interactions between Ca 2+ ions that also reflect these two length-scales, and that its well-known complexity emerges because the distances involved are in competition with one another.Our approach begins by obtaining a high quality measure of the Ca-pair correlation function in ACC.We do this by applying a hybrid Reverse Monte Carlo (HRMC) approach 15 to generate the first structural model of ACC that is simultaneously consistent with experiment and stable with respect to state-of-the-art potentials.The Ca-pair correlation function that emerges is then inverted using a recently-developed algorithm 16 to reveal the effective, carbonate/water-mediated Ca• • • Ca interaction potential.We find that this potential is closely related to the LJG formalism, with empirical parameters that are known to frustrate crystallisation.Monte Carlo (MC) simulations driven by our LJG model yield coarse-grained representations of ACC that capture key aspects of our fully-atomistic HRMC models.In this way, we explain the structural complexity of ACC in terms of geometric frustration of two competing energetically-favoured Ca• • • Ca separations.
Mapping the problem of ACC structure onto the phase behaviour of multi-well potentials is important not only because it establishes the first experimental system for which these potentials are relevant, but because it suggests how structural complexity in inorganic phases might be controllably targeted through suitable tuning of effective interactions.

Results
Structure of ACC.Our HRMC refinements made use of atomistic configurations containing 12 960 atoms (1620 CaCO 3 • H 2 O formula units), with simulation cell sizes of ∼5 nm.During refinement, atomic moves were proposed and then rejected or accepted according to a Metropolis Monte Carlo criterion, with a cost function that measured the quality of fit to X-ray total scattering data 17 and also the energetic stability evaluated using the state-of-the-art interatomic potential of Ref. 18.We favoured an HRMC approach over empirical potential structure refinement (EPSR) because calcium carbonate potentials are notoriously delicately balanced, 18 such that EPSR-derived potentials are unlikely to be physically meaningful for this system.To allow for direct comparison with previous studies, we used a Q-weighted variant of the experimental X-ray total scattering function (i.e.QF X (Q)) of synthetic ACC reported in Ref. 13, which has been shown to be indistinguishable from that obtained for biogenic ACC samples. 19For consistency, our HRMC configuration size and refinement constraints were also identical both goodness-of-fit to experimental data χ 2 and cohesive energy (E rel ), and therefore it finds a structure solution that is essentially as consistent with experiment as that obtained using RMC, whilst also as energetically sensible as that obtained using potentials alone (MD).c Representation of a converged structure of ACC obtained using HRMC.Water molecules connect to form filamentary strands (aquamarine strings) that separate calcium-carbonate-rich domains (Ca atoms as large beige spheres, carbonate ions shown as stick representation).
to those used in the reverse Monte Carlo (RMC) investigation of synthetic ACC in Ref. 14.The absence of key energetic considerations in that RMC study allowed the development of physically unreasonable charge separation to give a model of cationic Ca 2+ -rich domains separated by channels of carbonate anions and water molecules.Our expectation was that the explicit consideration of electrostatics in our HRMC implementation would guide refinement towards a solution that was equally consistent with experiment but also energetically sensible.
HRMC does indeed find a suitable compromise between experiment and theory.We show in give very similar high-quality fits to experiment-unsurprising, of course, since they have been refined against these data-whilst the MD simulation misses some aspects of the QF X (Q) function.By contrast, the HRMC and MD models give similar cohesive energies, whereas the RMC model is less stable than both by more than 800 kJ mol −1 per formula unit.These results are summarised in Fig. 1 The HRMC model itself is illustrated in Fig. 1(c).We observe that water is not homogeneously distributed throughout the configuration, but rather that the ACC structure consists of CaCO 3 -rich regions separated by a filamentary network of water: we term this a 'blue cheese' model.Qualitatively similar descriptions were obtained in MD simulations, 23 in EPSR refinements of combined neutron/Xray total scattering data, 22 and also inferred from solid-state NMR measurements. 24In all cases, the interpretation is that nearly all H 2 O molecules are bound to Ca 2+ (as implied by 1 H NMR; Ref. 25), but neighbouring water molecules are sufficiently close to form a network that percolates the ACC structure.
As anticipated, the charge separation that developed in the RMC model of Ref. 14 has now vanished: the water-rich filaments we observe do not contain any free carbonate and are significantly narrower than the nanopore channels reported previously.The fact that RMC and HRMC models give similar QF X (Q) fits despite their very different descriptions of the structure reflects the difficulty of discriminating between C and O atoms using X-ray scattering methods alone. 26cally, our HRMC model shows coordination environments that are consistent with the consensus of recent computational and experimental studies.For completeness, we show in Fig. 2  coordination geometries for Ca 2+ and CO 2− 3 are shown in Fig. 2(c) and (d).As anticipated, 14 we find that carbonate anions bridge calcium-ion pairs in two ways: either a pair of Ca 2+ ions share a common carbonate oxygen neighbour, or they bind distinct oxygens and so are connected formally by a Ca-O-C-O-Ca pathway.A full analysis of coordination environments, including bond-length and bond-angle distributions, is provided in the Supplementary Material.
Coarse graining.In seeking to improve our understanding of the structure of ACC, we focussed on the Ca-pair correlation function g Ca (r)-after all, this is the contribution to the pair distribution function that exhibits the strongest persistent well-defined oscillations.In Fig. 2(e), we show this function as extracted from our newly-generated HRMC configurations compared against that reported in the RMC study of Ref. 14.We also include the Fourier transform of the experimental F X (Q) function; this transform is an approximate (total) pair distribution function that emphasises contributions from Ca-Ca pairs as a consequence of the larger X-ray scattering cross-section for Ca relative to C, O, and H.While all three real-space correlation functions show maxima at r 4 and 6 Å, it is the HRMC result that resolves these features most clearly.By counting Ca-Ca pairs separated by Ca-O-Ca and Ca-O-C-O-Ca pathways in our HRMC configuration, we confirmed that the two maxima in g Ca (r) occur at the preferred Ca• • • Ca distances associated with these two different carbonate bridging motifs (see Supplementary Material).
Access to a smoothly-varying measure of g Ca (r), together with the HRMC configurations from which it is calculated, allowed us to exploit a recently-developed approach for the direct measurement of effective pair potentials from particle-coordinate data. 16The method works by equating the pair distribution functions measured directly, on the one hand, and calculated using a test-particle insertion approach, 29,30 on the other hand. 31Applying this methodology to the atomic coordinates in our HRMCrefined model, we extracted the effective Ca 2+ -ion pair potential u Ca (r) shown in Fig. 3(a).This function has two distinct potential wells, the minima of which are centred near distances for which g Ca (r) has its maxima.The physical meaning of u Ca (r) is that it captures the effective two-body interactions between Ca-Ca pairs, mediated by carbonate and water, as required to account for the observed g Ca (r).That this interaction energy is minimised when Ca-Ca pairs are separated by distances corresponding to either Ca-O-Ca or Ca-O-C-O-Ca pathways makes intuitive sense (of course); the longer (6 Å) separation appears energetically more favourable likely because it minimises the electrostatic repulsion between calcium ions bound to a common carbonate.
One assumption in using an isotropic effective pair potential to capture interactions mediated by anisotropic molecular species, such as H 2 O and CO 2− 3 , is that any local anisotropy is sufficiently shortranged.We checked this point by calculating from our original (all-atom) HRMC configurations the orientational correlation functions where P 2 (x) is the second-order Legendre polynomial, and the S are vectors parallel to a suitable local axis (e.g.C 3 , C 2 ) of molecules separated by distance r. 32 We find that φ(r) is essentially featureless-for both carbonate and water-for all but the very shortest intermolecular approaches, which occur for distances smaller than the nearest-neighbour Ca• • • Ca separation [inset to Fig. 3(a)].Hence the anisotropy of individual molecules is relevant only at a length-scale smaller than that that over which the effective interactions between Ca 2+ -ion pairs operate.
The double-well form of u Ca (r) is important for a number of reasons.First, it explains why sim-ple hard-sphere or LJ potentials fail to reproduce the experimental g Ca (r): if tuned to capture the 4 Å nearest-neighbour correlation, such potentials always predict a second maximum just below 8 Å that is not observed in experiment.Second, it allows us to rationalise the inability of earlier EPSR studies to determine the correct Ca distribution function; for example, the study of Ref. 20 included a LJ parameterisation of Ca• • • Ca interactions that then resulted in a g Ca (r) function with maxima at positions clearly inconsistent with the inverse Fourier transform of the experimental F X (Q).And, third, multi-well potentials are well known to support very complex phase behaviour, 9,11,12,[33][34][35] suggesting a qualitative explanation of the complexity of ACC.
LJG parameterisation.To test this hypothesis of a link between effective potential and structural complexity, we sought to map u Ca (r) onto a suitable multi-well potential for which the corresponding theory is already well established.We found that a double-well LJG interaction 8 described the observed functional form surprisingly well, so long as we included an additional broad, repulsive, Gaussian term that helps capture the local maxima between and beyond the two minima of the LJG function.The behaviour of the LJG potential is characterised by three parameters-, r 0 , σ-which describe, respectively, the depth, position, and width of the Gaussian well relative to the LJ component. 9,11  least-squares fit to u Ca (r) gave the values = 4.1, r 0 = 1.4 and σ = 0.14 [Fig.3(a)]; note that r 0 is particularly well defined because it is closely related to the ratio of the Ca• • • Ca separations resulting from the two different carbonate bridging motifs ( 6 Å / 4 Å).We will return to the importance of these empirical parameters in due course.We note in passing related work on the description of coarse-grained molecular systems with double-well potentials, where one minimum is explicitly assigned to enthalpic and one to entropic terms, 36 and on the 'learning' of effective pair potentials from simulation data by back-propagation. 37 a simple check that our combination of g(r) inversion and potential parameter fitting does indeed result in a meaningful effective pair-potential, we carried out direct MC simulations driven by our parameterised LJG potential.These simulations were performed using the experimentally-determined Ca

Discussion
From the perspective of soft-matter theory, much of the interest in the LJG potential lies in the complexity of phase behaviour which it supports. 9,11,12,38 Ths complexity arises because of competition between the structure-directing effects of LJ and Gaussian components, which operate at different length-scales.
When the positions of the LJ and Gaussian minima are related by a factor 1.2 ≤ r 0 ≤ 1.6, there is particularly strong geometric frustration that results in a general resistance to crystallisation 9 and the emergence of many competing low-symmetry ground-state structures. 11This behaviour extends across a wide variety of relative well depths ≥ 1 for the corresponding relative Gaussian width σ = 0.14. 11 Hence the empirical parameters we have determined to be relevant to effective Ca• • • Ca interactions in ACC-viz.r 0 = 1.4,= 4.1, and σ = 0.14-locate the system in one of the most complicated parts of the LJG phase diagram [Fig.3(c)]. 11In developing this mapping onto the LJG model, we assume that the additional broad Gaussian term we have included in our parameterisation does not strongly influence phase behaviour (note that the value of r 0 does not vary significantly if it is omitted from our fit).Likewise, we cannot be certain of the effects of finite temperature and fixed density on LJG phase behaviour, since these have not yet been studied from theory.It is nevertheless a general phenomenon in frustrated systems that regions of strong geometric frustration are characterised by the existence of many competing ground states with suppressed ordering temperatures, above which the system is disordered but far from random. 39,40 the context of ACC, the implications are twofold.First, the different favoured Ca• • • Ca separations for carbonate-bridged Ca pairs-viz.4 and 6 Å-are configurationally difficult to satisfy in a three-dimensionally periodic (i.e., crystalline) structure with the density of ACC.This rationalises, in general terms, why an amorphous form of CaCO 3 is so relatively (meta)stable.And, second, the sensitivity of the LJG potential around r 0 = 1.4 suggests-again in general terms-a qualitative explanation of why ACC might be directed to crystallise into a number of different polymorphs.Of course, as the water content of ACC varies during the ageing of the material, both the effective Ca• • • Ca interactions and the bulk density must change, which makes extrapolation from the LJG effective potential of ACC difficult.Increased material density also strengthens the orientational correlations in molecular components, and so one expects the isotropic effective pair potential description to break down as crystallisation is approached.Nevertheless, at the density of ACC and in the absence of orientational order, the relative stability of a disordered state is now more easily rationalised.Likewise, the incorporation of Mg 2+ or PO 3− 4 ions into ACC will introduce statistical variations into the effective potential governing cation arrangements, stabilising further the amorphous state-much as magnetic exchange disorder can stabilise spin-glasses. 41ereas the (hitherto unsuccessful) search for 'real-world' materials governed by isotropic multiwell potentials has focussed on the atomic (alloys) and mesoscopic (colloids) length-scales, 11,12 our study of ACC suggests that the intermediate domain of nanoscale materials may in fact be a more fertile source of relevant examples.By their very nature, molecular ions tend to exhibit a range of coordination modes, which must give rise to multiple specific preferred distances between the counter-ions they coordinate.Whenever these distances are geometrically frustrated, as in ACC, one expects that the corresponding effective potential must contain multiple wells so as to stabilise both separations at once.Through suitable choice of chemical components, one might hope to navigate the complicated phase space associated with such unconventional potentials.Doing so will provide an important test of theory, on the one hand, and also allow the targeted design of complex structures through self-assembly, on the other hand.
Our coarse-grained approach to understanding ACC structure may be applicable to other poorly ordered inorganic solids of particular scientific importance.Amorphous calcium phosphate (the precursor to bone) and calcium-silicate-hydrates (i.e.Portland cements) are obvious ever-contentious examples, 42,43 where one might expect the different bridging motifs-now of the phosphate or (poly)silicate anions-to again moderate an effective interaction potential with multiple minima.Better understanding the structures of long-studied materials is certainly one fruitful avenue for future research, but our work suggests also a new pathway towards complex materials design.By varying the inorganic cation radius for a given molecular counter-ion, or by changing the degree of hydration, one might hope to control the ratio of the distances at which minima in the effective potential occur: what emerges is a strategy of 'interaction engineering' to target a particular phase of interest, or to stabilise or destabilise an amorphous form of matter.

Figure 1 :
Figure 1: Hybrid reverse Monte Carlo yields a balanced model of ACC structure.a The experimental Q-weighted Xray total scattering function QF X (Q) is well fitted by both RMC (top) and HRMC (middle) refinements, but is meaningfully different to that calculated using the interatomic potential of Ref. 18 (bottom).Experimental data are shown as black lines, and fits or calculations as red lines.b Quality criteria for the different structural models.HRMC simultaneously optimises

Fig. 1 (
Fig. 1(a) the QF X (Q) function calculated from a representative configuration, and compare it to the equivalent functions predicted using pure RMC refinement, on the one hand, and unconstrained moleculardynamics (MD) simulations with the potentials of Ref. 18, on the other hand.Both RMC and HRMC (b), which captures the motivation for our use of HRMC as a suitable balancing act: it has, for the first time, allowed access to an atomistic representation of ACC structure that is consistent with experiment and gives sensible energies using established hydrated calcium carbonate potentials.We checked also whether the HRMC model could reproduce the neutron total scattering data of Refs.20-22 (it can), and whether it is stable in molecular dynamics (MD) simulations driven by the potentials of Ref. 18 (it is); both points are discussed in greater detail in the Supplementary Material.

Figure 2 :
Figure 2: Coordination environments and Ca-pair distributions in ACC. a Histogram of Ca 2+ coordination environments, decomposed into contributions from carbonate and water oxygen donors (O c and O w , respectively).b Histogram of CO 2− 3 coordination environments, now decomposed into contributions from Ca 2+ and water hydrogen donors (O w ).c Representative Ca 2+ coordination sphere for the modal coordination environment marked by a star in a. d Representative CO 2− 3 coordination sphere for the modal coordination environment marked by a star in b.Note that pairs of Ca 2+ ions within the same CO 2− 3 coordination sphere either share a common oxygen donor (e.g.red arrow) or are connected by Ca-O-C-O-Ca pathways (e.g.blue arrow).Ca, C, O, and H atoms shown in beige, grey, aquamarine, and white, respectively.e Ca-pair correlation functions g Ca (r) extracted from HRMC, RMC, and LJG configurations, compared against the normalised Fourier transform of the experimental X-ray total scattering function (which includes contributions from all atom-pairs, e.g. the Ca-O peak at 2.4 Å marked with an asterisk).The two principal peaks common to all functions, indicated by red and blue shading, can be assigned to the two types of Ca 2+ -ion pairs illustrated in d.

Figure 3 :
Figure 3: Effective Ca• • • Ca interactions in ACC. a The effective Ca• • • Ca interatomic potential extracted from our HRMC configurations 16 (open circles) and least-squares fit using a modified LJG model (line) as described in the text.The inset shows the orientational correlation functions (Eq.(1)), which vanish for distances relevant to the Ca• • • Ca separations.b Representative MC configuration of Ca atoms (beige spheres) generated by the LJG potential, parameterised by the fit shown in a.The Ca atoms are not uniformly distributed, but cluster to leave Ca-poor voids, shown as aquamarine surfaces.Note the qualitative similarity to the heterogeneous structure of ACC represented in Fig. 1(c).c Approximate location of ACC effective LJG parameters in the LJG phase space as reported in Ref. 11.
particle density in ACC.The corresponding pair correlation function matches closely our HRMC g Ca (r), as expected [Fig.2(e)], and the resulting structure resembles the Ca distribution in the fully atomistic model [Fig.3(b); cf.Fig. 1(c)].But the similarity between coarse-grained LJG-driven and HRMC Ca distributions turns out to extend beyond pair correlations, and we quantify this more extensive similarity in various forms (e.g.ring statistics, Voronoi volume distributions, and higher-order correlation functions) in the Supplementary Material.A particular point of interest is that the MC configurations contain Ca-rich/-poor regions that are qualitatively similar to those observed in the original atomistic model [Fig.3(b)]-the clear implication being that, at the density of ACC, an inhomogeneous Ca distribution is encoded within the effective Ca• • • Ca interaction itself.
(a) and (b) the distribution of calcium and carbonate coordination numbers.The average Ca 2+ coordination number of 7.0-defined for a cut-off distance of 2.8 Å-is similar to that reported in Refs.22,27,28.Whereas most Ca 2+ bind either five or six distinct carbonate anions, the carbonates tend to bind one fewer Ca 2+ ion each, reflecting a binding-mode distribution of ∼80:20 monodentate:bidentate.Representative modal