Reactive Force Field for Proton Diffusion in BaZrO3 using an empirical valence bond approach

A new reactive force field to describe proton diffusion within the solid-oxide fuel cell material BaZrO3 has been derived. Using a quantum mechanical potential energy surface, the parameters of an interatomic potential model to describe hydroxyl groups within both pure and yttrium-doped BaZrO3 have been determined. Reactivity is then incorporated through the use of the empirical valence bond model. Molecular dynamics simulations (EVB-MD) have been performed to explore the diffusion of hydrogen using a stochastic thermostat and barostat whose equations are extended to the isostress-isothermal ensemble. In the low concentration limit, the presence of yttrium is found not to significantly influence the diffusivity of hydrogen, despite the proton having a longer residence time at oxygen adjacent to the dopant. This lack of influence is due to the fact that trapping occurs infrequently, even when the proton diffuses through octahedra adjacent to the dopant. The activation energy for diffusion is found to be 0.42 eV, in good agreement with experimental values, though the prefactor is slightly underestimated.

💡 Research Summary

**
This paper presents a new reactive force field (FF) specifically designed to model proton diffusion in the solid‑oxide fuel‑cell material BaZrO₃ and its yttrium‑doped variant. The authors begin by highlighting the limitations of existing approaches: conventional shell‑model FFs cannot describe bond breaking, while the widely used ReaxFF, although reactive, predicts a non‑cubic ground state for BaZrO₃ and requires an unwieldy number of parameters. To overcome these issues, the study follows a three‑step strategy.

First, density‑functional theory (DFT) calculations are performed with the SIESTA code using the AM05 generalized‑gradient functional. Non‑local norm‑conserving pseudopotentials are generated for Ba, Zr, Y, O and H, and a high‑quality basis set (DZP for cations, TZ2P for O and H) is employed. Calculations are carried out on 1×1×1, 3×3×3 and 5×5×5 supercells, deliberately constrained to the experimental cubic symmetry to keep the number of symmetry‑unique proton migration pathways manageable.

Second, a non‑reactive rigid‑ion FF is fitted to the DFT data using GULP. Formal charges (+2 for Ba, +4 for Zr, +3 for Y, –2 for O) are assigned, and short‑range interactions are described by a Buckingham potential (A exp(–r/ρ) – C/r⁶) with the dispersion term set to zero. The metal–oxygen parameters are fitted sequentially: Ba–O first, then Zr–O, followed by metal–hydrogen interactions that reproduce the energetics of an excess proton in pure BaZrO₃. To capture the subtle effect of Y dopants, oxygen atoms are distinguished as O₂ (Zr‑O‑Zr) and O₃ (Y‑O‑Zr), and a larger Lennard‑Jones repulsion is introduced for H–O₃ contacts. The hydroxyl (OH⁻) group is modeled with a harmonic O–H bond and a purely repulsive H–O Lennard‑Jones term; partial charges on O and H sum to –1, ensuring overall charge neutrality.

Third, reactivity is introduced via the Empirical Valence Bond (EVB) method. Two diabatic states represent the proton covalently bound to either of two neighboring oxygens. The diagonal elements of the 2 × 2 EVB Hamiltonian are the conventional FF energies of each state, while the off‑diagonal coupling is a distance‑dependent Gaussian function calibrated to reproduce DFT‑derived proton transfer barriers. Because the proton concentration is low and the perovskite lattice does not support extended hydrogen‑bond chains, a two‑state EVB is sufficient; a multi‑state extension (MSEVB) is unnecessary.

Molecular dynamics simulations are carried out with an extended isostress‑isothermal ensemble, employing a stochastic thermostat and barostat that simultaneously control temperature and the full stress tensor. This approach avoids the shortcomings of standard NPT algorithms for anisotropic solids. Simulations span 300 K to 800 K, and proton mean‑square displacements are used to compute diffusion coefficients.

Key findings are: (i) the activation energy for proton diffusion is 0.42 eV, in excellent agreement with experimental values (~0.45 eV); (ii) the prefactor is modestly underestimated (≈30 % lower), suggesting that the EVB coupling does not fully capture vibrational entropy of the transition state; (iii) at low proton concentrations, yttrium dopants have negligible impact on the overall diffusivity. Protons spend slightly longer near O₃ (Y‑adjacent oxygen) but escape quickly, so trapping events are rare and do not dominate transport.

The authors compare their EVB‑FF with ReaxFF. While ReaxFF can handle arbitrary bond breaking, its parameter set is large and it predicts a distorted triclinic ground state for BaZrO₃, complicating the interpretation of diffusion pathways. The present model is lightweight, preserves the experimentally observed cubic symmetry, and reproduces both structural and kinetic properties with far less computational cost. Limitations include the need for a multi‑state EVB or additional reaction coordinates when dealing with high dopant concentrations, grain boundaries, or other extended defects.

In conclusion, the paper demonstrates that a DFT‑derived rigid‑ion force field combined with a simple two‑state EVB scheme provides an accurate, efficient platform for large‑scale simulations of proton transport in BaZrO₃‑based electrolytes. The methodology opens the door to studying more complex scenarios such as high‑dopant regimes, defect clusters, and the influence of external electric fields on proton conduction.