Method for locating low-energy solutions within $\text{DFT}+U$

The widely employed $\text{DFT}+U$ formalism is known to give rise to many self-consistent yet energetically distinct solutions in correlated systems, which can be highly problematic for reliably predicting the thermodynamic and physical properties of such materials. Here we study this phenomenon in the bulk materials ${\text{UO}}_2$, CoO, and NiO, and in a ${\text{CeO}}_2$ surface. We show that the following factors affect which self-consistent solution a $\text{DFT}+U$ calculation reaches: (i) the magnitude of $U$; (ii) initial correlated orbital occupations; (iii) lattice geometry; (iv) whether lattice symmetry is enforced on the charge density; and (v) even electronic mixing parameters. These various solutions may differ in total energy by hundreds of meV per atom, so identifying or approximating the ground state is critical in the $\text{DFT}+U$ scheme. We propose an efficient $U$-ramping method for locating low-energy solutions, which we validate in a range of test cases. We also suggest that this method may be applicable to hybrid functional calculations.

Figure 1

(Color online) The relative energies, as a function of $U$, of multiple solutions in $\text{DFT}+U$ for fluorite ${\text{UO}}_2$. The curves represent different initial correlated orbital occupation matrices that were applied to the calculations. These curves are plotted at two different unit-cell volumes: (a) corresponds to $a=1.02a_{expt}$ and (b) corresponds to $a=0.98a_{expt}$. The band gaps for the multiple solutions at the larger and smaller volumes are given in (c) and (d), respectively.

Figure 2

(Color online) The behavior of our $U$-ramping method at small values of $U$ in fluorite ${\text{UO}}_2$. (a) The number of fractionally occupied bands in calculations using the same enumerated occupation matrices as in Fig. 1, as compared to the iterative $U$-ramping scheme. (b) The relative energies of the calculations in (a); our method finds a low-energy solution at all $U$.

Figure 3

Schematic illustration of the principle underlying the proposed method. For simplicity, the horizontal axis denotes the occupation of a single localized orbital (out of the numerous localized orbitals present in the system). As the Hubbard $U$ parameter is gradually turned on, the energy surface (as a function of orbital occupation) goes from being convex (bottom curve, corresponding to conventional DFT) to being concave (top curve, corresponding to $\text{DFT}+U$). For small values of $U$, the $\text{DFT}+U$ global energy minimum is easy to find and gradually converges to the true $\text{DFT}+U$ global energy minimum (filled circles) as $U$ increases to its full correct value. The concavity of the $\text{DFT}+U$ curve (at the nominal value of $U$) would have otherwise necessitated a combinatorial search among the possible local minima (the true energy surface is a high-dimensional object with considerably more local minima).

Figure 4

(Color online) Results of our $U$-ramping scheme, as compared to manual enumeration of possible correlated orbital occupation matrices, for a variety of test cases using nominal values of $U$ and $J$ from the literature: (a) fluorite AFM ${\text{UO}}_2$; (b) fluorite AFM ${\text{UO}}_2$ with symmetry disabled; (c) rocksalt AFM CoO; (d) rhombohedrally distorted CoO; and (e) rocksalt AFM NiO with symmetry disabled.

Figure 5

(Color online) $\text{DFT}+U$ solutions for a ${\text{CeO}}_2$ (111) surface containing an O vacancy, obtained by two choices of electronic mixing parameters and by our $U$-ramping method. The difference between spin-up and spin-down charge densities is shown for the two present ${\text{Ce}}^{3+}$ ions. Internal degrees of freedom were relaxed.