First-principles calculations of the epsilon phase of solid oxygen

The crystal, electronic and magnetic structures of solid oxygen in the epsilon phase have been investigated using the strongly constrained appropriately normed (SCAN) + rVV10 method and the generalized gradient approximation (GGA) + vdW-D + U method. The spin-polarized SCAN + rVV10 method with an 8-atom primitive unit cell provides lattice parameters consistent with the experimental results over the entire pressure range, including the epsilon-zeta structural phase transition at high pressure, but does not provide accurate values of the intermolecular distances d1 and d2 at low pressure. The agreement between the intermolecular distances and the experimental values is greatly improved when a 16-atom conventional unit cell is used. Therefore, the SCAN + rVV10 method with a 16-atom unit cell can be considered the most suitable model for the epsilon phase of solid oxygen. The spin-polarized SCAN + rVV10 model predicts a magnetic phase at low pressure. Since the lattice parameters of the predicted magnetic structure are consistent with the experimental lattice parameters measured at room temperature, our results may suggest that the epsilon phase is magnetic even at room temperature. The GGA + vdW-D + U (with an ad hoc value of Ueff = 2 eV at low pressure instead of the first-principles value of Ulreff ~ 9 eV) and hybrid functional methods provide similar results to the SCAN + rVV10 method; however, they do not provide reasonable values for the intermolecular distances.


Lattice parameters calculated with different DFT functionals
Supplementary Figure S1. Lattice parameters a, b, and c and the angle β calculated with the LDA, PBE and BLYP functionals. Figure S2. Lattice parameters a, b, and c and the angle β calculated with the meta-GGA (M06L) 1 , GGA+semi-empirical van der Waals 2 (vdW-D) and van der Waals functional (vdW-DF1) 3 methods. Figure S3. Lattice parameters β, a, b, and c calculated with the GGA (PBE), vdW-DF1 functionals 3 , revised vdW-DF2 functionals 7 , and the semi-empirical GGA+vdW-D 2 method. The value of Ueff was changed from 0 eV to 9.6 eV.

Supplementary
In this Supplement, we compare the lattice parameters a, b, and c and the angle β calculated with different methods. Figure S1 shows the results from the LDA, PBE, and BLYP functionals. Figure S2 shows the data calculated with meta-GGA (M06L) 1 , the van der Waals functional 3 (vdW-DF1) and the GGA + semi-empirical van der Waals 2 (vdW-D) method. Figure S3 shows the lattice parameters calculated with the GGA (PBE), van der Waals functionals vdW-DF1 3 and revised van der Waals functionals vdW-DF2 7 , and the semi-empirical GGA+vdW-D 2 method. The value of Ueff was changed from 0 eV to 9.6 eV.

Enthalpy and magnetization of optimized structures calculated from the spin-polarized GGA+vdW-D+U method
Supplementary Figure S4. The difference in enthalpies between the spin-polarized and non-spin-polarized calculations and absolute magnetizations calculated for (a) group A and (b) group B at Ueff = 1, 2, 4, and 9.6 eV.
In Figure S4(a-b), we compare the enthalpies of optimized structures obtained from the non-spin-polarized calculations and spin-polarized calculations for groups A and B, respectively. The difference in enthalpy can be considered as the magnetic energy. Because we compare groups A and B with non-spin-polarized structures, the zero energies are the same for groups A and B. The absolute magnetizations are shown on the right-hand side. In groups A and B, the enthalpies of the magnetic structures are always lower than those for the non-magnetic structures before the epsilon-zeta transition occurs. The structures of group A are as stable as those of group B at low pressure. Moreover, the range of the low-pressure magnetic epsilon phase increases as Ueff increases. This means that the Hubbard U not only enhances the localization of the electron but also stabilizes the magnetic structure. Figure S5 shows the radial distribution function of the O-O distance g(r) simulated for the structure at 10 GPa and at 5 K and 300 K. The time step is 0.5 fs. We used the NpT ensemble with a Langevin thermostat, which was implemented in the VASP code 9 . The functional function was PBE without vdw-D and Hubbard U. The structures were simulated at a constant pressure of 10 GPa and constant temperatures of 5 K and 300 K. The structure consisted of 16 oxygen atoms in a conventional unit cell. Our non-spin-polarized PBE molecular dynamics (MD) simulations at 5 K show that d1 in the MD simulation is extended from 2 Å to 2.2 Å, which is closer to the measurement value at 2.35 Å. The MD simulation at 300 K shows a broadening distribution of d1 varying from 2 Å to 2.2 Å. This preliminary result suggests that there may be many meta-stable structures co-existing at low pressure and that the d1 values obtained from the measurements were an average value of these meta-stable structures. The MD simulation, thus, is more suitable to find a global minimum of solid oxygen than a conventional geometry optimization. Further spin-polarized molecular dynamics simulations will be carried out in the next study.

Molecular dynamics simulation
Supplementary Figure S5. Radial distribution function of the O-O distance g(r) simulated for the structure at P = 10 GPa with T = 5 K and 300 K. The time step is 0.5 fs.

More sophisticated methods
Finally, the poor description of DFT in many-body interactions 5,8 should be considered. In reference 8 , the authors used a multiconfigurational complete active space second-order perturbation (CASPT2) to estimate the singlet state of an (O2)4 cluster and obtained reasonable results for d1 and d2. However, a similar calculation for a periodic structure is still missing. Other studies suggest the necessity of a multi-reference wave function 5 or a van der Waals density functional with spin-polarized-dependent gradient correction 4,6 . Higher levels of calculation beyond the meta-GGA may be a key to solving this problem. Recent study of transition metal monoxides 10 shows that the SCAN+vdW+U with the self-consistent U calculated based on SCAN predicts good ground state of FeO. The SCAN+vdW+U may be a good choice for the calculation of solid oxygen's ground state.