Ab Initio bulk free energy surface of proper ferroelectrics

TL;DR

This work presents an end-to-end ab initio workflow to obtain the bulk free energy surface $\mathcal{F}(T,\boldsymbol{\mathcal{P}},\boldsymbol{\eta})$ for proper ferroelectrics by combining WT-MetaD sampling with neural-network potentials for the Hamiltonian and polarization, and a quadratic electrostriction framework to couple polarization and strain. Polarization thermodynamics are first recovered via WT-MetaD to yield $\mathcal{G}(T,\boldsymbol{\mathcal{P}})$, which is then denoised and represented by a sixth-order Landau-Devonshire form $\mathcal{G}_{\text{LD}}(T,\boldsymbol{\mathcal{P}})$ with temperature-dependent coefficients. Strain effects are incorporated through $\mathcal{F} = \mathcal{G} + \frac{N}{2} (\boldsymbol{\eta}-\hat{Z}^T\boldsymbol{Q})\hat{B}(\boldsymbol{\eta}-\hat{Z}\boldsymbol{Q})$, from which $l^{\text{ref}}(T)$, $\hat{Z}(T)$ and $\hat{B}(T)$ are extracted via least-squares and covariance analyses and then interpolated in $T$. Application to PbTiO$_3$ yields a phase-transition temperature around $T_c \approx 821$ K, with sub-meV per formula unit accuracy for the ab initio FES, and demonstrates a general methodology that can be extended to other ferroic materials and to the OpenFerro software toolkit.

Abstract

We report a systematic and accurate approach for deriving the bulk free energy surface (FES), a function of temperature, polarization, and strain, from the first-principles density functional theory (DFT) of proper ferroelectrics. The core of our approach is the metadynamics algorithm that extracts the polarization dependence of the FES from all-atom molecular dynamics simulations without an a priori ansatz. The rest of the FES is derived from the metadynamics trajectories that span the relevant phase space. We demonstrate our approach in the case of lead titanate. The errors across the phase transition, due to DFT numerics, all-atom molecular dynamics, and free energy evaluation by enhanced sampling, can be systematically controlled and are of the order of 1meV/atom. The accuracy of the resulting ab initio FES is only limited by the adopted functional approximation of DFT.

Figures (5)

• Figure 1: Schematic representation of the approach
• Figure 2: (a) Colored cross-sectional ($\mathcal{P}_z=0$) representation of $\mathcal{G} (T, \boldsymbol{\mathcal{P}})$. The contour lines are equipotential lines. Each panel in inset contains a sketch of the energy basins in $\mathcal{G} (T, \boldsymbol{\mathcal{P}})$ as blue shades. The arrows illustrate local dipole configurations. (b) Top: spontaneous polarization from unbiased MD simulations xie2024pto. Bottom: Optimal LD coefficients (scaled) as functions of temperature. Scaling factors are indicated in the legends. (c) Scatter plot of all the $\mathcal{G} (T, \boldsymbol{\mathcal{P}})$ data versus $\mathcal{G}_{\text{LD}}(T, \boldsymbol{\mathcal{P}})$. The inset shows the error distribution.
• Figure 3: ($\mathrm{a_1}$-$\mathrm{a_3}$) Optimal $l^{\text{ref}}$, $\hat{B}$, and $\hat{Z}$ as functions of temperature. Circles are from direct least squares fits. Solid lines are their linear or quadratic interpolation. Fading lines are phenomenological estimates assuming $T$ independence from Ref. haun1987thermodynamic. (b-c) 2D histograms of $\mathrm{p}(\mathcal{P}_x, \eta_x - Z_{12}(\mathcal{P}^2_y+\mathcal{P}^2_z))$ and $\mathrm{p}(\sqrt{\mathcal{P}_y^2+\mathcal{P}_z^2}, \eta_x - Z_{11}Q_x)$, obtained at $T=800$K. Red lines represent $\boldsymbol{\eta}=\hat{Z}(800\text{K})\boldsymbol{Q}$. (d) Comparison between $- k_BT \log \mathrm{p}_i(\zeta_i)$ (shifted by $N\epsilon_{0,i}(T)$ and represented by blue, green and orange markers) and the optimal theory $\frac{1}{2} N \zeta_i^2 B_{11}(T)$ (red dashed line).
• Figure 4: (a) $G(T=815\text{K},d_c)$ as a function of $d_c$ for different values of $L$. (b) $G(T,d_c)$ as a function of $d_c$ for different temperatures. (c) Free energy difference $\Delta G(T)$ computed with different supercell sizes. (d) Size dependence of the free energy barrier for the transition from the ferroelectric state to the paraelectric state at T=820K.
• Figure 5: (a) Heat and contour maps of the $\mathcal{P}_z=0$ section of the instantaneous free energy $\mathcal{G}^{[t]}( T, \boldsymbol{\mathcal{P}})$ at different instants of an 8ns-long WT-MetaD simulation. The white areas indicate unexplored regions. (b) The $\mathcal{P}_z=0$ sections of the $\mathcal{G}( T, \boldsymbol{\mathcal{P}})$ and $\mathcal{G}_{\text{LD}}( T, \boldsymbol{\mathcal{P}})$ at two different temperatures. The first and third panels are associated with $\mathcal{G}( T, \boldsymbol{\mathcal{P}})$ from WT-MetaD simulations. The second and fourth panels are associated with the optimal Landau-Devonshire (LD) model $\mathcal{G}_{\text{LD}}( T, \boldsymbol{\mathcal{P}})$.