High-throughput computation of electric polarization in solids via Berry flux diagonalization

Abstract

Electric polarization in the absence of an externally applied electric field is a key property of polar materials, but the standard interpolation-based ab initio approach to compute polarization differences within the modern theory of polarization presents challenges for automated high-throughput calculations. Berry flux diagonalization [J. Bonini et. al, Phys. Rev. B 102, 045141 (2020)] has been proposed as an efficient and reliable alternative, though it has yet to be widely deployed. Here, we assess Berry flux diagonalization using ab initio calculations of a large set of materials, introducing and validating heuristics that ensure branch alignment with a minimal number of intermediate interpolated structures. Our automated implementation of Berry flux diagonalization succeeds in cases where prior interpolation-based workflows fail due to band-gap closures or branch ambiguities. Benchmarking with ab initio calculations of 176 candidate ferroelectrics, we demonstrate the efficacy of the approach on a broad range of insulating materials and obtain accurate effective polarization values with fewer interpolated structures than prior automated interpolation-based workflows. Our real-space heuristics that can predict gauge stability a priori from ionic displacements enable a general automated framework for reliable polarization calculations and efficient high-throughput screening of chemically and structurally diverse polar insulators. These results establish Berry flux diagonalization as a robust and efficient method to compute the effective polarization of solids and to accelerate the data-driven discovery of functional polar materials.

Figures (6)

• Figure 1: Schematic of our automated workflow to compute effective polarization. First, a rigid translation is applied to the nonpolar reference cell to minimize maximal atomic displacement relative to the polar cell. If the maximal atomic displacement between the polar and translated nonpolar reference cell is larger than 0.3 Å, intermediate structures are created by linearly interpolating the atomic coordinates. Then, Bloch wavefunctions for a uniform $k$-grid are computed using DFT. The Berry flux diagonalization algorithm is directly applied to wavefunctions corresponding to the polar and nonpolar cells (and any intermediates) to obtain a value for polarization.
• Figure 2: Schematic of a Berry phase "plaquette" used in the BFD calculation of Berry phase differences $\phi_n^p$ along the $z$-direction. Each plaquette is defined by a closed $k$-point loop $\mathcal{P}$ involving adjacent $k$-points from both the polar and nonpolar structures. To compute $\phi$ for a single plaquette, we compute overlap matrices between Bloch wavefunctions along the path $\mathcal{P} = (1) \, \mathbf{k}_2^{\text{np}} \rightarrow (2) \, \mathbf{k}_2^{\text{pol}} \rightarrow (3) \, \mathbf{k}_3^{\text{pol}} \rightarrow (4) \, \mathbf{k}_3^{\text{np}} \rightarrow (1) \, \mathbf{k}_2^{\text{np}}$. That is, we compute overlap matrix $M^{\langle 1, 2 \rangle}$ corresponding to the overlap between Bloch wavefunctions at $k$-points labeled (1) $\mathbf{k}_2^{\text{np}}$ and (2) $\mathbf{k}_2^{\text{pol}}$, and then we compute overlap matrix $M^{\langle 2, 3 \rangle}$ corresponding to the overlap between Bloch wavefunctions at $k$-points labeled (2) $\mathbf{k}_2^{\text{pol}}$ and (3) $\mathbf{k}_3^{\text{pol}}$ in the polar cell, and likewise for $M^{\langle 3, 4 \rangle}$ and $M^{\langle 4, 1 \rangle}$. The Berry phase difference $\phi$ corresponding to this plaquette is computed via Eqs. \ref{['eq:overlap_matrix']}-\ref{['eq:berry_phase_eigs']}, and the Berry phase differences corresponding to all plaquettes defined by the $k$-points are summed to obtain the total Berry phase difference $\mathbf{\Phi}$.
• Figure 3: (a) Schematic of cubic perovskite BaTiO$_3$ test system with fixed Ba and O atoms and a Ti atom that is artificially displaced along the [001] direction (red arrow). The displacement of the Ti atom, $\Delta \mathbf{u}$, is the difference between its artificial position and its centrosymmetric position. (b) Minimal singular values of the overlap matrices $M^{\langle i, i+1 \rangle}$ and maximal eigenvalues of the unitary evolution matrices $U_{P}$ corresponding to the effective polarization calculated using Berry flux diagonalization for the BaTiO$_3$ with the artificially displaced Ti atom.
• Figure 4: Minimal singular values (top), maximal eigenvalues (middle), and effective polarization (bottom) computed from overlap matrices for standard ferroelectrics under different translations of the nonpolar structure relative to the polar structure. For BaTiO$_3$, KNbO$_3$, and PbTiO$_3$, the translations that minimize the maximal atomic displacement ensure sufficient overlap between wavefunctions at each $k$-point, with minimal singular values remaining above the threshold of 0.15 and maximal eigenvalues below $\pi$. These conditions guarantee proper branch alignment and accurate effective polarization calculations. For LiNbO$_3$, even the translation minimizing maximal atomic displacement results in insufficient overlap, with singular values approaching zero and eigenvalues near $\pi$, leading to improper resolution of polarization branches and inaccurate values of effective polarization. This can be addressed by adding interpolated structures, see main text. Note: polarization values for BaTiO$_3$ (45.8 $\mu C / \text{cm}^2$) are almost identical to those of KNbO$_3$ (47.3 $\mu C / \text{cm}^2$), so the green curve is obscured by the red one in the effective polarization plot.
• Figure 5: Calculated minimum singular values (top), maximal eigenvalues (middle) and effective polarization (bottom) as functions of the maximum atomic displacement for CrO$_3$, LiNbO$_3$, and HfO$_2$ computed from different numbers of interpolated structures. Results are presented for varying numbers of interpolations, where the points with the largest maximal atomic distance correspond to the polarizations computed directly from the polar and nonpolar reference structures (with the optimal translation applied to the nonpolar structure), and each successive point with a smaller maximal atomic displacement corresponds to one additional interpolated structure. The number of interpolations used in the calculation for each data point is labeled in the bottom panel. For points corresponding to calculations with more than one interpolation, the numbers plotted are the maxima of the maximum atomic displacements across the sets of adjacent interpolated structures used in the intermediate steps of the computation, and likewise for the maxima of the maximum eigenvalues and minima of the minium singular values.