The advances of machine-learned force fields have opened up molecular dynamics (MD) simulations for compounds for which ab initio MD is too resource intensive and phenomena for which classical force fields are insufficient. Here we describe a neural-network force field parametrized to reproduce the r²SCAN potential energy landscape of HfO₂. Based on an automatic differentiable implementation of the isothermal-isobaric (𝑁𝑃𝑇) ensemble with flexible cell fluctuations, we study the phase space of HfO₂. We find excellent predictive capabilities regarding the lattice constants and experimental x-ray diffraction data. The phase transition away from monoclinic is clearly visible at a temperature around 2000 K, in agreement with available experimental data and previous calculations. Another abrupt change in lattice constants occurs around 3000 K. While the resulting lattice constants are closer to cubic, they exhibit a small tetragonal distortion, and there is no associated change in volume. We show that this high-temperature structure is in agreement with the available high-temperature diffraction data.