Free energies of ferroelectric crystals from a microscopic approach

The free energy of barium titanate is computed around the Curie temperature as a function of polarization P→ from the first-principles derived Effective Hamiltonian of Zhong, Vanderbilt and Rabe [Phys. Rev. Lett. 73 (1994) 1861], through Molecular Dynamics simulations coupled to the method of the Thermodynamic Integration. The algorithms used to fix the temperature (Nosé–Hoover) and/or the pressure/stress (Parrinello–Rahman), combined with fixed-polarization molecular dynamics, allow to compute a Helmholtz free energy (fixed volume/strain) or a Gibbs free energy (fixed pressure/stress). The main feature of this approach is to calculate the gradient of the free energy in the 3-D space (PxPx, PyPy, PzPz) from the thermal averages of the forces acting on the local modes, that are obtained by Molecular Dynamics under the constraint of fixed P→. This work extends the method presented in [Phys. Rev. B 79 (2009) 064101] to the calculation of the Gibbs free energy and presents new features about the computation of the free energy of ferroelectric crystals from a microscopic approach. A careful analysis of the states of constrained polarization is performed at T=280 KT=280 K (≈15–17 K below TcTc) especially at low order parameter. These states are found reasonably homogeneous for small supercell size (L=12L=12 and L=16L=16), until inhomogeneous states are observed at low order parameter for large supercells (L=20L=20). The effect of this evolution towards multidomain configurations on the mean force and free energy curves is shown. However, for reasonable supercell sizes (L=12L=12), the free energy curves obtained are in very good agreement with phenomenological Landau potentials of the literature and the states of constrained polarization are homogeneous. Moreover, the free energy obtained is quite insensitive to the supercell size from L=12L=12 to L=16L=16 at T=280 KT=280 K, suggesting that interfacial contributions, if any, are negligible at these sizes around TcTc. The method allows a numerical estimation of the free energy barrier separating the paraelectric from the ferroelectric phase at TcTc (ΔG≈0.012–0.015 meV/5-atom cellΔG≈0.012–0.015 meV/5-atom cell). However, our tests evidence phase separation at low temperature and low order parameter, in agreement with the results of Tröster et al. [Phys. Rev. B 72 (2005) 094103]. Finally, the natural decomposition of the forces into onsite, short-range, dipole–dipole and elastic–local mode interaction allows to make the same decomposition of the free energy. Some parts of this decomposition can be directly calculated from the coefficients of the Effective Hamiltonian.