Thermodynamics and dielectric response of BaTiO₃ by data-driven modeling

Abstract

Modeling ferroelectric materials from first principles is one of the successes of density-functional theory and the driver of much development effort, requiring an accurate description of the electronic processes and the thermodynamic equilibrium that drive the spontaneous symmetry breaking and the emergence of macroscopic polarization. We demonstrate the development and application of an integrated machine learning model that describes on the same footing structural, energetic, and functional properties of barium titanate (BaTiO₃), a prototypical ferroelectric. The model uses ab initio calculations as a reference and achieves accurate yet inexpensive predictions of energy and polarization on time and length scales that are not accessible to direct ab initio modeling. These predictions allow us to assess the microscopic mechanism of the ferroelectric transition. The presence of an order-disorder transition for the Ti off-centered states is the main driver of the ferroelectric transition, even though the coupling between symmetry breaking and cell distortions determines the presence of intermediate, partly-ordered phases. Moreover, we thoroughly probe the static and dynamical behavior of BaTiO₃ across its phase diagram without the need to introduce a coarse-grained description of the ferroelectric transition. Finally, we apply the polarization model to calculate the dielectric response properties of the material in a full ab initio manner, again reproducing the correct qualitative experimental behavior.

Introduction

Ferroelectric materials possess a spontaneous electric polarization that can be switched with an external electric field. The discovery of ferroelectricity in barium titanate (BaTiO₃), the prototypical ferroelectric perovskite, changed the general understanding and perception of ferroelectrics due in large part to its relatively simple crystal structure¹. At low temperatures, BaTiO₃ is rhombohedral with polarization along the 〈111〉 direction; at higher temperatures, it undergoes three phase transitions, first to an orthorhombic phase with the polarization along the 〈110〉 direction at 183 K, then to a tetragonal phase with the polarization along 〈100〉 at 278 K, and finally, at 393 K, to a cubic, paraelectric phase². It has been long understood that the spontaneous polarization is a result of the titanium atom off-centering within the enclosing oxygen octahedron, but the detailed microscopic nature of the ferroelectric transition has been the subject of intense, ongoing research with a variety of experimental and theoretical techniques. The ferroelectric transitions were first described with a displacive model in which the Ti-displacements are driven by a transverse phonon instability³. Almost concurrently, an order-disorder model was proposed to explain the origin of the Ti-displacements along any one of the eight local 〈111〉 directions in the cubic phase, as driven by the pseudo Jahn-Teller effect ⁴, showing how these displacements order at lower temperatures in different ferroelectric phases^5,6. These models capture some of the phenomena experimentally observed in characterizing BaTiO₃, such as phonon softening at the transition temperatures^7,8 —consistent with the displacive model—and diffuse X-ray scattering in all phases except the rhombohedral one^9,10,11—consistent with the order-disorder model—leading also to approaches combining the two models^12,13,14. In this context, simulations – especially from first principles—can offer a precious microscopic understanding of the nature of the phase transitions.

A computer simulation of the ferroelectric phase transition of any given material requires three key ingredients: first, a model of the potential energy surface (PES) that describes the energetic response to atomic and structural distortions, second, the free-energy surface (FES) sampled at the relevant, finite-temperature thermodynamic conditions, and third, the polarization of individual configurations that determines, through averaging over samples, the macroscopic polarization.

Density-functional theory (DFT) calculations have long been used to explore the PES of BaTiO₃ as well as the soft phonons and their strong dependence on pressure^15,16,17,18. Further DFT investigations have found that Ti-displacements along local 〈111〉 directions can result in dynamically stable structures^19,20,21. The phase transitions and rhombohedral-orthorhombic-tetragonal-cubic (R-O-T-C) phase sequence of BaTiO₃ has been extensively studied and reproduced using effective Hamiltonians solved using both Monte-Carlo^22,23 and molecular dynamics (MD)^24,25,26; furthermore, similar studies have been carried out on other perovskite systems²⁷, including solid solutions²⁸. Despite their successes, effective models rely on the choice of an explicit parametrization of the Hamiltonian; therefore, in order to confidently make first-principles-accurate predictions of the thermodynamics, it is desirable to use an unbiased, agnostic approach without any prior assumption in the form of the PES.

To this aim, we introduce an integrated machine learning (ML) framework allowing us to carry out MD without the need to compromise on simulation size and time scales. This framework, based on a combination of an interatomic ML potential and a vector ML model for the polarization, is used to simultaneously predict the total energy, atomic forces, and polarization of a ferroelectric material in order to explore its complex, temperature-dependent phase diagram as well as to predict its functional properties. This approach allows us to compute macroscopic observables—chemical potentials and dielectric susceptibilities, specifically—with an accuracy equivalent to that of the level of theory of the underlying DFT calculations, but at a much smaller computational cost. Moreover, it is applicable with only minor changes to any perovskite or even any other type of ferroelectric material, including 2D ferroelectrics²⁹. Although we do not reach a quantitative agreement with the experimental R-O-T-C transition temperatures, we demonstrate that this limitation in accuracy stems from the DFT reference itself and not the approximation made in modeling the potential energy surface. Thus, we foresee clear, systematic pathways to improving the model potential, with only slight modifications of the ML methodology. Specifically, the generality of the framework and the relatively small size of the training dataset make it possible to improve the model accuracy by computing the reference structures with more advanced functionals such as Hubbard-corrected DFT^30,31, meta-GGAs,³² and hybrids³³.

The key advance underlying this work is a unified ML framework combining an interatomic potential based on the SOAP-GAP method³⁴ and a microscopic polarization model based on the symmetry-adapted Gaussian process regression (SA-GPR) method³⁵. The use of ML for materials modeling has gained considerable momentum in the past decade^{34,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56}. Specifically, the prediction of finite-temperature properties of materials as the ones we focus on in this paper relies on the construction of ML potential energy surfaces based on a set of reference structures computed with ab initio methods^53,57,58. Such potentials allow the simulation of molecules and complex solids with almost the same accuracy as the reference method used to generate the dataset. In this way, it is possible to investigate the meso- and macroscopic properties of materials at a considerably reduced computational effort compared to direct ab initio simulations. Notable successes of the machine learning potentials approach include the study of bulk and interfacial properties of metallic alloys from cryogenic temperatures up to the melting point⁴⁵; finite-temperature modeling of binary systems with variable concentration, such as GaAs⁴⁶; accurate calculations on the relative stability of competing phases of various compounds, such as sodium⁵⁸, carbon⁵⁹, water⁴⁷, iron⁴⁹, and silicon⁵⁰; as well as MD studies of polycrystalline phase-change materials⁵⁷ and hybrid perovskites⁶⁰.

Two important developments have enabled the great success of machine learning in condensed matter and chemical physics. First, appropriate regression schemes—such as kernel methods, typified by Gaussian approximation potentials (GAP)³⁴; neural networks (e.g., of the Behler–Parrinello type³⁶ or more recent graph convolutional approaches^43,61); or non-kernel-based linear fitting schemes (with appropriate representations^62,63,64,65)—have been designed and specialized for atomistic systems. The key to nearly all of these methods is the decomposition of a global (extensive) physical observable of the system into local contributions, each written as a function of the neighborhood of individual atoms. Note that this decomposition carries with it an implicit assumption of the locality of the potential energy surface, thus neglecting the effect of long-range electrostatic and dispersion forces. Several extensions have previously been proposed to include such forces within existing ML frameworks^66,67,68,69, but for the purpose of this work, we use an explicitly short-range model with an appropriately chosen cutoff.

The second advancement is the construction of suitable, physically-motivated representations to predict the target properties of interest^{42,70,71,72,73}. In particular, the representation of an atomic configuration should reflect all the physical symmetries of the target property. The framework built around the Smooth Overlap of Atomic Positions (SOAP) descriptor⁷⁴ and its covariant counterparts^35,73, which we call the atom-centered density-correlation framework, is well suited to the task of integrated machine learning modeling of multiple properties since it allows us to treat these properties within the same unified mathematical framework. We provide further details on the mathematical framework, as well as the construction of the unified ML model, including the definition of a polarization-derived collective variable, in the Methods section. Looking forward, the flexibility and extensibility of this framework will also allow us in future studies to include long-range interactions in a natural and general way, using a recent approach called LODE^75,76. This will allow us to address some of the observed disagreements with DFT benchmarks in the prediction of the phonon spectra, which likely derive from the neglect of long-range forces (see the subsection “Phonon dispersions” for additional details).

The modeling of multiple properties within a single ML framework is gaining importance as a way to extract richer information from simulations than the PES alone can provide. Such models combine the extensive, accurate, finite-temperature thermodynamic sampling afforded by an ML potential, as in a series of the previous works^{45,47,49,50,57,58,60,77,78}, with the expressiveness and utility of an ML property model. Particularly relevant are the studies using a potential energy surface combined with a dipole-moment model for studying the infrared spectra of isolated molecules^79,80. To date, such combined models have not yet been applied to ferroelectric materials; one important difficulty for ML modeling is the multi-valued character of the polarization in the condensed phase (although see ref. ⁸¹ and ref. ⁸² for applications in liquid water, where this difficulty is much less severe). We describe a method of overcoming this difficulty in a systematic and generalizable way in the Methods section and in Supplementary Note 3. In the following study, we show how a combined modeling study can advance the field of ferroelectrics by providing a rich array of experimentally relevant properties from one unified mathematical framework.

Results

Summary

In this section, we summarize the main results obtained via the integrated ML model described in the Methods section. In the subsection “Structural transition in BaTiO₃”, we recover the well-known sequence of R-O-T-C phases and highlight the key role of the ML-predicted polarization vector in distinguishing each of the phases. In the subsection “The microscopic mechanism of the ferroelectric phase transition”, we find that the Ti off-centering is the driving mechanism of ferroelectricity and not just a result of cell distorsions.

In the subsection “Thermodynamics of BaTiO₃”, we provide an explicit calculation of the phase transition temperatures. Finally, in the subsection “Dielectric response of BaTiO₃”, we compute its temperature- and frequency-dependent dielectric properties and compare these with the experiment.

Structural transitions in BaTiO₃

The detection of the ferroelectric transitions of BaTiO₃ in MD simulations is challenging due to the small lattice distortions and free-energy differences that differentiate the phases. To overcome these challenges, one has to choose a sufficiently large cell so as to make the transitions clearly visible while allowing for a well-converged statistical sampling of configurations across each of the coexistence regions. As a qualitative indicator of the phase transitions, we track the MD time evolution of the cell parameters (a, b, c, α, β, γ) and the histograms of the ML-predicted unit-cell polarization components for each phase. The unit-cell polarization correlates strongly with the magnitude and direction of the Ti-displacement (see SI).

Figure 1 shows that the high-temperature cubic phase consists of a collection of local minima arranged at the vertices of a cube, as proposed by the eight-site model^5,13,83,84. The presence of large thermal fluctuations, as compared to the energy barriers separating the minima, allows for diffusion of the polarization vector across the minima over a timescale of the order of a few ps in the MD trajectories. This makes the eight local Ti minima equally probable, yielding 〈P〉 = 0.

Fig. 1: The four phases of BaTiO₃.

a Time evolution of the lattice vector lengths (a, b, and c) and angles (α, β, γ) over fully flexible MD simulations at (from left to right) 250, 100, 50 and 15 K; b 3D histograms of the unit-cell polarization predicted by the SA-GPR polarization model, computed across the same simulations; c projections of those 3D histograms onto the three principal Cartesian planes (X-Y, X-Z, and Y-Z). The light gray boxes mark the range ±2e for easier comparison between the phases. The structure of the high-frequency (dark blue) regions characterizes the different phases of BaTiO₃, while the phase transitions are marked by a clear symmetry-breaking pattern that restricts the cell dipole to visit only a subset of the available off-centered sites.

A reduction of the temperature results in a structural first-order phase transition, both in agreement with ref. ²⁶ and with previous experimental⁸⁵ and theoretical works^22,86 showing a divergence of the latent heat at the Curie point. Such transition is characterized by a clear breaking of reflection symmetry of the cell dipoles across one Cartesian (X-Y, Y-Z, or X-Z) plane. The polarization vector can only visit four of the eight available cubic sites, marking the onset of the tetragonal phase. Any further decrease in the temperature further reduces the symmetry of the polarization histograms by successive freezing of the polarization components along a specific axis. At 50 K, the polarization densities show a single and broad minimum corresponding to an orthorhombic state, while at 15 K finally, the system is completely frozen in one minimum corresponding to the rhombohedral state. Note that each of the 6 tetragonal, 12 orthorhombic, and 8 rhombohedral minima that are trivially equivalent by symmetry can be reached depending on the initial configuration. These states are associated with different distortions of the lattice vectors (that are not symmetry invariants) and one can observe occasional transitions between them. Thus we can infer from Fig. 1 that the phenomenology of the ferroelectric-paraelectric transition agrees well with the eight-site model. Under this model, the breakdown of ferroelectricity is characterized by thermal fluctuations across cubic off-centered sites that restore, on average, the centrosymmetry of the Ti-displacements.

Furthermore, we see that the action of a large isotropic pressure of 30 GPa in the cubic state fully restores the isotropy of the polarization densities (see Supplementary Fig. 4) and generates a paraelectric BaTiO₃ phase down to 0 K. This is consistent with the experimentally observed loss of ferroelectricity in BaTiO₃ at high pressures⁸⁷, as well as with the flattening of the calculated PES ^15,16, the disappearance of all unstable phonon modes of cubic BaTiO₃, and the isotropic Ti-displacement distribution observed in Car-Parrinello MD simulations of cubic BaTiO₃ under pressure²⁰. This evidence strengthens the hypothesis that the fluctuations of the unit-cell polarization between preferential orientations act as a microscopic precursor of the macroscopic ferroelectricity of the material.

Moreover, the emergence of a ferroelectric state is facilitated by the presence of spatial and directional correlations across the structures, which have been proposed to explain X-ray diffraction results^10,88 and directly observed in ref. ⁸⁹, with further evidence from first-principles calculations^14,90,91.

Figure 2 shows the extent and directionality of the spatial correlations (specifically, component-wise Pearson correlation coefficients; see Supplementary Note 5 for the exact expression) of the unit-cell dipoles correlated against a central reference cell in the cubic phase at 250 K. The correlations are not only large and slowly decaying—they extend well up to the edges of the 5 × 5 × 5 supercell—but they are highly directional, with the slow decay taking place along the direction of the unit-cell dipole.

Fig. 2: Spatial correlations of the unit-cell dipoles computed on a 5 × 5 × 5 supercell simulated at 250 K.

a Slice at constant Z of the correlation of the X-component of the dipole of the unit-cell centered on one Ti-atom (circled in red) with the X-component of the dipole of every other Ti-centered cell. b 3D view of the dipole correlations for the lower half of the same 5 × 5 × 5 supercell; the full 3 × 3 correlation tensor of each dipole with the central Ti-atom (red sphere) is shown as an ellipsoid, where the elongation along the Cartesian axes shows the presence of highly directional, needle-like correlations⁹¹.

It has long been assumed that these correlations arise from a combination of an Ising-like nearest-neighbor interaction with a long-range dipole–dipole interaction, as typified e.g., in the model Hamiltonian of ref. ²². Indeed, the authors of that study observed that the Coulomb interaction was critical for reproducing the ferroelectric ground state in their model—when it was turned off, the ground state became antiferroelectric. However, our simulations show long-range correlations and a ferroelectric ground state even though the energy model itself is explicitly short-ranged—that is, the energy of an atom i is only sensitive to changes within a short-ranged, local environment of 5.5 Å of the atom (see Methods, subsection “Training the ML model for BaTiO₃”)—making the correlations observed in our simulations an emergent phenomenon, not relying on the existence of any explicit long-ranged interaction. This range is sufficient to capture short-range correlations between two neighboring Ti atoms, whose average distance in a typical MD run fluctuates about ≈4.0 Å. Furthermore, we observe these correlations even in the disordered, cubic phase, in contrast e.g., to ref. ⁹¹, where the strongest correlations were observed only in the ordered phase (albeit in a different material, and where the phase transition was triggered not by temperature but by including quantum nuclear effects).

Thus, our understanding of the nature of ferroelectricity in BaTiO₃ must take into account these emergent, long-range correlations. We will see, for instance, how they give rise to spontaneous ferroelectric states—even in the absence of lattice distortions—in the following subsection. On the other hand, these correlations also hamper the statistical and simulation-cell size convergence of various quantities computed from statistical averages of the total dipole moment, as will be discussed in the subsection “Dielectric response of BaTiO₃”.

The microscopic mechanism of the ferroelectric transition

A fundamental question that arises in relation to the structural transitions observed in perovskite ferroelectrics concerns the driving mechanism of the transitions. We have seen in the previous subsection that the presence of the ferroelectric behavior is accompanied by the onset of a macroscopic polarization, mostly driven by ordered displacements of the transition metal atoms and cell deformation. This indicates that reducing the temperature makes it energetically favorable to develop polarized states, even at the expense of an internal strain introduced by the subsequent cell deformation.

One might, however, question whether it is the Ti off-centering or the cell deformation that drives this sequence of transitions or whether these distinct mechanisms are equally present and competing. To this end, ambient-pressure simulations of a 4 × 4 × 4 cell over a wide range of temperatures (between 20 and 250 K) were carried out in a restricted cubic geometry. The geometric constraint inhibits the structural distortions, making it possible to investigate whether the displacements and ferroelectric states are still observable.

Two-dimensional histograms of the Ti-displacements for a series of representative trajectories at 25, 80, 130, and 250 K are provided in Fig. 3. The lowest-temperature trajectory is equivalent to a fast-freezing experiment, where the Ti atoms relax to the closest potential energy minimum. We still observe off-centered states according to the eight-site model picture, but no transitions between neighboring cubic sites take place due to the negligible thermal fluctuations.

Fig. 3: 2D displacement maps, mapped onto the three principal Cartesian planes (X-Y, X-Z, and Y-Z) in a series of NpT 4 × 4 × 4 simulations with restricted cubic geometry.

The 20 K trajectory shows fast fast-freezing of the Ti degrees of freedom to the nearest cubic sites (depending on the initial configuration) in the first few ps of trajectory. In this specific case, only states corresponding to positive displacements along the y-axis are sampled across the 4 × 4 × 4 supercell. At 80 K the presence of one single maximum in the 3D density marks the onset of a clear “rhombohedral-like” ferroelectric state (with a polarization parallel to the [−1, +1, +1] axis), showing how the GAP favors aligned displacements even in the presence of geometric constraints on the supercell. Simulations at higher temperatures (250 K and beyond) restore instead the eight-site structure of the displacement density.

By slightly increasing the temperature, the thermal fluctuations are still smaller than the energy barrier between neighboring cubic sites but are sufficient to induce rare jumps between them. The Ti atoms consequently freeze in a local minimum, but notably, they all collectively jump to a single off-centered state after a small transient of the order of a few ps. The system stays trapped in this state for the whole simulation time.

This gives us evidence that the GAP inherently favors correlated ferroelectric displacements, despite the short-ranged description of the interactions (in an Ising-like fashion) and indicates that a low-temperature ferroelectric state arises in fully flexible simulations and in experiments as a consequence of a dipolar ordering, which is suppressed at high-temperatures by thermal fluctuations.

At higher temperatures, thermal fluctuations enable transitions between cubic sites with a rate that increases with the temperature. Due to the restored cubic centrosymmetry of the displacements, states with a net polarization are no longer observed, provided sufficiently long MD runs are performed. We find this same sequence of states in NVT cubic simulations (see Supplementary Note 7), showing that constraining the volume does not affect the qualitative picture of Fig. 3. Additionally, we note that the intermediate tetragonal and orthorhombic states do not occur in these simulations, as opposed to the fully flexible ones.

In conclusion, the presence of a dipolar ordering is responsible for the emergence of a low-temperature ferroelectric state, as shown in Fig. 4 and in agreement with ref. ⁸⁸. At the same time, the absence of cell distortions considerably affects the shape of the FES, as the intermediate tetragonal and orthorhombic states do not occur with a fixed cell. This effect has also been reported in ref. ²², where Monte-Carlo simulation with no homogeneous strain showed the disappearance of such phases.

Fig. 4: Schematic of the ferroelectric transition.

2D sketch of the sequence of events marking the onset of a ferroelectric state starting from a perfect cubic configuration.

Thermodynamics of BaTiO₃

A challenge in modeling phase transitions such as the ones we focus on in this paper is that they are associated with small structural distortions that are comparable with the thermal fluctuations of individual atoms. A commonly-used strategy to improve the signal-to-noise ratio is to use collective variables (CVs), such as the lattice parameters⁶⁰, which are naturally averaged over multiple atomic environments and directly reflect the macroscopic observables associated with the transition. Cell vectors, however, are not symmetry-adapted so that multiple equivalent states are mapped to different values of the CVs. What is more, as we have seen in the previous subsection, cell distortions alone do not drive the different phase transitions of BaTiO₃, making them poor order parameters to distinguish these phases (see Supplementary Note 4 for a discussion of our metadynamics simulations that use a symmetrized combination of lattice parameters).

More effective characterization of the ferroelectric ordering can be obtained by explicitly using the predicted polarization P as an order parameter and, in particular, by building descriptors that show the orientation of the cell polarization relative to the atomic distortions. In Methods, subsection “Physically-inspired order parameters”, we provide the construction of a two-component CV, namely s = (s₁, s₂), that gives us an effective low-dimensional description of the phases of BaTiO₃.

In Fig. 5, we show 2D contour lines of s across fully flexible MD runs of a BaTiO₃ 4 × 4 × 4 supercell between 10 and 250 K. These represent molecular dynamics runs where the GAP predicts the coexistence of R-O, O-T, and T-C states respectively, with comparable probability. Four distinct phases are clearly visible, showing how a polarization-derived two-component CV can easily identify the subtle differences between the four phases. The relative positions of the clusters give additional (but only qualitative) physical insights: as the C-T-O clusters are maximally distinguishable by s₁ and the C center corresponds to the one with the lowest s₁ value, the first CV is clearly related to the average polarization magnitude. This is predicted to be exactly zero for paraelectric cubic BaTiO₃ in the thermodynamic limit while positive and increasingly large for the ferroelectric tetragonal and orthorhombic phases. The CV s₁ can then be used to discriminate ferroelectric and paraelectric BaTiO₃ states. Further evidence in this respect is provided in Supplementary Fig. 12, where we show how it is possible to reconstruct free-energy profiles as a function of s₁ across the T-C transition and thus capture their finite-temperature stability.

Fig. 5: Phases of BaTiO₃ in CV space.

2D contour plots of the two-component CV s = (s₁, s₂) plane across fully flexible MD trajectories, highlighting the presence of the R-O-T-C clusters.

On the other hand, s₂ maximizes the difference between the R-O-T ferroelectric states; thus, we can relate it physically to the polarization orientation, in agreement with our observations of the subsection “Structural transitions in BaTiO₃”. Additional evidence for this interpretation is provided in Supplementary Fig. 8, where we analyze the correlation between the CVs, constructed here, with a meaningful physical observable, namely the polarization magnitude, in a series of 5 × 5 × 5 fully flexible trajectories.

Based on 2D maps as the one shown in Fig. 5, it is possible to cluster the MD trajectories and compute directly the temperature-dependent FES by calculating the relative concentration of the R-O-T-C phases in sufficiently long MD runs so that many reversible transitions can be sampled. In practice, we never find more than two phases explored at each temperature. Fully flexible MD runs of a BaTiO₃4 × 4 × 4 cell, each with a total simulation time up to 1.6 ns, between 10 and 250 K, allow us to compute the relative Gibbs free energy ${{\Delta }}{G}^{k,k^{\prime} }(T)$ of the phase pair $(k,k^{\prime} )$ at temperature T and the corresponding chemical potential difference ${{\Delta }}{\mu }^{k,k^{\prime} }(T)$ using the following equations:

$$\begin{array}{ll}{{\Delta }}{G}^{k,k^{\prime} }(T)=-N{k}_{{{{\rm{B}}}}}T\ln \frac{{\sum }_{t}{P}_{k}(t)w(t)}{{\sum }_{t}{P}_{k^{\prime} }(t)w(t)}\\ {{\Delta }}{\mu }^{k,k^{\prime} }(T)=\frac{{{\Delta }}{G}^{k,k^{\prime} }(T)}{N}\end{array}$$

(1)

where N = 320 is the number of atoms, w(t) is the weight of the t-th structure, P_k(t) is the probability that the t-th structure belongs to the phase k, and k_B is the Boltzmann constant. For the O-T and T-C transition, unbiased MD simulations are used, i.e., with w(t) = 1, for every t. In the case of the R-O transition, w(t) represents the weight of the t-th structure computed via the iterative trajectory reweighting (ITRE) technique⁹², which is used to remove the time-dependent bias in the distribution of the microstates introduced by metadynamics.

The estimates of the critical temperatures at ambient pressure are computed by linear fits of the relative chemical potentials Δμ profiles, see Equation (1) and are reported for each pair of phases in Table 1. We note that our computed temperatures differ significantly from the experimentally observed transition temperatures, 393 K (T-C), 278 K (O-T), and 183 K (R-O)². This underestimation of the critical temperatures, stemming from an underestimation of the free-energy barriers between the phases, could come, in principle, from a variety of mechanisms, particularly the neglect of long-range electrostatics (and consequently of the LO-TO splitting in the training-set structures), as well as the presence of finite-size effects that could stabilize the high-temperature disordered phases in the MD. In fact, a significant size dependence of the finite-temperature properties of another perovskite, PbTiO₃, has been reported in a recent ML-driven study by ref. ⁹³.

Table 1 Critical temperatures of the T-C, O-T, and R-O transition at ambient pressure. R² represents the coefficient of determination of best linear fits of Δμ(T) across the coexistence regions, as shown in Fig. 6.

However, previous work based on effective Hamiltonian models already pointed out this same underestimation of the critical temperatures and connected them to a shortcoming of the underlying exchange-correlation functionals^22,26, which can be compensated by rescaling the potential energy surface or introducing an artificial negative pressure.

To confirm that pressure can significantly affect the transition temperature, in Fig. 9, we investigate the sensitivity of the Curie temperature to negative pressure (p = −2 GPa) within our ML framework. We observe a shift of the Curie temperature that is increased by a significant 82 K with respect to the ambient pressure estimate via NST simulations, while a very small variation in the lattice constant of the MD supercell (0.7% at 250 K) is seen in NpT. We note that the corresponding change in the volume (2.1%) is within the variation of calculated volumes of cubic BaTiO₃ with different DFT exchange-correlation functionals²⁰. Furthermore, experimental data⁹⁴ on the elastic properties of BaTiO₃ show that the bulk modulus of BaTiO₃ is in the range of ≈200 GPa, implying that the action of relatively high pressure would result in a small change in the cell volumes while completely modifying the free-energy landscape. This effect induces the Curie temperature shift. We note in this respect that the PBEsol functional that we used to compute energies and atomic forces of the training-set structures, predicts a slightly underestimated lattice constant −4.0 Å at 400 K from Car-Parrinello MD²⁰ as compared to the experimental value of 4.012 Å of cubic BaTiO₃⁹⁵. This is a minuscule underestimation that is, however, not negligible in the calculation of these tiny free-energy barriers. Moreover, we also rule out that the main source of discrepancy with the experimental transition temperatures might arise from an incorrect prediction of thermal expansion (see Supplementary Fig. 11), as previously shown in ref. ²³, where the underestimation of the critical temperatures by an effective Hamiltonian model was related to the approximations made in the construction of the PES.

This sensitivity of the relative free energies on the equilibrium volume shows how it is possible to tune the applied pressure to obtain a better agreement with the experiment. Moreover, while the use of external pressure is a common strategy to improve the accuracy of ab initio MD—including in the recent ML-driven study by ref. ⁹³, to correct the so-called supertetragonality problem—our strategy opens up other avenues for improvement. For instance, since our ML potential is trained on a relatively small set of 1458 self-consistent energy calculations, one could systematically test more accurate and demanding electronic structure approximations^30,31,32,33, by running them on the existing dataset to improve the quantitative agreement between simulations and experiment.

So far, we have shown the capability of the GAP of both qualitatively describing the emergence of ferroelectric states in BaTiO₃ and reproducing the correct phase sequence. Seeing, however, the substantial disagreement between the critical temperature predictions of the ML model with the experiments, we shall now assess its accuracy as compared to the underlying DFT method. A compelling test in this direction is provided by the free-energy perturbation (FEP) method. From the collected MD trajectories, we extract a validation set of 50 tetragonal and cubic structures just below (170 K) and above (194 K) the Curie point and recompute their energies with self-consistent DFT calculations. This allows us to compute how the error of the GAP-predicted energies on the test set propagates to the error of the chemical potential estimate at a given temperature.

The FEP on the chemical potentials is first computed as a correction on the Gibbs free energy G^k of phase k = T, C:

$${{\Delta }}{G}_{{{{\rm{FEP}}}}}^{k}(T)=-{k}_{{{{\rm{B}}}}}T\ln \left\langle \,\exp \left(-\frac{{E}_{{{{\rm{GAP}}}}}^{k}-{E}_{{{{\rm{DFT}}}}}^{k}}{{k}_{{{{\rm{B}}}}}T}\right)\,\right\rangle$$

(2)

where 〈⋅〉 represent the average over the test set structures, ${E}_{{{{\rm{GAP}}}}}^{k}-{E}_{{{{\rm{DFT}}}}}^{k}$ the deviation between GAP and DFT total energies in phase k, and T the temperature. The FEP on the Gibbs free energies can then be translated into a correction on the chemical potential differences as follows:

$${\mu }_{{{{\rm{FEP}}}}}^{k}=\frac{{G}_{{{{\rm{GAP}}}}}^{k}+{{\Delta }}{G}_{{{{\rm{FEP}}}}}^{k}}{N},$$

(3)

where N is the number of atoms. Equation (2) represents an average of Boltzmann factors: if the energy deviations between the DFT and GAP estimates are small compared to the thermal fluctuations at temperature T for both the tetragonal and cubic phases, the correction on the corresponding chemical potential is negligible, due to the exponential factors. This propagation of errors can, however become significant or even dominant if the energy deviations are of the same order of magnitude or larger than the thermal fluctuations.

Panel (f) of Fig. 6 shows the effect of the FEP correction on the estimate of the chemical potentials for the two selected temperatures. The GAP shows good performance in the prediction of both tetragonal and cubic structures (compared to k_BT) and the FEP correction is one order of magnitude smaller than the actual prediction of Δμ and is still well within the error bars computed with the MD runs. The correction is hence negligible and no shift in the Curie point is observed, providing strong numerical evidence of the DFT accuracy of the GAP in free-energy predictions.

Fig. 6: Phase diagram of BaTiO₃ at ambient pressure.

Panels a–c show the R-O, O-T, and T-C clusters in the (s₁, s₂) plane across the coexistence regions. Each configuration generated in the MD is a point in this plane, colored according to its probability P. The latter is computed via the probabilistic analysis of molecular motifs (PAMM)¹³⁹ algorithm (see Methods, subsection “Physically-inspired order parameters” for additional details). P is smoothly increasing in the [0, 1] interval while going from R to O (panel a), then from O to T (panel b), and finally from T to C (panel c). Panels d–f show the temperature-dependent chemical potential differences across the R-O (d), O-T (e), and T-C (f) transitions. Blue error bars represent standard deviations computed across multiple MD runs, the cyan line the best linear fit across the coexistence regions, and the purple error bars the propagated errors on the critical temperatures. In addition, the orange triangles in panel f show the free-energy perturbed chemical potentials, using 50 reference tetragonal and cubic structures just below (170 K) and above (194 K) the Curie point. 6× magnified insets corresponding to these two temperatures show how the FEP-corrected chemical potentials consistently fall within the error bars due to the MD sampling, confirming the DFT accuracy of the GAP (see subsection “Thermodynamics of BaTiO₃”).

Dielectric response of BaTiO₃

Let us now turn our attention to using the polarization model developed and described in Methods, subsection “Polarization model” to compute experimentally measurable quantities. As previously mentioned in ref. ⁹⁶ and elsewhere in the literature on the modern theory of polarization^97,98, the polarization of a condensed-phase system is well defined only modulo the quantum of polarization; however, we can still compute experimentally observable quantities as changes and fluctuations in its value.

The first of these experimentally relevant quantities is the static dielectric constant, which can be computed directly from the fluctuations of the system’s total dipole⁹⁶. In the cubic phase:

$${\varepsilon }_{r}=1+\frac{\left\langle {M}^{2}\right\rangle }{3{\varepsilon }_{0}V{k}_{{{{\rm{B}}}}}T}+{\varepsilon }_{\infty }$$

(4)

where both the total dipole magnitude M and the vacuum permittivity ε₀ are expressed in SI units and the average value of the cell dipole $\left\langle {{{\bf{M}}}}\right\rangle =0$ by symmetry. The optical (electronic) dielectric constant ε_∞ from both measurements and calculations^17,99,100 is in the range of 5 to 6, which is much smaller than the typical range of ε_r we calculate for this material, so this term will be neglected in the following analysis. In any case, the analyses below are nearly or completely insensitive to such a small constant shift. For non-cubic phases, we must modify the equation to subtract off the (now nonzero) average polarization, replacing $\left\langle {M}^{2}\right\rangle$ with $\left\langle {M}_{\alpha }{M}_{\beta }\right\rangle -\left\langle {M}_{\alpha }\right\rangle \left\langle {M}_{\beta }\right\rangle$, where α and β are Cartesian components of the total dipole vector¹⁰¹. Notably, since the non-cubic phases have an anisotropic structure, the dielectric tensor ε_r,αβ will also generally be anisotropic. Indeed, experimental measurements on single-domain crystals of BaTiO₃ have shown a pronounced dielectric anisotropy, especially in the tetragonal phase^2,102.

Comparison of our results with experiments is complicated by the dramatic variation in the measured value with temperature, composition, and grain size^103,104,105. We, therefore, study the temperature dependence explicitly, as shown in Fig. 7. The calculated values for the orthorhombic, tetragonal, and cubic phases agree qualitatively with the calculations of ref. ¹⁰⁰, which used a similar computational methodology but with a shell-model potential, as well as with measurements on single-domain crystals^2,102. In the tetragonal phase, we see the expected strong anisotropy between the components parallel (ε_∥) and perpendicular (ε_⊥) to the polarization axis—as we can already see from Fig. 1, the polarization fluctuations in the tetragonal phase are strongly suppressed along the polarization direction, which matches the much smaller value of ε_∥ seen here. In the orthorhombic phase, the experimental measurements are averages over different domains and thus do not show the same pattern of anisotropy—namely, the splitting into three separate principal components—seen here, but this splitting is present in ref. ¹⁰⁰.

Fig. 7: Static dielectric constant of BaTiO₃.

a Temperature dependence of the static dielectric constant computed for a 4 × 4 × 4 supercell across the orthorhombic, tetragonal, and cubic phases. The phase transition temperatures computed from the chemical potential (see Table 1 are shown as vertical dashed lines. b Supercell-size comparison of dielectric constants computed in the cubic phase, including fits to the Curie–Weiss law Equation (5)—note the log scale on the y-axis. The effective Curie temperature T_c,ε predicted by each fit is shown by the shaded vertical spans and the error bars towards the bottom left of the figure. Error bars show one standard deviation and are computed as described in Supplementary Note 6.

In the cubic phase, the expected temperature dependence follows a version of the Curie–Weiss law¹⁰³:

$${\varepsilon }_{r}=\frac{C}{T-{T}_{{{{\rm{c}}}},\varepsilon }}+{\varepsilon }_{T\to \infty }$$

(5)

where T_c,ε is the (dielectric) Curie temperature, which should—in the limit of infinite system size and statistical sampling—agree with the tetragonal-cubic phase transition temperature computed above, T_{c, C−T} = 182.4 K.

From the temperature dependence data in Fig. 7, we determine the best-fit parameters for the 4 × 4 × 4 cell data to be ε_T→∞ = 90 ± 18, T_c,ε = (167.4 ± 2.3) K, and C = (57100 ± 3600) K. The most important discrepancy to note here is that the Curie point predicted by this fit is still about 15 K lower than the thermodynamic phase transition temperature predicted for ambient pressure in the subsection “Thermodynamics of BaTiO₃”. This discrepancy is likely a result of finite-size effects due to the small supercell, which are known to broaden and shift critical points¹⁰⁶. The 5 × 5 × 5 fit, on the other hand, yields ε_T→∞ = 99 ± 18, T_c,ε = (171 ± 6) K, and C = (62200 ± 2500) K: the Curie temperature T_c,ε is slightly closer to the predicted phase transition temperature T_c,C−T, which is now within the 95% confidence interval of the fit parameters. However, even with the larger supercell, we still note a discrepancy from the parameters determined by fits to experimental data^103,107—namely, the Curie–Weiss constant C is under-predicted by a factor of about 2 with respect to the experiment. This difference could be due to approximations inherent in the underlying DFT functional, either directly or indirectly, due to the underestimation of the phase transition temperatures. We test this hypothesis in more detail below by investigating the negative-pressure simulations.

The equation for the static dielectric constant, Equation (4), is, in fact, only the zero-frequency limit of the whole frequency-dependent response function. We can compute the frequency-dependent susceptibility (and thus the relative dielectric constant) via linear response theory from the one-sided Fourier transform of the dipole–dipole autocorrelation function^108,109 (again for the cubic phase):

$$\begin{array}{l}\chi (\omega )=\frac{1}{3{\varepsilon }_{0}V{k}_{B}T}\left(\left\langle {M}^{2}\right\rangle \right.\\ \left.\quad\quad\,\,-\,{{{\rm{i}}}}\omega \int\nolimits_{0}^{\infty }\left\langle {{{\bf{M}}}}(\tau )\cdot {{{\bf{M}}}}(t+\tau )\right\rangle {e}^{-{{{\rm{i}}}}\omega t}{{{\rm{d}}}}\,t\right)\\ \quad\quad\,\,\,=\,({\varepsilon }_{r}-1)\left(1-{{{\rm{i}}}}\omega \int\nolimits_{0}^{\infty }{\tilde{C}}_{MM}(t){e}^{-{{{\rm{i}}}}\omega t}{{{\rm{d}}}}\,t\right)\end{array}$$

(6)

where ${\tilde{C}}_{MM}(t)=\frac{1}{\left\langle {M}^{2}\right\rangle }\left\langle {{{\bf{M}}}}(0)\cdot {{{\bf{M}}}}(t)\right\rangle$ is the normalized dipole–dipole autocorrelation function and ε_r is the static dielectric constant computed from Equation (4).

We show the frequency-dependent susceptibility for a 6 × 6 × 6 supercell trajectory at 250 K, computed using Equation (6), in Fig. 8. In general, we see the same structure as predicted for the high-temperature cubic phase by both theoretical effective-Hamiltonian MD calculations²⁵ and observed experimentally¹⁰³, namely, that of a large absorption peak corresponding to the soft-mode (TO₁) phonon frequency. Note the slight negative dip in the real dielectric constant is expected and seen in many previous observations^25,101,102. This does not imply that the real or imaginary part of the refractive index $n=\sqrt{{\epsilon }_{r}}$ is anywhere negative. It was previously proposed^110,111 that the “soft-mode” part of the absorption spectrum of BaTiO₃ could be described with a single, strongly damped harmonic oscillator of the form

$$\chi (\omega )=\frac{A{\omega }_{0}^{2}}{{\omega }_{0}^{2}-{\omega }^{2}+{{{\rm{i}}}}\gamma \omega }$$

(7)

with amplitude A, damping constant γ, and resonant frequency ω₀ (which is always larger than the actual apparent peak frequency). However, a later study¹¹² uncovered possible inadequacies of this single-oscillator model, especially in the high-frequency range (ω₀ ≈ 100 cm⁻¹), and suggested a two-oscillator model as a possible replacement, though it was not yet justified by the available experimental data.

Fig. 8: Frequency-dependent dielectric response of cubic BaTiO₃ at 250 K computed with a 6 × 6 × 6 supercell.

The imaginary part of the spectrum is fitted to a sum of two simple harmonic oscillator response functions, as in ref. ²⁵; the two separate oscillator responses are shown in thin light blue lines, while the fitted sum is shown as a thick light blue line.

More recently, ref. ²⁵ both measured high-accuracy infrared spectra and computed theoretical spectra from MD simulations of the effective Hamiltonian model of ref. ²⁸, and they found strong evidence that the spectrum indeed is best modeled by two harmonic oscillators. The computed spectrum from our model at 250 K, shown in Fig. 8, further supports this picture: we also find that the imaginary part of the spectrum could only be satisfactorily described with two oscillators, although with different parameters from those calculated in ref. ²⁵: we find one oscillator with fundamental frequency ω₁ = 86 cm⁻¹ and damping ratio γ₁/ω₁ = 3.0, and another with fundamental frequency ω₂ = 187 cm⁻¹ and damping ratio γ₂/ω₂ = 1.1. Comparing these parameters with those calculated in ref. ²⁵, we find both frequencies to be rather high, so our agreement with their results remains mostly qualitative for now.

On the one hand, the differences we observe could be due to the large oscillations and lack of resolution at high frequencies due to the limited sampling time imposed by the relatively large computational cost of our model. However, it is more likely that both these discrepancies have the same origin as the underestimation of the phase transition temperatures discussed above – either inaccuracies in the underlying DFT model or some other effect not yet accounted for. As noted in the subsection “Thermodynamics of BaTiO₃”, the phase transition temperatures can be compensated by applying negative pressure. Indeed, ref. ²⁵ associate the higher-frequency mode with short-range correlations between (mostly) neighboring unit-cell dipoles, so it is likely that this frequency shift has the same origin as the pressure effect.

To investigate this discrepancy further and to assess the effect of negative pressure on the dielectric response, we compute frequency-dependent susceptibility spectra for all the negative-pressure simulations previously run for subsection “Thermodynamics of BaTiO₃” (specifically, Fig. 9), where the material remained in the cubic phase. The spectra are also compared to those derived from ambient-pressure simulations, specifically those used to compute the temperature dependence of the dielectric constant in Fig. 7. The comparison is shown in Fig. 10. On the one hand, we see the main peak shifting towards higher frequencies as the temperature increases, as expected from previous theoretical and experimental studies^25,111. On the other hand, we also see the peak shifting towards lower frequencies when negative pressure is applied at any given temperature. While the peak frequencies for the negative-pressure simulations still do not match experimental data for the same temperatures, the shifts are in the right direction.

Fig. 9: Results of the MD simulations with an externally applied negative pressure p = −2 GPa.

Panel a shows the chemical potential differences across the T-C transition, obtained via NST fully flexible simulations and with an estimated Curie temperature of T_c = (264 ± 1) K, as opposed to (182.4 ± 0.7) K calculated at ambient pressure (see subsection “Thermodynamics of BaTiO₃”). Panel b shows instead a comparison between the time evolution of the lattice constant in NpT simulations at 250 K with negative pressure and ambient pressure. These show how the effect of a negative pressure slightly increases the average unit-cell lattice parameter (by only 0.7%) due to the large bulk modulus of BaTiO₃. This has, however, important consequences on the relative stability of T and C phases, with a shift of the Curie point of 82 K with respect to the ambient pressure estimate. Error bars have the same meaning as in Fig. 6.

Fig. 10: Temperature and pressure dependence of the imaginary part of the dielectric response spectrum, all computed in the cubic phase on a 4 × 4 × 4 supercell.

For each two pressures, spectra computed at different temperatures are displayed on the same graph, vertically offset from one another for clarity. The units of the y-axis are arbitrary, though all spectra on the same plot have the same scale. While the increase in temperature broadens the main peak and shifts it to higher frequencies, the negative pressure instead shifts the main peak to lower frequencies for a given temperature (i.e., closer to the experiment, as well as the theoretical calculations of ref. ²⁵).

Furthermore, all simulations show a small narrow peak or edge at around 340 cm⁻¹, independent of the temperature. The frequency of the mode does depend on pressure, but due to the large bulk modulus of BaTiO₃, the mode shifts very little: only about 7 cm⁻¹ under −2 GPa of pressure.

Although this peak likely represents a feature of our model and not just a simulation artifact, we do not yet have enough information to confidently identify this peak with known vibrational modes of BaTiO₃^103,110,113.

In fact, the difficulties we encounter here in reproducing the results of simpler, experimentally accurate—but empirically adjusted—models are reminiscent of the difficulties encountered previously, e.g., in ref. ⁶⁸, in applying more accurate (in the sense of reproducing the quantum PES) ML potentials that must, in turn, account for more accurate physics, such as many-body dispersion and quantum nuclear effects, in order to arrive at the right predictions for the right reasons. Rather than being a deficiency in the machine learning simulation approach, we see this as an opportunity to discover interesting physical behaviors and mechanisms that were overlooked before.

The calculations presented here are a promising first step towards using the ML-PES and polarization framework as a generally applicable tool to predict experimentally relevant response properties. This tool will be a valuable future asset for investigating new candidate ferroelectric materials or gaining more insight into the underlying behavior of existing ones.

Discussion

In this work, we introduce a modern, general ML framework to describe at once the finite-temperature and functional properties (dielectric response) of perovskite ferroelectrics and apply it specifically to model barium titanate (BaTiO₃). This framework matches the accuracy of the underlying DFT method and does not require to preselect a given effective Hamiltonian model^22,114. The simulations made possible by this framework recover the correct R-O-T-C phase ordering in fully flexible simulations and allow to investigate of the emergence of Ti off-centerings. In particular, we highlight the driving mechanism of the ferroelectric transition, showing how the presence of these off-centered displacements gives rise to a low-temperature ferroelectric phase due to a long-range dipolar ordering. Moreover, the interplay between the displacements and the cell deformations leads to the emergence of intermediate tetragonal and orthorhombic phases. We further proceed to reconstruct the thermodynamics of BaTiO₃ (see subsection “Thermodynamics of BaTiO₃”) by means of a two-component, polarization-derived CV.

Finally, we apply the ML polarization model to calculate dielectric response properties of experimental interest, including the static and frequency-dependent dielectric constants, and investigate their dependence on temperature. While we do not reach a quantitative agreement with experimental measurements for many of the properties computed here, we see several clear, systematic pathways to improving the model potential and its predictions, such as including long-range electrostatic effects, simulating larger system sizes, as well as addressing the possible deficiencies in the underlying DFT model for both energies and polarizations. We expect that such improvements will allow us to reach a quantitative agreement with the experiments. Our results obtained with negative pressure calculations and the FEP show how this discrepancy can, in fact, be traced back in part to the sensitivity of the transition temperatures to cell volume combined with the deviation of the DFT cell volume from the experimental one. This effect suggests that a more in-depth investigation of the effects of pressure—which is well known to influence the onset of ferroelectricity—could provide further insights into the deviation from experiments. A closer agreement could also be obtained by combining, as recently proposed, different DFT schemes to describe simultaneously the energy, structure, and electronic density of perovskite oxides¹¹⁵.

We also plan to make improvements to the model performance by means of the feature sparsification technique, as detailed in ref. ⁷⁰. The latter has proven to reduce the computational cost (in energies and force predictions) by a factor of 3 or 4 for realistic systems and, in combination with larger-scale parallelization techniques, will allow us to treat larger, more complex systems.

Importantly, this ML framework automates the construction of a model of the PES and the polarization and can then be used to investigate finite-temperature properties in detail and with first-principles accuracy. Since the ML-PES was made with no explicit assumption on the functional form of the underlying PES and no prior definition of the relevant degrees of freedom of the system, this strategy is generalizable to other materials to study, e.g., 2D ferroelectrics²⁹ and solid solutions with variable stoichiometries, that are known to possess different and more complex ferroelectric states. For instance, Ba_xSr_1−xTi_yZr_1−yO₃ is known to display a rich phase diagram, depending on composition, and shows both ferroelectric and relaxor ferroelectric phases¹¹⁶. Furthermore, the framework developed is easily applicable to study the role of nuclear quantum effects, for instance, in incipient ferroelectrics such as SrTiO₃ and KTaO₃, where quantum fluctuations appear to suppress the ferroelectric state^91,117. Further extensions of this framework include the investigation of the role of a finite electric field in the MD and its effect on polarization. This will allow us to simulate, for instance, hysteresis loops, which are key to measure the energy storage of ferroelectric devices.

In conclusion, we have shown how a comprehensive, data-driven modeling framework for a perovskite ferroelectric material, based on DFT reference data, can capture the mechanisms of the ferroelectric transition, as well as make predictions of thermodynamic and functional properties with first-principles accuracy. The work opens the door for a new avenue of fruitful research into the understanding and characterization of known ferroelectric materials, as well as the discovery and design of new candidate compounds with improved industrially relevant properties.

Methods

In this section, we first summarize the construction and properties of the symmetry-adapted features used to train the ML models; a more thorough discussion of this family of features and an introduction to the notation we use here is given in Section 3 of ref. ⁷³. With these features defined, we detail how the potential energy surface and the polarization models are constructed. Turning our attention to the specifics of modeling BaTiO₃, we report the training and validation of our ML model. Furthermore, we develop physically-inspired order parameters, which we use to characterize and interpret our results from subsection “Thermodynamics of BaTiO₃”. Finally, we report the computational details of the ML-MD simulations.

Symmetry-adapted features

To construct the family of features that are relevant for this paper, we make use of the atom-centered density-correlation framework¹¹⁸. The starting point is the definition of a set of features, namely 〈anlm∣A; ρ_i〉, from an expansion of the atomic density for an environment i of structure A, as in Equation (31) of ref. ¹¹⁸. The different indices in the bra identify the chemical species (a), radial function (n), and angular momentum (l, m), the latter being especially important to track the symmetry of the features.

Symmetry-adapted descriptors can be obtained as a symmetrized average (referred to by an overline decoration) of the tensor product of ν sets of expansion coefficients, resulting in density-correlation features $\langle q| \overline{{\rho }_{i}^{\otimes \nu };\lambda \mu }\rangle$. While the generic index q only enumerates the features, the other indices encode the physical meaning of these descriptors. There are two fundamental parameters: (a) the body-order exponent ν, which indicates that the features describe the relative position of ν neighbors of the central atom and (b) the λ, μ coefficients, which determine how the descriptor transforms under rotations—namely as spherical harmonics ${Y}_{\lambda }^{\mu }$. This framework allows us to build features that are not only invariant to rotations but also explicitly covariant (more generally called equivariant) features of any tensor order. Such equivariant features were first introduced by ref. ¹¹⁹, for vector features, and in ref. ³⁵ for tensors of arbitrary order. Equivariant features are now gaining considerable popularity, especially for graph convolutional neural networks to predict scalar and tensor properties^{61,120,121,122,123}. In this work we only deal with spherical invariants or SOAP descriptors⁷⁴—corresponding to λ = 0, μ = 0 and λ = 1, μ = (−1, 0, +1) features, representing spherical equivariants of order 1. For instance, SOAP power spectrum features, which are invariant under rotations, are obtained from the contraction of two sets of coefficients (ν = 2):

$$\langle an;a^{\prime} n^{\prime} ;l| A;\overline{{\rho }_{i}^{\otimes 2};00}\rangle \equiv \mathop{\sum}\limits_{m}\langle anlm| A;{\rho }_{i}\rangle \langle A;{\rho }_{i}| a^{\prime} n^{\prime} lm\rangle .$$

(8)

These features can thus be written as $\langle q| A;\overline{{\rho }_{i}^{\otimes 2}}\rangle$.

Similarly, the simplest example of equivariant features only encodes information on the radial distribution of neighbors. They are equivalent to the density coefficients themselves:

$$\langle q\,=\,(an)| A;\overline{{\rho }_{i}^{\otimes 1};\lambda \mu }\rangle =\langle an\lambda \mu | A;{\rho }_{i}\rangle ,$$

(9)

An extension of this construction allows one to build symmetry-adapted tensors of arbitrary rank and body order⁶⁵.

Given that, in order to learn dipole moments and polarizations, we only need the special case of vector-valued features, we find it convenient to exploit the relationship between real-valued spherical harmonics of order λ = 1 and the Cartesian coordinates α = (x, y, z) to define Cartesian equivariants

$$\langle q| A;\overline{{\rho }_{i}^{\otimes \nu };\alpha \,=\,(x,y,z)}\rangle \equiv \langle q| A;\overline{{\rho }_{i}^{\otimes \nu };\lambda \,=\,1\,\mu \,=\,(1,-1,0)}\rangle ,$$

(10)

The Cartesian equivariants of Equation (10) now explicitly transform as a 3-vector under rotations:

$$\langle q| \hat{R}A;\overline{{\rho }_{i}^{\otimes \nu };\alpha }\rangle =\mathop{\sum}\limits_{\alpha ^{\prime} }{R}_{\alpha \alpha ^{\prime} }\langle q| A;\overline{{\rho }_{i}^{\nu };\alpha ^{\prime} }\rangle .$$

(11)

$\hat{R}A$ indicates an arbitrary rotation of a structure A, while ${R}_{\alpha \alpha ^{\prime} }$ is its representation as a 3 × 3 Cartesian matrix. We use this family of features to model the polarization of a BaTiO₃ structure and to build an order parameter to distinguish the R-O-T-C phases (see Results, subsection “Thermodynamics of BaTiO₃”). We refer the reader to refs. ^70,73 and the documentation of librascal¹²⁴ for implementation details.

Potential energy surface

A Gaussian approximation potential (GAP) is constructed by linear regression of energies E and atomic force components ${\left\{{{{{\bf{f}}}}}_{r}\right\}}_{r = 1}^{N}$, where N is the number of atoms, in the space of the kernels of these descriptors, representing the degree of correlation between the structures.

In order to control the computational cost of the calculation of energies and forces, we also construct a sparse set of representative atomic environments J that are used to define a basis of kernels k(⋅, J_j) in order to approximate the structure-energy relation. This is discussed further in subsection “Training the ML model for BaTiO₃”.

We write the target properties as a sum of kernel contributions:

$$\begin{array}{ll}E=\mathop{\sum}\limits_{i\in A}\mathop{\sum}\limits_{j\in J}{b}_{j}k({A}_{i},{J}_{j})\\ {{{{\bf{f}}}}}_{r}=-{\nabla }_{r}E,\end{array}$$

(12)

where the kernel is built as a function of a set of atom-centered invariant features 〈q∣A_i〉, the index j runs over all environments J_j in the sparse set J and b_j are the weights on each sparse environment to be determined via ridge regression. Here we use SOAP power spectrum features, $\langle q| A;\overline{{\rho }_{i}^{\otimes 2}}\rangle$, and we compute the kernel between atomic environments as a scalar product raised to an integer power $k({A}_{i},{A}_{i^{\prime} }^{\prime})={({\sum }_{q}\langle {A}_{i}| q\rangle \langle q| {A}_{i^{\prime} }^{\prime}\rangle )}^{\zeta }$, using ζ = 4 here, to introduce nonlinear behavior.

Polarization model

Besides this potential energy model, we construct a fully flexible, conformationally sensitive dipole moment surface for the material by employing the symmetry-adapted Gaussian process regression (SA-GPR) framework³⁵, previously benchmarked in the context of molecules in ref. ¹²⁵ and proven to extend to the condensed phase in ref. ⁸¹. Even though the cell polarization (or dipole) is not uniquely defined in periodic boundary conditions^97,98, we can still make a model for only a single branch of this polarization manifold with suitable pre-processing of the training data, detailed in Supplementary Note 3. This branch choice is essentially equivalent to fixing the polarization to be a single continuous function whose linearization about P = 0 is the product of Born effective charges and displacement from some non-polar reference structure, in the spirit of ref. ¹²⁶.

As with existing SA-GPR approaches, the total dipole of the cell is decomposed into vector-valued atomic contributions. In analogy to Equation (12), we express the total dipole M and polarization P of a structure A as:

$${{{\bf{M}}}}(A)=\mathop{\sum}\limits_{i\in A}\mathop{\sum}\limits_{j\in J^{\prime} }{b}_{j}{{{\bf{k}}}}({A}_{i},{J}_{j}),$$

(13)

$${{{\bf{P}}}}(A)=\frac{{{{\bf{M}}}}(A)}{{V}_{A}}.$$

(14)

Our model works with total dipoles rather than polarizations as only the former are size extensive. A key advantage of this model is that we represent the dipole moment as a sum of atom-centered contributions (effectively, “partial dipoles”, in analogy to partial charges), giving us a spatially resolved, atomistic picture of how the different parts of the system contribute to the total polarization. Note that in contrast to the model described in ref. ¹²⁵, we do not define an additional partial-charge model, since such a model would depend on the choice of the unit-cell and be incompatible with the modern theory of polarization. The only situation in which we use nonzero partial charges is in the linearized effective-charge model used to shift the training-set polarizations to the same branch; these effective charges are not used in the production model. As already remarked in ref. ⁹⁶ and later in ref. ¹²⁵, this information can give us a much deeper insight into the physics of the system than predicting the total dipole alone. In this study, we use this information to define Ti-centered unit-cell dipoles by an appropriate sum of atomic partial dipoles. The dipole of the Ti-atom is added to the dipoles from neighboring O and Ba atoms, with the neighboring contributions weighted (by 1/2 for O and 1/8 for Ba) so that the sum of the unit-cell dipoles is still equal to the total cell dipole. These unit-cell dipoles were used to make Figs. 1 and 2.

The model is trained on the same set of structures as the potential energy surface from subsection “Potential energy surface” but uses different training data and, generally, a different sparse set $J^{\prime}$ each with a different set of weights {b}. These weights take the form of 3-component vectors, corresponding to the kernel $k({A}_{i},{A}_{i^{\prime} }^{\prime})$, which is now a rank-2 Cartesian tensor (i.e., a 3 × 3 matrix) for any pair of environments. This kernel is computed, as in the scalar case, as an inner product of symmetry-adapted features $\langle q| A;\overline{{\rho }_{i}^{\otimes 2};\alpha }\rangle$

$${k}_{\alpha \alpha ^{\prime} }({A}_{i},{A}_{i^{\prime} }^{\prime})=\mathop{\sum}\limits_{q}\langle A;\overline{{\rho }_{i}^{\otimes 2};\alpha }| q\rangle \langle q| A^{\prime} ;\overline{{\rho }_{i}^{\otimes 2};\alpha ^{\prime} }\rangle .$$

(15)

Training the ML model for BaTiO₃

As pointed out in the subsection “Potential energy surface”, constructing a GAP model requires defining a representative set of environments to control the computational cost in evaluating energies and atomic forces of structures. The representative environments should be ideally as diverse as possible so as to provide a good extrapolation across all the phases of interest. Specifically, for our case study of BaTiO₃, we use Farthest-Point Sampling (FPS) to select a total of 250 environments centered around barium and titanium atoms and 500 around oxygen atoms from the initial training dataset obtained via DFT optimizations (additional details are given below).

A second crucial parameter is the radial cutoff in the neighbor density ρ_i(x), defined in subsection “Symmetry-adapted features”. This defines the size of the atomic environment, centered around atom i. Choosing large cutoff radii means including more neighbors in the density expansion and allows, in general, a more accurate representation of the environment. This happens, however, at the expense of increasing the computational complexity. For the purpose of constructing a GAP for BaTiO₃, we choose a radial cutoff of 5.5 Å around each center which is larger than the average separation of the first nearest Ti neighbors (≈4.0 Å). This cutoff allows us to capture the short-ranged Ti–Ti interactions that ultimately result in long-range emergent dipole correlations, a distinctive feature of polarized states in BaTiO₃, as seen in Results, subsection “Structural transitions in BaTiO₃”.

The training dataset is constructed in an iterative fashion, which also means it can be systematically extended. Energies and forces are calculated using DFT as implemented in Quantum ESPRESSO^127,128 with the PBEsol¹²⁹ functional, and managed with AiiDA^130,131,132; further details can be found in Supplementary Note 2. An initial training set of N₀ = 518 cubic structures (obtained from DFT optimizations with the PBEsol functional) is used to train a preliminary GAP. Molecular dynamics simulations with i-PI¹³³ are then performed in all the R-O-T-C geometries and for a total simulation time of up to 500 ps. Among all uncorrelated structures thus generated with MD—the correlation being computed via the time-dependent autocorrelation function of the total energy—only the most diverse according to their SOAP descriptors are then selected via FPS and recomputed with DFT self-consistent calculations. These are then used to extend the training dataset and refit the GAP, thus restarting the loop and obtaining an increasingly accurate description of the PES. The final dataset built with this procedure has a total of 1458 structures, with an adequate sampling of all the phases of BaTiO₃. Specifically, on top of the initial training set of 518 structures, we added 100 structures coming from the first round of replica-exchange molecular dynamics (REMD) simulations in the NVT ensemble and additional 840 structures coming from a sampling of each of the R-O-T-C phases (210 per phase) in the second round of NpT REMD calculations.

The learning curve of the GAP, trained on a total of 1200 training structures, is shown in Fig. 11a, with 258 randomly selected structures used as a validation set. The root mean square error (RMSE) decreases significantly with an increasing number of training points and the final accuracy of the potential in energy estimations is about 6 meV per formula unit (f.u.). This level of accuracy is sufficient to capture several interesting features of the physics of BaTiO₃, including the structural R-O-T-C phases, the presence of needle-like correlations even in the high-temperature paraelectric phase, and to enable predictions of the free-energy surface, that have the same degree of accuracy as the underlying DFT method (see Results, subsection “Thermodynamics of BaTiO₃”).

Fig. 11: Learning curves of the ML models.

a The GAP model for energies and atomic force components, trained on a total of 1200 training structures with 258 randomly selected structures used as the validation set and b the SA-GPR polarization model, trained on a total of 640 structures with 200 validation structures.

The polarization model, in contrast to the GAP, is trained only on the set of N_T,pol = 840 structures sampled from the NpT REMD calculations described above, with 210 structures coming from each of the four phases. A total of 200 randomly selected structures are withheld for testing; the largest model has therefore been trained with N_max,pol = 640 structures. The learning curve of the polarization model, shown in Fig. 11b, shows good performance; the largest model (N = 640 structures) achieves an accuracy of 3% of the intrinsic variation of the total dipoles in the training set, corresponding to an RMSE of 0.013 a₀ per atom, or 0.07 a₀ per unit-cell—which is still small compared to the scale of unit-cell polarizations shown, for example, in Fig. 1.

Phonon dispersions

A crucial test to evaluate the performance of the GAP is to compute phonon spectra and the corresponding density of states (DOS) and compare them with the DFT phonon spectra. In Fig. 12, we directly compare the outcome on a 4 × 4 × 4q-mesh, taking two representative structures as reference: the 5-atom cubic structure and the rhombohedral ground state, optimized via variable-cell DFT calculations. The calculations were carried out via the finite difference method using the atomic simulation environment¹³⁴ (ASE) for the GAP calculations and phonopy¹³⁵ in conjunction with Quantum ESPRESSO for the DFT calculations. Since no explicit correction for the long-range electrostatics was explicitly taken into account in constructing the ML model, we compare the GAP predictions with the DFT calculations without such contributions. We stress, however, that this contribution due to long-range electrostatic interactions should be included to recover, e.g., the correct LO and TO mode splitting at Γ and to stabilize the TA mode of the rhombohedral structure along the T-Γ and Γ-F paths (see the Supplementary Fig. 10 for the DFT dispersion with LO-TO splitting). It has been shown in the work of ref. ¹³⁶ that short-ranged potentials in polar materials can capture the correct phonon dispersions if the appropriate long-range dielectric model is subtracted before fitting the short-ranged potential and then added back analytically—in analogy to what is done to Fourier interpolate phonon dispersions¹³⁷. We also show, in Supplementary Fig. 10, the full phonon spectra once these dielectric contributions are considered. The spectra show an overall good agreement, especially for the low-frequency acoustic modes, with the most apparent discrepancies occurring for the highest LO mode. These discrepancies are likely to be caused by two main effects: (a) the training set construction and (b) the locality of the GAP. First, we recall that the interatomic potential is only trained on 2 × 2 × 2 structures so that long-wavelength modes that correspond to the periodicity of a 4 × 4 × 4 cell lie in the extrapolative regime of the potential. Second, the GAP is only sensitive to atomic displacements within the chosen radial cutoff, so phonon modes with a small momentum q, and thus involving long-wavelength excitations outside this radial cutoff, are not guaranteed to be well reproduced. These effects are likely the root of the disagreement between modes that lie along the Γ-X and Γ-L paths, like $(0,0,\frac{1}{4})$. Additional studies in this direction to investigate the role of the long-range electrostatic contribution on top of the GAP will shed light on this discrepancy and likely offer a better agreement with the reference DFT calculations. Furthermore, the inclusion of the LO-TO splitting will allow us to perform a finite-temperature study of the phonon dispersion across the T-C transition, to be compared with a recent study by ref. ¹³⁸. As we have seen, however, long-range electrostatic contributions are not essential to model the thermodynamics and phase transitions of BaTiO₃.

Fig. 12: Phonon dispersion and density of states (DOS) of the stable (ω > 0) phonons for BaTiO₃ calculated using finite differences with a 4 × 4 × 4 q-mesh for.

a The 5-atom rhombohedral ground state and b the 5-atom cubic structure (high-symmetry k-point labels from ref. ¹⁴⁶). We compare the GAP predictions with the DFT calculations without long-range electrostatic contributions, i.e., without LO-TO splitting. The vertical dashed lines indicate the explicitly calculated q-points.

Validation with local dipole rotations

As a further test, we evaluate the accuracy of the GAP by modeling some of the distortions associated with the ferroelectric transition. In particular, the states associated with the presence of off-centered Ti atoms relative to the O cage and the energy barrier separating them is key. As we have discussed in Results, subsection “Structural transitions in BaTiO₃”, the long-range ordering of these displacements is the fundamental driver of ferroelectricity in BaTiO₃.

To test the performance of the GAP in reproducing these states, we construct two paths, representing a local dipole rotation, across the phase space of a 2 × 2 × 2 cubic supercell with a lattice parameter of 8 Å. We start with the DFT-optimized structure with all Ti displaced by 0.082 Å along the 〈111〉 direction resulting in local dipoles, as depicted by the arrows in panel a of Fig. 13. This is a rhombohedral structure—spacegroup R3m (160)—with Ba and Ti occupying the 1a position (z_Ba = −0.0004 and z_Ti = 0.51116), and O occupying the 3b position (x_O = 0.48823, z_O = −0.01872). For reference, the cubic structure with no dipole moment would have z_Ba = 0, z_Ti = 0.5, x_O = 0.5, and z_O = 0, resulting in a cubic structure with spacegroup Pm$\bar{3}$m (221). One dipole, depicted in cyan, is then rotated about the barycenter of the enclosing oxygen octahedron to align with $\langle \bar{1}\bar{1}\bar{1}\rangle$ while keeping the magnitude of the Ti-displacement constant and all other atoms fixed. The two paths, shown as the insets in panel b, have the same endpoints but visit different vertices of the cube centered at the barycenter of the octahedron with the $\langle \bar{1}\bar{1}\bar{1}\rangle$ and 〈111〉 displacements defining a diagonal.

Fig. 13: Validation of the GAP.

Panel a shows the starting configuration of the two selected paths: the Ti dipoles are aligned along the 〈111〉 state and marked by black arrows. They are kept frozen along the two paths shown as insets in panel (b). The dipole that is rotated is instead marked by a cyan arrow. Panel b shows a comparison between GAP and DFT predictions for the two paths that connect the 〈111〉 and $\langle \bar{1}\bar{1}\bar{1}\rangle$ states.

Physically, these paths represent the energy cost due to a relative rotation of one local dipole starting from a perfect ferroelectric state. A comparison between the GAP and the DFT energy variations across these paths (see panel b of Fig. 13) shows that the GAP correctly reproduces the energy profile and favors states that correspond to aligned Ti-displacements, a feature that we have also seen in low-temperature MD simulations (see Results, subsection “The microscopic mechanism of the ferroelectric transition”). From a quantitative perspective, the GAP overestimates the energy barriers by some nonnegligible, but still reasonable, 18% for both paths. We stress, however, that these paths lie within the extrapolative regime of the potential, as they are constructed artificially and no MD simulation visits configurations that are close to them, except for the starting, completely ordered structure that is visited at low temperature (see Results, subsection “Structural transitions in BaTiO₃”).

Physically-inspired order parameters

As mentioned in Results, subsection “Thermodynamics of BaTiO₃”, the construction of a CV that can effectively distinguish the structural phases of BaTiO₃ is key for the prediction of its phase diagram. In this section, we provide the construction of a two-component CV, namely s = (s₁, s₂), by explicitly using the predicted polarization P as an order parameter. As we shall see, we will build a set of invariant descriptors that correspond, for each structure, to a scalar product of vectors. These are constructed using the equivariant features $\langle an| A;\overline{{\rho }_{i}^{\otimes \nu };\alpha }\rangle$ defined in Equation (10), averaged over Ti-centered environments. Physically, they will carry information about the orientation of P relative to the ’mean’ atomic distortion, which we call Q (see also Fig. 14).

Fig. 14: Sketch of the i-th Ti-center of structure A, along with its enclosing oxygen cage and the two vectors ${\tilde{{{{\bf{P}}}}}}_{i}(A)$ and ${\tilde{Q}}_{i}^{an}(A)$ associated with it.

The final CV is computed as a scalar product of these quantities after first summing over all Ti-centers.

Firstly, in order to compute the CV efficiently for long MD runs, we need to define an easy-to-compute proxy for P, which we will denote as $\tilde{{{{\bf{P}}}}}$. In practice, we find that some of the neighbor density coefficients 〈anlm∣A; ρ_i〉 introduced in the subsection “Symmetry-adapted features” correlate strongly with P (see the correlation plots in Supplementary Fig. 7). We can then define $\tilde{{{{\bf{P}}}}}$ by restricting ourselves to Ti-centered environments, as follows:

$${\tilde{P}}_{\alpha }(A)=\mathop{\sum}\limits_{i\in A,{{{\rm{Ti}}}}}{\tilde{P}}_{\alpha ,i}(A)=\mathop{\sum}\limits_{i\in A,{{{\rm{Ti}}}}}\langle a\,=\,{{{\rm{O}}}}\,n\,=\,6| \overline{{\rho }_{i}^{\otimes 1};\alpha }\rangle$$

(16)

where a = O represents the atomic species (the oxygen) onto which we project the Ti-centered density. Note that here we use the expression for the Cartesian equivariants defined in subsection “Symmetry-adapted features” so that α = (x, y, z) and $\tilde{{{{\bf{P}}}}}$ transforms like a vector under rotations. It represents, in fact a sum of vectors $\tilde{{P}_{i}}$, each assigned to one Ti-center, as shown in Fig. 14. Similarly, we average the full neighbor density coefficient $\langle an| \overline{{\rho }_{i}^{\otimes 1};\alpha }\rangle$ over all Ti-centers to obtain a measure of the mean structural deformations:

$${\tilde{Q}}_{\alpha }^{an}(A)=\mathop{\sum}\limits_{i\in A,{{{\rm{Ti}}}}}{\tilde{Q}}_{\alpha ,i}^{an}(A)=\mathop{\sum}\limits_{i\in A,{{{\rm{Ti}}}}}\langle an| \overline{{\rho }_{i}^{\otimes 1};\alpha }\rangle .$$

(17)

Finally, we compute the scalar product of $\tilde{{{{\bf{P}}}}}$ and $\tilde{Q}$ to construct a set of invariants:

$${O}^{an}(A)=\tilde{P}(A)\cdot {\tilde{Q}}^{an}(A)$$

(18)

and perform a principal component analysis (PCA) on the scalar descriptors O^an(A) to obtain two physically-motivated and symmetry-invariant order parameters. This step allows us to obtain the scalar components that mostly contribute to the observed variance of the O^an invariants across a dataset of structures. In particular, by performing a PCA analysis over the entirety of the MD trajectories as a function of all simulated temperatures, we find that the first two PCs, corresponding to s₁ and s₂ can neatly separate all four phases (see Results, subsection “Thermodynamics of BaTiO₃”).

At each temperature, we then perform a separate clustering using the probabilistic analysis of molecular motifs (PAMM)¹³⁹ algorithm, that determines a Gaussian mixture model in which each cluster corresponds to a different phase. Using the posterior probabilities associated with the mixture model (named probabilistic motif identifiers in ref. ¹³⁹), we can associate with each MD frame a smooth probability P_k(t), based on the corresponding values of the CVs (s₁(t), s₂(t)), that represents the probability that the corresponding structure at time t belongs to the cluster k = (R, O, T, C). These probabilities are then used to determine the relative stability of the different phases. The advantage of this technique, as compared to perhaps simpler methodologies, such as tracking the temperature evolution of the lattice parameters, is the fact that it is fully automatized, rotationally invariant, and makes direct use of the polarization vector—the key ingredient to physically describe the onset of ferroelectricity.

ML-MD Computational details

All the machine learning data that we have generated to investigate the physics of BaTiO₃ combines the use of molecular dynamics simulations performed with i-PI¹⁴⁰—the MD integrator—and librascal^70,124—the engine to compute the total energy, atomic force components, and stress tensor of a BaTiO₃ structure. In all cases, we choose the smallest simulation-cell size that provides converged results; this is to optimize the tradeoff between adequate sampling in time and adequate sampling in system size that is possible under a given computational budget.

In particular, the results of subsection “Structural transitions in BaTiO₃” correspond to NST fully flexible simulations of a 5 × 5 × 5 cell, i.e., with an external constant stress tensor σ = diag(p, p, p) with p = 1 atm. The full flexibility of the cell allows the system to relax the off-diagonal components of the MD computed stress tensor as the system undergoes the structural R-O-T-C phase transitions as a function of the temperature. In this case, we choose the simulation size so as to show well-separated structural minima as a function of the temperature while maintaining the simulations computationally inexpensive.

The results of Results, subsection “The microscopic mechanism of the ferroelectric transition“ correspond instead to isotropic NpT simulations of a 4 × 4 × 4 cell over a wide range of temperatures (between 20 and 250 K) with a restricted cubic geometry. This supercell size is sufficient to identify the Ti off-centering as the physical mechanism governing the emergence of ferroelectricity.

Fully flexible MD runs of a BaTiO₃ 4 × 4 × 4 cell with a total simulation time up to 1.6 ns between 10 and 250 K are performed for quantitative estimation of the temperature-dependent free energies (see Results, subsection “Thermodynamics of BaTiO₃”). In particular, unbiased MD is used to generate trajectories across the coexistence regions of the O-T and T-C transitions (between 40 and 250 K), while well-tempered metadynamics¹⁴¹ runs across the R-O transition are needed to enable collective jumps between R and O states within times that are affordable by classical MD runs. Additional details on the metadynamics runs are given in Supplementary Note 4. The relatively small supercell size, in this case, allows both efficient sampling of the structural transitions and simulation times, on the order of nanoseconds, that are required to converge the chemical potential estimates.

The spatial correlations shown in Results, subsection “Structural transitions in BaTiO₃” are calculated on a 5 × 5 × 5 supercell trajectory of length 400 ps, while the static and frequency-dependent dielectric constant in Results, subsection “Dielectric response of BaTiO₃” were calculated on a 6 × 6 × 6 supercell trajectory of length 200 ps in order to ensure supercell-size convergence of the static value. The temperature dependence of the dielectric constant, being a more expensive calculation requiring multiple trajectories, instead used both a 4 × 4 × 4 and a 5 × 5 × 5 supercell, simulated for 250 ps each, to explicitly assess the rate of supercell-size convergence.

All the NpT/NST simulations were carried out with an isotropic/anisotropic barostat, leaving the cell volume/vectors free to equilibrate at finite temperature. Thermalization of the cell degrees of freedom is achieved by means of a generalized Langevin equation (GLE) thermostat¹⁴², while thermalization of the atomic velocity distribution is realized via stochastic velocity rescaling (SVR)¹⁴³. This combination of thermostats allows for an optimal equilibration of the system’s relevant degrees of freedom on a timescale of the order of picoseconds without significantly interfering with the dynamical properties of the system, especially the polarization vectors. The characteristic times of the barostat, the SVR thermostat, and the MD timestep are 1 ps, 10 fs, and 2 fs, respectively.