Dirac's equation and its implications for density functional theory based calculations of materials containing heavy elements

Electronic structure calculations based on density functional theory (DFT) give quantitatively accurate predictions of properties of most materials containing light elements. For heavy materials, and in particular for $f$-electron systems, DFT based methods can fail both qualitatively and quantitatively for two distinct reasons: their failure to describe confinement effects arising from localized $f$-electron behavior and their incomplete or approximate treatment of relativity. In addition, different methods for incorporating relativistic effects, which give identical results in most light materials, can give different predictions in heavy elements. In order to develop a quantitative capability for calculating the properties of these materials, it is essential to separate the predictions of the underlying equations from the uncertainty introduced in approximations used in computation. Working toward that goal, we have developed a code, called dirac-fp, which is based directly on solving the Dirac-Kohn-Sham equations and uses the full-potential linear muffin-tin orbital (FP-LMTO) approach to electronic structure. In order to assess the performance of dirac-fp, we perform calculations on three different face-centered cubic materials using different approximate treatments of relativity: the scalar relativistic (SR) approach commonly used in most solid-state DFT codes, the scalar relativistic plus spin-orbit coupling corrections (SR+SO) approach which includes spin-orbit coupling self-consistently using the SR states inside the muffin tins, and the Dirac-Kohn-Sham (Dirac) approach implemented in dirac-fp. Performing calculations on thorium, in which relativistic effects should be strong, aluminum, in which relativistic effects should be negligible, and gold, in which relativistic effects play an intermediate role, we find that the Dirac approach is able to provide theoretically consistent results in the electronic structure and ground-state properties across all three materials.

Figure 1

Schematic of the SR, $\mathrm{SR}+\mathrm{SO}$, and Dirac methods within the FP-LMTO methodology. Energy ($E$) vs radial distance from the center of a muffin-tin sphere ($r$) is shown, with a representation of full effective potential $v_{\mathrm{eff}}\left(\mathbf{r}\right)$ in the solid black line. All core states are treated using the Dirac equation. The core states and muffin-tin radii are chosen so that the core states do not extend spatially outside of the muffin tins and therefore do not interact with other core states. The valence states, which can correspond to states other than the valence electrons, are split spatially between a region inside the muffin tins and outside the muffin tins (interstitial region). The SR and $\mathrm{SR}+\mathrm{SO}$ methods use the SR basis for both spatial regions, while the SO term is only applied within the muffin tins where it is possible to associate the orbital angular momentum with the state. The Dirac method solves the Dirac equation throughout all regions of the cell.

Figure 2

Energy-volume curves for thorium using the SR (first row), $\mathrm{SR}+\mathrm{SO}$ (second row), and Dirac (third row) methods and the PBE exchange-correlation functional. The first column shows results using a double basis and the second column shows results using a full basis. In each subplot, the fixed-volume-fraction (FF) curve is shown as a solid line and the fixed-radius (FR) curve is shown as a dashed line. The zero of energy $E_0$ is different for each subplot, and is taken to be the minimum energy among both curves in the subplot, so that the curves can be compared in a more straightforward manner. The experimental value of $V_0=220.0{\mathrm{bohrs}}^3$ (Ref. [73]) is used to normalize the $x$ axis. The bulk modulus for thorium is 54 GPa (Ref. [74]).

Figure 3

Electronic density of states of thorium for SR, $\mathrm{SR}+\mathrm{SO}$, and Dirac methods. The energy-level splitting of the $6p$ states appears only for the $\mathrm{SR}+\mathrm{SO}$ and Dirac methods. All results are computed using PBE at the equilibrium volume and FF ratio $S=0.85$.

Figure 4

Electronic density of states for thorium using $\mathrm{SR}+\mathrm{SO}$ and Dirac, but with different FF radii. The $S=0.50$ results are shown as dashed lines, while the solid lines use $S=0.85$. The lower-energy $\mathrm{SR}+\mathrm{SO}$ $6p_{1/2}$ states exhibit a large shift in energy upon changing $S$, while the Dirac method shows no such shift. All curves use PBE.

Figure 5

Energy-volume curves for aluminum for the SR, $\mathrm{SR}+\mathrm{SO}$, and Dirac methods. The zero of energy is the same for all curves, so the curves overlap in shape and also in value. For aluminum, $V_0=109.5{\mathrm{bohrs}}^3$ [73].

Figure 6

Energy-volume curves for gold using the SR, $\mathrm{SR}+\mathrm{SO}$, and Dirac methods with different choices for which states are chosen as valence states. In each subplot, the solid line shows the case where the $5s^{2}4f^{14}5p^6$ semicore and valence electrons are treated in the valence, while the darker dashed line shows the case where the $5s^{2}4f^{14}5p^6$ semicore electrons are treated in the core and only the valence electrons are treated as valence states. Each subplot has a unique value of $E_0$. For the SR and $\mathrm{SR}+\mathrm{SO}$ methods, the $\mathrm{semicore}+\mathrm{valence}$ treatment is significantly higher in energy than the valence-only results, requiring the values to be plotted on two separate $y$ axes (valence-only values on the left, $\mathrm{semicore}+\mathrm{valence}$ values on the right). Arrows point from each curve to the appropriate axis. In the Dirac case, the same $y$ axis is used for both curves since the energy values are nearly identical. The range of all $y$ axes are the same (104 mRy). Notably, the predicted $V_0$ for valence-only vs $\mathrm{semicore}+\mathrm{valence}$ are the same for the SR and Dirac methods, while the $\mathrm{SR}+\mathrm{SO}$ method predicts very different values of $V_0$, depending on which states are treated in the valence. All curves are computed using PBE and FF ratio $S=0.85$. For gold, $V_0=113.0{\mathrm{bohrs}}^3$ (Ref. [73], computed with a zero-point anharmonic expansion correction).

Figure 7

Electronic densities of $5p$ states for gold at different unit-cell volumes using the $\mathrm{semicore}+\mathrm{valence}$ treatment and the PBE exchange-correlation functional. SR and Dirac FF curves are shown in blue and green, respectively, while the $\mathrm{SR}+\mathrm{SO}$ method shows curves for the FF (orange lines) and FR (brown dashed lines) treatments. The muffin-tin radii at $V/V_0=0.8$ are the same for the FR and FF methods, leading identical densities of states. The SR $5p$ states are split into a $5p_{1/2}$ (below $-65$ eV) and $5p_{3/2}$ (above $-55$ eV) when adding the SO term, though the $\mathrm{SR}+\mathrm{SO}$ $5p_{1/2}$ states do not match the Dirac $5p_{1/2}$ states in energy, owing primarily to the different $r\to 0$ behavior of the SR and Dirac $5p_{1/2}$ radial functions. For gold, $V_0=113.0{\mathrm{bohrs}}^3$ (Ref. [73], computed with a zero-point anharmonic expansion correction).