curie temperature of barium titanate experiment

Dissemination of IT for the Promotion of Materials Science (DoITPoMS)

Barium titanate and phase changes

The temperature at which the spontaneous polarisation disappears is called the Curie temperature, T C .

Above 120°C, barium titanate has a cubic structure. It is therefore centro-symmetric and possesses no spontaneous dipole. With no spontaneous dipole the material behaves like a simple dielectric, such that its polarisation varies linearly with field. T C for barium titanate is 120°C.

Below 120°C, it changes to a tetragonal phase, with an accompanying movement of the atoms. The movement of Ti atoms inside the O 6 octahedra may be considered to be significantly responsible for the dipole moment:

Cooling through 120°C causes the cubic phase of barium titanate to transform to a tetragonal phase with the lengthening of the c lattice parameter (and a corresponding reduction in a and b). The dipole moment may be considered to arise primarily due to the movement of Ti atoms with respect to the O atoms in the same plane, but the movement of the other O atoms (i.e. those O atoms above and below Ti atoms) and the Ba atoms is also relevant.

The diagram below shows the BaTiO 3 structure with an O 6 octahedron surrounding the Ti atom.

The switching to a cubic structure is the reason for the polarisation spontaneously disappearing above 120°C. Barium titanate has two other phase transitions on cooling further, each of which enhances the dipole moment:

The phase which is reached after cooling to ~ 0°C from tetragonal is orthorhombic.

And then rhombohedral below -90°C:

All of these ferroelectric phases have a spontaneous polarisation based to a significant extent on movement of the Ti atom in the O 6 octahedra in the following way (using pseudo-cubic notation):

Physical Review Journals Archive

Published by the american physical society, surface space-charge layers in barium titanate, a. g. chynoweth, phys. rev. 102 , 705 – published 1 may 1956.

• Citing Articles (348)

Above the Curie temperature, pyroelectric currents can be produced in single crystals of barium titanate even though there is no electric field applied. The polarization that remains at these temperatures is ascribed to space-charge fields in the crystal. From studies of the wave forms of the pyroelectric current signals it is concluded, tentatively, that space-charge layers of up to 10 − 5 cm in thickness reside at the crystal surface and that these space charges also produce a field through the interior of the crystal. Further evidence for space-charge fields is provided by the occurrence of an associated photovoltaic effect and by asymmetrical hysteresis loops. The space-charge fields vary considerably in magnitude from crystal to crystal. In general, they can be modified by suitable heat treatment but return to their original condition when fields are applied above the Curie point.

The space charge fields apparently influence the direction in which the domains polarize when the crystal is cooled through the Curie point. They will also affect capacity measurements above the transition and will influence the actual temperature of the transition. It is quite possible that the fields play an important role in the process of domain nucleation.

• Received 23 November 1955

DOI: https://doi.org/10.1103/PhysRev.102.705

©1956 American Physical Society

Authors & Affiliations

• Bell Telephone Laboratories, Murray Hill, New Jersey

References (Subscription Required)

Vol. 102, Iss. 3 — May 1956

Access Options

• Buy Article »
• Log in with individual APS Journal Account »
• Log in with a username/password provided by your institution »
• Get access through a U.S. public or high school library »

Authorization Required

Other options.

• Buy Article »
• Find an Institution with the Article »

Download & Share

Sign up to receive regular email alerts from Physical Review Journals Archive

• Forgot your username/password?
• Create an account

Article Lookup

Paste a citation or doi, enter a citation.

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

• View all journals
• Explore content
• About the journal
• Publish with us
• Sign up for alerts
• Open access
• Published: 29 September 2022

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

• Lorenzo Gigli   ORCID: orcid.org/0000-0002-6935-657X 1 ,
• Max Veit   ORCID: orcid.org/0000-0001-7813-4015 1 ,
• Michele Kotiuga   ORCID: orcid.org/0000-0002-2592-7342 2 ,
• Giovanni Pizzi   ORCID: orcid.org/0000-0002-3583-4377 2 ,
• Nicola Marzari   ORCID: orcid.org/0000-0002-9764-0199 2 &
• Michele Ceriotti   ORCID: orcid.org/0000-0003-2571-2832 1  

npj Computational Materials volume  8 , Article number:  209 ( 2022 ) Cite this article

7909 Accesses

23 Citations

20 Altmetric

Metrics details

• Computational methods
• Ferroelectrics and multiferroics

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 3 ), 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 3 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.

Similar content being viewed by others

An automatically curated first-principles database of ferroelectrics

Two-dimensional ferroelectrics from high throughput computational screening

Experimental discovery of structure–property relationships in ferroelectric materials via active learning

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 3 ), the prototypical ferroelectric perovskite, changed the general understanding and perception of ferroelectrics due in large part to its relatively simple crystal structure 1 . At low temperatures, BaTiO 3 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 2 . 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 3 . 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 4 , 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 3 , 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 3 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 3 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 27 , including solid solutions 28 . 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 29 . 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, 32 and hybrids 33 .

The key advance underlying this work is a unified ML framework combining an interatomic potential based on the SOAP-GAP method 34 and a microscopic polarization model based on the symmetry-adapted Gaussian process regression (SA-GPR) method 35 . 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 45 ; finite-temperature modeling of binary systems with variable concentration, such as GaAs 46 ; accurate calculations on the relative stability of competing phases of various compounds, such as sodium 58 , carbon 59 , water 47 , iron 49 , and silicon 50 ; as well as MD studies of polycrystalline phase-change materials 57 and hybrid perovskites 60 .

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) 34 ; neural networks (e.g., of the Behler–Parrinello type 36 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 74 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. 81 and ref. 82 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.

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 3 ”, 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 3 ”, we provide an explicit calculation of the phase transition temperatures. Finally, in the subsection “Dielectric response of BaTiO 3 ”, we compute its temperature- and frequency-dependent dielectric properties and compare these with the experiment.

Structural transitions in BaTiO 3

The detection of the ferroelectric transitions of BaTiO 3 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.

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 ±2 e for easier comparison between the phases. The structure of the high-frequency (dark blue) regions characterizes the different phases of BaTiO 3 , 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. 26 and with previous experimental 85 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 3 phase down to 0 K. This is consistent with the experimentally observed loss of ferroelectricity in BaTiO 3 at high pressures 87 , as well as with the flattening of the calculated PES 15 , 16 , the disappearance of all unstable phonon modes of cubic BaTiO 3 , and the isotropic Ti-displacement distribution observed in Car-Parrinello MD simulations of cubic BaTiO 3 under pressure 20 . 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. 89 , 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.

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 91 .

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. 22 . 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 3 ”)—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. 91 , 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 3 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 3 ”.

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.

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. 88 . 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. 22 , where Monte-Carlo simulation with no homogeneous strain showed the disappearance of such phases.

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

Thermodynamics of BaTiO 3

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 60 , 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 3 , 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 1 , s 2 ), that gives us an effective low-dimensional description of the phases of BaTiO 3 .

In Fig. 5 , we show 2D contour lines of s across fully flexible MD runs of a BaTiO 3 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 1 and the C center corresponds to the one with the lowest s 1 value, the first CV is clearly related to the average polarization magnitude. This is predicted to be exactly zero for paraelectric cubic BaTiO 3 in the thermodynamic limit while positive and increasingly large for the ferroelectric tetragonal and orthorhombic phases. The CV s 1 can then be used to discriminate ferroelectric and paraelectric BaTiO 3 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 1 across the T-C transition and thus capture their finite-temperature stability.

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

On the other hand, s 2 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 3 ”. 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 3 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:

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 92 , 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) 2 . 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 3 , has been reported in a recent ML-driven study by ref. 93 .

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 3 with different DFT exchange-correlation functionals 20 . Furthermore, experimental data 94 on the elastic properties of BaTiO 3 show that the bulk modulus of BaTiO 3 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 20 as compared to the experimental value of 4.012 Å of cubic BaTiO 3 95 . 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. 23 , 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. 93 , 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 3 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:

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:

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 B T ) 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.

Panels a – c show the R-O, O-T, and T-C clusters in the ( s 1 , s 2 ) 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) 139 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 3 ”).

Dielectric response of BaTiO 3

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. 96 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 96 . In the cubic phase:

where both the total dipole magnitude M and the vacuum permittivity ε 0 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 101 . Notably, since the non-cubic phases have an anisotropic structure, the dielectric t ensor ε r , α β will also generally be anisotropic. Indeed, experimental measurements on single-domain crystals of BaTiO 3 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. 100 , 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. 100 .

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 103 :

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 3 ”. This discrepancy is likely a result of finite-size effects due to the small supercell, which are known to broaden and shift critical points 106 . 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):

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 25 and observed experimentally 103 , namely, that of a large absorption peak corresponding to the soft-mode (TO 1 ) 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 3 could be described with a single, strongly damped harmonic oscillator of the form

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

The imaginary part of the spectrum is fitted to a sum of two simple harmonic oscillator response functions, as in ref. 25 ; 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. 25 both measured high-accuracy infrared spectra and computed theoretical spectra from MD simulations of the effective Hamiltonian model of ref. 28 , 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. 25 : we find one oscillator with fundamental frequency ω 1  = 86 cm −1 and damping ratio γ 1 / ω 1  = 3.0, and another with fundamental frequency ω 2  = 187 cm −1 and damping ratio γ 2 / ω 2  = 1.1. Comparing these parameters with those calculated in ref. 25 , 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 3 ”, the phase transition temperatures can be compensated by applying negative pressure. Indeed, ref. 25 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 3 ” (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.

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 3 ”). 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 3 . 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 .

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. 25 ).

Furthermore, all simulations show a small narrow peak or edge at around 340 cm −1 , independent of the temperature. The frequency of the mode does depend on pressure, but due to the large bulk modulus of BaTiO 3 , the mode shifts very little: only about 7 cm −1 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 3 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. 68 , 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.

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 3 ). 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 3 (see subsection “Thermodynamics of BaTiO 3 ”) 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 115 .

We also plan to make improvements to the model performance by means of the feature sparsification technique, as detailed in ref. 70 . 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 29 and solid solutions with variable stoichiometries, that are known to possess different and more complex ferroelectric states. For instance, Ba x Sr 1− x Ti y Zr 1− y O 3 is known to display a rich phase diagram, depending on composition, and shows both ferroelectric and relaxor ferroelectric phases 116 . Furthermore, the framework developed is easily applicable to study the role of nuclear quantum effects, for instance, in incipient ferroelectrics such as SrTiO 3 and KTaO 3 , 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.

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. 73 . 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 3 , 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 3 ”. 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 118 . The starting point is the definition of a set of features, namely 〈 a n l m ∣ A ; ρ i 〉, from an expansion of the atomic density for an environment i of structure A , as in Equation (31) of ref. 118 . 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. 119 , for vector features, and in ref. 35 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 74 —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):

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:

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

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

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

\(\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 3 structure and to build an order parameter to distinguish the R-O-T-C phases (see Results, subsection “Thermodynamics of BaTiO 3 ”). We refer the reader to refs. 70 , 73 and the documentation of librascal 124 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 3 ”.

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

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 35 , previously benchmarked in the context of molecules in ref. 125 and proven to extend to the condensed phase in ref. 81 . 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. 126 .

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:

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. 125 , 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. 96 and later in ref. 125 , 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\)

Training the ML model for BaTiO 3

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 3 , 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 3 , 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 3 , as seen in Results, subsection “Structural transitions in BaTiO 3 ”.

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 129 functional, and managed with AiiDA 130 , 131 , 132 ; further details can be found in Supplementary Note 2 . An initial training set of N 0  = 518 cubic structures (obtained from DFT optimizations with the PBEsol functional) is used to train a preliminary GAP. Molecular dynamics simulations with i-PI 133 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 3 . 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. 11 a, 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 3 , 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 3 ”).

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. 11 b, 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 0 per atom, or 0.07  a 0 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 × 4 q -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 134 (ASE) for the GAP calculations and phonopy 135 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. 136 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 137 . 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. 138 . As we have seen, however, long-range electrostatic contributions are not essential to model the thermodynamics and phase transitions of BaTiO 3 .

a The 5-atom rhombohedral ground state and b the 5-atom cubic structure (high-symmetry k -point labels from ref. 146 ). 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 3 ”, the long-range ordering of these displacements is the fundamental driver of ferroelectricity in BaTiO 3 .

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.

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 3 ”).

Physically-inspired order parameters

As mentioned in Results, subsection “Thermodynamics of BaTiO 3 ”, the construction of a CV that can effectively distinguish the structural phases of BaTiO 3 is key for the prediction of its phase diagram. In this section, we provide the construction of a two-component CV, namely s  = ( s 1 , s 2 ), 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 ).

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 〈 a n l m ∣ 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:

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:

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

and perform a principal component analysis (PCA) on the scalar descriptors O a n ( 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 a n 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 1 and s 2 can neatly separate all four phases (see Results, subsection “Thermodynamics of BaTiO 3 ”).

At each temperature, we then perform a separate clustering using the probabilistic analysis of molecular motifs (PAMM) 139 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. 139 ), we can associate with each MD frame a smooth probability P k ( t ), based on the corresponding values of the CVs ( s 1 ( t ), s 2 ( 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 3 combines the use of molecular dynamics simulations performed with i-PI 140 —the MD integrator—and librascal 70 , 124 —the engine to compute the total energy, atomic force components, and stress tensor of a BaTiO 3 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 3 ” 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 3 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 3 ”). 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 141 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 3 ” 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 3 ” 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 142 , while thermalization of the atomic velocity distribution is realized via stochastic velocity rescaling (SVR) 143 . 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.

Data availability

All numerical data supporting the results of this paper and allowing to reproduce the results are openly available on the Materials Cloud Archive 144 .

Code availability

In order to generate the data needed for this paper, we made use of the librascal 124 , i-PI 133 , and TenSOAP 145 codes. These codes are all publicly available on Github. Some additional scripts necessary for data analysis and processing beyond that provided in these codes is provided along with the research data in the Materials Cloud Archive 144 .

Jona, F. & Shriane, G. Ferroelectric Crystals (Dover, 1962).

Merz, W. J. The electric and optical behavior of BaTiO 3 single-domain crystals. Phys. Rev. 76 , 1221–1225 (1949).

Article   CAS   Google Scholar  

Cochran, W. Crystal stability and the theory of ferroelectricity. Adv. Phys. 9 , 387–423 (1960).

Bersuker, I. B. Pseudo-Jahn-Teller effect-a two-state paradigm in formation, deformation, and transformation of molecular systems and solids. Chem. Rev. 113 , 1351–1390 (2013).

Bersuker, I. B. On the origin of ferroelectricity in perovskite-type crystals. Phys. Lett. 20 , 589–590 (1966).

Chaves, A. S., Barreto, F. C. S., Nogueira, R. A. & Zéks, B. Thermodynamics of an eight-site order-disorder model for ferroelectrics. Phys. Rev. B 13 , 207–212 (1976).

Yamada, Y., Shirane, G. & Linz, A. Study of critical fluctuations in BaTiO 3 by neutron scattering. Phys. Rev. 177 , 848–857 (1969).

Vogt, H., Sanjurjo, J. A. & Rossbroich, G. Soft-mode spectroscopy in cubic BaTiO 3 by hyper-Raman scattering. Phys. Rev. B 26 , 5904–5910 (1982).

Comès, R., Lambert, M. & Guinier, A. The chain structure of BaTiO 3 and KNbO 3 . Solid State Commun. 6 , 715–719 (1968).

Article   Google Scholar  

Comès, R., Lambert, M. & Guinier, A. Désordre linéaire dans les cristaux (cas du silicium, du quartz, et des pérovskites ferroélectriques). Acta Crystallogr. Sect. A 26 , 244–254 (1970).

Paściak, M., Welberry, T., Kulda, J., Leoni, S. & Hlinka, J. Dynamic displacement disorder of cubic BaTiO 3 . Phys. Rev. Lett. 120 , 167601 (2018).

Girshberg, Y. & Yacoby, Y. Ferroelectric phase transitions and off-centre displacements in systems with strong electronphonon interaction. J. Phys. Condens. Matter 11 , 9807–9822 (1999).

Pirc, R. & Blinc, R. Off-center Ti model of barium titanate. Phys. Rev. B 70 , 134107 (2004).

Paściak, M., Boulfelfel, S. E. & Leoni, S. Polarized cluster dynamics at the paraelectric to ferroelectric phase transition in BaTiO 3 . J. Phys. Chem. B 114 , 16465–16470 (2010).

Cohen, R. E. & Krakauer, H. Lattice dynamics and origin of ferroelectricity in BaTiO 3 : linearized-augmented-plane-wave total-energy calculations. Phys. Rev. B 42 , 6416–6423 (1990).

Cohen, R. E. Origin of ferroelectricity in perovskite oxides. Nature 358 , 136–138 (1992).

Ghosez, P., Gonze, X. & Michenaud, J.-P. Lattice dynamics and ferroelectric instability of barium titanate. Ferroelectrics 194 , 39–54 (1997).

Ghosez, P. H. S. H., Gonze, X. & Michenaud, J. P. Ab initio phonon dispersion curves and interatomic force constants of barium titanate. Ferroelectrics 206 , 205–217 (1998).

Zhang, Q., Cagin, T. & Goddard, W. A. The ferroelectric and cubic phases in BaTiO 3 ferroelectrics are also antiferroelectric. Proc. Natl Acad. Sci. USA 103 , 14695–14700 (2006).

Kotiuga, M. et al. Microscopic picture of paraelectric perovskites from structural prototypes. Phys. Rev. Res. 4 , L012042 (2022).

Zhao, X.-G., Malyi, O. I., Billinge, S. J. L. & Zunger, A. Intrinsic local symmetry breaking in nominally cubic paraelectric BaTiO 3 . Phys. Rev. B 105 , 224108 (2022).

Zhong, W., Vanderbilt, D. & Rabe, K. M. First-principles theory of ferroelectric phase transitions for perovskites: the case of BaTiO 3 . Phys. Rev. B 52 , 6301–6312 (1995).

Tinte, S., Íñiguez, J., Rabe, K. M. & Vanderbilt, D. Quantitative analysis of the first-principles effective Hamiltonian approach to ferroelectric perovskites. Phys. Rev. B 67 , 064106 (2003).

Tinte, S., Stachiotti, M. G., Sepliarsky, M., Migoni, R. L. & Rodriguez, C. O. Atomistic modelling of BaTiO 3 based on first-principles calculations. J. Phys. Condens. Matter 11 , 9679–9690 (1999).

Ponomareva, I., Bellaiche, L., Ostapchuk, T., Hlinka, J. & Petzelt, J. Terahertz dielectric response of cubic BaTiO 3 . Phys. Rev. B 77 , 012102 (2008).

Qi, Y., Liu, S., Grinberg, I. & Rappe, A. M. Atomistic description for temperature-driven phase transitions in BaTiO 3 . Phys. Rev. B 94 , 134308 (2016).

Krakauer, H., Yu, R., Wang, C.-Z., Rabe, K. M. & Waghmare, U. V. Dynamic local distortions in KNbO 3 . J. Phys. Condens. Matter 11 , 3779–3787 (1999).

Walizer, L., Lisenkov, S. & Bellaiche, L. Finite-temperature properties of (Ba, Sr)TiO 3 systems from atomistic simulations. Phys. Rev. B 73 , 144105 (2006).

Zhang, J., Wei, D., Zhang, F., Chen, X. & Wang, D. Structural phase transition of two dimensional single-layer SnTe from artificial neural network. Preprint at https://arxiv.org/abs/2012.11137 (2020).

Liechtenstein, A. I., Anisimov, V. I. & Zaanen, J. Density-functional theory and strong interactions: orbital ordering in Mott-Hubbard insulators. Phys. Rev. B 52 , 5467 (1995).

Dudarev, S. L., Botton, G. A., Savrasov, S. Y., Humphreys, C. J. & Sutton, A. P. Electron-energy-loss spectra and the structural stability of nickel oxide: an LSDA+U study. Phys. Rev. B 57 , 1505 (1998).

Perdew, J. P. Jacob’s ladder of density functional approximations for the exchange-correlation energy. In AIP Conference Proceedings 1–20 (AIP, 2001).

Becke, A. D. A new mixing of hartree-fock and local density-functional theories. J. Chem. Phys. 98 , 1372 (1993).

Bartók, A. P., Payne, M. C., Kondor, R. & Csányi, G. Gaussian approximation potentials: the accuracy of quantum mechanics, without the electrons. Phys. Rev. Lett. 104 , 136403 (2010).

Grisafi, A., Wilkins, D. M., Csányi, G. & Ceriotti, M. Symmetry-adapted machine learning for tensorial properties of atomistic systems. Phys. Rev. Lett. 120 , 036002 (2018).

Behler, J. & Parrinello, M. Generalized neural-network representation of high-dimensional potential-energy surfaces. Phys. Rev. Lett. 98 , 146401 (2007).

Rupp, M., Tkatchenko, A., Müller, K.-R. & von Lilienfeld, O. A. Fast and accurate modeling of molecular atomization energies with machine learning. Phys. Rev. Lett. 108 , 058301 (2012).

Montavon, G. et al. Machine learning of molecular electronic properties in chemical compound space. N. J. Phys. 15 , 095003 (2013).

Behler, J. First principles neural network potentials for reactive simulations of large molecular and condensed systems. Angew. Chem. Int. Ed. 56 , 12828–12840 (2017).

Deringer, V. L., Caro, M. A. & Csányi, G. Machine learning interatomic potentials as emerging tools for materials science. Adv. Mater. 31 , 1902765 (2019).

Noé, F., Tkatchenko, A., Müller, K.-R. & Clementi, C. Machine learning for molecular simulation. Annu. Rev. Phys. Chem. 71 , 361–390 (2020).

Butler, K. T., Davies, D. W., Cartwright, H., Isayev, O. & Walsh, A. Machine learning for molecular and materials science. Nature 559 , 547–555 (2018).

Schütt, K. T., Sauceda, H. E., Kindermans, P.-J., Tkatchenko, A. & Müller, K.-R. SchNet - A deep learning architecture for molecules and materials. J. Chem. Phys. 148 , 241722 (2018).

Friederich, P., Häse, F., Proppe, J. & Aspuru-Guzik, A. Machine-learned potentials for next-generation matter simulations. Nat. Mater. 20 , 750–761 (2021).

Lopanitsyna, N., Ben Mahmoud, C. & Ceriotti, M. Finite-temperature materials modeling from the quantum nuclei to the hot electron regime. Phys. Rev. Mater. 5 , 043802 (2021).

Imbalzano, G. et al. Uncertainty estimation for molecular dynamics and sampling. J. Chem. Phys. 154 , 074102 (2021).

Cheng, B., Engel, E. A., Behler, J., Dellago, C. & Ceriotti, M. Ab initio thermodynamics of liquid and solid water. Proc. Natl Acad. Sci. USA 116 , 1110–1115 (2019).

Carleo, G. et al. Machine learning and the physical sciences. Rev. Mod. Phys. 91 , 045002 (2019).

Dragoni, D., Daff, T. D., Csányi, G. & Marzari, N. Achieving DFT accuracy with a machine-learning interatomic potential: thermomechanics and defects in bcc ferromagnetic iron. Phys. Rev. Mater. 2 , 013808 (2018).

Bartók, A. P., Kermode, J., Bernstein, N. & Csányi, G. Machine learning a general-purpose interatomic potential for silicon. Phys. Rev. X 8 , 041048 (2018).

Google Scholar  

Isayev, O. et al. Materials cartography: representing and mining materials space using structural and electronic fingerprints. Chem. Mater. 27 , 735–743 (2015).

Sanchez-Lengeling, B. & Aspuru-Guzik, A. Inverse molecular design using machine learning: generative models for matter engineering. Science 361 , 360–365 (2018).

Szlachta, W. J., Bartók, A. P. & Csányi, G. Accuracy and transferability of Gaussian approximation potential models for tungsten. Phys. Rev. B 90 , 104108 (2014).

Deringer, V. L. & Csányi, G. Machine learning based interatomic potential for amorphous carbon. Phys. Rev. B 95 , 094203 (2017).

Morawietz, T., Singraber, A., Dellago, C. & Behler, J. How van der waals interactions determine the unique properties of water. Proc. Natl Acad. Sci. USA 113 , 8368–8373 (2016).

Caro, M. A., Deringer, V. L., Koskinen, J., Laurila, T. & Csányi, G. Growth mechanism and origin of high sp 3 content in tetrahedral amorphous carbon. Phys. Rev. Lett. 120 , 166101 (2018).

Sosso, G. C., Miceli, G., Caravati, S., Behler, J. & Bernasconi, M. Neural network interatomic potential for the phase change material GeTe. Phys. Rev. B 85 , 174103 (2012).

Eshet, H., Khaliullin, R. Z., Kühne, T. D., Behler, J. & Parrinello, M. Microscopic origins of the anomalous melting behavior of sodium under high pressure. Phys. Rev. Lett. 108 , 115701 (2012).

Khaliullin, R. Z., Eshet, H., Kühne, T. D., Behler, J. & Parrinello, M. Graphite-diamond phase coexistence study employing a neural-network mapping of the ab initio potential energy surface. Phys. Rev. B 81 , 100103 (2010).

Jinnouchi, R., Lahnsteiner, J., Karsai, F., Kresse, G. & Bokdam, M. Phase transitions of hybrid perovskites simulated by machine-learning force fields trained on the fly with Bayesian inference. Phys. Rev. Lett. 122 , 225701 (2019).

Park, C. W. et al. Accurate and scalable graph neural network force field and molecular dynamics with direct force architecture. npj Comput. Mater. 7 , 1–9 (2021).

Thompson, A., Swiler, L., Trott, C., Foiles, S. & Tucker, G. Spectral neighbor analysis method for automated generation of quantum-accurate interatomic potentials. J. Comput. Phys. 285 , 316–330 (2015).

Shapeev, A. V. Moment tensor potentials: a class of systematically improvable interatomic potentials. Multiscale Model. Simul. 14 , 1153–1173 (2016).

van der Oord, C., Dusson, G., Csányi, G. & Ortner, C. Regularised atomic body-ordered permutation-invariant polynomials for the construction of interatomic potentials. Mach. Learn. Sci. Technol. 1 , 015004 (2020).

Nigam, J., Pozdnyakov, S. & Ceriotti, M. Recursive evaluation and iterative contraction of N -body equivariant features. J. Chem. Phys. 153 , 121101 (2020).

Artrith, N., Morawietz, T. & Behler, J. High-dimensional neural-network potentials for multicomponent systems: applications to zinc oxide. Phys. Rev. B 83 , 153101 (2011).

Bereau, T., DiStasio, R. A., Tkatchenko, A. & von Lilienfeld, O. A. Non-covalent interactions across organic and biological subsets of chemical space: physics-based potentials parametrized from machine learning. J. Chem. Phys. 148 , 241706 (2018).

Veit, M. et al. Equation of state of fluid methane from first principles with machine learning potentials. J. Chem. Theory Comput. 15 , 2574–2586 (2019).

Ko, T. W., Finkler, J. A., Goedecker, S. & Behler, J. A fourth-generation high-dimensional neural network potential with accurate electrostatics including non-local charge transfer. Nat. Commun. 12 , 398 (2021).

Musil, F. et al. Efficient implementation of atom-density representations. J. Chem. Phys. 154 , 114109 (2021).

Himanen, L., Geurts, A., Foster, A. S. & Rinke, P. Data-driven materials science: status, challenges, and perspectives. Adv. Sci. 6 , 1900808 (2019).

Csányi, G., Willatt, M. J. & Ceriotti, M. in Machine Learning Meets Quantum Physics (eds Schütt, K. T. et al.) Ch. 6 (Springer International Publishing, 2020).

Musil, F. et al. Physics-inspired structural representations for molecules and materials. Chem. Rev. 121 , 9759–9815 (2021).

Bartók, A. P., Kondor, R. & Csányi, G. On representing chemical environments. Phys. Rev. B 87 , 184115 (2013).

Grisafi, A. & Ceriotti, M. Incorporating long-range physics in atomic-scale machine learning. J. Chem. Phys. 151 , 204105 (2019).

Grisafi, A., Nigam, J. & Ceriotti, M. Multi-scale approach for the prediction of atomic scale properties. Chem. Sci. 12 , 2078–2090 (2021).

Imbalzano, G. & Ceriotti, M. Modeling the Ga/As binary system across temperatures and compositions from first principles. Phys. Rev. Mater. 5 , 063804 (2021).

Deringer, V. L. et al. Origins of structural and electronic transitions in disordered silicon. Nature 589 , 59–64 (2021).

Gastegger, M., Behler, J. & Marquetand, P. Machine learning molecular dynamics for the simulation of infrared spectra. Chem. Sci. 8 , 6924–6935 (2017).

Laurens, G., Rabary, M., Lam, J., Peláez, D. & Allouche, A.-R. Infrared spectra of neutral polycyclic aromatic hydrocarbons based on machine learning potential energy surface and dipole mapping. Theor. Chem. Acc. 140 , 66 (2021).

Kapil, V., Wilkins, D. M., Lan, J. & Ceriotti, M. Inexpensive modeling of quantum dynamics using path integral generalized Langevin equation thermostats. J. Chem. Phys. 152 , 124104 (2020).

Zhang, L. et al. Deep neural network for the dielectric response of insulators. Phys. Rev. B 102 , 041121 (2020).

Chaves, A. S., Barreto, F. C. S., Nogueira, R. A. & Zẽks, B. Thermodynamics of an eight-site order-disorder model for ferroelectrics. Phys. Rev. B 13 , 207–212 (1976).

Comes, R., Lambert, M. & Guinier, A. The chain structure of BaTiO 3 and KNbO 3 . Solid State Commun. 6 , 715–719 (1968).

Roberts, S. Adiabatic study of the 128° C transition in barium titanate. Phys. Rev. 85 , 925–926 (1952).

Zhong, W., Vanderbilt, D. & Rabe, K. M. Phase transitions in BaTiO 3 from first principles. Phys. Rev. Lett. 73 , 1861–1864 (1994).

Decker, D. L. & Zhao, Y. X. Dielectric and polarization measurements on BaTiO 3 at high pressures to the tricritical point. Phys. Rev. B 39 , 2432–2438 (1989).

Senn, M., Keen, D., Lucas, T., Hriljac, J. & Goodwin, A. Emergence of long-range order in batio 3 from local symmetry-breaking distortions. Phys. Rev. Lett. 116 , 207602 (2016).

Bencan, A. et al. Atomic scale symmetry and polar nanoclusters in the paraelectric phase of ferroelectric materials. Nat. Commun. 12 , 3509 (2021).

Vanderbilt, D. & Zhong, W. First-principles theory of structural phase transitions for perovskites: competing instabilities. Ferroelectrics 206 , 181–204 (1998).

Akbarzadeh, A. R., Bellaiche, L., Leung, K., Íñiguez, J. & Vanderbilt, D. Atomistic simulations of the incipient ferroelectric KTaO 3 . Phys. Rev. B 70 , 054103 (2004).

Giberti, F., Cheng, B., Tribello, G. A. & Ceriotti, M. Iterative unbiasing of quasi-equilibrium sampling. J. Chem. Theory Comput. 16 , 100–107 (2020).

Xie, P., Chen, Y., E, W. & Car, R. Ab initio multi-scale modeling of ferroelectrics: the case of PbTiO 3 Preprint at https://arxiv.org/abs/2205.11839 (2022).

Fischer, G. J., Wang, Z. & Karato, S.-i Elasticity of CaTiO 3 , SrTiO 3 and BaTiO 3 perovskites up to 3.0 Gpa: the effect of crystallographic structure. Phys. Chem. Miner. 20 , 97–103 (1993).

Kay, H. & Vousden, P. XCV. Symmetry changes in barium titanate at low temperatures and their relation to its ferroelectric properties. Lond. Edinb. Dublin Philos. Mag. J. Sci. 40 , 1019–1040 (1949).

Sharma, M., Resta, R. & Car, R. Dipolar correlations and the dielectric permittivity of water. Phys. Rev. Lett. 98 , 247401 (2007).

Resta, R. & Vanderbilt, D. In Physics of Ferroelectrics: A Modern Perspective (eds Rabe, K. M., Ahn, C. H. & Triscone, J.-M.) pp. 31–68 (Springer, 2007).

Spaldin, N. A. A beginner’s guide to the modern theory of polarization. J. Solid State Chem. 195 , 2–10 (2012).

Hashimoto, T. & Moriwake, H. Dielectric properties of BaTiO 3 by molecular dynamics simulations using a shell model. Mol. Simul. 41 , 1074–1080 (2015).

MacDowell, L. G. & Vega, C. Dielectric constant of ice Ih and ice V: a computer simulation study. J. Phys. Chem. B 114 , 6089–6098 (2010).

Li, Z., Grimsditch, M., Foster, C. M. & Chan, S. K. Dielectric and elastic properties of ferroelectric materials at elevated temperature. J. Phys. Chem. Solids 57 , 1433–1438 (1996).

Ostapchuk, T., Petzelt, J., Savinov, M., Buscaglia, V. & Mitoseriu, L. Grain-size effect in BaTiO 3 ceramics: study by far infrared spectroscopy. Phase Transit. 79 , 361–373 (2006).

Davis, L. & Rubin, L. G. Some dielectric properties of barium-strontium titanate ceramics at 3000 megacycles. J. Appl. Phys. 24 , 1194–1197 (1953).

Chu, F., Sun, H.-T., Zhang, L.-Y. & Yao, X. Temperature dependence of ultra-low-frequency dielectric relaxation of barium titanate ceramic. J. Am. Ceram. Soc. 75 , 2939–2944 (1992).

Binder, K. Finite size effects on phase transitions. Ferroelectrics 73 , 43–67 (1987).

Rupprecht, G. & Bell, R. O. Dielectric constant in paraelectric perovskites. Phys. Rev. 135 , A748–A752 (1964).

Löffler, G., Schreiber, H. & Steinhauser, O. The frequency-dependent conductivity of a saturated solution of ZnBr 2 in water: a molecular dynamics simulation. J. Chem. Phys. 107 , 3135–3143 (1997).

Frenkel, D. Understanding Molecular Simulation: From Algorithms to Applications 2nd edn (Academic Press, 2002).

Luspin, Y., Servoin, J. L. & Gervais, F. Soft mode spectroscopy in barium titanate. J. Phys. C Solid State Phys. 13 , 3761–3773 (1980).

Vogt, H., Sanjurjo, J. A. & Rossbroich, G. Soft-mode spectroscopy in cubic BaTiO 3 by hyper-raman scattering. Phys. Rev. B 26 , 5904–5910 (1982).

Presting, H., Sanjurjo, J. A. & Vogt, H. Mode softening in cubic BaTiO 3 and the problem of its adequate description. Phys. Rev. B 28 , 6097–6099 (1983).

Hlinka, J., Petzelt, J., Kamba, S., Noujni, D. & Ostapchuk, T. Infrared dielectric response of relaxor ferroelectrics. Phase Transit. 79 , 41–78 (2006).

García, A. & Vanderbilt, D. Temperature-dependent dielectric response of BaTiO 3 from first principles. AIP Conf. Proc. 436 , 53–60 (1998).

Williams, K., Wagner, L. K., Cazorla, C. & Gould, T. Combining density functional theories to correctly describe the energy, lattice structure and electronic density of functional oxide perovskites. Preprint at https://arxiv.org/abs/2005.03792 (2020).

Maiti, T., Guo, R. & Bhalla, A. S. Structure-property phase diagram of BaZr x Ti 1− x O 3 system. J. Am. Ceram. Soc. 91 , 1769–1780 (2008).

Zhong, W. & Vanderbilt, D. Effect of quantum fluctuations on structural phase transitions in SrTiO 3 and BaTiO 3 . Phys. Rev. B Condens. Matter Mater. Phys. 53 , 5047–5050 (1996).

Willatt, M. J., Musil, F. & Ceriotti, M. Atom-density representations for machine learning. J. Chem. Phys. 150 , 154110 (2019).

Glielmo, A., Sollich, P. & De Vita, A. Accurate interatomic force fields via machine learning with covariant kernels. Phys. Rev. B 95 , 214302 (2017).

Anderson, B., Hy, T. S. & Kondor, R. Cormorant: covariant molecular neural networks. In 33rd Conference on Neural Information Processing Systems (eds Wallach, H. et al.) (Curran Associates, Inc., 2019).

Batzner, S. et al. E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. Nat. Commun. 13 , 2453 (2022).

Schütt, K., Unke, O. & Gastegger, M. Equivariant message passing for the prediction of tensorial properties and molecular spectra. in Proc. 38th International Conference on Machine Learning (eds Meila, M. & Zhang, T.) 9377–9388 (PMLR, 2021).

Qiao, Z. et al. Informing geometric deep learning with electronic interactions to accelerate quantum chemistry. Proc. Natl Acad. Sci. USA 119, e2205221119 (2022).

Musil, F. et al. librascal. https://github.com/cosmo-epfl/librascal (2020).

Veit, M., Wilkins, D. M., Yang, Y., DiStasio, R. A. & Ceriotti, M. Predicting molecular dipole moments by combining atomic partial charges and atomic dipoles. J. Chem. Phys. 153 , 024113 (2020).

Zhong, W., King-Smith, R. D. & Vanderbilt, D. Giant LO-TO splittings in perovskite ferroelectrics. Phys. Rev. Lett. 72 , 3618–3621 (1994).

Giannozzi, P. et al. Quantum espresso: a modular and open-source software project for quantum simulations of materials. J. Phys. Condens. Matter 21 , 395502 (2009).

Giannozzi, P. et al. Advanced capabilities for materials modelling with quantum espresso. J. Phys. Condens. Matter 29 , 465901 (2017).

Perdew, J. P. et al. Restoring the density-gradient expansion for exchange in solids and surfaces. Phys. Rev. Lett. 100 , 136406 (2008).

Pizzi, G., Cepellotti, A., Sabatini, R., Marzari, N. & Kozinsky, B. Aiida: automated interactive infrastructure and database for computational science. Comput. Mater. Sci. 111 , 218–230 (2016).

Huber, S. P. et al. Aiida 1.0, a scalable computational infrastructure for automated reproducible workflows and data provenance. Sci. Data 7 , 300 (2020).

Uhrin, M., Huber, S. P., Yu, J., Marzari, N. & Pizzi, G. Workflows in aiida: engineering a high-throughput, event-based engine for robust and modular computational workflows. Comput. Mater. Sci. 187 , 110086 (2021).

Kapil, V. et al. I-PI Software. http://ipi-code.org (2018).

Larsen, A. H. et al. The atomic simulation environment-a Python library for working with atoms. J. Phys. Condens. Matter 29 , 273002 (2017).

Togo, A. & Tanaka, I. First principles phonon calculations in materials science. Scr. Mater. 108 , 1–5 (2015).

Libbi, F., Bonini, N. & Marzari, N. Thermomechanical properties of honeycomb lattices from internal-coordinates potentials: the case of graphene and hexagonal boron nitride. 2D Mater. 8 , 015026 (2020).

Giannozzi, P., de Gironcoli, S., Pavone, P. & Baroni, S. Ab initio calculation of phonon dispersions in semiconductors. Phys. Rev. B 43 , 7231–7242 (1991).

Zhang, X., Zhang, C., Zhang, C., Zhang, P. & Kang, W. Finite-temperature phonon dispersion and vibrational dynamics of BaTiO 3 from first-principles molecular dynamics. Phys. Rev. B 105 , 014304 (2022).

Gasparotto, P., Meißner, R. H. & Ceriotti, M. Recognizing local and global structural motifs at the atomic scale. J. Chem. Theory Comput. 14 , 486–498 (2018).

Kapil, V. et al. I-PI 2.0: a universal force engine for advanced molecular simulations. Comput. Phys. Commun. 236 , 214–223 (2019).

Barducci, A., Bussi, G. & Parrinello, M. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. Phys. Rev. Lett. 100 , 020603 (2008).

Ceriotti, M., Manolopoulos, D. E. & Parrinello, M. Accelerating the convergence of path integral dynamics with a generalized Langevin equation. J. Chem. Phys. 134 , 84104 (2011).

Bussi, G., Donadio, D. & Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 126 , 14101 (2007).

Gigli, L. et al. Thermodynamics and dielectric response of BaTiO 3 by data-driven modeling. Materials Cloud Archive 2022.88, https://doi.org/10.24435/materialscloud:9g-k6 (2022).

Wilkins, D. M. & Grisafi, A. TENSOAP repository. https://github.com/dilkins/TENSOAP (2021).

Hinuma, Y., Pizzi, G., Kumagai, Y., Oba, F. & Tanaka, I. Band structure diagram paths based on crystallography. Comput. Mater. Sci. 128 , 140–184 (2017).

Download references

Acknowledgements

We thank Federico Grasselli for his insightful suggestions and critical reading of the manuscript. L.G., M.K. and M.C. were supported by the Samsung Advanced Institute of Technology (SAIT). M.V., G.P., N.M. and M.C. acknowledge support from the MARVEL National Centre of Competence in Research (NCCR), funded by the Swiss National Science Foundation (grant agreement ID 51NF40-182892). G.P. acknowledges the swissuniversities “Materials Cloud” project (number 201-003). G.P. and N.M. acknowledge support from the European Centre of Excellence MaX “Materials design at the Exascale” (824143). This work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project IDs mr0 and s1073.

Author information

Authors and affiliations.

Laboratory of Computational Science and Modeling (COSMO), Institute of Materials, École Polytechnique Fédérale de Lausanne, CH-1015, Lausanne, Switzerland

Lorenzo Gigli, Max Veit & Michele Ceriotti

Theory and Simulation of Materials (THEOS) and National Centre for Computational Design and Discovery of Novel Materials (MARVEL), École Polytechnique Fédérale de Lausanne, CH-1015, Lausanne, Switzerland

Michele Kotiuga, Giovanni Pizzi & Nicola Marzari

You can also search for this author in PubMed   Google Scholar

Contributions

G.P., N.M. and M.C. jointly supervised the project. L.G. and M.V. jointly developed, trained, and benchmarked the ML framework. L.G. ran the MD simulations. L.G. and M.V. analyzed the results of the MD simulations. M.K. performed the DFT calculations. All authors contributed to the discussion and writing of the paper.

Corresponding author

Correspondence to Lorenzo Gigli .

Ethics declarations

Competing interests.

The authors declare no competing interests.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Supplementary information, rights and permissions.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/ .

Reprints and permissions

About this article

Cite this article.

Gigli, L., Veit, M., Kotiuga, M. et al. Thermodynamics and dielectric response of BaTiO 3 by data-driven modeling. npj Comput Mater 8 , 209 (2022). https://doi.org/10.1038/s41524-022-00845-0

Download citation

Received : 16 February 2022

Accepted : 06 July 2022

Published : 29 September 2022

DOI : https://doi.org/10.1038/s41524-022-00845-0

Share this article

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

This article is cited by

Understanding of dielectric properties of cellulose.

• Mathias Boström
• Oleksandr I. Malyi

Cellulose (2024)

Dynamics of lattice disorder in perovskite materials, polarization nanoclusters and ferroelectric domain wall structures

• Jan Očenášek
• Jorge Alcalá

npj Computational Materials (2023)

Quick links

• Explore articles by subject
• Guide to authors
• Editorial policies

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

• < Previous

Home > MSE > Faculty Works > 491

Materials Science and Engineering Faculty Research & Creative Works

Dielectric constant of barium titanate powders near curie temperature.

Vladimir Petrovsky Tatiana Petrovsky Swetha Kamlapurkar Fatih Dogan , Missouri University of Science and Technology Follow

Impedance spectroscopy is an effective method to investigate the dielectric properties of powders by dispersing particles in proper liquids. It allows reliable calculation of the dielectric constant by analysis of impedance spectra for a variety of particulate materials. The temperature dependence of the dielectric constant near the Curie temperature was investigated for hydrothermal barium titanate powders suspended in ethylene glycol or butoxyethanol. Room-temperature measurements revealed that the dielectric constant of hydrothermal BaTiO3 powders was ɛ∼500. After annealing of powders at 1250°C for 1 h, the dielectric constant increased to ɛ∼3000. The temperature dependence of the impedance spectra was determined by conducting measurements within the temperature range from 20° to 150°C. It was shown that hydrothermal BaTiO3 powders with a primarily cubic crystal structure had a suppressed temperature dependence of the dielectric constant near the Curie temperature. The dielectric constant of annealed powders with a tetragonal crystal structure was significantly higher at the Curie temperature similar to the values obtained for sintered bulk barium titanate.

Recommended Citation

V. Petrovsky et al., "Dielectric Constant of Barium Titanate Powders Near Curie Temperature," Journal of the American Ceramic Society , John Wiley & Sons, Sep 2008.

The definitive version is available at https://doi.org/10.1111/j.1551-2916.2008.02693.x

Department(s)

Materials Science and Engineering

National Science Foundation (U.S.) United States. Office of Naval Research

Keywords and Phrases

Dielectric Constant; Impedance spectroscopy; Powders

International Standard Serial Number (ISSN)

0002-7820; 1551-2916

Document Type

Article - Journal

Document Version

Language(s).

© 2008 John Wiley & Sons, All rights reserved.

Publication Date

01 Sep 2008

Since June 03, 2016

Advanced Search

• Notify me via email or RSS
• Collections
• Disciplines
• All Authors
• Faculty Authors

Faculty Gallery

• View Gallery

Author Corner

• About Faculty Author Profiles

Related Content

• Departmental Web Page

Useful Links

• Library Resources

Article Locations

• View articles on map
• View articles in Google Earth

Home | About | FAQ | My Account | Accessibility Statement

Privacy Copyright

• DOI: 10.1038/160058A0
• Corpus ID: 4104585

Curie Point of Barium Titanate

• M. Harwood , P. Popper , D. F. RUSHMAN
• Published in Nature 12 July 1947
• Materials Science, Physics

64 Citations

On the phase transition in barium-lead titanate (ii), anomalous temperature coefficient of permittivity in barium titanate, dielectric residual effects in titanates, some factors influencing the dielectric properties of barium titanates, on the electrical conductivity of some alkaline earth titanates, xcvi. theory of barium titanate, specific heat and thermal expansion of batio 3, ferroelectrics and antiferroeletrics, the influence of suspended nanoparticles on the frederiks threshold of the nematic host, strain engineering of piezoelectric properties of strontium titanate thin films, 5 references, effect of temperature on the permittivity of barium titanate, corrigendum: crystal structure of double oxides of the perovskite type, high permittivity crystalline aggregates, dielectric constants of some titanates, the permittivity of polycrystals of the perovskite type, related papers.

Showing 1 through 3 of 0 Related Papers

Advertisement

Structure analyses and ferroelectric behaviour of barium titanate-doped glass–ceramic nanocrystals for energy storage applications

• Published: 16 February 2023
• Volume 129 , article number  196 , ( 2023 )

Cite this article

• M. M. El-Desoky 1 ,
• Ibrahim Morad   ORCID: orcid.org/0000-0003-0359-8805 1 , 2 , 3 , 4 ,
• H. Elhosiny Ali 2 , 3 &
• F. A. Ibrahim 4  

277 Accesses

4 Citations

Explore all metrics

New glass–ceramic (GC) nanocrystals of x BaTiO 3 –(80– x )V 2 O 5 –20PbO glasses (where x  = 5, 10, 15, 20 and 25 mol%) were synthesized via heat treatment at crystallization peak temperature ( Tp ) according to DSC thermograms. XRD together with dielectric measurements and E-P hysteresis loop were used to evaluate the microstructural and ferroelectric characteristics. Combining these methods made it feasible to improve the conditions for the production of the obtained nanomaterial and to identify correspondences among its nanostructure and ferroelectric features. The ability of appropriate heat treatment to transform glasses into nanocrystalline materials with crystallites smaller than 60 nm embedded in the glassy matrix was demonstrated by XRD measurements. The present glasses’ fulfilled dielectric constant values do not show any ferroelectric behavior. Nevertheless, by thermal treatment of the glass system at T p , GC nanocrystals exhibited an average broad peak of around 330 K in the dielectric constant. The Curie temperature of BaTiO 3 with particle size smaller than 100 nm is extremely close to the average Curie temperature of 338 K measured in the current glass system. By properly adjusting teat-treatment time and BaTiO 3 content, this finding of these samples can be employed to manage BaTiO 3 crystal size and, consequently, transition temperature. As a result, the glass–ceramic samples segregated with nanocrystalline BaTiO 3 are supported by this result’s dipolar direction and phase transition. A GC nanocrystal has an intentional energy storage density of 104 mJ cm −3 . These findings indicate that the current glass–ceramic nanocrystals are a promising material for creating energy storage devices.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Subscribe and save.

• Get 10 units per month
• Download Article/Chapter or eBook
• 1 Unit = 1 Article or 1 Chapter
• Cancel anytime

Price includes VAT (Russian Federation)

Instant access to the full article PDF.

Rent this article via DeepDyve

Institutional subscriptions

Similar content being viewed by others

Relaxor ferroelectric-like behavior in barium titanate-doped glass via formation of polar clusters, grain size effects on the transport properties of li 3 v 2 (po 4 ) 3 glass–ceramic nanocomposites for lithium cathode batteries.

The Role of Nanocrystallization for the Enhancement of Structural, Electrical, and Transport Properties of BaTiO 3 -V 2 O 5 -PbO Glasses

Availability of data and materials.

All data generated or analyzed during this study are included in this article. Requests for material should be made to the corresponding author (M. M. El-Desoky).

H.C. Zeng, K. Tanaka, K. Hirao, N. Soga, J. Non-Cryst. Solids. 209 , 112–121 (1997)

ADS   Google Scholar  

M.V. Shankar, K.B.R. Varma, J. Non-Cryst. Solids. 226 , 145–154 (1998)

M.S. Al-Assiri, M.M. El-Desoky, A. Al-Hajry, A. Al-Shahrani, A.M. Al-Mogeeth, A.A. Bahgat, Phys. B. 404 , 1437–1445 (2009)

A.E. Harby, A.E. Hannora, M.S. Al-Assiri, M.M. El-Desoky, J. Mater. Sci. 27 , 8446–8454 (2016)

Google Scholar  

L. Lin, J. Huang, W. Yu, L. Zhu, H. Tao, P. Wang, Y. Xu, Z. Zhang, J. Magn. Magn. Mater. 500 , 166380 (2020)

N. Syam Prasad, K.B.R. Varma, J. Non-Cryst. Solids. 351 , 1455–1465 (2005)

Y. Hu, C.L. Huang, J. Non-Cryst. Solids. 278 , 170–177 (2000)

Montagne, L., L. Cormier, and D. Caurant, 25. Bibliographie , in Du verre au cristal . EDP Sciences, 501–578 (2021)

B.C. Babu, B.V. Rao, M. Ravi, S. Babu, J. Mol. Struct. 1127 , 6–14 (2017)

M.M. El-Desoky, F.A. Ibrahim, A.G. Mostafa, M.Y. Hassaan, Mater. Res. Bull. 45 , 1122–1126 (2010)

M.M. El-Desoky, Mater. Chem. Phys. 119 , 389–394 (2010)

J. Garbarczyk, P. Jozwiak, M. Wasiucionek, J. Nowinski, J. Power Sour. 173 , 743–747 (2007)

N. SyamPrasad, K.R. Varma, Y. Takahashi, Y. Benino, T. Fujiwara, T. Komatsu, J. Solid State Chem. 173 , 209–215 (2003)

M.M. El-Desoky, Phys. Stat. Sol. (A). 195 , 422–428 (2003)

K. Tanaka, K. Kashima, K. Hirao, N. Soga, A. Mito, H. Nasu, J. Non-Cryst. Solids. 185 , 123–126 (1995)

Y. Ohta, M. Kitayama, K. Kaneko, S. Toh, F. Shimizu, K. Morinaga, J. Am. Ceramic Soc. 88 , 1634–1636 (2005)

L.A. Thomas, Ferroelectrics 3 , 231–238 (1972)

A.S. Aricò, P. Bruce, B. Scrosati, J.-M. Tarascon, W. van Schalkwijk, Nat. Mater. 4 , 366–377 (2005)

X. Hao, J. Adv. Dielect. 03 , 1330001–1330009 (2013)

L. Yang, X. Kong, F. Li, H. Hao, Z. Cheng, H. Liu, J.-F. Li, S. Zhang, Progress Mater. Sci. 102 , 72–108 (2019)

J.E. Shelby, Introduction to glass science and technology (Royal Society of Chemistry, 2020)

E.P. Gorzkowski, M.J. Pan, B. Bender, C.C.M. Wu, J. Electroceramics. 18 , 269–276 (2007)

A.M. Ali, A.E. Hannora, E. El-Falaky, M.M. El-Desoky, J. Non-Cryst. Solids. 584 , 121382 (2022)

M. Al-Assiri, S. Salem, M. El-Desoky, J. Phys. Chem. Solids. 67 , 1873–1881 (2006)

M. El-Desoky, M. Al-Assiri, A. Bahgat, J. Alloys Compd. 590 , 572–578 (2014)

M.M. El-Desoky, M.S. Al-Assiri, Mater. Sci. Eng. B. 137 , 237–246 (2007)

N.K. Wally, E. Sheha, B.M. Kamal, A.E. Hannora, M.M. El-Desoky, J. Alloys Compd. 895 , 162644 (2022)

S. Sakka, J. Mackenzie, J. Non-Cryst. Solids 6 , 145–162 (1971)

M. El-Desoky, F. Ibrahim, A. Mostafa, M. Hassaan, J Mater Res Bull. 45 , 1122–1126 (2010)

J.E. Garbarczyk, P. Jozwiak, M. Wasiucionek, J.L. Nowinski, Solid State Ionics 177 , 2585–2588 (2006)

I. Morad, H.E. Ali, M. Wasfy, A. Mansour, M. El-Desoky, Vacuum 181 , 109735 (2020)

M. Salah, I. Morad, H.E. Ali, M.M. Mostafa, M.M. El-Desoky, J. Inorg. Organomet. Polym. Mater. 31 , 3700–3710 (2021)

M. Sadhukhan, D. Modak, B. Chaudhuri, J. Appl. Phys. 85 , 3477–3487 (1999)

W. Xiaoyong, F. Yujun, Y. Xi, Appl. Phys. Lett. 83 , 2031–2033 (2003)

Y.M. Poplavko, V. Bovtun, N.N. Krainik, GAe. Smolenskii, J. Fizika Tverdogo Tela. 27 , 3161–3163 (1985)

A.A. Bokov, M. Maglione, A. Simon, Z.G. Ye, Ferroelectrics 337 , 169–178 (2006)

M.E. Lines, A.M. Glass, Principles and applications of ferroelectrics and related materials (Oxford University Press, 2001)

T. Hoshina, H. Kakemoto, T. Tsurumi, S. Wada, M. Yashima, J. Appl. Phys. 99 , 054311 (2006)

M. Al-Assiri, M. El-Desoky, J. Non-cryst. Solids. 358 , 1605–1610 (2012)

K.A. Shore, Contemp. Phys. 55 , 337–337 (2014)

D.J. Thouless, Science 207 , 1196–1197 (1980)

J. Wang, N. Sun, Y. Li, Q. Zhang, X. Hao, X. Chou, Ceram. Int. 43 , 7804–7809 (2017)

Download references

Acknowledgements

The authors extend their appreciation to the Ministry of Education in KSA for funding this research work through the project number KKU-IFP2-DA-5.

Author information

Authors and affiliations.

Department of Physics, Faculty of Science, Suez University, Suez, 43518, Egypt

M. M. El-Desoky & Ibrahim Morad

Department of Physics, Faculty of Science, King Khalid University, P. O. Box 9004, Abha, Saudi Arabia

Ibrahim Morad & H. Elhosiny Ali

Department of Physics, Faculty of Science, Zagazig University, Zagazig, Egypt

Department of Physics, Faculty of Science, El-Arish University, Al-Arish, Egypt

Ibrahim Morad & F. A. Ibrahim

You can also search for this author in PubMed   Google Scholar

Corresponding author

Correspondence to M. M. El-Desoky .

Ethics declarations

Conflict of interest.

The authors declare that there is no conflict of interest in the current article.

Additional information

Publisher's note.

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions

About this article

El-Desoky, M.M., Morad, I., Ali, H.E. et al. Structure analyses and ferroelectric behaviour of barium titanate-doped glass–ceramic nanocrystals for energy storage applications. Appl. Phys. A 129 , 196 (2023). https://doi.org/10.1007/s00339-023-06474-8

Download citation

Received : 28 October 2022

Accepted : 06 February 2023

Published : 16 February 2023

DOI : https://doi.org/10.1007/s00339-023-06474-8

Share this article

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

• Ferroelectricity
• Barium titanate
• Glass–ceramic nanocrystals
• Dielectric constant
• Find a journal
• Publish with us
• Track your research

Structural, optical, and electronic properties of barium titanate: experiment characterisation and first-principles study

• University of Miskolc

• Delft University of Technology

• University of Leeds

Discover the world's research

• 25+ million members
• 160+ million publication pages
• 2.3+ billion citations
• Muneerah Alomar

• OPT QUANT ELECTRON

• Hongjuan Zheng
• Jingjin Liu

• Rajashree Khatua

• T.P. Vishaul

• CHINESE J PHYS

• J PHYS-CONDENS MAT

• CRYST RES TECHNOL

• Huijing Yang

• Hanxing Liu
• SENSORS-BASEL

• J Phys Conf

• J AM CERAM SOC

• Andy H. Choi

• Recruit researchers
• Join for free
• Login Email Tip: Most researchers use their institutional email address as their ResearchGate login Password Forgot password? Keep me logged in Log in or Continue with Google Welcome back! Please log in. Email · Hint Tip: Most researchers use their institutional email address as their ResearchGate login Password Forgot password? Keep me logged in Log in or Continue with Google No account? Sign up