## Computational Methods  
### Electronic Structure Method  
The electronic structure results we present here are obtained from density-functional theory. This theory is in principle exact, but it relies on practical simplifications and the foundational one is that for the electron interactions. For the actinides [16; 17] the most robust formulation is the so-called generalized gradient approximation (GGA) [18]. Regarding the DFT implementation, we utilize an all-electron full-potential linear muffin-tin orbital (FPLMTO) method that is well known [19] and detailed [20]. FPLMTO makes no simplifications for the core electrons that occupy deeper energy levels than the valence states, as opposed to the pseudopotential method. The core approximation made in pseudopotential methods is computationally efficient but less accurate. Basis functions, electron densities, and potentials are not subject to geometrical approximations and they are expanded in spherical harmonics inside non-overlapping (muffin-tin (MT)) spheres centered at each atomic site and in Fourier series in the region between these MTs. There is a need to determine the MT sphere radius (r\({}_{\rm MT}\)) and here it is selected so that r\({}_{\rm MT}\sim 0.7\)r\({}_{\rm WS}\), where r\({}_{\rm WS}\) is the Wigner-Seitz (atomic-sphere) radius. The radial part of the basis functions inside the MT spheres is calculated from a wave equation for the \(l=0\) component of the potential that includes all relativistic corrections, including spin-orbit coupling for d and f states but not for the p states, as explained in [17]. The orbital-orbital interaction (orbital polarization (OP)) is only operating on the uranium f states [21]. This interaction is not explicit in conventional density-functional theory but has been shown to be important for some f-electron systems, particularly plutonium [17; 22]. The crystal structure of the orthorhombic UI\({}_{3}\) compound is optimized with respect to the total energy. In other words, all degrees of freedom, given the symmetry of the crystal, have been varied in order to produce the lowest (ground-state) total energy. This procedure is particularly important for an accurate calculation of the energy of formation of the compound. For this we also similarly relax the orthorhombic structure of \(\alpha\)-uranium. Typically, structural relaxation from electronic structure calculations are guided by the forces on the atoms. However, in the presence of strong spin-orbit interaction and orbital polarization these forces are difficult to compute accurately [23]. The relaxation is therefore here accomplished without an explicit need to compute forces on the atoms but by separately moving them small amounts from their equilibrium positions, followed, iteratively, by optimizing the unit cell (_b/a_ and _c/a_ axial ratios). Although a somewhat tedious and demanding scheme, it is very robust for relaxing crystal structures and it was developed for our study of \(\epsilon\)-uranium phonons [23]. The obtained lattice constants for UI\({}_{3}\) are \(a\) = 4.48 A, \(b\) = 15.94 A, and \(c\) = 10.26 A. Both \(a\) and \(c\) are good while the \(b\) lattice constant is overestimated relative to the experimental data (4.33, 14.01, and 10.0 A) [2]. We should mention that atomic forces can be evaluated resulting from this structural-relaxation procedure (they approach zero at the equilibrium). Such forces were used for the calculation of temperature-dependent phonons for uranium mononitride recently [15], but for the low-symmetry orthorhombic phase relevant here, the applied phonon method requires further development. For determining the formation enthalpy, we calculate the energies of UI\({}_{3}\), \(\alpha\)-uranium, and the iodine atom. For these energies to be consistent, the calculations need to be set up as similarly as possible, using the same computational tool, and for the iodine atom we thus construct a simple-cubic cell with the atoms far apart.