First-principles based Landau-Devonshire potential for BiFeO$_3$

TL;DR

This paper introduces a first-principles computational method to accurately derive the Landau-Devonshire potential for BiFeO$_3$, focusing on the energy surface near the ground state to better understand its phase transitions.

Contribution

It develops a novel approach using eigenvector-based configuration space to determine the Landau-Devonshire potential from first principles, improving the analysis of structural phase transitions.

Findings

Successfully applied to BiFeO$_3$ to compute energy contributions

Probes energy surface near the ground state more effectively

Provides detailed insights into strain, polarization, and tilt degrees of freedom

Abstract

The work describes a first-principles-based computational strategy for studying structural phase transitions, and in particular, for determination of the so-called Landau-Devonshire potential - the classical zero-temperature limit of the Gibbs energy, expanded in terms of order parameters. It exploits the configuration space attached to the eigenvectors of the modes frozen in the ground state, rather than the space spanned by the unstable modes of the high-symmetry phase, as done usually. This allows us to carefully probe the part of the energy surface in the vicinity of the ground state, which is most relevant for the properties of the ordered phase. We apply this procedure to BiFeO$_3$ and perform ab-initio calculations in order to determine potential energy contributions associated with strain, polarization and oxygen octahedra tilt degrees of freedom, compatible with its two-formula unit cell periodic boundary conditions.

Figures

Tables

Table 1. Table 1: Paths in the order-parameter space, sampled by a path parameter ξ 𝜉 \xi , for which the first-principles energies were evaluated and plotted in Fig.  1 . The last column corresponds to the position of the calculated point on the horizontal axis in the Fig  1 .

	$P_i/{\left(P_s\right)}_i$	$A_i/{\left(A_s\right)}_i$	$e_{ij}/{\left(e_s\right)}_{ij}$	x
1	($\xi$,$\xi$,$\xi$,	0,0,0,	0,0,0,0,0,0)	$\xi$
2	(0,0,0,	$\xi$,$\xi$,$\xi$,	0,0,0,0,0,0)	$\xi$
3	(0,0,0,	0,0,0,	$\xi$,$\xi$,$\xi$,$\xi$,$\xi$,$\xi$)	$\xi$
4	(1,1,1,	$\xi$,$\xi$,$\xi$,	0,0,0,0,0,0)	1+$\xi$
5	(1,1,1,	0,0,0,	$\xi$,$\xi$,$\xi$,$\xi$,$\xi$,$\xi$)	1+$\xi$
6	($\xi$,$\xi$,$\xi$,	1,1,1,	0,0,0,0,0,0)	1+$\xi$
7	(0,0,0,	1,1,1,	$\xi$,$\xi$,$\xi$,$\xi$,$\xi$,$\xi$)	1+$\xi$
8	($\xi$,$\xi$,$\xi$,	0,0,0,	1,1,1,1,1,1)	1+$\xi$
9	(0,0,0,	$\xi$,$\xi$,$\xi$,	1,1,1,1,1,1)	1+$\xi$
10	($\xi$,$\xi$,$\xi$,	1,1,1,	1,1,1,1,1,1)	2+$\xi$
11	(1,1,1,	$\xi$,$\xi$,$\xi$,	1,1,1,1,1,1)	2+$\xi$
12	(1,1,1,	1,1,1,	$\xi$,$\xi$,$\xi$,$\xi$,$\xi$,$\xi$)	2+$\xi$
13	($\xi$,$\xi$,$\xi$,	$\xi$,$\xi$,$\xi$,	$\xi$,$\xi$,$\xi$,$\xi$,$\xi$,$\xi$)	3$\xi$

Table 2. Table 2: Coefficients of the LD potential for BiFeO 3 . Second column provide one representant of symmetry-equivalent terms associated with the respective coefficient, third column gives the numerical value of the coefficient or its simple multiple, indicated in the first column.

	Type			SI	Unit
$a_1$	$P_x^2$	$-3.362$	$\times$	${10}^9$	J m C-2
$a_{11}^{\left(e\right)}$	$P_x^4$	$2.646$	$\times$	${10}^9$	J m C-2
$2a_{12}^{\left(e\right)}$	$P_x^2{P}_y^2$	$3.274$	$\times$	${10}^9$	J m C-2
$a_{111}$	$P_x^6$	$-5.960$	$\times$	${10}^8$	J m C-2
$3a_{112}$	$P_x^4{P}_y^2$	$2.634$	$\times$	${10}^8$	J m C-2
$6a_{123}$	$P_x^2{P}_y^2{P}_z^2$	$-7.132$	$\times$	${10}^9$	J m C-2
$a_{1111}$	$P_x^8$	$9.043$	$\times$	${10}^7$	J m C-2
$4a_{1112}$	$P_x^6{P}_y^2$	$-2.284$	$\times$	${10}^8$	J m C-2
$6a_{1122}$	$P_x^4{P}_y^4$	$4.636$	$\times$	${10}^8$	J m C-2
$12a_{1123}$	$P_x^4{P}_y^2{P}_z^2$	$1.493$	$\times$	${10}^9$	J m C-2
$b_1$	$A_x^2$	$-1.585$	$\times$	${10}^7$	J m-3 deg.-2
$b_{11}^{\left(e\right)}$	$A_x^4$	$5.396$	$\times$	${10}^4$	J m-3 deg.-2
$2b_{12}^{\left(e\right)}$	$A_x^2{A}_y^2$	$6.314$	$\times$	${10}^4$	J m-3 deg.-2
$b_{111}$	$A_x^6$	$-6.598$	$\times$	${10}^1$	J m-3 deg.-2
$3b_{112}$	$A_x^4{A}_y^2$	$-5.203$	$\times$	${10}^1$	J m-3 deg.-2
$6b_{123}$	$A_x^2{A}_y^2{A}_z^2$	$-5.910$	$\times$	${10}^1$	J m-3 deg.-2
$b_{1111}$	$A_x^8$	$4.890$	$\times$	${10}^{-2}$	J m-3 deg.-2
$4b_{1112}$	$A_x^6{A}_y^2$	$1.598$	$\times$	${10}^{-2}$	J m-3 deg.-2
$6b_{1122}$	$A_x^4{A}_y^4$	$1.194$	$\times$	${10}^{-1}$	J m-3 deg.-2
$12b_{1123}$	$A_x^4{A}_y^2{A}_z^2$	$1.002$	$\times$	${10}^{-1}$	J m-3 deg.-2
$t_{1111}^{\left(e\right)}$	$P_x^2{A}_x^2$	$1.720$	$\times$	${10}^7$	J m C-2 deg.-2
$t_{1122}^{\left(e\right)}$	$P_x^2{A}_y^2$	$2.273$	$\times$	${10}^7$	J m C-2 deg.-2
$4t_{1212}^{\left(e\right)}$	$P_x{P}_y{A}_x{A}_y$	$-2.844$	$\times$	${10}^7$	J m C-2 deg.-2
$t_{111111}^{42}$	$P_x^4{A}_x^2$	$2.371$	$\times$	${10}^6$	J m5 C-4 deg.-2
$t_{111111}^{24}$	$P_x^2{A}_x^4$	$-5.689$	$\times$	${10}^4$	J m C-2 deg.-4
$t_{111122}^{42}$	$P_x^4{A}_y^2$	$-9.069$	$\times$	${10}^6$	J m5 C-4 deg.-2
$t_{112222}^{24}$	$P_x^2{A}_y^4$	$-4.608$	$\times$	${10}^4$	J m C-2 deg.-4
$6t_{112233}^{42}$	$P_x^2{P}_y^2{A}_z^2$	$-8.438$	$\times$	${10}^6$	J m5 C-4 deg.-2
$6t_{112233}^{24}$	$P_x^2{A}_y^2{A}_z^2$	$-2.421$	$\times$	${10}^4$	J m C-2 deg.-4
$6t_{112211}^{42}$	$P_x^2{P}_y^2{A}_x^2$	$6.805$	$\times$	${10}^6$	J m5 C-4 deg.-2
$6t_{111122}^{24}$	$P_x^2{A}_x^2{A}_y^2$	$2.594$	$\times$	${10}^4$	J m C-2 deg.-4
$8t_{111212}^{42}$	$P_x^3{P}_y{A}_x{A}_y^2$	$-1.954$	$\times$	${10}^7$	J m5 C-4 deg.-2
$24t_{123312}^{42}$	$P_x{P}_y{P}_z^2{A}_x{A}_y$	$2.660$	$\times$	${10}^7$	J m5 C-4 deg.-2
$8t_{121112}^{24}$	$P_x{P}_y{A}_x^3{A}_y$	$-8.314$	$\times$	${10}^3$	J m C-2 deg.-4
$24t_{121233}^{24}$	$P_x{P}_y{A}_x{A}_y{A}_z^2$	$-5.457$	$\times$	${10}^4$	J m C-2 deg.-4
$t_{11111111}^{62}$	$P_x^6{A}_x^2$	$-2.573$	$\times$	${10}^6$	J m9 C-6 deg.-2
$t_{11111111}^{26}$	$P_x^2{A}_x^6$	$1.132$	$\times$	${10}^2$	J m C-2 deg.-6
$t_{11111122}^{62}$	$P_x^6{A}_y^2$	$1.370$	$\times$	${10}^6$	J m9 C-6 deg.-2
$t_{11222222}^{26}$	$P_x^2{A}_y^6$	$5.255$	$\times$	${10}^1$	J m C-2 deg.-6
$t_{11111111}^{44}$	$P_x^4{A}_x^4$	$4.747$	$\times$	${10}^4$	J m5 C-4 deg.-4
$t_{11112222}^{44}$	$P_x^4{A}_y^4$	$1.888$	$\times$	${10}^4$	J m5 C-4 deg.-4
$q_{1111}$	$e_{\mathrm{xx}}{P}_x^2$	$1.447$	$\times$	${10}^{10}$	J m C-2
$q_{1122}$	$e_{\mathrm{xx}}{P}_y^2$	$4.776$	$\times$	${10}^9$	J m C-2
$2q_{1212}$	$e_{\mathrm{xy}}{P}_x{P}_y$	$7.186$	$\times$	${10}^9$	J m C-2
$r_{1111}$	$e_{\mathrm{xx}}{A}_x^2$	$2.319$	$\times$	${10}^7$	J m-3 deg.-2
$r_{1122}$	$e_{\mathrm{xx}}{A}_y^2$	$-4.886$	$\times$	${10}^6$	J m-3 deg.-2
$2r_{1212}$	$e_{\mathrm{xy}}{A}_x{A}_y$	$-2.526$	$\times$	${10}^7$	J m-3 deg.-2
$C_{1111}$	$e_{\mathrm{xx}}^2$	$2.666$	$\times$	${10}^{11}$	Pa
$C_{1122}$	$e_{\mathrm{xx}}{e}_{\mathrm{yy}}$	$1.435$	$\times$	${10}^{11}$	Pa
$C_{1212}$	$e_{\mathrm{xy}}^2$	$9.548$	$\times$	${10}^{10}$	Pa

Table 3. Table 3: Elastic stiffness constants of rhombohedral BiFeO 3 ground state (sign of C 14 subscript 𝐶 14 C_{14} depends on adopted coordinate system).

	$C_{11}$	$C_{12}$	$C_{13}$	$C_{14}$	$C_{33}$	$C_{44}$	$C_{66}$
LD potential (this work)	278	122	95	-22	228	57	78
LDA calculations (this work)	264	147	63	-16	132	53	54
Borissenko et al.Borissenko et al. (2013), LDA	249	151	75	9	160	44	49
Shang et al.Shang et al. (2009), GGA	222	110	50	16	150	49	56

Equations

$$f=$$

$$f_{\mathrm{a}}^{(\mathrm{e})}=a_i{P}_i^2+a_{ij}^{(\mathrm{e})}{P}_i^2{P}_j^2+a_{ijk}{P}_i^2{P}_j^2{P}_k^2+a_{ijkl}{P}_i^2{P}_j^2{P}_k^2{P}_l^2,$$

$$f_{\mathrm{t}}^{(\mathrm{e})}=t_{ijkl}^{(\mathrm{e})}{P}_i{P}_j{A}_k{A}_l+t_{ijklmn}^{42}{P}_i{P}_j{P}_k{P}_l{A}_m{A}_n+t_{ijklmn}^{24}{P}_i{P}_j{A}_k{A}_l{A}_m{A}_n+t_{ijklmnpq}^{44}{P}_i{P}_j{P}_k{P}_l{A}_m{A}_n{A}_p{A}_q$$

$$+t_{ijklmnpq}^{62}{P}_i{P}_j{P}_k{P}_l{P}_m{P}_n{A}_p{A}_q+t_{ijklmnpq}^{26}{P}_i{P}_j{A}_k{A}_l{A}_m{A}_n{A}_p{A}_q,$$

$$f_{\mathrm{est}}=-q_{ijkl}{e}_{ij}{P}_k{P}_l,f_{\mathrm{rot}}=-r_{ijkl}{e}_{ij}{A}_k{A}_l,f_{\mathrm{ela}}=\frac{1}{2}{C}_{ijkl}{e}_{ij}{e}_{kl}.$$

Full text

BiFeO3, Landau potential, first-principles calculations, single-domain properties

Phenomenological models, taking into account important structural order parameters and coupling between them, can greatly help to grasp the physical mechanisms involved in various crystal structure based phenomena, such as piezoelectricity, ferroelectricity, electrostriction etc. Typically, when simple models for structural phase transitions are derived from first-principles calculations, the microscopic order parameters relevant for structural phase transition are selected from unstable phonon modes of the high-symmetry reference phase. Then, the configurational space attached to the Landau-Devonshire (LD) potential (zero-temperature limit of the Gibbs free energy functional) is defined by perturbation of the high-symmetry state along the coordinates associated with these unstable modes.Zhong (1994); Rabe (2007); Paul (2017); Nakhmanson (2010); Olsen (2016)

Here, we present an alternative approach consisting in the adjustment of the LD potential landscape in the vicinity of the ground state configuration, which is the most relevant region of the order parameter space when the ordered phase itself is in the focus of interest. This is particularly important when the disordered phase and ordered ground state have rather distinct atomic and electronic structures. It actually happens for many materials showing phase transitions at high temperatures.

An interesting material, where such an approach is desirable, is the ferroelectric BiFeO3 – an insulating material with a metallic paraelectric phase, known by its large spontaneous polarization and a large $T_{\rm C}$.Catalan07 ; Sando14 ; Neaton et al. (2005) The symmetry-breaking order parameters of its ferroelectric $Pm\bar{3}m>R3c$ phase transition, namely the ferroelectric polarization vector $\bf P$ and the oxygen-octahedron tilt vector $\bf A$, are in our approach associated with the atomic displacement patterns frozen in the fully relaxed ground state. Therefore, the eigenvectors of the high-symmetry phase dynamical matrix are not used here. Actually, in this particular case, the displacement patterns associated with the $\bf P$ and $\bf A$ vectors together with the spontaneous deformation tensor $\bf e$ define the difference between the disordered reference and the ordered ground state completely.

The energy profiles along selected paths connecting the reference paraelectric state with (${\bf P}=0$, ${\bf A}=0$ and ${\bf e}=0$) to the ab-initio ground state (${\bf P}={\bf P}_{\rm s}$, ${\bf A}={\bf A}_{\rm s}$ and ${\bf e}={\bf e}_{\rm s}$) are shown in Fig. 1. We believe that the numerical values of the LD potential parameters given in this work can be readily used for a variety of purposes such as estimation of nonlinear electromechanical properties of BiFeO3, influence of the epitaxial strain or in a range of phase-field modeling tasks.

Analytic form of the potential. The LD potential (zero-temperature potential density) $f$ considered here is expanded around the reference paraelectric state around in the usual form of the symmetry-restricted Taylor expansion in terms of 12 independent variables, covering the selected set of order parameters ${\bf P}=(P_{\rm x},P_{\rm y},P_{\rm z})$, ${\bf A}=(A_{\rm x},A_{\rm y},A_{\rm z})$ and ${\bf e}=(e_{\rm xx},e_{\rm yy},e_{\rm zz},e_{\rm yz},e_{\rm xz},e_{\rm xy})$. Resulting expression can be expressed as a sum

[TABLE]

The first three terms describe the energy expansion in term of the ferroelectric polarization and oxygen octahedron tilt only. They have been expanded till 8-th order because stopping at the 6-th order did not reproduced the ab-initio calculated potential landscape satisfactorily. The ferroelectric part of the Landau energy, $f_{\rm a}^{\rm({\rm e})}$, contains all symmetry allowed terms

[TABLE]

and the same holds for the $f_{\rm b}^{\rm({\rm e})}$, which has the same form but in terms of angles $\{A_{i}\}$ and coefficients $b$. In case of the coupling between polarization and oxygen octahedron tilt coupling terms

[TABLE]

we have included all symmetry allowed terms up to the 6-th order, while the 8-th order terms were limited to the pair interaction between single $P_{i}$ and single $A_{i}$ components only. On the other hand, since the strain contribution to the ground state LD energy are rather small, only the lowest order coupling and self-energy terms in strain were introduced. Their role is to describe electrostriction, rotostriction (introduced in an obvious analogy to electrostriction), and elastic energy contributions, respectively:

[TABLE]

Let us note that the term including all three order parameters is not included and that the adopted potential form allows a straightforward analytical elimination of the strain degree of freedom using linear equations of mechanical equilibrium,Nambu and Sagala (1994); Hlinka and Marton (2006) which facilitates considerations about a mechanically free crystal.

The superscript (e) marks the terms and parameters which are renormalized upon such strain elimination, and it emphasizes that they are related to material clamped to the reference cubic shape.

First-principles-calculations details. The total-energy calculations are based on density functional theory (DFT) within the local spin-density approximation (LSDA) method as implemented in the Vienna ab-initio Simulation Package (VASP).Kresse and Hafner (1993); Kresse and Furthmüller (1996) The projector-augmented plane-wave method was used.Blöchl (1994) There were 15 explicitly treated electrons for Bi (5d10 6s2 6p3), 14 for Fe (3p6 4s1 3d7) and six for oxygen (2s2 2p4). Energy cutoff for plane-waves was set to 500 eV. The Brillouin-zone integrations were carried out using $3\,\times\,3\,\times\,3$ Monkhorst-Pack $k$-point mesh.Monkhorst and Pack (1976) Gaussian broadening of 0.1 eV was applied.Fu and Ho (1983) The antiferromagnetic (AFM) G-type order on Fe atoms has been taken into account. The LSDA+U method is applied with the Hubbard termAnisimov et al. (1997) added to the iron $d$-orbitals. In the used DudarevDudarev et al. (1998) approach the difference $U-J$ was set to 3.0 eV.

All calculations presented here were performed using a 10-atom supercell of BiFeO3 compatible with the G-type AFM order. The two perovskite cells in the supercell host two oppositely rotated oxygen octahedra. The tilt vector ${\bf A}$ discussed in the text always refers to the axial rotation vector of the oxygen octahedron in the first cell, while the rotation in the other cell is automatically opposite (-${\bf A}$). Polarization vector ${\bf P}$ is the same in both perovskite cells.

Microscopic content of the selected order parameters. As a starting point in the presented procedure the atomic structure of the reference cubic and the rhombohedral ground states of the BiFeO3 were determined. They are fully optimized using first-principles calculations within cubic and $R3c$ symmetries, respectively, leading to configurations in good agreement with the available theoreticalGraf et al. (2014); Neaton et al. (2005) and experimentalArnold et al. (2009) data.

The difference between the ground-state and the reference configurations can be specified by differences in fractional atomic coordinates $\Delta{\bf u}_{\rm s}$ and by the change of the supercell lattice vectors. Evaluation of the strain tensor is a straightforward procedure, because the lattice vectors of any slightly distorted unit cell described by the tensor of deformation ${\bf e}$ can be obtained by applying the ${\bf 1}+{\bf e}$ matrix multiplication operation to the three reference-state lattice vectors. However, the $\Delta{\bf u}_{\rm s}$ (which is a 30-component vector for the ten-atom BiFeO3 supercell) requires more attention. It can be further decomposed into two parts, $\Delta{\bf u}_{\rm s}^{(a)}+\Delta{\bf u}_{\rm s}^{(b)}$ . The first part, which transforms as the (111) component of the $F_{1u}$ Brillouin zone center polar mode of the parent cubic phase, defines displacements $\Delta{\bf u}_{\rm s}^{(a)}$ related to the spontaneous ferroelectric polarization ${\bf P}_{\rm s}\parallel(111)$, while the remaining part, $\Delta{\bf u}_{\rm s}^{(b)}$, transforming as the Brillouin zone corner irreducible representation, defines the spontaneous oxygen octahedron tilt ${\bf A}_{\rm s}\parallel(111)$. This decomposition is thus unique.

Furthermore, we have assumed that the space of vectors ${\bf P}$ are attached to a single zone-center mode, i.e. to a mode with a fixed atomic pattern. In other words, the atomic displacements corresponding to the polarization ${\bf P}=\xi{\bf P}_{\rm s}$ are given by $\Delta{\bf u}^{(a)}=\xi\Delta{\bf u}_{\rm s}^{(a)}$ and the atomic displacements corresponding to the polarization ${\bf P}$ equal to the spontaneous one but directed along, say (100) Cartesian direction, are obtained by rotation of each individual atomic displacement vector comprised in $\Delta{\bf u}_{\rm s}^{(a)}$ by the same proper rotation that turns (111) into (100) direction in the space of the attached polarization vectors ${\bf P}$. Similar procedure is applied to relate $\Delta{\bf u}_{\rm s}^{(b)}$ with ${\bf A}$. This construction defines a consistent linear subspace of the atomic coordinates compatible with the $Z=2$ supercell boundary conditions.

Fitting parameters of the Landau-Devonshire potential. The described procedure establishes a one-to-one correspondence between point in the order parameters space and atomic structure. Using this link, it is possible to sample the parameter space, use first-principles calculations to determine energies, and fit the energy surface in order to obtain parameters of the potential.

The individual terms in (1) are fitted separately. For each term, a suitable set of paths was chosen. For example, to fit the parameters $\alpha$’s in the $f_{\rm a}^{\rm({\rm e})}(\{P_{i}\})$ several paths along the high symmetry directions in the polarization space were chosen, while keeping the oxygen octahedra tilt and strain zero. Angular dependence of energy on polarization vector was probed employing circular paths around the paraelectric reference state with several different diameters and orientations of the circles. Similar procedure was then adopted for $f_{\rm b}^{\rm({\rm e})}(\{A_{i}\})$. In evaluation of the coupling terms, e.g. between polarization and tilt, we subtract the already determined contribution of polarization and tilt and fit only the energy difference corresponding to the coupling energy. The so obtained set of parameters is consequently utilized as an initial condition for fitting of all parameters together under an additional constraint ensuring that the position of the global minimum of the potential corresponds exactly to the ground state obtained from first-principles. This fitting procedure has been conveniently accomplished in dimensionless variables normalized to the ground state values of $\Delta{\bf u}_{\rm s}^{(a)}$, $\Delta{\bf u}_{\rm s}^{(b)}$ and ${\bf e}_{\rm s}$. The final set of the numerical values of defining the LD potential parameters of BiFeO3 were rescaled with spontaneous values $({\bf e}_{\rm s})_{\rm xx}=0.136$, $({\bf e}_{\rm s})_{\rm xy}=0.0012$, ${A}_{\rm s}=|{\bf A}_{\rm s}|=14.355$  degree (our ab-initio data) and ${P}_{\rm s}=|{\bf P}_{\rm s}|=0.91$ C/m2 (adopted from Ref. Diéguez et al., 2011). Resulting Landau-Devonshire model for BiFeO3 is given in Tab. 2.

The energy profiles of the LD potential along several important directions in the order-parameter space (path listed in Tab. 1) are displayed in Fig. 1). The direct ab-initio calculated energies, presented as dots in the figure, exhibit visually perfect agreement with the predictions of the fitted LD potential. Therefore, it can be expected that the present potential also describes well the various linear response properties, in particular within its ground state. For example, the elastic tensor of BiFeO3 predicted by second order strain derivatives of the present LD potential in its rhombohedral groundstate compares fairly well with other literature data (see Tab. 3). As an even more representative performance test, we have generated 300 random configurations with the 12 order parameter components falling within $\pm 20$ percent around their spontanous values and calculated their total energies from LD potential as well as from the first principles. The agreement is also very satisfactory (see Fig. 2). Moreover, since the present potential goes well behind the so far adopted quartic anharmonicity in the polarization and tilt degrees of freedomXue et al. (2014); Kornev et al. (2007); Daraktchiev et al. (2010); Kulagin et al. (2015); Kornev et al. (2006), one can expect that the this LD potential will be more appropriate when dealing with nonlinear responses.

To conclude, we present a comprehensive and efficient procedure for extraction of the Landau-Devonshire-type potential from quantum-mechanical calculations. We believe that this scheme, based on the microscopic content of order parameters derived from the full amplitude distortions of the low-symmetry phase ground states, will enable real methodological progress for systematic development of models for a large family of materials with structural phase transitions. In the specific case of the prototypical multiferroic material BiFeO3, we present a carefully engineered Landau-Devonshire potential, which can be readily used for any analytical or computer simulations requiring its realistic phenomenological LD model. We believe that the potential and its applications will improve understanding of BiFeO3 as well as open an avenue for computationally supported engineering of BiFeO3-based functional structures.

This work was supported by the Czech Science Foundation (project No. 15-04121S). Access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum, provided under the program ”Projects of Large Research, Development, and Innovations Infrastructures” (CESNET LM2015042), is greatly appreciated.