The diffusion of oxygen vacancies in BaTiO ₃ doped with Mg and Ca has been investigated by employing molecular dynamics (MD). The supercell consisted of 10 ​× ​10 ​× ​10 unit cells. The Mg ions may occupy the B-site sublattice, serving as acceptors which are compensated by oxygen vacancies. The Mg – content was equal to 1%, corresponding to 10 oxygen vacancies in the simulation box. The isovalent Ca ions were added to the A-site sublattice with concentrations ranging from 5 to 15%. The diffusivity of oxygen vacancies was obtained from the mean square displacement of oxygen at temperatures varying from 1273 to 2500 ​K. The self-diffusion coefficient of oxygen vacancies is found to be almost independent of the Ca-content. The activation energies extracted from Arrhenius plots decrease with increasing temperature with typical values between 1.3 ​eV and 1.6 ​eV. The MD results for oxygen diffusion in BaTiO ₃ are interpreted in detail in terms of radial pair distribution functions for cation – oxygen correlations as well as migration barriers for various different oxygen vacancy hops derived from nudged elastic band (NEB) calculations. The coordination number for Mg–O (first coordination shell) is close to 5 ​at ​T ​< ​1500 ​K, almost irrespective of the Ca-concentration, and approaching 5.7 ​at 2500 ​K, indicating the formation of Mg – oxygen vacancy associates which gives rise to the enhancement of the activation energy with decreasing temperature. Furthermore, the energy barriers for oxygen vacancy hopping, involving free oxygen vacancies as well as Mg - oxygen vacancy associates, have been found to range from 0.82 to 1.95 ​eV.