Ab initio calculations suggest that the robust ferroelectricity of HfO2 at reduced dimensions arises from a unique mode coupling in which an antipolar mode stabilizes a polar distortion. Based on these insights, Landau–Devonshire energy models have been proposed using these lattice modes as order parameters. However, most existing models remain limited to simplified one-dimensional descriptions because of the high computational cost of ab initio calculations and the limitations of conventional Landau polynomials. In this study, we constructed a three-dimensional Landau–Devonshire potential for HfO2 by employing tetragonal, antipolar, and polar modes as coupled order parameters using machine-learning techniques. We generated a large-scale dataset of energies over a three-dimensional structural space, with the computational cost drastically reduced through the use of machine-learning interatomic potentials. Using this dataset, we trained a neural network (NN) to learn the complex relationship between the order parameters and the energy. The trained NN captured the characteristic coupling behavior in which antipolar modes induce the polar mode. Furthermore, we extended this NN-based Landau potential into a position-dependent functional for phase-field modeling. Using this framework, we showed that, in thin films, the critical strain required for the onset of spontaneous polarization increases, while the coercive field for polarization switching decreases, as the polarization magnitude is reduced by surface effects. This study presents the on-demand construction of phase-field models using modern machine-learning techniques, thereby enabling multiscale analysis of complex ferroelectric phenomena.