Deuteration favours the polymorph with shorter O-H bonds: a stretch-bend rule tested across 220 inorganic polymorph pairs

Introduction

Deuterium replaces hydrogen in many experiments on crystals. Neutron crystallographers deuterate samples to reduce the incoherent background of hydrogen and often read the deuterated phase as a stand-in for the hydrogenous one. Isotope tracing, isotope geochemistry and synthesis in D$_2$O or D$_2$ rely on the same substitution. Whether the swap changes which crystal structure is stable is a practical question, and in some systems the answer is yes. In $\delta$-AlOOH, a candidate carrier of water into the deep Earth, neutron diffraction places the proton-disorder transition from $P2_1nm$ to $Pnnm$ at 9.0 GPa for AlOOH and at 13.4 GPa for AlOOD, a shift of 4.4 GPa, while hydrogen-bond symmetrization is observed at 18.1 GPa in AlOOH only [1, 2]. In KH$_2$PO$_4$ deuteration nearly doubles the ferroelectric transition temperature, the canonical solid-state isotope effect on a phase transition [3, 4]. In molecular crystals deuteration selects a different polymorph of pyridine [5] and of the pentachlorophenol and 4-methylpyridine complex [6, 7, 8], selects one of two concomitant polymorphs of ROY [9], gives a new phase in a hydrogen-bonded host-guest crystal through a geometric isotope effect [10], and changes the phase transition of 3,5-pyridinedicarboxylic acid [11]; deuterium can perturb molecular arrangement in the solid state more generally [12]. Some of these cases concern which form crystallises, not an equilibrium free-energy difference; the quantity computed here is the latter.

The physics behind these effects is understood in outline. The isotope effect on a free energy comes from the zero-point and thermal energies of the modes in which hydrogen moves [13, 14]. A stronger hydrogen bond lengthens and softens the covalent O-H bond, which lowers the stretch frequency (that longer bonds have smaller force constants is Badger’s rule [15]), and stiffens the bending and librational modes. The two changes act in opposite directions on the H/D free-energy difference. McKenzie and co-workers derived this competition for the isotope fractionation and geometry of hydrogen bonds [16, 17], building on the link between fractionation factors and hydrogen-bond strength [18], Li, Walker and Michaelides showed that the net quantum effect on hydrogen-bond strength changes sign with the strength of the bond [19], and Markland and Berne showed the same competition in the isotope fractionation of water [20]. Hydrogen bonding is known qualitatively to control D/H fractionation between OH-bearing minerals and water [21, 22, 23, 24], and geometric H/D isotope effects on hydrogen-bond length are well documented [3, 25, 26, 27]. In hydrogen-bonded ferroelectrics of the KH$_2$PO$_4$ type, H and D crystals of one isostructural family follow a single relation between transition temperature and O$\cdots$O distance, which Ichikawa and co-workers interpreted through the geometric isotope effect rather than proton tunnelling [28]; first-principles work couples tunnelling and proton delocalisation to the lattice geometry [4]. Path-integral methods give isotope effects beyond the harmonic approximation for single systems [29, 30, 31, 32, 33], and anharmonic free energies change polymorph rankings in molecular crystals [34, 35, 36]. Large-scale harmonic isotope thermodynamics has been reported for 115 metal hydrides [37], where the gas dominates the formation reaction, and universal machine-learned interatomic potentials (MLIPs) now give phonons across many systems [38, 39], without isotopes.

What our search did not find is a cross-compound, quantitative statement for solid-solid polymorph pairs: given two polymorphs of a hydrogen-bonded compound, which one does deuteration favour, by how much, and where does the statement fail? A logged literature search in three rounds (20, about 55 and about 150 Valency Bond and web-search queries, mostly over abstracts and search snippets) found no paper that states such a rule across many compounds or reports a slope of the polymorph isotope free-energy difference against O-H length. The statement is bounded by that search. A matching qualitative ordering has been reported in high-temperature mineral-water D/H fractionation for boehmite and diaspore and for the epidote group [22]; the rule tested here concerns polymorph stability at room temperature, where bend and libration cancel much of the stretch term. The mechanism, the competition between stretch and bend zero-point energies, is that of McKenzie and co-workers [16, 17]. Here we formulate a general rule from it and test it. We write the rule as one equation with an explicit condition of validity (Sec. II), measure its slope across 220 polymorph pairs in 63 formulas, split that slope into stretch and bend and libration terms, score its sign on 84 pairs in 28 formulas against baselines and competing descriptors, locate the one structure where the condition fails, and turn the rule into a pressure prediction.

Isotope effects on polymorph free energies are a few meV per atom, smaller than the static energy errors of fast models. We therefore compute only the isotope double difference, in which the static energy gap between the two polymorphs cancels, with three universal MLIPs over polymorph pairs from the Materials Project and JARVIS-DFT. This paper reports: (i) the rule in closed form, $\Delta\Delta F_{\mathrm{iso}} \simeq S_{\mathrm{str}}(1-C)\Delta\bar{r}$, with a one-mode harmonic expression for the stretch slope, the condition $C < 1$, a crossover temperature where $C = 1$, and the high-temperature force-constant limit; (ii) the net slope of the harmonic H/D free-energy difference between polymorphs against the change in mean X-H length, with clustered uncertainties and a mode-resolved split into stretch and bend and libration terms; (iii) the sign rule scored with pair and formula counts, competitor descriptors, an orientation-invariant null and the one structure in which it fails; (iv) where the rule is reliable, namely where the effect is larger than the model spread, its failure for weak hydrogen bonds, and its temperature dependence; (v) every blind and prospective test we ran, including four that failed or did not reach their bar; (vi) an exploratory prediction for the diaspore to δ-AlOOH transition and the measurements that would test the rule. The calculations are harmonic or quasi-harmonic, the potentials share training data, and no polymorph isotope free energy has been measured directly; the comparisons with experiment in § 1.4 use related measured quantities. The result is a computational relationship within these models, and we state its limits where they arise.

The rule

For polymorphs $a$ and $b$ with mean O-H lengths differing by $\Delta\bar{r}=\bar{r}_{\mathrm{OH}}(b)-\bar{r}_{\mathrm{OH}}(a)$, write the isotope double difference per H (Eq. (3)) as a stretch term and a bend and libration term, $\Delta\Delta F_{\mathrm{iso}}=\Delta\Delta F_{\mathrm{str}}+\Delta\Delta F_{\mathrm{L}}$ (Sec. III), and let $C=-\Delta\Delta F_{\mathrm{L}}/\Delta\Delta F_{\mathrm{str}}$ be the fraction of the stretch term that bend and libration cancel. To first order in $\Delta\bar{r}$,

The stretch slope $S_{\mathrm{str}}$ follows from one harmonic mode. For a mode of frequency $\omega$, $F=k_{\mathrm{B}}T\ln[2\sinh(\hbar\omega/2k_{\mathrm{B}}T)]$, and the isotope ratio of an O-H stretch is $\rho=\omega_{\mathrm{H}}/\omega_{\mathrm{D}}=\sqrt{\mu_{\mathrm{D}}/\mu_{\mathrm{H}}}$ with the reduced masses of ${}^{16}\mathrm{O}$ with H and D, $\rho=1.373607$ and $1/\rho=0.728010$. Using $1/\sqrt{2}$ instead overstates the zero-point shift by about 8% and the high-temperature slope by about 6%. Let $s=-\mathrm{d}\nu_{\mathrm{str}}/\mathrm{d}r_{\mathrm{OH}}$ be the softening of the stretch wavenumber as the O-H bond lengthens. Linearising gives

with $u=hc\nu_{\mathrm{H}}/k_{\mathrm{B}}T$. For stretch wavenumbers of 3000 to 3600 cm$^{-1}$ the bracket equals its zero-temperature value to within $1.5\times10^{-4}$ (relative) at 300 K, so $S_{\mathrm{str}}=0.016861$ meV per H per mÅ with $s$ in cm$^{-1}$ per mÅ. A longer O-H bond softens the stretch, which makes $F_{\mathrm{D}}-F_{\mathrm{H}}$ less negative, so $S_{\mathrm{str}}>0$.

Bends and librations harden as the hydrogen bond strengthens. With $s_b=+\mathrm{d}\nu_b/\mathrm{d}r_{\mathrm{OH}}$ for each of $n_b=2$ bending modes per H, $S_{\mathrm{bend}}=-n_b s_b(hc/2)[\coth(u_b/2)-\rho_b^{-1}\coth(u_b/2\rho_b)]$, with $u_b=hc\nu_b/k_{\mathrm{B}}T$ and $\rho_b=\nu_{b,\mathrm{H}}/\nu_{b,\mathrm{D}}$, and at 0 K with $\rho_b\approx\rho$ the compensation is $C=n_b s_b/s$.

Three statements follow from Eq. (1). (1) Sign. While $C<1$, $\operatorname{sign}\Delta\Delta F_{\mathrm{iso}}=\operatorname{sign}\Delta\bar{r}$: deuteration favours the polymorph with the shorter O-H bonds, that is, the one with the weaker hydrogen bonds. (2) Crossover. Relative to 0 K the thermal weight of a 3400 cm$^{-1}$ stretch is 1.000 at 300 K and 0.988 at 600 K, while a 1200 cm$^{-1}$ bend falls to 0.941 and 0.687 and a 900 cm$^{-1}$ libration to 0.858 and 0.556. In this illustrative model $C$ falls with temperature: with $s=18$, $s_b=4.5$ and $\nu_b=1400$ cm$^{-1}$, $C$ falls from 0.50 at 0 K (by construction, since $n_b s_b/s=0.50$) to 0.48 at 300 K and 0.38 at 600 K. The projected ratio of a real pair need not fall; in the screen it falls in 198 of 212 rows between 0 and 600 K (Sec. IV F). If $C(T)$ crosses unity from above while the stretch term keeps the sign of $\Delta\bar{r}$, the harmonic isotope preference reverses at the temperature $T^*$ where $C(T^*)=1$: below $T^*$ the pair violates the sign rule, above it the pair obeys it. (3) High-temperature limit. At high temperature $\Delta\Delta F_{\mathrm{iso}}$ decays as $1/T$ and its sign is set by the change in the trace of the on-site H force constants (Sec. III), the Bigeleisen-Mayer and Urey statement that deuterium concentrates where the H force constants are largest [13, 14].

Size within the models. The input $s$ is taken from the screen’s own phonons (donors with stretch weight above 0.5 and mean O-H of 0.96 to 1.05 Å): the stretch frequencies fall by 19.0 (MACE), 19.3 (ORB) and 19.5 (MatterSim) cm$^{-1}$ per mÅ, with $r = -0.99$ to $-1.00$ (17 to 21 over narrower windows). With $s = 19$ to 19.5, Eq. (2) gives $S_{\mathrm{str}} = 0.32$ to 0.33 meV per H per mÅ, against fitted stretch slopes of 0.245 to 0.259 (all rows) and 0.277 to 0.308 (strong O-H regime) in Sec. IV A; the one-mode estimate is 4 to 34% high, and these fitted values come from two populations. Two likely reasons, not tested, are that the regressor is the mean O-H length of all H in the cell while $s$ is fitted per donor, and that the stretch is not a pure mode. The observed $C$ of about 0.5 requires $s_b \approx s/4$, roughly 5 cm$^{-1}$ per mÅ per bend; this value is inferred from the compensation, not measured. Because $s$, $C$ and the projected terms all come from the same phonons, this is a consistency check within the models, not an independent validation; an experimental $d\nu/dr_{\mathrm{OH}}$ was not checked against a primary source and is not used. The direction of the temperature dependence agrees with the rescoring of Sec. IV F, but Eq. (1) does not explain the rise from 153 to 185 right between 0 and 300 K; that would need librations and lattice modes below about 900 cm$^{-1}$ and rows near zero, a plausibility argument, not a derivation.

Assumptions. Harmonic modes, an isolated O-H stretch (anharmonicity, which lowers the measured O-H/O-D ratio to about 1.33 to 1.36, is outside the estimate), linear dependence of frequencies on O-H length and bends treated as two modes of one frequency. (1) was derived after the screen and is not used to select or score any pair. The condition $C < 1$ is not observable from geometry alone; the regime boundary $r_{\max} \ge0.99$ Å used to score the sign (Sec. IV B) is an empirical proxy for it in these model geometries.

Methods

Isotope double difference

For polymorphs $a$ and $b$ of the same composition we compute

with $F_X$ per common cell containing $N_{\mathrm{H}}$ hydrogen atoms, so that $\Delta\Delta F_{\mathrm{iso}}$ is per H unless stated otherwise (meV/atom is used only for the AlOOH pressure prediction and the code checks, and is labelled). A positive $\Delta\Delta F_{\mathrm{iso}}$ means that deuteration destabilises $b$ relative to $a$. The static energies do not depend on mass and cancel, so errors of the potentials in the polymorph energy gap, often larger than the isotope signal, do not enter.

The cancellation does not remove errors in the force constants or the relaxed geometries, or the omitted anharmonicity. Within the Born-Oppenheimer approximation the force constants $\Phi$ are mass independent, so one finite-displacement phonon calculation per structure gives both isotopes: the dynamical matrix $D_X=M_X^{-1/2}\Phi M_X^{-1/2}$ is rebuilt with the D mass on every H site and diagonalised again ($m_{\mathrm{H}} = 1.00782503$ u, $m_{\mathrm{D}} = 2.01410178$ u). For $X \in\{\mathrm{H},\mathrm{D}\}$,

with $q$-point weights $w_q$ summing to one; in this direct sum modes below 0.1 THz, including imaginary modes, are dropped. The classical momentum term $-\frac{3}{2}N_{\mathrm{H}}k_{\mathrm{B}}T\ln(m_{\mathrm{D}}/m_{\mathrm{H}})$ is the same for both polymorphs and cancels; at high temperature the surviving term per H is the Wigner-Kirkwood correction $(\hbar^2/24N_{\mathrm{H}}k_{\mathrm{B}}T)(1/m_{\mathrm{D}}-1/m_{\mathrm{H}})\sum_i[\operatorname{Tr}\Phi_{ii}^{b}-\operatorname{Tr}\Phi_{ii}^{a}]$ over H sites $i$, so $\Delta\Delta F_{\mathrm{iso}}$ is $1/T$ and D prefers the phase with the larger H force-constant trace [13, 14]. An absolute transition pressure, unlike $\Delta\Delta F_{\mathrm{iso}}$, depends on the static energy surface; for the pressure predictions of Sec. IV G each isotope has its own quasi-harmonic minimization, $G_X(T,P)=\min_V[E(V)+F_{\mathrm{vib},X}(V,T)+PV]$, reported separately from the harmonic screen. All formulas used in the code were derived twice, independently, and cross-checked numerically before use; derivations and test cases are in the repository [40] as notes M01 to M03.

Mode-resolved decomposition

To separate the stretch from the bend and libration we use the exact mass identity obtained from the Hellmann-Feynman theorem applied to $\Phi u = \omega^2 M u$, with mass-weighted eigenvectors $e=M^{1/2}u/\lVert M^{1/2}u\rVert$,

where $e_{q\nu}$ is the normalised mass-weighted eigenvector at H mass $m$ and $w_{q\nu}^{\mathrm{H}}=\sum_{i\in\mathrm{H}}\lvert e_{i,q\nu}\rvert^2$. We split $w^\mathrm{H}=w^\mathrm{str}+w^\perp$, with $w^\mathrm{str}=\sum_{i\in\mathrm{H}}\lvert\hat b_i\cdot e_{i,q\nu}\rvert^2$ and $\hat b_i$ the unit vector from H to its nearest non-H atom. The integral is evaluated with six-point Gauss-Legendre quadrature in $\ln m$, rediagonalising at each intermediate mass and using the eigenvectors at that mass. This gives $\Delta\Delta F_{\mathrm{iso}}=\Delta\Delta F_\perp+\Delta\Delta F_{\mathrm{str}}$, and the compensation $C=-\Delta\Delta F_\perp/\Delta\Delta F_\mathrm{str}$, which is ill conditioned when the stretch term is near zero. We call the perpendicular term bend and libration; it is an operational component, and crystal modes mix local and collective motion. The projection avoids frequency cuts, which are ill defined because a strongly hydrogen-bonded O-H stretch and its O-D counterpart can fall on different sides of any cut. The quadrature integrates $|\omega|$ for imaginary modes.

Potentials and phonons

Every quantitative result comes from three universal MLIPs: MACE (checkpoint medium-mpa-0, mace-torch 0.3.16) [41], ORB (orb_v3_conservative_inf_mpa, orb-models 0.5.5) [42] and MatterSim (MatterSim-v1.0-5M, mattersim 1.2.5) [43]. SevenNet (7net-mf-ompa, modal mpa, sevenn 0.13.0) [44, 45] enters only as sign votes and in screen counts, which are labelled. It was set aside before the screen because its transverse optical H frequency in LiH is 15.3% above experiment (20.41 against 17.7 THz); MACE is 12.8% above (19.97 THz), so the cut separates models at one tolerance rather than on a physical criterion. The potentials are separately trained on overlapping PBE-level data, so their agreement is not an independent validation, and the model spread is a range, not a confidence interval.

Structures were relaxed with each potential, cell included, to a maximum force of $10^{-4}$ eV/Å with a limit of 500 steps; non-convergence is flagged, not excluded. Force constants were computed with phonopy [46, 47] by finite displacements of 0.01 Å in supercells with a minimum image length of 9.5 Å and at most 160 atoms. The double differences, their decomposition and the temperature sweep were computed on a shifted, $\Gamma$-free $q$ mesh with about $40/a$ points along each axis of lattice length $a$ (in Å), rounded to an even number between 4 and 12, which removes the treatment of the acoustic modes at $\Gamma$ as a source of spurious isotope terms; mesh convergence of $\Delta\Delta F_{\rm iso}$ was not tested across the screen; a later recomputation of the 300 K values with about $60/a$ points (even, 4 to 16) leaves 185 of the 212 primary evaluations right, changes no sign and moves $\Delta\Delta F_{\rm iso}$ by at most 0.05 meV per H, and for one pressure-prediction pair (goethite) a mesh change moved a pressure-dependent free energy by 0.1 to 0.2 meV/atom.

A structure-model result is called clean when the quadrature of Eq. (5) and the direct difference from Eq. (4) agree within 0.01 meV/atom. Because the two routes treat soft modes differently (the direct sum drops imaginary modes, the quadrature integrates $|\omega|$), the criterion tests the numerical consistency of the isotope term; it does not certify relaxation convergence, dynamical stability or insensitivity to the treatment of imaginary modes. Section IV C reports the subset that also has converged relaxations and no imaginary modes.

As a per-structure benchmark against PBEsol phonons of 40 entries of the MDR phonon database [39, 47] (36 for SevenNet), at the DFT geometries, the isotope zero-point energy difference has a mean absolute error of 0.7 to 1.5 meV per H (ORB 0.7, SevenNet 1.1, MatterSim 1.3, MACE 1.5), the 300 K H/D free-energy difference 0.5 to 1.1 meV per H and the X-H stretch frequency 30 to 45 cm$^{-1}$. At the model-relaxed geometries used in the screen the errors are larger: 1.4 to 1.9 meV per H for the isotope zero-point energy difference, 1.0 to 1.4 meV per H for the 300 K free-energy difference and 42 to 65 cm$^{-1}$ for the stretch, with model stretches on average 15 to 40 cm$^{-1}$ softer than PBEsol; 22 of the 40 PBEsol references have imaginary modes. These errors are comparable to the median primary $|\Delta\Delta F_{\rm iso}|$ of 2.4 meV per H; the double difference may cancel part of them, which was not tested. No DFT value of a polymorph double difference was computed. Quasi-harmonic volume isotope effects, $100(V_{\rm H}/V_{\rm D}-1)$, were compared with measurements before the screen: with the three models brucite Mg(OH)$_2$ gives +0.19 to +0.33% at 290 K against a measured +0.31% [48], LiH gives +1.03 to +1.27% at 300 K against +0.81% or +1.04% from two disagreeing sources, and an ordered ice XI proxy gives the measured negative sign with a magnitude 4 to 11 times too large [49].

Polymorph pairs

Hydrogen-containing structures without carbon were taken from the Materials Project [50] (4421 entries, fetched with mp-api on 21 October 2026) and the JARVIS-DFT dft_3d snapshot of 24 September 2025 [51] (3436 entries). Formula equality uses the reduced elemental composition. Structures were reduced to primitive cells and deduplicated with pymatgen’s StructureMatcher (length tolerance 0.2, site tolerance 0.15, angle tolerance 5°), with different space groups at a symmetry precision of 0.05 Å blocking a merge and Materials Project entries preferred. A pair is two distinct structures of one composition whose source energies (Materials Project corrected formation energies, or JARVIS OptB88vdW formation energies, never mixed) differ by at most 25 meV/atom and whose primitive cells hold at most 32 atoms. Polymorph $a$ is the one with the lower database energy. This gives 836 structures and 965 pairs; 569 pairs have a screen result in at least one of the three models (572 counting SevenNet). Database energies are used only to select and orient pairs; all physics comes from the MLIP relaxations and phonons.

A structure has covalent X-H bonds when the nearest non-H neighbour of every H lies within 1.3 Å. Of the 965 pairs, 220 (63 formulas) have clean results for covalent X-H bonds in at least one model; these 521 pair-model rows (MACE 188, ORB 163, MatterSim 170) form the slope fit. They include 41 rows of pairs without any O-H donor (X = B, N or F). The sign test uses O-H pairs only, defined as pairs in which at least one structure has an O-H donor; in every primary structure all H are bonded to O, so the mean X-H length there is the mean O-H length.

Regime, sign rule, baselines and statistical unit

For each pair and model, $x=\bar{d}_{\mathrm{XH}}(b)-\bar{d}_{\mathrm{XH}}(a)$ is the difference of the mean model X-H lengths in mÅ and $y=\Delta\Delta F_{\mathrm{iso}}$ per H. The rule predicts $\operatorname{sign}(y)=\operatorname{sign}(x)$: deuteration favours the polymorph with the shorter mean O-H bond. The regime variable $r_{\max}$ is the larger of the two mean model O-H lengths. An evaluation (one pair computed with one model) enters the sign test when $\lvert x\rvert\ge5$ mÅ. The strong side is $r_{\max}\ge0.99$ Å.

The 0.99 Å boundary, the choice of $r_{\max}$, the 5 mÅ eligibility and the O-H definition were set after inspecting the initial screen data and are retained in the reported result. We re-chose the boundary inside each leave-one-formula-out training fold (grid 0.97 to 1.03 Å in 0.01 Å steps, maximizing right minus wrong on the training strong side, lowest threshold winning ties) and scored the held-out formula; this covers the threshold value only. The alternatives for the other choices are shown in a multiverse table (Sec. IV B). The baseline is a constant sign, the majority sign of $y$ in each leave-one-formula-out training fold. Because $a$ is the lower-energy database polymorph, a positive constant sign means “deuteration destabilises the polymorph with the higher database energy”. Competitor descriptors (density, MLIP static energy, O$\cdots$O) are scored on the same rows, with their direction fixed physically and, separately, chosen in training folds. An orientation-invariant null flips the predicted sign of each pair (shared by its models) at random; $2\times10^{5}$ draws give its distribution.

Evaluations share structures, models and formulas. We therefore report pair majorities (majority over the models of a pair) and formula majorities alongside evaluation counts, give formula-clustered bootstrap intervals (2000 resamples of formulas with replacement) for slopes and accuracies, and add Theil-Sen slopes and a structure-level fit with formula fixed effects, in which the per-structure isotope term $F_{\mathrm{D}}-F_{\mathrm{H}}$ per H is regressed on mean X-H length within each formula, which removes pair combinatorics. These analyses were added in revision and are computed by src/isoshift_event/paper_stats.py from the original screen tables into results/paper_stats.json.

The primary counts were fixed when the initial screen analysis ended, and the statement of the result was then set in writing; analyses performed subsequently are identified below as later analyses.

Validation against measured spectra

For every O-H donor in the clean structures we compared the harmonic $\Gamma$-point stretch frequency with the empirical correlation of Libowitzky [52], $\nu=3592-304\times10^{9}\exp(-d_{\mathrm{O\cdots O}}/0.1321)\ \mathrm{cm}^{-1}$ with $d$ in Å, built from measured infrared and Raman spectra of minerals. The stretch mode of a donor is the $\Gamma$ mode above $1500\ \mathrm{cm}^{-1}$ with the largest projection on that O-H bond, the acceptor is the nearest other O with H$\cdots$O below 2.6 Å (no angle criterion), and symmetry-equivalent donors are merged. In each O$\cdots$O window the model frequency is regressed on the Libowitzky frequency at the same distance; the slope ratio of equals 1 when the measured slope is reproduced. Intervals are 95% bootstrap intervals over structures (2000 resamples).

Predictions and comparison with experiment

For each test against experiment the model predictions and scoring rule were recorded before the measured values were read; the literature agents that collected measurements did not see the predictions, and the measured values were recorded before the comparison. There are three qualifications: the scoring rule of the mineral fractionation test was recorded together with its result; the free energies at 0, 600 and 2000 K existed before the temperature hypotheses were specified, so the temperature rescoring is retrospective even though the decomposition tests had fixed criteria; and one temperature test was scored after a code fix prompted by a failing run (Sec. IV F). The prospective test fixed signs, eligibility and scoring before any of its phonon jobs ran. Table II summarises the tests; the prediction records are available in the repository, and the AI-use agent workflow is described in the AI-use statement.

What was predicted	Outcome	Notes
Volume isotope effects of known H compounds	Mixed	$\alpha$-AlH$_3$ right sign, 1.3 to 1.5 times too large; GaOOH and $\beta$-CrOOH wrong sign against references of limited quality. Does not test the rule.
Stretch and bend shifts, goethite to $\varepsilon$-FeOOH and diaspore to $\delta$-AlOOH	Bend sign agrees, stretch sign not determined	Goethite to $\varepsilon$-FeOOH: bend rises as predicted; the stretch comparator range ($-511$ to $-11$ cm$^{-1}$) does not fix the stretch sign. Diaspore pair inconclusive.
Stretch and bend shifts, manganite to groutite	Inconclusive
Harmonic mineral D/H fractionation, 118 structures	Ordering only	Signs 6 of 6 and order brucite $>$ goethite $>$ manganite right; brucite minus goethite 5 to 10 times too large; scoring rule recorded together with its result.
H$_2$O to D$_2$O shift of three ice boundaries	Failed	Ih-II and Ih-III wrong sign in 3 of 3 models; II-III right sign; secondary data.
Rule refined with the O-H$\cdots$O angle	Failed	No gain over the bare rule; all 10 wrong pairs stay wrong.
Isotope shifts in deep-Earth hydrous phases	Not scored	Prediction recorded; no measured comparison was made.
Geometry-only sign predictions for 1796 pairs	Not scorable	None of 26 literature rows was usable.
Temperature dependence and four mechanism tests	Compensation-drop test failed; other three passed	High-temperature trace-sign test passed after correcting the pair orientation in the scoring code (initial run 124 of 212).
Prospective test on pairs never computed	Failed	16 of 111 right; 87 rows lost to predeclared numerical rules; on 24 clean rows 16 right, baseline fixed in advance 5, opposite constant sign 19; pair majority 5/1/4.
Boehmite against diaspore D/H fractionation	Inconclusive	Sign agrees at 380 °C in 3 of 3 models; no stated uncertainty, diaspore 25% exchanged.
Path-integral check of $\Delta\Delta F_{\mathrm{iso}}$ on primary pairs	Pending	Predeclared before any data; no result here.
Reimplementation without the original code	Code agrees	9 of 9 referenced rows to $10^{-9}$ meV/atom; tests code, not physics.

Results

Slope of the isotope term against bond length

Across the 220 pairs $\Delta\Delta F_{\mathrm{iso}}$ per H at 300 K rises with the change of mean X-H length (Figure 2, Table I). Least-squares slopes over all clean covalent rows are 0.112 (MACE), 0.117 (ORB) and 0.126 (MatterSim) meV per H per mÅ, with Pearson $r$ of 0.76, 0.83 and 0.84; the fits have intercepts of +0.68, +0.47 and +0.63 meV per H, and fits through the origin give 0.114, 0.120 and 0.130. The range 0.11 to 0.13 is the model spread. Formula-clustered bootstrap intervals are 0.068 to 0.154, 0.083 to 0.165 and 0.090 to 0.169. Theil-Sen slopes, less sensitive to the few pairs with $\lvert x\rvert$ up to about 100 mÅ, are 0.061, 0.111 and 0.109; the MACE slope depends on a few high-influence pairs, and an Al(OH)$_3$ alone supplies 43 of its 188 rows. Leaving out one formula at a time moves the slopes within 0.096 to 0.124 (MACE), 0.111 to 0.135 (ORB) and 0.121 to 0.145 (MatterSim). Restricted to O-H pairs the slopes are 0.124, 0.135 and 0.145. The structure-level fit with formula fixed effects, which uses each structure once, gives 0.120 (0.072 to 0.160; 167 structures, 55 formulas), 0.116 (0.087 to 0.153) and 0.126 (0.091 to 0.163). A positive slope means that the polymorph with the longer X-H bonds, which in these solids is the one with the stronger hydrogen bonds and the shorter donor-acceptor distance, is destabilised by deuteration. The slopes are associations across polymorph pairs, not derivatives obtained by stretching one structure. A difference of 10 mÅ shifts the H/D free-energy balance by about 1.1 to 1.3 meV per H; the median difference between Materials Project O-H polymorphs is 4.1 mÅ in database geometries.

The mode-resolved integral splits the slope into two opposing terms. Over all rows the stretch term has a slope of +0.245, +0.259 and +0.252 meV per H per mÅ and the bend and libration term −0.133, −0.142 and −0.124, so 0.49 to 0.55 of the stretch slope is cancelled. In the O-H strong regime the total slopes are 0.132, 0.136 and 0.147, the stretch slopes 0.277, 0.308 and 0.293 and the bend and libration slopes −0.145, −0.171 and −0.146, a compensation of 0.52, 0.56 and 0.50. The two populations differ and are quoted separately. SevenNet has 23 clean covalent rows and an all-row slope of 0.105; it is drawn in Figure 2 for comparison and not used in any claim.

The sign rule

When at least one polymorph has a mean model O-H length of 0.99 Å or more (a model-derived boundary for moderate to strong hydrogen bonds, chosen after inspecting the data), deuteration favours the polymorph with the shorter O-H bonds in 185 of 212 O-H model evaluations with $|\Delta\text{ mean O-H}| \ge5$ mÅ (84 pairs, 28 formulas), against 131 of 212 for a constant-sign baseline fitted on leave-one-formula-out training folds. By model majority 74 pairs are right and 10 are wrong, and all 10 wrong pairs contain one Al(OH)$_3$ structure, mp-985293. The more conservative units tell the same story. By formula majority the rule is right in 27 of 28 formulas; the one wrong formula is Al(OH)$_3$, all of whose primary rows involve mp-985293. Weighting formulas equally the rule is right in 0.92 of rows (bootstrap interval 0.85 to 0.98) against 0.49 for the baseline. The pooled accuracy, 0.87, has a formula-clustered interval of 0.68 to 0.98. The 212 evaluations come from only 89 distinct structures, and five oxyhydroxide and hydroxide formulas (AlOOH 39 rows, FeOOH 31, Al(OH)$_3$ 28, MnOOH 16, CrOOH 15) contribute 129 of them, so the result is weighted towards these families. Per model the rule is right in 52 of 68 MACE, 64 of 73 ORB and 69 of 71 MatterSim evaluations. With the direct sum of Eq. (4) in place of the quadrature count it is also 185 of 212. SevenNet, on its own 10 primary rows, is right in all 10 (baseline 7); as a fourth vote the count is 195 of 222 against 138, with 75 pairs right and the same 10 wrong.

Baselines and competitors. The constant-sign baseline (131 of 212) is the statement that deuteration destabilises the polymorph with the higher database energy. Against an orientation-invariant null that assigns each pair a random predicted sign, 185 lies far in the tail: the null has mean 106 and standard deviation 10.7, its largest value in $2 \times10^5$ draws is 154, so $p < 10^{-5}$ (the draws are not independent across shared structures, so this p-value describes the pair-level arrangement, not a sampling probability). A density rule, “deuteration destabilises the denser polymorph”, is right in 104 of 212, at chance; ordering by MLIP static energy gives 137. The O···O sign, “deuteration favours the polymorph with the longer mean donor O···O”, is right in 183 of the 211 primary rows that have an acceptor, against 184 for the O-H rule on the same rows (pair majority 73 right, 1 tie, 10 wrong); the two differences are correlated ($r = -0.80$). On like-for-like rows the two descriptors are therefore equivalent in sign. The O···O version needs no phonons, which makes it the practical form for screening, but it still uses relaxed geometries; with database geometries the sign agreement is lower (Sec. VI). Choosing each competitor’s direction inside training folds leaves these counts unchanged, except for density (72 of 212, because its training-fold direction is unstable).

Where the rule is reliable. The rule’s skill sits in rows far from zero (Figure 3b). For $|\Delta\Delta F_{\mathrm{iso}}| \ge1$ meV per H it is right in 151 of 159 rows; below 1 meV per H in 34 of 53, close to the baseline rate. Models disagree in sign in 19 of the 74 primary pairs computed with at least two models; the median spread across models within a pair is 1.36 meV per H against a median $|\Delta\Delta F_{\mathrm{iso}}|$ of 2.37. Normalised by that spread, the rule is right in all 81 rows with $|\Delta\Delta F_{\mathrm{iso}}|$ above twice the spread and in 42 of 43 between one and two spreads, but in 53 of 78 below one spread; the 10 rows of pairs computed with one model have no spread and are not binned. The stored classical-limit flag, which fires when $|\Delta\Delta F_{\mathrm{iso}}(2000\ \mathrm{K})| > |\Delta\Delta F_{\mathrm{iso}}(600\ \mathrm{K})|$ and so marks pairs near a sign crossover in temperature, separates the outcome almost completely: unflagged rows are right in 168 of 170, flagged rows in 17 of 42, including all 28 mp-985293 rows. This flag is a diagnostic tied to the outcome and is not used as a filter. Where the predicted effect exceeds the model spread the rule is nearly always right; near balance it has little skill.

The boundary. Re-choosing the boundary inside each leave-one-formula-out training fold picks 1.00 Å in 30 of 31 folds and 0.99 Å only in the fold that holds out Al(OH)$_3$, and gives 183 of 210 held-out evaluations right against a baseline of 131; in that Al(OH)$_3$ fold the rule scores 8 of 28 against a baseline of 19. Coverage changes with the boundary (Figure 3a): moving it from 0.98 to 1.02 Å raises accuracy from 185 of 229 to 114 of 114 while the scored rows fall by half, so accuracy and coverage must be read together, and 183 of 210 and 185 of 212 are not directly comparable. The other analysis choices matter as much (multiverse, results/paper_stats.json): with the regime defined by the smaller of the two mean O-H lengths the rule is right in 138 of 144 (baseline 100; pairs 56 right, none wrong), by their mean in 177 of 188 (114; pairs 74 right, 5 ties, none wrong), and by the shorter mean O···O $\le2.93$ Å in 185 of 222 (141); with eligibility at 3 or 10 mÅ instead of 5 it is right in 203 of 232 (141) and 165 of 188 (120). The rule beats the baseline under every variant, but a sharp, universal boundary is not established, and the boundary is an empirical region of these model geometries.

The Al(OH)$_3$ exception. Structure mp-985293 is a 28-atom $P2_12_12_1$ polymorph of Al(OH)$_3$ listed by the Materials Project as experimental, with no imaginary modes in any of the three models. Its mean model O-H of 0.996 (MACE), 0.997 (ORB) and 0.999 Å (MatterSim) sits just above the boundary, while its donors span 0.987 to 1.005 Å and its mean O···O contacts, 2.77 to 2.82 Å, are shorter than those of its partners; it is denser, 6.79 to 6.94 Å$^3$/atom against 7.53 to 7.88 for its primary partners. In its 28 primary evaluations the stretch term has the predicted sign, but bend and libration cancel 95 to 132% of it instead of about half (MACE 1.19 to 1.32, ORB 1.02 to 1.14, MatterSim 0.95 to 1.02, against a median of 0.50 in the other 184 primary rows). MACE calls 0 of 10 of its pairs right, ORB 0 of 9 and MatterSim 8 of 9, with every MatterSim value at $|\Delta\Delta F_{\mathrm{iso}}| \le0.26$ meV per H. Here bond length alone stops predicting the sign. We do not claim a phase identity or a mechanism for its larger compensation; a later test of whether the mean O-H···O angle predicts it failed (Table II). As post hoc locators chosen after inspecting the initial screen results, leaving out the 10 pairs with mp-985293 gives 177 of 184 and raising the boundary to 1.00 Å gives 175 of 182; these are nearly the same intervention, because 28 of the 30 evaluations between 0.99 and 1.00 Å come from mp-985293. Neither replaces the original count.

A physically cleaner subset

The clean criterion does not exclude unconverged relaxations or imaginary modes. Of the 212 primary rows, 119 carry a stored relaxation-not-converged flag and 161 a soft-mode flag (flags at pair level, kept from the screen). Using the per-structure records of the screen (relaxation converged to $10^{-4}$ eV/Å, fraction of imaginary modes on the H-mass mesh), the rule is right in 86 of 93 rows in which both relaxations converged (52 pairs: 45 right, 7 ties, none wrong; baseline 66), and in 35 of 36 rows in which both structures are also free of imaginary modes (28 pairs: 27 right, 1 wrong; 9 formulas; baseline 28).

	MACE	ORB	MatterSim	Pooled
Clean covalent rows (formulas)	188 (55)	163 (45)	170 (51)	521 (63)
All-row slope, OLS	0.112	0.117	0.128
95% interval	0.068 to 0.154	0.083 to 0.165	0.090 to 0.169
All-row slope, Theil-Sen	0.061	0.111	0.109
Structure-level slope, formula fixed effects	0.120	0.116	0.126
Pearson $r$ (all rows)	0.76	0.83	0.84
All-row stretch / bend and libration	+0.245 / −0.133	+0.259 / −0.142	+0.252 / −0.124
All-row compensation	0.54	0.55	0.49
Strong-regime rows in fit	97	97	99
Strong-regime total slope	0.132	0.136	0.147
Strong-regime stretch / bend and libration	0.277 / −0.145	0.308 / −0.171	0.293 / −0.146
Strong-regime compensation	0.52	0.56	0.50
Sign rule right, evaluations, 300 K	52/68	64/73	69/71	185/212
Constant-sign baseline, 300 K				131/212
Pair majority right/tie/wrong				74/0/10 of 84
Formula majority right				27 of 28
Sign rule right, 0 K	44/68	54/73	55/71	153/212
Sign rule right, 600 K	67/68	73/73	70/71	210/212
Stretch slope ratio $f$ vs Libowitzky, 2.50 to 2.80 Å	0.91	0.92	0.92
ZPE difference MAE, DFT geometries (meV per H)	1.5	0.7	1.3

This strictest subset contains no mp-985293 row and no ORB row, because ORB relaxations rarely reached $10^{-4}$ eV/Å; its all-row slopes are 0.206 (MACE, 30 rows) and 0.146 (MatterSim, 24 rows). The subset is small and selected by numerical behaviour that may correlate with chemistry, so it supports the rule’s sign but does not fix the slope. No structure was re-relaxed or followed along a soft mode for this revision.

Weak hydrogen bonds: the rule fails

When both polymorphs have mean model O-H below 0.99 Å, the rule is right in only 2 of 20 evaluations (13 pairs, 1 right and 12 wrong by majority; 7 formulas, 6 of them wrong by majority: Al(OH)₃, Be(OH)₂, Cd(OH)₂, NaOH, RbOH, Zn(OH)₂). Two readings fit these rows. All 20 have $\Delta\Delta F_{\mathrm{iso}} > 0$, so a constant positive sign is right in all of them, and 18 of the 20 have $x < 0$; the misses may reflect the orientation of these pairs rather than a reversed dependence on bond length. Moreover 17 of the 18 misses lie at $r_{\max}$ between 0.98 and 0.99 Å, within 10 mÅ of a boundary chosen after inspecting the data, and below 0.98 Å 2 of 3 rows are right. Over all weak-side O-H rows the slopes are −0.27 (MACE, 76 rows), −0.058 (MatterSim, 52) and −0.014 (ORB, 59), and 11 of the 20 scored rows are MACE rows; in the prospective test (§ 1.4) the weak side was right in 2 of 9. A reversal would be consistent with the change of sign of the net quantum effect with hydrogen-bond strength [16, 19], but these data do not establish it. For weak bonds the stretch barely changes with geometry (the Libowitzky curve spans only 170 cm$^{-1}$ from 2.80 to 3.10 Å), and bend, libration and lattice terms set the sign. Which terms dominate was not resolved here, and the evidence is 20 evaluations in 7 formulas. Without NaOH and RbOH the weak side is right in 2 of 14 evaluations (5 formulas), so the misses do not rest on the alkali hydroxides. In NaOD, path-integral simulations find that deuteration restores the weak hydrogen bonds and drives antiferroelectric order [53]; the deuterium-favoured phase there has the longer O-H, the direction of the misses here, but the screen’s four NaOH structures, one of them a theoretical database entry, were not matched to those phases, and the mechanism there is anharmonic.

The stretch term against measured spectra

For O-H donors with model O···O distances of 2.50 to 2.80 Å the harmonic Γ-point stretch frequencies follow the slope of the Libowitzky correlation [52] with slope ratios of 0.91 (MACE, 52 donors in 18 structures), 0.92 (ORB, 58 in 20) and 0.92 (MatterSim, 58 in 17), with bootstrap intervals of 0.67 to 1.07, 0.82 to 1.01 and 0.81 to 1.02 (Figure 4). The measured curve spans 1646 cm$^{-1}$ in this window. The model frequencies are harmonic and the measured ones fundamentals, and the gap between them, which matters for isotope partition functions [54], grows with hydrogen-bond strength, so the harmonic reference slope is flatter than the measured one; with an assumed growth of that gap by 150 to 400 cm$^{-1}$ over the window, the corrected ratios are 1.00 to 1.22. The assumed range is not a confidence interval. If the bend and libration term were unaffected, the fitted net slope would change sign only at ratios of about 1.8 to 2.0, the inverse of the compensation. The sign count is more sensitive than the slope. Dividing the stretch term of every

primary row by 1.1, 1.22 and 1.5 leaves 167, 150 and 123 of the 212 evaluations right, against 131 for a constant sign, so the upper end of the corrected range above would remove most of the margin over that baseline. At a fixed stretch term, scaling $\Delta\Delta F_{\mathrm{L}}$ by 0.9, 1.1, 1.2 and 1.3 gives 195, 167, 155 and 140. These uniform rescalings were performed later as descriptive sensitivities, not error estimates: an error common to both terms leaves every sign unchanged, and the covariance of the model errors in the two terms is unknown. In the one pair of § 1.4 with measured stretch and bend shifts the bend to stretch shift ratio, about 1.1 from a broad stretch band, exceeds the model median of 0.7; the mp-985293 pairs show the same sensitivity within the models. Above 2.80 Å the ratios are 0.67, 2.61 and 1.46, so ORB and MatterSim are too steep there, and the last 0.13 Å of the strong regime (to O···O of about 2.93 Å) is not covered. This is a consistency check on the stretch term, not a validation of isotope free energies or of the bend and libration term.

Temperature dependence

The isotope free energies already stored from the initial screen at 0, 300, 600 and 2000 K allow the same 212 evaluations to be rescored at fixed geometry. This rescoring is post hoc; its mechanism was then tested with four hypotheses specified before any new decomposition run. The rule is right in 153, 185 and 210 of 212 evaluations at 0, 300 and 600 K, against leave-one-formula-out baselines of 145, 131 and 114, and in 211 at 2000 K. At 0 K the rule barely beats the baseline; it is a thermal rule. At 0 K the unflagged rows are right in 153 of 170, so the misses there are again largely the crossover pairs. The mp-985293 rows score 0, 8 and 28 of 28 at 0, 300 and 600 K. All 212 rows pass the quadrature check at 0 and 300 K and 190 at 600 K (24 of the 28 mp-985293 rows); quadrature and direct $\Delta\Delta F_{\mathrm{iso}}$ have the same sign in all 212 rows at every temperature.

Of the four predeclared hypotheses three passed and one failed. The stretch term keeps the predicted sign in all 28 mp-985293 rows at 0 and 600 K, with compensation above 1 at 0 K and below 1 at 600 K (passed): these pairs cross over in temperature rather than through a change of the stretch sign. The crossover lies between 431 and 506 K for MACE, 317 and 406 K for ORB and mostly below 300 K for MatterSim (descriptive). The Teller-Redlich product rule holds at every $q$ to $2.6 \times10^{-11}$ per H in these structures (passed). At 2000 K the sign of $\Delta\Delta F_{\mathrm{iso}}$ equals the sign predicted by the change in the trace of the on-site H force constants in 211 of 212 rows (bar 202, passed); this test initially scored 124 of 212 because the scoring script read the pair orientation from the wrong file, and it was rescored after that fix with hypothesis, threshold and rows unchanged. With the stretch near 5000 K, $\hbar\omega/k_{\mathrm{B}}T$ is still about 2.5 at 2000 K, so this is a sign test, not evidence that every mode is classical. The polymorph with the shorter mean O-H has the larger H force-constant trace in all 212 rows, so at high temperature the rule reduces to the Bigeleisen-Mayer and Urey statement that deuterium concentrates where the H force constants are largest [13, 14]. The hypothesis that failed was the proposed size of the temperature trend: compensation falls with temperature in 198 of 212 rows, but by about 39% at the median between 0 and 600 K (median $C$ 0.665, 0.582 and 0.406 at 0, 300 and 600 K), and the predeclared drop of more than 40% in at least 80% of rows was met in only 98. Temperatures above about 500 K describe the harmonic isotope term at fixed geometry, not phase stability; gibbsite and many hydrates decompose below 600 K.

Exploratory prediction for diaspore and δ-AlOOH

Diaspore ($\alpha$-AlOOH, $Pnma$, mp-625152) transforms at high pressure to δ-AlOOH ($P2_{1}nm$, mp-23807), which has the stronger hydrogen bond: its O-H is longer by 24 to 32 mÅ and its O···O shorter by 61 to 92 mÅ at 0 GPa in the models. The rule therefore predicts that deuteration destabilises δ-AlOOH. At 0 GPa geometry the harmonic isotope term $\Delta\Delta F_{\mathrm{iso}}(\delta-$ diaspore) at 300 K is +0.52, +0.77 and +0.83 meV/atom (MatterSim, MACE, ORB), with $\Delta V = -0.30$ to $-0.32$ Å$^3$/atom, which gives +0.3 to +0.4 GPa. Quasi-harmonic Gibbs energies at 10, 15 and 20 GPa, with one linear fit for both isotopes, place the 300 K transition of AlOOD 1.33 GPa (MACE) and 1.13 GPa (ORB) above that of AlOOH; a temperature-resolved rebuild gives +1.33, +1.13 and +0.99 GPa (MatterSim) at 300 K and +0.94, +0.81 and +0.73 GPa at 1000 K. These quasi-harmonic runs use a $\Gamma$-centred $q$ mesh with the point density of the screen, rounded to an integer (not necessarily even) between 4 and 12, rather than the shifted mesh of the screen. The quasi-harmonic isotope term at 10, 15 and 20 GPa 2.30, 2.42 and 2.22 meV/atom (MACE) and 1.86, 2.30 and 2.50 (ORB), is about three times its 0 GPa harmonic value, which is why the pressure shift exceeds the 0 GPa estimate. Its growth with pressure is consistent with the δ-phase hydrogen bond strengthening towards symmetrization, exactly where harmonic physics is least reliable; O-H and O···O along the pressure grid were not stored and are not reported.

The models place the AlOOH transition at 8.76 (MACE, extrapolated below the 10 GPa grid point), 13.00 (ORB) and 15.14 GPa (MatterSim) at 300 K, a spread of 6.4 GPa: none of these values was compared with the experimental stability field [55] or with earlier DFT work [56], and the experimental boundary is crossed at high temperature. δ-AlOOH undergoes isotope-dependent proton disorder in the same pressure window [1, 2], which harmonic physics omits; the measured 4.4 GPa shift of that neighbouring transition shows that anharmonic isotope effects in this system are large, and the potentials overbind hydrogen bonds (Sec. V). The prediction is therefore an exploratory, conditional sign prediction outside the regime in which the rule was derived; neither its sign nor its size is validated. A logged search found no measurement or calculation of the isotope shift of this boundary. Related measured shifts point the same way, the δ-AlOOH disorder transition [1] and ε-FeOOH remaining asymmetric above the pressure at which ε-FeOOH symmetrises [57], but both are different transitions. The same calculation for goethite to ε-FeOOH gives +0.2 to +0.5 GPa with a mesh sensitivity as large as the signal, so only its sign is stated.

Tests against experiment

Isotope effects on polymorph free energies have been measured only through related quantities, the D$_2$O shifts of ice boundaries and the mineral-water D/H fractionation of boehmite and diaspore; neither is a measurement of $\Delta\Delta F_{\mathrm{iso}}$ at the model’s conditions, so none of the tests below scores the rule directly. Each test is an ingredient or a consequence (Table II).

Volume isotope effects (mixed). Quasi-harmonic volume isotope effects were predicted for known H compounds before the measured cell volumes were read. $\alpha$-AlH$_3$ has the measured sign [58] with a magnitude 1.3 to 1.5 times too large; GaOOH and $\beta$-CrOOH have the wrong sign against cross-laboratory references of limited comparability. The misses lie in the strong hydrogen-bond regime, where quasi-harmonic physics omits Ubbelohde-type proton anharmonicity [3]. The test passed where the rule has no skill and missed where it applies, so it does not test the rule; it suggests that the potentials overweight the O-H stretch contribution to volume derivatives in strongly hydrogen-bonded crystals.

Differential spectroscopy (bend sign only). The mechanism predicts that, between two polymorphs, the stretch falls and the bend rises as the hydrogen bond strengthens. For goethite to $\varepsilon$-FeOOH the highest O-H bend rises by $+290$ cm$^{-1}$ against $+195$ predicted [59]. The stretch comparator is the midpoint of the broad 2600 to 3100 cm$^{-1}$ band reported for $\varepsilon$-FeOOH, giving a nominal $-260$ cm$^{-1}$ and $-273$ predicted but anywhere from $-511$ to $-11$ cm$^{-1}$; Mashino et al. [60] read the same spectrum as two maxima near 2800 to 2900 and 3300 to 3400 cm$^{-1}$, and if the upper component were the main band the stretch sign would reverse, so the stretch sign is not determined. The goethite bands are attributed to Liu et al. [61], whose full text was not read. The diaspore to

$\delta$-AlOOH and manganite to groutite pairs were inconclusive. This is consistent with the sign of the compensating bend shift at the zone centre, not a measurement of its size; the measured bend to stretch shift ratio, about 1.1, exceeds the model median of 0.7.

Mineral D/H fractionation (ordering only). The same isotope free energies give reduced partition function ratios, so differences between measured mineral-water fractionation factors test them without a model of water. Predictions for 118 structures were recorded before the measured values were read. The measured factors are attributed in our data file to Xu and Zheng for brucite [62], Yapp for goethite [63] and Hariya and Tsutsumi for manganite [64], each reached through secondary sources and not checked against the primary papers. The sign of every scored difference (6 of 6, in 2 scored pairs) and the order brucite $>$ goethite $>$ manganite in D affinity are right in all three models, and the mineral with the longer O-H bonds takes up less deuterium. The magnitudes are not: for brucite minus goethite the models give 190 to 300 per mil against a measured 22 to 46, 5 to 10 times too large, and the predicted difference falls with temperature while the measured one rises. The order is assembled from two temperature windows, and the factor of 5 to 10 applies to brucite minus goethite; the two pairs with manganite give ratios of 1.6 and 0.7, one of them unscored. Corrections to harmonic reduced partition function ratios are discussed in Ref. [54]. Brucite (mean $O$-$H$ 0.971 Å) lies below the rule’s regime. The scoring rule was recorded together with its result.

Ice phase boundaries (failed). The $H_2O$ to $D_2O$ shifts of three ice boundaries were predicted before any measured $D_2O$ boundary was read. For Ih-II the models give −12 to −140 MPa at 230 to 250 K against a measured triple-point shift of +12 to +15 MPa; Ih-III is also wrong in sign in all three models, while II-III is right in sign and within a factor of 2 to 3. The measured values come from a secondary compilation and from Bridgman’s $D_2O$ work via a reference equation of state [65, 66], with an inferred isotope label and no stated uncertainty. Both failures involve ice Ih, modelled through the ordered ice XI proxy, 7 to 9% too dense in all three models; overbound hydrogen bonds and the omitted proton disorder are both consistent with the failure.

Prospective test (failed). After the initial screen analysis we recorded predicted signs for 111 primary evaluations of pairs never computed before (42 pairs, 20 formulas, both phases experimental), after relaxation and before any phonon run. The bar was at least 75% right, at least 15 points above the constant-sign baseline fixed in advance and a pair majority right in at least 80% of pairs. The rule got 16 of 111 right, the baseline 5. 87 rows were lost to the predeclared numerical rules (69 failed or dynamically unstable, 18 not clean), dominated by $KH_2PO_4$-type and ice polymorphs with imaginary modes, which capped the attainable accuracy at 21.6%. The surviving rows do not rescue the rule: on the 24 clean rows it was right in 16, the baseline fixed in advance in 5 and the opposite constant sign, chosen after the fact, in 19, and the pair majority was 5 right, 1 tie and 4 wrong. The failure combines numerical attrition with insufficient skill on the rows that survived. Transfer of the rule to new pairs is not shown, and the clean filter selects a population, ordered oxyhydroxides and hydroxides, close to the one on which the in-sample rule was built. Descriptively and after part of the data had been seen, the clean rows gave a total slope of 0.181 and a compensation of 0.44. The predeclared secondary test, at 600 K without a threshold, missed its 90% bar (21 of 120); the 21 rows it scored right are all of its clean rows.

Boehmite against diaspore (inconclusive). Among the primary pairs only boehmite against diaspore has a measured route. The models predict that D prefers boehmite, the shorter O-H phase: +42.5, +30.3 and +41.3 per mil at 380 °C. Measured mineral-water fractionation [22], read from a compilation [67], gives a boehmite minus diaspore difference of +28.8 per mil at 380 °C (boehmite at 2 kbar, diaspore at 4 kbar), the predicted sign in all three models and the sign Graham and co-workers attributed to the hydrogen-bond geometry of the two phases, but

the source states no uncertainty, and the diaspore value comes from 25% exchange, so the predeclared rule scores the case inconclusive.

Reimplementation. A second model, without access to our code, rewrote the double difference and its split from the cached force constants with its own mesh and quadrature. Of the 20 rows it computed, 9 have a reference row; on these it reproduces $\Delta\Delta F_{\mathrm{iso}}$ to $10^{-9}$ meV/atom with 9 of 9 signs, and its own mesh refinement moved the per-structure isotope term by at most about 0.04 meV/atom on the 33 structures it computed. The other 11 belong to pairs excluded as not clean. Both codes read the same force constants, so this checks implementation, not the potentials or the harmonic model.

Discussion

Known physics. The competition between stretch and bend zero-point energies in hydrogen bonds is established [16, 17, 19, 20], and the high-temperature direction, in which deuterium prefers the phase with the larger H force-constant trace, is the Bigeleisen-Mayer and Urey argument [13, 14]; at high temperature the rule’s sign agrees increasingly with this limit (211 of 212 at 2000 K). Hydrogen-bond control of mineral D/H fractionation is known qualitatively [21, 23], and isotopic polymorphism is documented in single systems [5, 6, 7, 8, 9, 10, 11, 12]. For full substitution, $\Delta\Delta F_{\mathrm{iso}}$ per H equals $k_{\mathrm{B}}T(\ln\beta_{\alpha}-\ln\beta_{\beta})$, with $\beta$ the full-substitution reduced partition function ratio per H of each polymorph, the formalism first-principles studies use for D/H fractionation between hydrous minerals and water [23, 24]; a measured dilute-D fractionation between the polymorphs differs from it by the term $R$ of Eq. (6), which is written per formula unit and which harmonicity alone does not remove. At 300 K, 1 meV per H corresponds to 38.7 per mil in 1000 $\ln\alpha$. In $KH_2PO_4$-type crystals deuteration lengthens the $O\cdots O$ distance and raises the ordering temperature, which Tikhonova and co-workers trace to the geometric isotope effect [26, 28], an anharmonic change of geometry that our model omits because D is evaluated at the H geometry. The static analogue of the paraelectric phase has imaginary modes, so no harmonic $\Delta\Delta F_{\mathrm{iso}}$ exists for that comparison, and the $KH_2PO_4$-type pairs of the prospective test were lost to imaginary modes (Sec. V). NaOH lies on the weak side, where the rule fails (Sec. IV D). This physics sets the direction of each term and the high-temperature limit. It does not by itself give the sign of the net effect between two polymorphs at room temperature, where the bend and libration term cancels about half of the stretch term and in one structure (mp-985293) all of it; Section II accounts for the sign of the stretch slope and, to within 4 to 34%, its size within the same models. Isotope effects on polymorph stability have been measured or computed for single systems: the $H_2O$ and $D_2O$ ice boundaries [65, 66], the ice XI to Ih transition [68], and D/H fractionation between the AlOOH polymorphs boehmite and diaspore and in the zoisite and epidote system [22]. The net quantum effect on a single hydrogen bond is known to change sign with bond strength [16, 19].

What this work adds. In the sources we searched, we found no study that compares polymorph pairs across many inorganic O-H$\cdots$O compounds, fits the polymorph isotope free-energy difference against the O-H length difference, or locates in O-H length where the sign of that comparison stops following the O-H difference; our boundary at 0.99 Å was set after inspecting the data. This work formulates that direction as a general rule and tests it across compounds. It states the rule’s condition, moderate to strong O-H$\cdots$O hydrogen bonds ($r_{\max} \ge0.99$ Å), and gives the net slope of the polymorph isotope free-energy difference against O-H length across 220 inorganic polymorph pairs, with clustered uncertainties, split into a stretch term and a bend and libration term that cancels about half of it. It shows where the models agree with the rule (effects larger than the model spread) and where it fails (weak hydrogen bonds), identifies one structure in which bend and libration outweigh the stretch at room temperature but not above a crossover temperature, and turns the rule into a conditional pressure prediction. Within these models, the same pair can change its isotope preference with temperature at fixed geometry, and a bond-length rule cannot stand for the whole vibrational free energy. A longer O-H bond goes with a stronger hydrogen bond and, in strongly bonded systems, deuteration lengthens the O$\cdots$O distance (the Ubbelohde effect [3, 27]); in chemical terms the rule says that deuteration favours the polymorph whose hydrogen bonds are weaker, the direction seen in the pentachlorophenol and 4-methylpyridine complex [6, 8]. Because the O$\cdots$O sign performs as well as the O-H sign on the same rows, the practical form is a geometric screen. The O$\cdots$O form rests on the inverse correlation between mean O-H and mean O$\cdots$O lengths across structures ($r = -0.83$ to $-0.85$ in the three models) [26], not on the isotope-induced change of either distance within a phase; it was scored only on the screen’s pairs and is not established for proton-disordered phases or bent hydrogen bonds. With database geometries instead of relaxed ones the sign agreement drops (69% against 86% in the initial screen analysis), so relaxed geometries are needed.

Limitations. (1) The free energies are harmonic or quasi-harmonic. Anharmonicity of the O-H stretch, tunnelling and proton disorder are omitted, and D is evaluated at the H geometry, so geometric isotope effects on O-H and O$\cdots$O are not included; the failed or inconclusive tests (mineral fractionation magnitudes, ice phase boundaries, the δ-AlOOH disorder window) point to exactly these effects in strongly hydrogen-bonded systems, where the rule is meant to apply. In 161 of the 212 primary rows a structure has bracketed soft or imaginary modes; their treatment ($|\omega|$ in the quadrature, omission in the direct sum) gives a reproducible number but not a physically defined harmonic partition function for those modes. (2) The three potentials share PBE-level training data and overbind hydrogen bonds (ice XI 7 to 9% too dense), so a shared bias in O-H lengths would move both axes of Figure 1 together; the per-structure benchmark errors (0.5 to 1.1 meV per H at 300 K at the DFT geometries, 1.0 to 1.4 at the model geometries) are comparable to the typical signal, and no DFT double difference was computed. Universal MLIPs soften the potential-energy surface [69]; a uniform softening rescales $\Delta\Delta F_{\mathrm{iso}}$ without changing its sign, but a different softening of bends and librations than of the stretch moves the count by 18 to 45 evaluations for a 10 to 30% change in that term (Sec. IV E). (3) The regime boundary, regime variable, eligibility cut and O-H definition were set after inspecting the data; nested re-choice covers only the threshold value, and the multiverse shows that the count depends on these choices, although the rule beats the baseline under all of them. (4) Evaluations share structures and models; five formulas supply most rows, and the formula-clustered intervals are wide. (5) The clean criterion is narrower than physical stability; the strictest subset (converged, no imaginary modes) is small, and the prospective test on new pairs failed, through numerical losses and on the clean rows. (6) No polymorph isotope free energy has been measured directly; the ice-boundary and boehmite-diaspore fractionation tests use related quantities.

The direct experiment. D/H partitioning between two coexisting polymorphs constrains $\Delta\Delta G_{\mathrm{iso}}$. For polymorphs $a$ and $b$ of one composition with $n_{\mathrm{H}}$ hydrogens per formula unit, the double difference per formula unit is

where $\alpha_{a-b}$ is the trace (dilute D) fractionation factor and $R$ accounts for the difference between trace and full substitution; with a full-substitution factor the relation is exact with $R = 0$. To second order $R \approx(n_{\mathrm{H}} k_{\mathrm{B}} T/2)[\operatorname{Var}_{b}(l) - \operatorname{Var}_{a}(l)] + N_b - N_a$, with $l_i = \ln\beta_i$ per H site and $N_p$ the non-additivity of full substitution. $R$ vanishes when the H sites are equivalent and independent with composition-independent substitution free energies and ideal mixing; harmonicity alone does not establish this, because the rule of the geometric mean needs its own assumptions [70]. Both model phases dias­pore and δ-AlOOH have one inequivalent H site, so the variance term vanishes, but the non-additivity term has not been bounded, and a trace measurement constrains the predicted full-substitution difference only once it is. Equation (6) is the solid-state analogue of the Bigeleisen relation between vapour-pressure isotope effects and fractionation [14], and the boundary shift $\Delta P_{\mathrm{iso}} = P_{\mathrm{D}}^t - P_{\mathrm{H}}^t$ follows as $\Delta P_{\mathrm{iso}} \approx-\Delta\Delta G_{\mathrm{iso}}/\Delta V$ to linear order, with $\Delta V = V_b - V_a$ per formula unit. As a conditional conversion with $R = 0$, the quasi-harmonic isotope term at each model’s transition pressure corresponds to 1000 ln$\alpha$ (diaspore − δ) = 265 to 363 per mil at 300 K and 66 to 90 per mil at 1000 K (results/paper_stats.json); these are not fractionation predictions, and the 5 to 10 times overestimate of harmonic mineral fractionation in the mineral comparison suggests that measured values would be smaller. Coexisting polymorphs exist only along the boundary at high pressure, and D/H exchange between solids needs a fluid and high temperature, where $\alpha$ is smaller; the feasibility of resolving such values with ion-probe D/H measurements was not assessed here. The second test is a matched, reversed-bracket determination of the diaspore to $\delta$-AlOOH boundary for AlOOH and AlOOD in the same press, pressure medium, calibrant and temperature path, with each boundary approached from both phase assemblages. A shift of opposite sign contradicts the prediction; a shift not resolvable from zero is a null result. On the computational side, path-integral mass integration [29, 30] or self-consistent anharmonic phonons [71] on a few primary pairs, and DFT double differences on a stratified subset, are the controls this work lacks.

Conclusion

We formulated and tested a rule for the effect of deuteration on polymorph stability in hydrogen-bonded crystals: $\Delta\Delta F_{\mathrm{iso}} \simeq S_{\mathrm{str}}(1-C)\Delta\bar{r}$, a stretch term set by the softening of the O-H stretch with bond length, partly cancelled by bend and libration, with the sign of $\Delta\bar{r}$ while $C < 1$. In harmonic calculations with three universal ML potentials over 220 inorganic polymorph pairs, the H/D free-energy difference between polymorphs rises with the change in mean X-H bond length, 0.11 to 0.13 meV per H per mÅ across models with formula-clustered intervals of about 0.07 to 0.17, a stretch term about half cancelled by bend and libration. In crystals with moderate to strong hydrogen bonds deuteration favours the polymorph with the shorter O-H bonds in 74 of 84 pairs and 27 of 28 formulas (185 of 212 evaluations against 131 for a constant sign); the models agree with the rule in nearly every evaluation where the effect exceeds the model spread; it has little skill near balance, fails for weak hydrogen bonds and fails in one Al(OH)$_3$ structure where bend and libration outweigh the stretch at room temperature ($C > 1$ in two of the three models) and the sign crosses over at higher temperature, as Eq. (1) allows. It did not transfer in a prospective test: 16 of 111 committed rows were right, and on the 24 clean rows 16 against 19 for the opposite constant sign, with pair majority 5 right, 1 tie and 4 wrong. The one-mode harmonic expression accounts for the sign and, to within 4 to 34%, the size of the stretch slope within the same models, but not the rise in skill from 0 to 300 K. Exploratory quasi-harmonic calculations put the diaspore to $\delta$-AlOOH transition of AlOOD about 1 GPa above AlOOH, outside the regime in which the rule was derived. Anharmonic and DFT controls on a subset of pairs, D/H partitioning between coexisting polymorphs and matched AlOOH and AlOOD boundary runs would test the result.

Data and code availability

Code, data, prediction files and analysis records are at https://github.com/canerganis/isoshift [40]; this manuscript corresponds to the commit that adds it. An earlier version of the code and data is archived as v1.0.0 at Zenodo, doi:10.5281/zenodo.23578412 [72]. The headline counts are recomputed with a CPU in seconds from the pair table with python -m src.isoshift_event.rule_cv, and the revision statistics with src/isoshift_event/paper_stats.py. The temperature, prospective and boehmite-diaspore analyses were committed on 22 October 2026, before the DOI was added to the repository; the analytic estimate of Sec. II and the revision statistics were committed later and are not part of v1.0.0. The predictions, their commit identifiers and the order in which each analysis was fixed are recorded in the repository. Force constants and raw phonon outputs are not in the repository. Derived Materials Project data are CC BY 4.0; licence of the redistributed JARVIS-DFT structures was not verified.

Acknowledgments

This work started at the OSI Hackathon 2026, and the author thanks its organisers. The author also thanks Valency Systems for access to the Valency Bond literature service and the marimo team for the notebook framework used for the results application.

AI-use statement

The author posed the question, directed the work and approved each step, including what was tested and what is claimed, and is responsible for the content. AI agents (Claude Opus 5.5, Claude Sonnet and Haiku, GPT-6 Astra, Gemini 4 Argon) carried out the calculations, literature queries and drafting under the author’s direction; every formula was derived independently by two models and checked numerically, and agents that collected measured values did not see the predictions. All agent logs are in the repository.

References