## 2 Computational Methods  
We show results derived from density functional theory that was developed in the mid-1960s [42; 43]. Density functional theory is actually correct, but it relies on practical approximations. The central assumption relates to the electrons exchange and correlation interactions. For the actinide metals, we employ the generalized gradient approximation that is more proper than the previously suggested local density approximation [6]. Specifically, we are applying the Perdew-Burke-Ernzerhof (PBE) parameterization of the generalized gradient approximation (GGA) [44]. To our knowledge, there is no established, parameter-free, functional that results in larger atomic volumes in the actinides than that of the GGA. Therefore, there is no choice to replace the GGA to expand the volumes, while a +\(U\) approach can accomplish this objective and has often been applied for the actinides, as discussed in the introduction. The actinide metals exhibit interesting crystal structures, beginning with an innocent face-centered cubic for Th with a gradual escalation of complexity to tetragonal (Pa), orthorhombic (U and Np), and finally culminating with a monoclinic structure (\(\alpha\)-plutonium). These structures are described in the compilation by [PERSON] [45]. For best possible accuracy, when calculating the equilibrium volumes, we have optimized these crystal structures (not the cubic phase) so that the model gives the lowest total energy (structural relaxation). This relaxation has generally a rather weak effect on the actinide volumes, except for \(\alpha\)-plutonium where it is significant. In this instance, a structural relaxation contracts the volume substantially relative to a calculation where the experimental (room-temperature) structure is assumed [40; 46]. It has been reported that orbital-orbital interactions, absent in conventional DFT, is important for plutonium [14; 40; 47]. This interaction is efficiently included in an orbital-polarization scheme for some of the calculations presented in this study. In terms of the technical details of the computations, we are applying three techniques. The bulk of the calculations are performed with the pseudopotential plane-wave Vienna ab initio simulation package (VASP) code. In VASP, the computational scheme relies on the GGA with PBE exchange and correlation functional [44] and the projector-augmented-wave (PAW) method with plane-wave basis set as implemented in VASP [48; 49; 50]. The computational parameters consist of an energy cut-off of 450 eV, energy convergence of 100 eV, and a Monkhorst-Pack scheme \(k\)-point mesh of typically \(12\times 12\times 12\) for a 1-atom cell and correspondingly less for multi-atom cells. In the cases where we perform DFT + \(U\) computations, we employ the technique suggested by [PERSON] et al. [51] for the 5f electrons with an effective \(U\) = 0.4 eV. The spin-orbit interaction and orbital polarization are included, allowing for non-collinear spins. For the spin-orbit coupling, we have introduced a correction that retains the scalar-relativistic treatment for the p states. Orbital polarization has been implemented similarly to the scheme proposed by [PERSON] et al. [52]. In this approach, the DFT total-energy functional, EDFT, is amended by a quadratic term favoring orbital polarization,  
\[E_{DFT\ +\ OP}\ =\ E_{DFT}\ -\ \frac{1}{2}\sum_{I}E_{I}^{3}[\langle L_{I,z}^{f }\rangle]^{2} \tag{1}\]  
where \(E_{I}^{3}\) is the Racah parameter that is a linear combination of the radial Slater integrals F\({}^{2}\), F\({}^{4}\), and F\({}^{6}\) over f-electron orbitals within a muffin-tin sphere centered on site I. For consistency and simplicity, we use a constant Racah parameter that is an average of spin-decomposed values obtained from all-electron calculations (see below). \(L_{I,z}^{f}\) is the \(z\) component of the angular-momentum operator within the same sphere operating on the f states (5f in this case).