Atomistic simulation methods based on pair-wise interatomic potentials and energy minimization have been applied to elucidate the energetics of cation vacancies and the incorporation of 13 trivalent M3+ cations (Cr3+, Ga3+, Fe3+, Lu3+, Yb3+, Er3+, Y3+, Tb3+, Gd3+, Eu3+, Sm3+, Nd3+, La3+) in gamma-Al2O3. Calculations have been carried out using Al64O96 defect spinel supercells containing eight aluminum vacancies. The lowest energy configurations correspond to a random distribution of tetrahedral and octahedral vacancies. The energy gain in comparison with exclusive tetrahedral or octahedral vacancies is rather small (0.03 and 0.09 eV/Al2O3, respectively). Unit cell volume, density, and lattice properties of optimized structures are in good agreement with the experimental values or the results of high-quality density functional theory calculations. The trends observed for the solution energy of the M2O3 oxides in the supercell with minimum energy indicate the preferential incorporation of the foreign ions at the tetrahedral site and an increase of the solubility of M2O3 in the defect spinel in comparison with alpha-Al2O3. Configurations with the lowest energy have negative solution energies and, consequently, incorporation of trivalent ions can improve the thermodynamic stability of gamma-Al2O3 in comparison with alpha-Al2O3 and increase the gamma ->alpha transition temperature.