Prediction of equilibrium Li isotope fractionation between minerals and aqueous solutions at high P and T: an efficient ab initio approach

The mass-dependent equilibrium stable isotope fractionation between different materials is an important geochemical process. Here we present an efficient method to compute the isotope fractionation between complex minerals and fluids at high pressure, P, and temperature, T, representative for the Earth’s crust and mantle. The method is tested by computation of the equilibrium fractionation of lithium isotopes between aqueous fluids and various Li bearing minerals such as staurolite, spodumene and mica. We are able to correctly predict the direction of the isotope fractionation as observed in the experiments. On the quantitative level the computed fractionation factors agree within 1.0 permil with the experimental values indicating predictive power of ab initio methods. We show that with ab initio methods we are able to investigate the underlying mechanisms driving the equilibrium isotope fractionation process, such as coordination of the fractionating elements, their bond strengths to the neighboring atoms, compression of fluids and thermal expansion of solids. This gives valuable insight into the processes governing the isotope fractionation mechanisms on the atomic scale. The method is applicable to any state and does not require different treatment of crystals and fluids.

💡 Research Summary

The paper presents a computationally efficient ab initio approach for predicting equilibrium stable isotope fractionation of lithium (⁶Li/⁷Li) between complex silicate minerals and aqueous solutions under high‑pressure, high‑temperature (HP‑HT) conditions typical of the Earth’s crust and mantle. Traditional methods for calculating isotope fractionation rely on full normal‑mode analysis of the vibrational spectrum of the entire system or on cluster models that treat the fluid as a static hydration shell. Both strategies become prohibitively expensive for systems containing dozens to hundreds of atoms, especially when periodic boundary conditions are required to model fluids realistically.

The authors adopt the single‑atom approximation originally derived by Bigeleisen and Mayer. In the high‑temperature limit (T > 600 K) and for vibrational frequencies satisfying ħω < 2kBT (u < 2), the reduced partition function ratio (β‑factor) can be expressed solely in terms of the three force constants acting on the isotopic atom along the Cartesian axes. The resulting formula (Eq. 3) reduces the computational effort from O(N) to O(1) with respect to the number of atoms N, because only the Hessian elements associated with the isotope of interest need to be evaluated.

To test the method, the authors performed density‑functional theory (DFT) calculations using the CPMD code with the BLYP exchange‑correlation functional and Goedecker pseudopotentials. Crystalline phases (staurolite, spodumene, and four polytypes of mica) were modeled with periodic supercells containing 40–80 atoms, while the aqueous solution was represented by a periodic box of 64 water molecules, one Li⁺ ion, and one F⁻ counter‑ion. Geometry optimizations employed a 70 Ryd cutoff, and vibrational analyses used a higher 140 Ryd cutoff to ensure convergence of force constants. Ab‑initio molecular dynamics (AIMD) simulations in the NVT ensemble (Car‑Parrinello scheme with Nosé–Hoover chain thermostat) generated statistically independent snapshots of the fluid; for each snapshot the force constants on Li were extracted and averaged to obtain the fluid β‑factor.

Two routes to β were compared: (i) the conventional full‑spectrum calculation (Eq. 1) using all normal modes, and (ii) the single‑atom approximation (Eq. 3). Across the temperature range 600–1200 K the two approaches yielded virtually identical β‑values; at 1000 K the discrepancy was less than 0.01 ‰. Fractionation factors α = β_min/β_fluid were then computed for each mineral–fluid pair. The predicted α values reproduced the direction of fractionation observed experimentally (⁷Li preferentially entering the fluid) and matched measured fractionation factors within 1 ‰, demonstrating that the simplified approach retains quantitative accuracy.

Beyond validation, the study provides mechanistic insight. By examining the computed force constants, the authors show that (a) Li in four‑fold coordination with short Li–O bonds exhibits larger force constants, favoring the heavier isotope; (b) increasing pressure shortens Li–O bonds, raising force constants and thus enhancing ⁷Li enrichment; (c) at high temperature the hydration shell of Li⁺ is highly dynamic, with coordination numbers fluctuating between 4, 5, and 6, which reduces average force constants and promotes ⁶Li incorporation. These atom‑scale observations rationalize experimental trends linking coordination environment, bond length, and pressure to isotope partitioning.

Importantly, the method treats solids and fluids within the same periodic framework, eliminating the need for separate cluster or lattice‑dynamics treatments. This unified approach is readily extensible to more complex multi‑component systems, such as mixed‑valence oxides, fluid‑rock interfaces, or high‑pressure melts, where traditional normal‑mode calculations would be infeasible. The authors argue that the combination of reduced computational cost, high accuracy, and mechanistic transparency makes the single‑atom force‑constant approach a powerful tool for geochemical modeling of isotope systems under mantle‑relevant conditions.

In summary, the paper delivers (1) a theoretically sound, temperature‑dependent simplification of the β‑factor based on force constants; (2) a robust implementation using DFT‑AIMD that accurately reproduces Li isotope fractionation between several silicates and aqueous solutions; (3) a clear demonstration of how pressure, temperature, and coordination affect fractionation at the atomic level; and (4) a scalable computational protocol applicable to a broad range of geologically relevant isotope systems.