<p>Topological constraints such as the Bernal-Fowler ice rules govern atomic arrangements in proton-disordered crystals. Machine learning force fields (MLFFs) provide a computationally efficient route to exploring such disordered energy landscapes using on-the-fly sampling without explicit topological supervision. Whether such a data-driven approach can autonomously recover ice-rule constraints remains an open question. Here, we show that a data-efficient GAP-SOAP force field trained on 212 DFT calculations achieves a energy RMSE of 19 meV/H<sub>2</sub>O on Ice VII configurations and reproduces phonon dispersions and infrared spectra in good agreement with DFT. A representative trained model relaxes 8 × 8 × 8 supercells (3072 atoms) with randomized proton orientations into structures consistent with ice-rule constraints, recovering a body-centred cubic oxygen lattice at 30 GPa. Under compression, the model captures the qualitative features of the Ice VII → Ice X transformation and remains stable over the 30–90 GPa range for the tested configurations. While extended training does not guarantee monotonic improvement in all properties, the model enables efficient large-scale molecular dynamics, reaching approximately one frame per second on a single CPU core. These results demonstrate that local environment descriptors can encode key aspects of hydrogen-bond topology, providing a data-efficient route to structural relaxation in disordered molecular crystals.</p>
