Geometrically frustrated interactions drive structural complexity in amorphous calcium carbonate

Calcium carbonate is relatively unusual among simple inorganic salts in that it precipitates from aqueous solution in a metastable hydrated, amorphous form. Amorphous calcium carbonate (ACC)—with a nominal composition of CaCO₃·xH₂O (x ≈ 1)—can be stabilized for weeks by incorporating dopants such as Mg²⁺ or PO₄³⁻, or alternatively directed to crystallize into different polymorphs by varying pH or temperature. Nature exploits this complex phase behaviour in biomineralization to control the development of shells and other skeletal structures. The amorphous nature of biogenic ACC allows transformation to different crystalline CaCO₃ polymorphs and supports hierarchical morphologies important to biomineral architectures. In developing bio-inspired approaches for synthetic control over phase and morphology selection, there is a need to understand why such a chemically simple system exhibits complex phase behaviour. Theoretical studies of soft-matter systems show that multi-well pair potentials can drive structural complexity through geometric frustration, unlike single-well isotropic potentials that favour simple crystals. Double-well Lennard-Jones-Gauss (LJG) potentials can stabilize quasicrystals or complex crystals and in some cases frustrate crystallization. For ACC, experiments have hinted at two preferred Ca…Ca distances dominating medium-range order, corresponding to distinct carbonate bridging motifs (Ca–O–Ca and Ca–O–C–O–Ca). This work explores whether effective Ca²⁺–Ca²⁺ interactions reflecting these two length scales govern ACC structure and whether competition between them underlies ACC’s structural complexity and metastability. The approach is to obtain a high-quality measure of the Ca-pair correlation function in ACC via hybrid reverse Monte Carlo (HRMC), invert it to derive the effective Ca…Ca interaction potential, map it to an LJG form, and test via simulations whether such interactions reproduce the observed structural features and resistance to crystallization.

Prior studies established that ACC is a metastable hydrated phase formed from aqueous solution and can be stabilized by dopants such as Mg²⁺ or PO₄³⁻, or transformed into various crystalline polymorphs by environmental conditions. Experimental X-ray and neutron total scattering and NMR have characterized ACC, with observations of two dominant Ca…Ca separations in synthetic ACC attributed to different carbonate bridging modes. Earlier reverse Monte Carlo (RMC) models fitted scattering data but, lacking energetic constraints, produced unphysical charge-separated structures with Ca²⁺-rich domains and carbonate/water channels. Empirical potential structure refinement (EPSR) for calcium carbonate is challenging due to delicately balanced potentials; LJ-type Ca…Ca interactions used in some EPSR studies led to Ca–Ca pair distributions inconsistent with experimental transforms. In soft-matter theory, isotropic multi-well pair potentials, especially the double-well Lennard-Jones-Gauss (LJG) model, are known to generate complex and frustrated phase behaviour, including quasicrystals and large-unit-cell crystals, and can frustrate crystallization for certain parameter ranges. These insights motivate examining whether ACC’s two Ca–Ca preferred separations correspond to an effective double-well interaction that could rationalize its amorphous character and structural heterogeneity.

HRMC refinements were performed on atomistic configurations containing 12,960 atoms (1,620 CaCO₃·H₂O units) in simulation boxes ~5 nm, using a custom Python code interfaced with LAMMPS. Moves (atomic displacements up to 0.1 Å; rigid-body translations 0.1 Å and rotations 1° for water and carbonate) were accepted by a Metropolis criterion with a cost function combining fit to X-ray total scattering data and energetic stability. The fit term used the weighted X-ray total scattering function QF_x(Q) from experimental data (1.2 ≤ Q ≤ 25 Å⁻¹), computed by weighted Fourier transforms of partial pair correlation functions with Δr = 0.02 Å. The energetic term used the Raiteri et al. (2010) rigid-body force field for aqueous calcium carbonate via LAMMPS at T = 300 K. The relative weighting parameter σ was determined empirically (σ = 0.057) to balance goodness-of-fit and energy. Starting configurations were generated by placing 1,620 Ca atoms randomly in a box (52.8 × 54.8 × 45.2 ų), applying short-range repulsive interactions to avoid Ca–Ca < 3.0 Å, and adding stoichiometric rigid carbonate (r_CO = 1.284 Å) and water molecules (TIP4P-Ew; r_OH = 0.9572 Å, θ_HOH = 104.52°), followed by minimization and brief NVT MD to distribute particles homogeneously. No closest-approach constraints were needed due to explicit potentials in HRMC. Convergence was assessed by stabilization of χ² and energy; stability was also tested by 1 ns NVT MD at 300 K using LAMMPS (time step 1 fs, Nosé–Hoover thermostat). The HRMC configuration decreased in energy by ~−0.3 eV per Ca during equilibration, and NPT MD indicated densification from 2.43 to 2.91 g cm⁻³ under pressure. Pair distribution function inversion: The effective Ca–Ca pair potential u_Ca(r) was extracted from the HRMC configurations using a test-particle insertion based inversion method adapted to 3D. Ca coordinates and averaged g(r) up to 12 Å (Δr = 0.2 Å) from the final 12 HRMC frames were input, using ~10,000 test-particle insertions per configuration to achieve convergence. Orientational correlation functions φ(r) = ⟨P₂(S·S)⟩ were computed for carbonate and water local axes to assess anisotropy; φ(r) was featureless beyond very short ranges, supporting an isotropic effective potential description for Ca–Ca at relevant separations. LJG parameterization and coarse-grained simulations: The extracted u_Ca(r) showed two wells and was fit by a double-well Lennard-Jones-Gauss (LJG) potential with an added broad repulsive Gaussian term to reproduce local maxima; least-squares fit yielded ε = 4.1, r₀ = 1.4, σ = 0.14. Coarse-grained Monte Carlo (MC) and NVT molecular dynamics (MD) simulations of Ca particles at the experimental ACC Ca density were performed using the parameterized LJG potential in LAMMPS. Twelve independent starting configurations (1,620 Ca atoms per box with the same cell dimensions as HRMC) were generated randomly. MC used maximum displacements of 0.1 Å and was run until 4 million accepted moves with energy convergence. MD runs employed a 0.5 fs timestep, 50 fs thermostat damping, and 1 ns duration after initial energy minimization. Structural analyses included pair correlation functions, ring statistics, Voronoi volumes, and higher-order correlations to compare with HRMC.

• HRMC provided atomistic ACC structures that simultaneously fit X-ray total scattering data and are energetically reasonable with state-of-the-art aqueous CaCO₃ potentials. HRMC and RMC fit QF_x(Q) similarly well, but the RMC model is >800 kJ mol⁻¹ per formula unit less stable than HRMC and MD-based models. - The HRMC structure exhibits heterogeneous distribution: water molecules form filamentary, percolating strands separating CaCO₃-rich regions (“blue cheese”), without unphysical charge-separated channels seen in earlier RMC-only models. - Local coordination is consistent with prior reports: average Ca²⁺ coordination number ≈ 7.0 (cutoff 2.8 Å). Most Ca²⁺ bind five or six distinct carbonates; carbonates bind one fewer Ca²⁺ on average, with ~80:20 monodentate:bidentate binding. - The Ca–Ca pair correlation function g_Ca(r) shows two principal peaks at ~4 Å and ~6 Å. These correspond to two carbonate-bridged Ca–Ca motifs: sharing a common oxygen (shorter separation) and Ca–O–C–O–Ca bridging (longer separation). - Inversion of g_Ca(r) yields an effective Ca···Ca pair potential u_Ca(r) with two minima located near the two preferred Ca–Ca separations; the ~6 Å minimum is deeper, consistent with reduced electrostatic repulsion for Ca bound to the same carbonate. - Orientational correlations of water and carbonate vanish at distances relevant to Ca–Ca separations, validating an isotropic effective pair potential for Ca–Ca at ACC density. - u_Ca(r) maps well to a double-well LJG model with fitted parameters ε = 4.1, r₀ = 1.4, σ = 0.14 (plus a broad repulsive Gaussian to reproduce local maxima). - Coarse-grained MC/MD simulations using this LJG potential reproduce g_Ca(r), higher-order structural features (ring statistics, Voronoi volumes), and the emergence of Ca-rich/poor regions similar to the atomistic HRMC model. - The fitted LJG parameters place ACC in a strongly frustrated region of LJG phase space (1.2 ≤ r₀ ≤ 1.6; σ ≈ 0.14; ε ≥ 1), known to support complex structures and resist crystallization. This provides a mechanistic explanation for ACC’s structural complexity and metastability.

The findings address the central question of why ACC, a chemically simple system, exhibits structural complexity and resists crystallization. The effective Ca···Ca interaction possesses two competing minima at ~4 and ~6 Å, reflecting distinct carbonate-bridging modes. These incompatible preferred separations introduce geometric frustration that disfavors periodic, three-dimensional crystalline order at ACC’s density when interactions are effectively isotropic. Mapping u_Ca(r) to an LJG potential with parameters in the highly frustrated regime rationalizes the persistence of amorphous structure and the emergence of heterogeneous Ca-rich/poor regions. This framework explains why simple single-well potentials cannot reproduce the observed g_Ca(r) and why earlier EPSR approaches mislocated Ca–Ca peak positions. It also suggests that modest changes in interaction parameters (e.g., via hydration level, density, or dopant incorporation such as Mg²⁺ or PO₄³⁻) can tune the relative depths and positions of the wells, directing ACC towards different crystalline polymorphs or further stabilizing the amorphous state. As ACC densifies and orientational correlations strengthen during dehydration/crystallization, the isotropic effective pair description is expected to break down, consistent with a transition to anisotropy-governed crystalline phases. The study provides a first experimental mapping of an inorganic system to multi-well isotropic potentials and suggests broader applicability to complex inorganic amorphous materials.

This work delivers the first atomistically informed, energetically sensible ACC structural model that reproduces X-ray scattering and reveals the effective Ca···Ca interaction via pair-distribution inversion. The effective interaction has a double-well form corresponding to two carbonate-bridging motifs and maps quantitatively to an LJG potential with parameters in a strongly frustrated regime. Coarse-grained simulations with this potential replicate key structural features, explaining ACC’s complexity and resistance to crystallization through geometric frustration of competing Ca–Ca length scales. The framework suggests an ‘interaction engineering’ strategy for designing complex or amorphous inorganic materials by tuning cation size, hydration, and composition to control the ratio and depths of preferred separations. Future work should explore temperature and density effects on LJG phase behaviour relevant to ACC, investigate the impact of hydration changes and dopants on effective potentials, and extend the approach to other poorly ordered inorganic solids such as amorphous calcium phosphate and calcium-silicate-hydrates.

• The effective Ca···Ca interaction is modeled as an isotropic pair potential; while orientational correlations are negligible at ACC density for relevant separations, anisotropy is expected to become important at higher densities or near crystallization. - The LJG fit includes an additional broad repulsive Gaussian term to reproduce local features; its influence on phase behaviour is assumed to be minor but not rigorously established. - The mapping to LJG phase space is based on ground-state studies; finite-temperature and fixed-density effects for the specific parameters used here have not been comprehensively charted. - Changes in water content and bulk density during ageing/dehydration of ACC will alter effective interactions and may invalidate the isotropic coarse-grained description as crystallization is approached, limiting direct extrapolation. - EPSR limitations and RMC artifacts underscore sensitivity to potential parameterization; while HRMC balances data and energetics, residual model dependence remains.