non-parabolicity of the Birch–Murnaghan curve away from its minimum. If we
consider a very small volume range, the curve will be very close to parabolic,
and we therefore expect a large error on $B_{1}$ since the fit has very little
information on the non-parabolic behavior. For larger volume ranges, the curve
starts to deviate significantly from a parabola, thus providing more
information to the fitting algorithm on the actual value of $B_{1}$, in turn
resulting into a smaller relative error on $B_{1}$ vs. $B_{0}$ or $V_{0}$.

  * •

The stability of the fit, especially on $B_{1}$, is significantly affected by
the choice of fitting algorithm. For instance, we realized that if one uses
the optimize.curve_fit subroutine of SciPy (https://www.scipy.org), which is
_not_ the algorithm used in this work, the choice of the fitting starting
point is very important, and we also observe that iterating the procedure a
few times (using the results of the previous step as starting points for the
next fit) improves the stability. Instead, the function used in this work
(that is the same also used in Ref. 29, 30) is a non-iterative fitting
algorithm that proves to be much more robust.

Figure S3.1: Histograms of the typical error ratios of $B_{0}$ vs. $V_{0}$
(${\bar{\eta}_{B_{0}}}/{\bar{\eta}_{V_{0}}}$, top row) and of $B_{1}$ vs.
$V_{0}$ (${\bar{\eta}_{B_{1}}}/{\bar{\eta}_{V_{0}}}$, bottom row) for the
unaries set (left column) and the oxides set (right column). The simulations
were run for $N_{s}=100$ random samples and a standard deviation on the error
of the energy on the datapoints of $n_{\sigma}=10^{-5}$ eV. The histograms
indicate that the error of $B_{0}$ ($B_{1}$) is approximately 20 (400) times
larger than the error on $V_{0}$, justifying our choice of weights for the
metric $\nu$.

##  S4 Reference all-electron results for $V_{0}$, $B_{0}$ and $B_{1}$

This section reports the complete reference dataset of the EOS parameters
obtained with the two all-electron codes FLEUR and WIEN2k, and the absolute
value of their percentage difference (that we indicate with $\eta$). Moreover,
it reports the averaged parameters among them, that constitutes our reference
average dataset presented in this manuscript. Data is divided in 10 tables,
one for each crystal structure (4 unaries and 6 oxides). The agreement for
$V_{0}$ is within 0.3% for all materials except Cs2O5 (0.323%), Fr2O5
(0.645%), Ra2O5 (0.333%), NeO3 (0.302%) and RbO3 (0.343%). Not surprisingly,
these are 5 crystals with very small bulk moduli B0: it has been shown in Ref.
98 that the error in the volume scales inversely with the value of the bulk
modulus.

Parameters are expressed per formula unit (see also SI Sec. S1). Note that,
for X2O3 and for X2O5, the primitive cell has twice the number of atoms (10
and 14, respectively) than the number of atoms in the formula unit (5 and 7,
respectively). This is reflected in a factor of 0.5 in the volumes reported in
this table with respect to the volume of the unit cells in the input files
available in Ref. 62.

##  S5 Simulation parameters for FLEUR and WIEN2k

The AE codes used for producing the reference results in this work implement
the (linearized) augmented-plane-wave + local orbitals ((L)APW+LO) method.
This method is named after a family of different, but related basis sets, each
having different characteristics and demands on cutoff parameters. It is based
on a partitioning of the unit cell into muffin-tin (MT) spheres at the atom
positions and an interstitial region in between the spheres. The MT sphere
radii hereby have to be adapted to each structure to avoid overlapping
spheres, but are typically chosen to be nearly touching, because large spheres
reduce demands on the LAPW basis set cutoff parameter. The exact choice of the
radii typically differs between (L)APW+LO codes.

The choice of the core/valence electron separation is related to the MT sphere
radii. To avoid instabilities in the calculations and to obtain precise
results, the extent of the core-electron states beyond the MT sphere boundary
has to be limited. The AE codes employed here share the same core/valence
separation for many structures, but differ in this choice for a considerable
amount of other structures. Taking into account the overall excellent
agreement between the results from the two codes it can be deduced that there
is no significant dependence of the results on the differing description of
physics for core and valence electrons, as long as the core/valence separation
is within a reasonable range.

The valence electrons are represented by the (L)APW+LO basis chosen by the
respective code. The different possible choices of basis sets nevertheless
share a common set of parameters with which they are specified. These are
cutoff parameters for the basis set size and the angular momentum expansion in
the MT spheres, as well as energy parameters defining the linearization
centers for each atom. Additionally, the local orbital setup has to be defined
for each MT sphere. The two AE codes employed in this work make use of
different kinds of (L)APW+LO basis sets and their choices for the setup
aspects discussed here differ strongly. A sketch on the parameter setup
recipes is provided in the following two subsections. A detailed list of the
setup parameter choices for each structure is available in the supplementary
data, as referenced in the respective sections below.

###  S5.1 FLEUR

The FLEUR code[31, 32] is an open-source implementation of the all-electron
full-potential linearized augmented-plane-wave (FLAPW) method[99, 100]. The
calculations with this code make use of a conventional LAPW basis in
combination with local orbitals (LOs) of different types to describe semicore
states[101] and to eliminate the linearization error for the valence
states[102, 103, 104]. The employed setup profile defines global parameters,
identical for all calculations and element-specific parameters. Further
parameters are automatically adapted to each investigated structure.

In the context of this profile, global parameters are the reciprocal LAPW
basis set cutoff parameter $K_{\text{max}}=5.0\leavevmode\nobreak\ a_{0}^{-1}$
and the plane-wave cutoff parameters
$G_{\text{max}}=G_{\text{max,XC}}=25.0\leavevmode\nobreak\ a_{0}^{-1}$, where
$a_{0}$ is the Bohr radius. The latter cutoff parameter covers the expansions
of the plane-wave part of the density and the potential, the interstitial-
region indicator function, and the exchange-correlation contribution to the
potential. The Fermi–Dirac smearing and the k-point density are also fixed to
the common choices of this work. Element-specific setup aspects include the
core-valence separation of the electron states and the LO setup. The radii of
the MT spheres centered on atom $\alpha$, $R_{\text{MT}}^{\alpha}$, are
adapted to the smallest unit cell within an equation-of-states (EOS) workflow.
For this, element-specific initial MT sphere radii are expanded to cover up to
about $92\%$ of the distance between the atoms, with a limit of a maximal MT
radius of $2.66\leavevmode\nobreak\ a_{0}$. The procedure also implies an
adaption of the angular momentum cutoffs in the spheres to
$l_{\text{max}}^{\alpha}\approx K_{\text{max}}R_{\text{MT}}^{\alpha}$. A
detailed description of the different parameters is available in the FLEUR
user guide[31, 32].

The parameter profile is the result of an iterative refinement process aiming
at precision and stability on the basis of the structures investigated in this
work, i.e., structures with a wide range of neighboring atom distances.
However, it is not designed to provide an absolute convergence of the total
energy, which may still be affected by structure-adapted numerical parameters
like the MT radii. The determination of these parameters for each structure on
the basis of the smallest unit cell solves this issue for the EOS workflow.
Comparing the total-energy results from the EOS workflows for different
structures may also yield reasonable numbers, but these do not reflect the
precision capabilities of the code when used in an adequate way to perform
such a comparison, see also discussion in SI Sec. S18.

With the exception of the k-point integration-mesh generation, the used
parameter profile is implemented in the openly available releases of FLEUR
starting with the MaX-R6.1 release. The results presented here have been
obtained using the development version as of 2022/03/31. The profile is
employed by invoking the FLEUR input generator with the command line option
"-profile oxides_validation". Beyond using this profile, the AiiDA common
workflows package[36] protocol "verification-pbe-v1" in combination with
AiiDA-FLEUR[105, 106] also sets the k-point integration mesh.

The resulting parametrization for each structure is discussed in the all-
electron-setups folder of the supplementary data available in Ref. 62, in
particular in the files

  * •

setup-oxides-verification-PBE-v1-fleur.json

  * •

setup-unaries-verification-PBE-v1-fleur.json

that are documented in the file all-electron-data.md.

In an ongoing effort the data from this work is also used to define further
FLEUR parameter profiles for different precision levels and computational
effort. This is done by reducing in a controlled way cutoff parameters and
changing other aspects of the setup and relating the corresponding results to
those presented here.

###  S5.2 WIEN2k

The WIEN2k calculations[33, 34] employ the (linearized) augmented plane wave +
local orbitals ((L)APW+lo) method[107] with additional local orbitals[101]
(LOs) for states with an energy 1 Ry below the Fermi energy (semicore states)
and high-derivative LOs (HDLOs)[108] for all “chemical” $l$ values (except
when there is already a semicore LO). The APW+lo basis set is used for all
“chemical” $l$ values ($s$, $p$, $d$ or $f$ \- depending on the atom), while
LAPW is used for higher angular momentum up to $l_{\text{max}}$ = 10.

The WIEN2k calculations have been initialized (see the WIEN2k users guide[33])
for the smallest volume of each case using:

init_lapw -b -prec 3 -nokshift -fermits 0.0045 -red 3 -numk -1 0.0317506

For subsequent volumes we use:

init_lapw -b -prec 3 -nokshift -fermits 0.0045 -red Element:RMT -numk 0 kx ky
kz -fft ix iy iz

where Element:RMT, kx ky kz, and ix iy iz are inserted from the output of the
first volume to ensure identical parameters.

This high-precision setup limits the maximal atomic sphere radius
$R_{\text{MT}}^{\alpha}$ to $2.35\leavevmode\nobreak\ a_{0}$, but otherwise
sets the sphere sizes automatically depending on nearest neighbor distances
and type of atom (largest for $f$ elements, intermediate for $d$ elements and
smallest for $sp$ elements). Also the plane wave cutoff $R_{\text{MT,
min}}^{\alpha}K_{\text{max}}$ is set automatically depending on the type of
atom and the smallest atomic sphere radius $R_{\text{MT, min}}^{\alpha}$ and
varies from 7.08 (H2O5) to 11. All states with an energy above $-6$ Ry or with
a charge density of more than 0.01 $e^{-}$ outside the atomic sphere are
considered as valence states and treated scalar-relativistically, while lower
energy states are considered as core and solved numerically with a radial
symmetric Dirac equation. Note that with this choice the definition of core
states for an element may change depending on its $R_{\text{MT}}^{\alpha}$ and
this makes the calculation of formation energies in certain cases unrealistic
as discussed also in SI Sec. S18. An SCF cycle was considered converged when
both the change in total energy was less than $10^{-6}$ Ry and the change in
the electron charge density within $R_{\text{MT}}^{\alpha}$ was less than
$10^{-6}\leavevmode\nobreak\ e$.

The charge density and potential inside spheres is expanded into lattice
harmonics up to $L_{\text{max}}$ = 6 and for the non-spherical Hamiltonian
matrix elements the angular momentum of the wave functions is restricted to
$l^{ns}_{\text{max}}$ = 8. In the interstitial region the density/potential is
expanded into a Fourier series with cutoff parameter
$G_{\text{max}}\leavevmode\nobreak\ =\leavevmode\nobreak\
25\leavevmode\nobreak\ a_{0}^{-1}$ (except for alkali metals, noble gases, and
Hg, where $G_{\text{max}}\leavevmode\nobreak\ =\leavevmode\nobreak\
40\leavevmode\nobreak\ a_{0}^{-1}$).

The resulting parametrization for each structure is discussed in the all-
electron-setups folder of the supplementary data available in Ref. 62, in
particular in the files

  * •

setup-oxides-verification-PBE-v1-wien2k.json

  * •

setup-unaries-verification-PBE-v1-wien2k.json

that are documented in the file all-electron-data.md.

##  S6 Dependence of the metrics on the size of the simulation cell and on
bond stiffness

The equation of state is often expressed in terms of the absolute energy $E$
and volume $V$ of the simulation cell, but other times quantities "per-
formula-unit" or "per-atom" are considered. We want here to demonstrate that
the new metrics $\varepsilon$ and $\nu$, introduced here, are intrinsic
quantities, i.e., they are independent of the simulation cell size (while the
original $\Delta$ metric is extensive).

Dependence on the number of atoms in the simulation cell: We first show that
the original metric $\Delta$ is an extensive quantity that depends linearly on
the number of atoms (or, equivalently, on the volume) of the simulation cell.

Let us consider a supercell where the number of atoms is increased by a factor
of $C$. The volume and the total energy (assuming $E_{0}=0$) will scale
accordingly as

$E^{\prime}=CE,\quad V^{\prime}=CV.$ (S5)

The integrand in Eq. (3) for the $\Delta$ metric will scale as

$[E^{\prime}_{a}(V^{\prime})-E^{\prime}_{b}(V^{\prime})]^{2}\leavevmode\nobreak\
dV^{\prime}=C^{3}[E_{a}(V)-E_{b}(V)]^{2}\leavevmode\nobreak\ dV$ (S6)

and the denominator as

$V^{\prime}_{M}-V^{\prime}_{m}=C(V_{M}-V_{m}).$ (S7)

The effect of increasing the number of atoms on the $\Delta$ metric is

$\displaystyle\Delta^{\prime}(a,b)$
$\displaystyle=\sqrt{\frac{1}{V^{\prime}_{M}-V^{\prime}_{m}}\int_{V^{\prime}_{m}}^{V^{\prime}_{M}}[E^{\prime}_{a}(V^{\prime})-E^{\prime}_{b}(V^{\prime})]^{2}\leavevmode\nobreak\
dV^{\prime}}=\sqrt{\frac{C^{3}}{C(V_{M}-V_{m})}\int_{V_{m}}^{V_{M}}[E_{a}(V)-E_{b}(V)]^{2}\leavevmode\nobreak\
dV}$ (S8)
$\displaystyle=C\sqrt{\frac{1}{V_{M}-V_{m}}\int_{V_{m}}^{V_{M}}[E_{a}(V)-E_{b}(V)]^{2}\leavevmode\nobreak\
dV}=C\Delta(a,b).$

Thus, $\Delta$ scales linearly with the number of atoms in the simulation
cell. This shortcoming is typically addressed by computing the $\Delta$ metric
renormalized per atom.

For the $\varepsilon$ metric we proceed in the same manner to demonstrate,
instead, its independence from the simulation cell size.

$\displaystyle\varepsilon^{\prime}(a,b)$
$\displaystyle=\sqrt{\frac{\int_{V^{\prime}_{m}}^{V^{\prime}_{M}}[E^{\prime}_{a}(V^{\prime})-E^{\prime}_{b}(V^{\prime})]^{2}\leavevmode\nobreak\
dV^{\prime}}{\sqrt{\left(\int_{V^{\prime}_{m}}^{V^{\prime}_{M}}[E^{\prime}_{a}(V^{\prime})-\langle
E^{\prime}_{a}\rangle]^{2}\leavevmode\nobreak\
dV^{\prime}\right)\leavevmode\nobreak\
\left(\int_{V^{\prime}_{m}}^{V^{\prime}_{M}}[E^{\prime}_{b}(V^{\prime})-\langle
E^{\prime}_{b}\rangle]^{2}\leavevmode\nobreak\ dV^{\prime}\right)}}}$ (S9)
$\displaystyle=\sqrt{\frac{C^{3}\int_{V_{m}}^{V_{M}}[E_{a}(V)-E_{b}(V)]^{2}\leavevmode\nobreak\
dV}{\sqrt{\left(C^{3}\int_{V_{m}}^{V_{M}}[E_{a}(V)-\langle
E_{a}\rangle]^{2}\leavevmode\nobreak\ dV\right)\leavevmode\nobreak\
\left(C^{3}\int_{V_{m}}^{V_{M}}[E_{b}(V)-\langle
E_{b}\rangle]^{2}\leavevmode\nobreak\ dV\right)}}}$
$\displaystyle=\sqrt{\frac{\int_{V_{m}}^{V_{M}}[E_{a}(V)-E_{b}(V)]^{2}\leavevmode\nobreak\
dV}{\sqrt{\left(\int_{V_{m}}^{V_{M}}[E_{a}(V)-\langle
E_{a}\rangle]^{2}\leavevmode\nobreak\ dV\right)\leavevmode\nobreak\
\left(\int_{V_{m}}^{V_{M}}[E_{b}(V)-\langle
E_{b}\rangle]^{2}\leavevmode\nobreak\ dV\right)}}}=\varepsilon(a,b).$

Since $C$ cancels out, we have proven that $\varepsilon$ is independent of the
number of atoms in a simulation cell considered.

Analogously, it is also easy to see that $\nu$ is also an intrinsic quantity,
independent of the number of atoms in the simulation cell: indeed, $\nu$ is
defined as a function of the relative errors of the parameters $V_{0}$,
$B_{0}$ and $B_{1}$, that are all intrinsic quantities.

We stress that this means that, while for $\Delta$ is recommended to report it
normalized (e.g., per atom, or per formula unit), $\varepsilon$ and $\nu$
should _not_ be normalized.

Sensitivity to the value of the bulk modulus: Let us compare the metrics
obtained comparing results for two different materials. We assume that the
first material has bulk modulus $B_{0}$ and the second is identical except for
its bulk modulus, that is scaled by factor of $C$, i.e.,
$B_{0}^{\prime}=CB_{0}$. We assume that there is no other difference between
the two materials (same $V_{0}$ and $B_{1}$) and we are considering the same
simulation volume (or number of atoms in the simulation cell). The total
energy of the second material will then scale as

$E^{\prime}=CE$ (S10)

according to Eq. (1) and assuming $E_{0}=0$ (minimum energy of both materials
have been shifted to zero). The integrand in Eq. (3) for the $\Delta$ metric
scales as

$[E^{\prime}_{a}(V)-E^{\prime}_{b}(V)]^{2}\leavevmode\nobreak\
dV=C^{2}[E_{a}(V)-E_{b}(V)]^{2}\leavevmode\nobreak\ dV.$ (S11)

Following similar steps as for the number of atoms above, we arrive to the
conclusion that $\Delta$ scales also linearly with the bulk modulus, while the
$\varepsilon$ metric is insensitive to the it. Similarly, also $\nu$ is
insensitive to the reference value of the bulk modulus of the material, since
it only depends on the relative change of the bulk modulus $B_{0}$ between the
two systems, that does not depend on the factor $C$.

We stress that this fact does not mean that $\varepsilon$ or $\nu$ will not
capture a difference between bulk moduli when comparing two computational
approaches $a$ and $b$. It means, instead, that two datasets with a similar
discrepancy in the bulk moduli (say 2%) will result in the same $\varepsilon$
or $\nu$ irrespective of overall stiffness of their chemical bonds (i.e.,
their bulk modulus). We highlight that this shortcoming of the $\Delta$ metric
was already recognized in the literature and addressed by defining a modified
metric $\Delta_{1}$ in Ref. 63.

##  S7 Sensitivity of $\Delta$, $\varepsilon$ and $\nu$ to perturbations of
the EOS parameters and choice of thresholds for excellent and good agreement

Any of the three metrics $\Delta$, $\varepsilon$ or $\nu$ expresses the
difference between two EOS curves by a single number. It is not _a priori_
obvious, however, when those numbers can be considered small or large, and
which features of the EOS have the largest impact on the value. In this
section, we address these points.

In SI Fig. S7.1 we compare the EOS of a hypothetical material with the EOS
obtained after four different perturbations of the material (see caption for
details). The values for $\Delta$, $\varepsilon$ and $\nu$ that express the
difference between the original and perturbed EOS are listed for each case.
Analyzing these results allows to associate typical orders of magnitude to
each of the metrics, and quantify their variation with respect to changes in
$V_{0}$, $B_{0}$, and $B_{1}$.

Figure S7.1: EOS for a hypothetical material with $V_{0}=50.61\
\text{\AA}^{3}$ per formula unit, $B_{0}=0.71\ \text{eV}/\text{\AA}^{3}$, and
$B_{1}=4.67$ (“reference” curve, dashed blue), compared with the EOS where a
perturbation has been applied to some of the EOS-defining parameters
(“perturbed” curve, solid red). Each panel reports the resulting value for the
three metrics $\varepsilon$, $\nu$ and $\Delta$ obtained comparing the two
curves. The perturbed parameters and the magnitude of the perturbation are
also indicated in each panel. In particular, in panel (a) a perturbation of
0.12% is applied to $V_{0}$, of 0.7% to $B_{0}$ and of 4% to $B_{1}$. Panels
(b), (c) and (d) present the cases when a perturbation is applied
independently to $V_{0}$ (0.4%), $B_{0}$ (8%) or $B_{1}$ (160%), respectively,
while the other parameters are kept unchanged.

The hypothetical reference material that is used in SI Fig. S7.1 has
$V_{0}=50.61$ Å3 per formula unit, $B_{0}=0.71$ eV/Å3, and $B_{1}=4.67$. These
values are obtained as averages over the entire crystals set and thus
represent a hypothetical “average” EOS.

In panel (a), a perturbation is applied to all three parameters, namely 0.12%
to $V_{0}$, of 0.7% to $B_{0}$ and of 4% to $B_{1}$. These values are twice
the standard deviations of the discrepancies between the two AE codes in our
reference dataset, see Fig. 1 in the main text. Since the two EOS curves are
almost undistinguishable, this result highlights the high level of agreement
between our two AE codes. Based on this observation, we define a qualitative
range of $\varepsilon\lesssim 0.06$ or $\nu\lesssim 0.1$ for which we consider
two codes display excellent agreement. The threshold of $\varepsilon=0.06$
(approximately) corresponds to a determination coefficient $R^{2}\approx
1-\varepsilon^{2}=0.9964$ when one EOS curve is considered as a fit to the
other.

In panels (b), (c) and (d), instead, a larger perturbation is applied to only
one of the three parameters $V_{0}$, $B_{0}$ or $B_{1}$, respectively. The
perturbation to $V_{0}$ is chosen as $0.4\%$, and the magnitude of the
perturbations for $B_{0}$ and $B_{1}$ is scaled following the inverse ratios
1/20/400 of the weights for the $\nu$ metric (see SI Sec. S3). This results in
visually similar discrepancies between the two curves. This is an expected
result, and is another way to interpret the results discussed in SI Sec. S3
for the weights of $\nu$. Indeed, those weights were obtained by inferring the
error propagated on the fitted parameters from a given amount of random noise
on the datasets; the inverse weights can be conversely interpreted,
intuitively, as the relative magnitude of the perturbation to each of the
parameters required to induce similar changes to the EOS curve. These panels
help us make a number of observations:

  * •

All metrics ($\Delta$ and $\epsilon$ intrinsically, and $\nu$ by explicit
definition of the weights) give a much stronger importance to changes of the
equilibrium volume $V_{0}$ than to changes of the other parameters. This is a
positive feature of the metrics, as the EOS shape is mostly sensitive to
$V_{0}$ as well. The metric $\nu$ has the additional advantage, as already
discussed, that weights can be tuned to give more importance to other
parameters, if an application requires it.

  * •

With our definitions of $\epsilon$ and $\nu$, the two metric often return
similar values for a given pair of EOS curves, with $\varepsilon$ typically
slightly smaller (in the SI Sec. S8, we actually identify an approximate
proportionality ratio between the two, valid for small values of the metrics).

  * •

The values of the metrics on these last three panels allow us to define
threshold values for noticeable, but still relatively small, changes between
two EOS curves: $\varepsilon\lesssim 0.2$ or $\nu\lesssim 0.33$. Therefore, we
take the ranges of $0.06<\varepsilon\lesssim 0.20$ or $0.10<\nu\lesssim 0.33$
as the signature for good (but not excellent) agreement between two codes. The
upper end of this range $\varepsilon=0.20$ (approximately) corresponds to a
determination coefficient $R^{2}\approx 1-\varepsilon^{2}=0.96$ when one EOS
curve is considered as a fit to the other.

  * •

We can assign an intuitive meaning to the metric $\nu$. If two EOS only differ
in the equilibrium volume, its numerical value corresponds to the percentage
error on the equilibrium volume between the two curves. If also $B_{0}$ and
$B_{1}$ change, then $\nu$ will also take into account the discrepancies on
these two parameters, rescaled so that similar contributions to $\nu$ result
to similar quantitative changes to the EOS curve (in the volume range of
interest, $\pm 6\%$ in this work).

  * •

We note that in panels (b), (c) and (d) the value of $\nu$ is not exactly 0.4
as one might naively expect, because the perturbation that we apply refers to
the reference curve, but the $\nu$ metric is defined in a symmetric way, with
percentage differences with respect to the average of the two curves.

  * •

SI Fig. S7.2, finally, illustrates the clear disagreement between EOS curves
when $\varepsilon\geq 1.0$ or $\nu\geq 1.65$ (these values are used as upper
limit for the colorbars of the figures in SI Sec. S9). As a note,
$\varepsilon=1$ is an estimator for the situation where the coefficient of
determination $R^{2}\approx 1-\varepsilon^{2}$ starts to be negative (even if
the approximation $R^{2}\approx 1-\varepsilon^{2}$ does not hold exactly
anymore for such large values of $\varepsilon$). We highlight that a negative
$R^{2}$ value indicates that a horizontal line at the average value of the
data provides a better fit than the actual fit function. It is clear from SI
Fig. S7.2 that in such cases there is no agreement at all between the results
of two codes. Therefore, when $\varepsilon>1.0$ or $\nu>1.65$, two codes are
said to be clearly different in SI Sec. S9.

Figure S7.2: EOS for a hypothetical material with $V_{0}=50.61\
\text{\AA}^{3}$ per formula unit, $B_{0}=0.71\ \text{eV}/\text{\AA}^{3}$, and
$B_{1}=4.67$ (“reference” curve, dashed blue), compared with the EOS where a
perturbation has been applied to some of the EOS-defining parameters
(“perturbed” curve, solid red). With respect to SI Fig. S7.1, we highlight
here typical values for the metrics for large changes of the parameters that
make the EOS curves clearly different. Each panel reports the resulting value
for the three metrics $\varepsilon$, $\nu$ and $\Delta$ obtained comparing the
two curves. The perturbed parameters and the magnitude of the perturbation are
also indicated in each panel. In particular, in panel (a) a perturbation is
applied to $V_{0}$ to obtain a value of $\varepsilon=1$; panels (b) and (c)
apply a (positive and negative, respectively) perturbation to $B_{0}$
resulting in $\varepsilon=1$, and finally in panel (d) an even larger
perturbation is applied to $V_{0}$ to obtain a value of $\varepsilon=1.5$.

##  S8 Mutual correlation between the metrics

Figure S8.1: Cross-correlation between the $\Delta$, $\varepsilon$, and $\nu$
metrics for the entire data (unaries and oxides) presented in the main text.
Black lines are helpers to indicate a linear relation between the quantities
(slope of 1 in a log-log plot). Dashed lines on $\varepsilon$ and $\nu$ axes
show the “excellent” (green) and “good” (red) agreement thresholds recommended
in the main text. From the plots, these correspond to values of $\Delta$ of
approximately 0.3 and 0.95 meV/atom, respectively.

In order to assess the correlation between the various metrics, we present in
SI Fig. S8.1 for all codes mentioned in the main text their $\Delta$,
$\varepsilon$ or $\nu$ with respect to the reference EOS, for each of the
unaries and oxides. These $\Delta$, $\varepsilon$ or $\nu$ are plotted against
one of the other metrics. The results support the presence of an approximately
linear correlation between the three metrics ($\Delta$, $\varepsilon$, and
$\nu$). However, there is more scattering in the correlation between
$\varepsilon$ and $\Delta$ (and, similarly, between $\nu$ and $\Delta$), while
$\varepsilon$ and $\nu$ agree more consistently on a global scale.

We first note that the scattering between $\varepsilon$ and $\Delta$ can be
reduced by normalizing $\Delta$ by the bulk modulus, similar to the
$\Delta_{1}$ metric introduced in Ref. 63, as shown in SI Fig. S8.1(c). We
then observe that the new metrics $\varepsilon$ and $\nu$ are instead almost
linearly correlated when their values are $\lesssim 1$; for larger
discrepancies, the values of $\nu$ tend to grow faster than the values of
$\varepsilon$, i.e., $\varepsilon$ becomes relatively less sensitive to
further small changes to the EOS curves if they are already significantly
different.

This almost linear correlation can be justified with some approximations. Let
us consider the simple case of two parabolic EOS curves with the same $B_{0}$,
differing only in the equilibrium volume $V_{0}$. This is a valid
approximation, since we discussed above that both metrics are mostly sensitive
to changes of $V_{0}$ rather than $B_{0}$ or $B_{1}$. If we call
$2V_{R}=V_{M}-V_{m}$ the volume range for the integration in $\varepsilon$
(see also Eq. (4) in the main text for the definition of $V_{m}$ and $V_{M}$)
and $\tilde{V}$ the average volume of the two curves (with actual minima for
$V_{0}=\tilde{V}\pm\Delta V$), we obtain $\nu=200\frac{\Delta V}{\tilde{V}}$.
We now consider the limit in which $\Delta V\ll V_{R}$ (i.e., of a small
discrepancy of the two curves in the volume range of interest, corresponding
to the regime of small $\varepsilon$ and $\nu$ in which our data show an
almost linear relation between the two metrics).

It is then straightforward to show that $\varepsilon\approx
2\sqrt{15}\frac{\Delta V}{V_{R}}$. Indeed, writing the two curves as
$E_{1,2}(V)=A(V-\tilde{V}\pm\Delta V)^{2}$, with $A$ an appropriate
coefficient (the same for both curves with our assumptions of same bulk
modulus), we get that the integral in the numerator of $\varepsilon$,
$\langle[E_{1}(V)-E_{2}(V)]^{2}\rangle$, is given by:

$\displaystyle\frac{1}{2V_{R}}\int_{\tilde{V}-V_{R}}^{\tilde{V}+V_{R}}[A(V-\tilde{V}-\Delta
V)^{2}-A(V-\tilde{V}+\Delta V)^{2}]^{2}dV$
$\displaystyle=\frac{1}{V_{R}}\int_{\tilde{V}}^{\tilde{V}+V_{R}}[A(V-\tilde{V}-\Delta
V)^{2}-A(V-\tilde{V}+\Delta V)^{2}]^{2}dV$ (S12)
$\displaystyle=\frac{1}{V_{R}}\int_{\tilde{V}}^{\tilde{V}+V_{R}}[4A(V-\tilde{V})\Delta
V]^{2}dV=\frac{16A^{2}\Delta V^{2}V_{R}^{2}}{3}.$ (S13)

Similarly, we can obtain (using our assumption $\Delta V\ll V_{R}$) that

$\displaystyle\langle
E_{1}(V)\rangle=\frac{1}{2V_{R}}\int_{\tilde{V}-V_{R}}^{\tilde{V}+V_{R}}A(V-\tilde{V}-\Delta
V)^{2}dV\approx\frac{1}{V_{R}}\int_{\tilde{V}}^{\tilde{V}+V_{R}}A(V-\tilde{V})^{2}dV=\frac{AV_{R}^{2}}{3},$
(S14)

and $\langle E_{2}(V)\rangle=\langle E_{1}(V)\rangle$. Using similar steps,
one can also obtain $\langle[E_{1}(V)-\langle
E_{1}\rangle]^{2}\rangle=\langle[E_{2}(V)-\langle
E_{2}\rangle]^{2}\rangle\approx\frac{4}{45}A^{2}V_{R}^{4}$. Putting all
results together, we obtain the final result $\varepsilon\approx
2\sqrt{15}\frac{\Delta V}{V_{R}}$.

Finally, considering our choice of a $\pm 6\%$ volume range ($V_{R}\approx
0.06\tilde{V}$) gives $\epsilon\approx\frac{100}{3}\sqrt{15}\frac{\Delta
V}{\tilde{V}}$. Therefore, the two metrics are linearly dependent with the
ratio of

$\left.\frac{\nu}{\varepsilon}\right|_{\Delta V}=\frac{6}{\sqrt{15}}\approx
1.55,$ (S15)

which also matches well with the ratio of $\nu$ and $\epsilon$ calculated for
a specific value of $V_{0}$ perturbation shown in SI Fig. S7.1(b).

A similar analysis can be performed for two parabolic EOS curves with the same
$V_{0}$, but differing in $B_{0}$. In this case, $\nu=10\frac{\Delta
B_{0}}{\tilde{B}_{0}}$, where $\tilde{B}_{0}$ is the average $B_{0}$ of the
two curves and $\pm\Delta B_{0}$ is the difference from the average value for
the two curves. It can be shown, similarly to the case of differing $V_{0}$,
that in this case $\varepsilon\approx 3\frac{\Delta B_{0}}{\tilde{B}_{0}}$.
The resulting ratio of the two metrics now becomes

$\left.\frac{\nu}{\varepsilon}\right|_{\Delta B_{0}}\approx\frac{10}{3}\approx
3.33,$ (S16)

which also matches well with the ratio of $\nu$ and $\epsilon$ in SI Fig.
S7.1(c). The difference by a factor of $2$ between (S15) and (S16) shows that,
compared to $\varepsilon$, $\nu$ is $2$ times more sensitive to $B_{0}$
variation relative to $V_{0}$ variation.

To assess the typical ratio of $\nu$ and $\varepsilon$ in our dataset, in SI
Fig. S8.2 we show a cross-correlation plot between $\nu$ and $\epsilon$ for
the entire dataset of calculated crystals. One can see that for smaller values
of the metrics ($\nu<1$), the relation is approximately linear, and a
numerical fit gives the slope of $\frac{\nu}{\varepsilon}\approx 1.65$, which
is close to 1.55 found in Eq. (S15). This is expected, as both metrics give
larger weight to $V_{0}$ errors compared to $B_{0}$ or $B_{1}$ errors, and in
our dataset errors on $V_{0}$ (once rescaled with these weights) dominate over
$B_{0}$ and $B_{1}$ errors. This is also visible, for instance, from the
histograms of Fig. 1 in the main text, where the ratios of the standard
deviation of the histograms on $V_{0}$, $B_{0}$ and $B_{1}$ do not follow the
1, $\frac{1}{20}$, $\frac{1}{400}$ ratio of $\nu$.

Figure S8.2: Cross-correlation between the $\varepsilon$ and $\nu$ metrics on
a linear-axes plot. Left panel: all data; right panel: zoom on the data with
$\nu<1$ where the relation is approximately linear, together with a linear fit
of the data (red) resulting in a slope of $\approx 1.65$, and a line with a
slope $\approx 1.55$ (yellow curve) taken from the parabolic model.

The cross correlation plots allow to establish a data-driven relation between
different metrics. For instance, the two boundaries $\varepsilon=0.06$ and 0.2
selected in this project as “excellent” and “good” agreement between two EOS
s, and the corresponding thresholds for $\nu$ (0.10 and 0.33) have been chosen
according to these cross correlations. In addition, we can see from SI Fig.
S8.1a that these thresholds translate, for $\Delta$, approximately
$\Delta\approx 0.3$ and 0.95 meV/atom. This result is comparable with the
average $\Delta$ across AE codes $\langle\Delta\rangle=0.5-0.9$ meV/atom
obtained in the earlier benchmark[29, 30] for monoelemental solids, and is
consistent with the conclusion obtained there that $\Delta=1$ meV/atom is a
threshold under which one can speak about good agreement (for materials with
bulk moduli that are not particularly small). This consistency between former
and present benchmarks is about the metric; we refer to SI Sec. S15 for a
demonstration of consistency regarding the crystal set, and an illustration of
the added value of the present benchmark study.

##  S9 Detailed results for all computational approaches

In this section, we report the comparison of each of the computational
approaches considered in the main text with the average all-electron reference
dataset, using both metrics $\varepsilon$ and $\nu$. For each metric, the same
colorbar is used for all approaches, based upon the ranges of agreement
identified in SI Sec. S7 (in addition, the ratio of the threshold values for
$\varepsilon$ and $\nu$ is in agreement with their approximate linear
relationship, see SI Sec. S8):

  * •

“excellent agreement” ($\varepsilon\leq 0.06$, $\nu\leq 0.10$): a very dark
shade of blue (not evolving very much over this narrow interval);

  * •

“good agreement” ($0.06<\varepsilon\leq 0.20$, $0.10<\nu\leq 0.33$): color
evolving from a dark shade of blue to yellow as the values of $\varepsilon$ or
$\nu$ increase;

  * •

threshold for good agreement ($\varepsilon=0.20$, $\nu=0.33$): yellow;

  * •

“noticeably different” ($0.20<\varepsilon\leq 1.0$, $0.33<\nu\leq 1.65$):
color evolving from yellow to red as the values of $\varepsilon$ or $\nu$
increase;

  * •

“clearly different” ($\varepsilon>1.0$, $\nu>1.65$): one uniform darker shade
of red, regardless of the value.

Crystals that were not computed are left in white. The caption of every plot
mentions the number of crystals belonging to each of these categories. The
results for all the codes are shown in SI Figs. LABEL:fig:si-ptables-ABINIT to
LABEL:fig:si-ptables-WIEN2k.

<EMAIL_ADDRESS>ABINIT, 720, 232, 377, 111, 0, 244, 378, 98, 0

Figure S9.1: Value of the comparison metrics $\varepsilon$ (top) and $\nu$
(bottom) for [1] with respect to the average all-electron reference dataset.
Left panels: unaries; right panels: oxides. [3] out of 960 crystals were
calculated. The number of crystals that land in the excellent, good,
noticeably different, and clearly different agreement ranges for the
$\varepsilon$ metric are [4], [5], [6], [7], respectively. For the $\nu$
metric, they are [8], [9], [10], [11], respectively.

, BigDFT@DW|HGH-K(Valence), BigDFT, 402, 45, 97, 173, 87, 29, 106, 173, 94

Figure S9.2: Value of the comparison metrics $\varepsilon$ (top) and $\nu$
(bottom) for [1] with respect to the average all-electron reference dataset.
Left panels: unaries; right panels: oxides. [3] out of 960 crystals were
calculated. The number of crystals that land in the excellent, good,
noticeably different, and clearly different agreement ranges for the
$\varepsilon$ metric are [4], [5], [6], [7], respectively. For the $\nu$
metric, they are [8], [9], [10], [11], respectively.

, CP2K/Quickstep@TZV2P|GTH, CP2K-quickstep, 709, 57, 171, 317, 164, 55, 169,
302, 183

Figure S9.3: Value of the comparison metrics $\varepsilon$ (top) and $\nu$
(bottom) for [1] with respect to the average all-electron reference dataset.
Left panels: unaries; right panels: oxides. [3] out of 960 crystals were
calculated. The number of crystals that land in the excellent, good,
noticeably different, and clearly different agreement ranges for the
$\varepsilon$ metric are [4], [5], [6], [7], respectively. For the $\nu$
metric, they are [8], [9], [10], [11], respectively.

, FLEUR@LAPW+LO, FLEUR, 960, 936, 23, 1, 0, 938, 22, 0, 0

Figure S9.4: Value of the comparison metrics $\varepsilon$ (top) and $\nu$
(bottom) for [1] with respect to the average all-electron reference dataset.
Left panels: unaries; right panels: oxides. [3] out of 960 crystals were
calculated. The number of crystals that land in the excellent, good,
noticeably different, and clearly different agreement ranges for the
$\varepsilon$ metric are [4], [5], [6], [7], respectively. For the $\nu$
metric, they are [8], [9], [10], [11], respectively.

<EMAIL_ADDRESS>GPAW, 670, 130, 156, 350, 34, 128, 155, 347, 40

Figure S9.5: Value of the comparison metrics $\varepsilon$ (top) and $\nu$
(bottom) for [1] with respect to the average all-electron reference dataset.
Left panels: unaries; right panels: oxides. [3] out of 960 crystals were
calculated. The number of crystals that land in the excellent, good,
noticeably different, and clearly different agreement ranges for the
$\varepsilon$ metric are [4], [5], [6], [7], respectively. For the $\nu$
metric, they are [8], [9], [10], [11], respectively.

, CASTEP@PW|C19MK2, CASTEP, 960, 197, 410, 277, 76, 206, 399, 267, 88

Figure S9.6: Value of the comparison metrics $\varepsilon$ (top) and $\nu$
(bottom) for [1] with respect to the average all-electron reference dataset.
Left panels: unaries; right panels: oxides. [3] out of 960 crystals were
calculated. The number of crystals that land in the excellent, good,
noticeably different, and clearly different agreement ranges for the
$\varepsilon$ metric are [4], [5], [6], [7], respectively. For the $\nu$
metric, they are [8], [9], [10], [11], respectively.

, Quantum<EMAIL_ADDRESS>QE, 960, 388, 300, 199, 73, 395, 300,
184, 81

Figure S9.7: Value of the comparison metrics $\varepsilon$ (top) and $\nu$
(bottom) for [1] with respect to the average all-electron reference dataset.
Left panels: unaries; right panels: oxides. [3] out of 960 crystals were
calculated. The number of crystals that land in the excellent, good,
noticeably different, and clearly different agreement ranges for the
$\varepsilon$ metric are [4], [5], [6], [7], respectively. For the $\nu$
metric, they are [8], [9], [10], [11], respectively.

<EMAIL_ADDRESS>SIESTA, 698, 30, 117, 444, 107, 18,
137, 424, 119

Figure S9.8: Value of the comparison metrics $\varepsilon$ (top) and $\nu$
(bottom) for [1] with respect to the average all-electron reference dataset.
Left panels: unaries; right panels: oxides. [3] out of 960 crystals were
calculated. The number of crystals that land in the excellent, good,
noticeably different, and clearly different agreement ranges for the
$\varepsilon$ metric are [4], [5], [6], [7], respectively. For the $\nu$
metric, they are [8], [9], [10], [11], respectively.

<EMAIL_ADDRESS>SIRIUS-CP2K, 700, 363, 251, 81, 5, 374, 247,
72, 7

Figure S9.9: Value of the comparison metrics $\varepsilon$ (top) and $\nu$
(bottom) for [1] with respect to the average all-electron reference dataset.
Left panels: unaries; right panels: oxides. [3] out of 960 crystals were
calculated. The number of crystals that land in the excellent, good,
noticeably different, and clearly different agreement ranges for the
$\varepsilon$ metric are [4], [5], [6], [7], respectively. For the $\nu$
metric, they are [8], [9], [10], [11], respectively.

, VASP@PW|GW-PAW54*, VASP, 960, 403, 348, 200, 9, 419, 341, 189, 11

Figure S9.10: Value of the comparison metrics $\varepsilon$ (top) and $\nu$
(bottom) for [1] with respect to the average all-electron reference dataset.
Left panels: unaries; right panels: oxides. [3] out of 960 crystals were
calculated. The number of crystals that land in the excellent, good,
noticeably different, and clearly different agreement ranges for the
$\varepsilon$ metric are [4], [5], [6], [7], respectively. For the $\nu$
metric, they are [8], [9], [10], [11], respectively.

, WIEN2k@(L)APW+lo+LO, WIEN2k, 960, 936, 23, 1, 0, 938, 22, 0, 0

Figure S9.11: Value of the comparison metrics $\varepsilon$ (top) and $\nu$
(bottom) for [1] with respect to the average all-electron reference dataset.
Left panels: unaries; right panels: oxides. [3] out of 960 crystals were
calculated. The number of crystals that land in the excellent, good,
noticeably different, and clearly different agreement ranges for the
$\varepsilon$ metric are [4], [5], [6], [7], respectively. For the $\nu$
metric, they are [8], [9], [10], [11], respectively.

##  S10 Smearing and k-point convergence

For the study presented in this paper, a fixed choice of k-point integration
mesh and smearing has been implemented. In particular, the k-point mesh is a
uniform regular grid including the $\Gamma$ point, that guarantees a linear
spacing of 0.06 Å-1 in each of the three reciprocal-space directions, and the
smearing is a Fermi–Dirac type with broadening of 0.0045 Ry. This choice of
parameters is essential in order to compare with the reference dataset
presented in this manuscript, as explained in Box 3 of our recommendations. In
this section, we present the reasoning for our choice of parameters.

The Fermi–Dirac smearing is widely used in the community and it is implemented
in all the codes that participate in the study. The choice of the broadening
has been made according to the recommendation of Ref. 66: the smearing
parameter should not be too small to avoid sampling errors, nor too large to
prevent systematic deviations due to the dependency of the total energy on the
smearing broadening. The latter problem is explained in details in Ref. 66,
that shows the quadratic dependence of the free energy with respect to the
Fermi–Dirac smearing temperature. The problem of sampling error is instead
demonstrated in SI Fig. S10.1, where we analyze a FCC aluminum crystal
(conventional cell with 4 atoms) with an atom displaced by 0.1 Å with respect
to its equilibrium position. The figure reports the magnitude of the force on
the displaced atom as a function of the k-point integration mesh and smearing
broadening. The instability of the force for very small broadening is clearly
visible. Approaching the zero limit of the smearing, it becomes more and more
necessary to have a dense k-point integration mesh in order to maintain 0.001
eV/Å convergence on the forces. The sampling error is not a peculiar feature
of the Fermi–Dirac smearing; any other smearing type suffers from this
drawback, as demonstrated in Ref. 66. Our choice of smearing broadening
(0.0045 Ry $=\approx$ 61 meV) lays on the extreme right of SI Fig. S10.1 and
in this region a k-point integration mesh of $36\times 36\times 36$ is
sufficient to obtain a converged value of the force within 0.001 eV/Å. At the
same time, 0.0045 Ry is a small enough value to reduce to a minimum issues due
to the dependency of the total energy on the smearing broadening. In any case,
for the goal of verification, the exact value is not so important, as long as
all codes perform the very same choice of smearing and k-point integration
mesh.

Figure S10.1: Total force on a displaced atom as a function of the smearing
temperature and the k-point integration mesh. The system under investigation
is an FCC Al with an atom displaced by 0.1 Å in the $x$ direction. The volume
per atom is 16.47 Å3, for which a distance of 0.06 Å-1 between k-points
correspond to a $26\times 26\times 26$ mesh. The smearing is a Fermi–Dirac
smearing. Calculations of the forces are made with Siesta and DZP (double-zeta
polarized) basis set.

With our choice of the smearing broadening, we expect a rapid convergence of
the electronic-structure properties with respect to the the k-point density.
We test this assumption looking at the effect of the k-point integration mesh
on the estimation of the EOS parameters. Using WIEN2k results, we compare
calculations with k-point distance of 0.06 Å-1 and 0.045 Å-1 for all materials
in the study. The comparison is reported in SI Fig. S10.2. This figure shows
an overall discrepancy that is significantly smaller with respect to the
FLEUR-WIEN2k comparison presented in Fig. 1 in the main text (note that the
x-axis range is half of the one in Fig. 1). Looking at the histograms, we can
estimate that the overall agreement is at least a factor of 2 better for
$B_{1}$ with respect to the FLEUR-WIEN2k comparison, more than a factor of 4
better for $B_{0}$ and even one order of magnitude better for $V_{0}$. The
same conclusion cannot be drawn looking at the standard deviations reported in
the figure, due to the presence of two important outliers: RbO3 (3.7%
difference in $V_{0}$, -142.29% difference in $B_{0}$ and 149.35% difference
in $B_{1}$) and HeO (0.16% difference in $V_{0}$, 3.33% difference in $B_{0}$
and -7.90% difference in $B_{1}$). RbO3 and HeO are the only two materials
that are not converged with a k-point distance of 0.06 Å-1. All other
materials are converged within 0.07% of $V_{0}$. It is interesting to notice
that, even though RbO3 and HeO are not converged, their discrepancy is not so
dramatic when comparing FLEUR and WIEN2k (Fig. 1 in the main text). This, once
more, justifies our recommendation of adopting the same k-point integration
mesh for all computational approaches.

For completeness, we mention that the k-point mesh comparison of SI Fig. S10.2
has been performed on the crystal-structure set described in the main text for
the unaries. For the oxides set, instead, we present results of a previous
iteration of the volume refinement. Therefore the structures used for the
oxides calculation in SI Fig. S10.2 have central volumes that slightly differs
with respect to the ones used in the main text.

Figure S10.2: Histograms of the percentage difference between the results
obtained with a k-point integration mesh with linear density 0.06 Å-1 and
0.045 Å-1 for the three EOS parameters $V_{0}$, $B_{0}$, and $B_{1}$. Results
are obtained with WIEN2k code. Numbers close to the arrows indicate outliers
beyond the x-axis range.

##  S11 Band structure of erbium in the diamond crystal structure

In SI Fig. S11.1 we show the band structure of erbium in the diamond crystal
structure, obtained using the initial crystal structure (see SI Table
LABEL:table-starting-volume) with conventional cubic lattice parameter of
8.6296 Å. The simulation has been run with the Quantum ESPRESSO code (version
7.0) using the PBE pseudopotential from the SSSP PBE Precision 1.2 library.
78, and the recommended cutoffs of 40 Ry and 320 Ry for the wavefunctions and
charge density, respectively. A $9\times 9\times 9$ k-point integration mesh
is chosen (note: this is less dense than the mesh recommended in the main
text, but sufficient to demonstrate qualitatively the key features of the band
structure, that is the goal of this section), and a Fermi–Dirac smearing with
0.0045 Ry of broadening. Input and output files are available in Ref. 62.

The band structure clearly displays a set of “almost flat” bands (with a
dispersion of $\sim 0.1$ eV) very close to the Fermi energy in the range
between $-0.2$ eV and $0$ eV, originating from the $f$ states of erbium. For
this specific calculation and volume, these $f$ bands are just below the Fermi
level, but their position can shift with volume and cross the Fermi level,
significantly affecting the nature of the occupied states in the material.
This explains the unconventional shape of the EOS displayed in the main text
in Fig. 3.

An interesting note is that the location of the $f$ bands, determining the
lowest-energy minimum of the EOS, also depends on the value $Z$, i.e., on the
column of the periodic table. For our choice of parameters (k-point
integration mesh and smearing) erbium is at the boundary between elements
favoring the minimum at lower volume (for smaller $Z$) and elements favoring
the minimum at higher volume (for higher $Z$). This is clearly visible in SI
Fig. 3(a) as a jump of the first-neighbor distance for the diamond structures
for elements before and after erbium.

Figure S11.1: Band structure of erbium in the diamond crystal structure. The
zero of energy is set at the Fermi level.

##  S12 Total energy versus free energy

Here we want to show a couple of examples where the exact choice of the energy
to be computed in the EOS, namely the internal energy $E$, the free energy
$E-TS$ (where $-TS$ is the smearing contribution) or the approximation
$E-TS/2$, can significantly affect the EOS curves.

Figure S12.1: Comparison of EOS curves for two systems computed using the
FLEUR code, using different quantities on the energy axis. Left panel: Er in
the diamond structure; right panel: CsO2.

While in many cases the choice of one of these three quantities does not have
any effect on the curves (e.g., in the case of large-gap insulators), in SI
Fig. S12.1 we show two EOS curves that we selected as they show a significant
deviation (even larger deviations exist for other systems in our dataset). The
first case is erbium in the diamond structure. We already discussed its band
structure in SI Sec. S11: we expect that there are significant flat $f$ bands
(thus with high density of states) crossing the Fermi level as a function of
volume; therefore, the $-TS$ contribution will also depend on volume, since it
originates from the contribution of bands within a small energy range
(comparable to the chosen smearing broadening) from the Fermi level.

Figure S12.2: Band structure of CsO2 for two different volumes computed with
the Quantum ESPRESSO code. The shaded red area indicates an energy range of
$\pm 0.0045$ Ry, that is the smearing value recommended in this work.

One does not need, however, to consider such pathological cases: even other
systems might show important discrepancies between the various curves if the
density of states at the Fermi level changes significantly as a function of
volume. This is illustrated for instance by CsO2 (right panel in SI Fig.
S12.1), whose band structure for two different volumes is reported in SI Fig.
S12.2. Here, we see that bands with saddle points or almost-flat bands can be
found within $\pm\sigma$ from the Fermi level (with $\sigma$ being the value
of smearing broadening recommended here). These bands shift with volume and
therefore their contribution to $-TS$ will be a function of volume.

##  S13 Code specific parameters for the pseudopotential codes

In this section, we discuss the technical choices adopted for the verification
work and implemented, for each code, in the verification-PBE-v1 protocol. We
remind the reader that a number of parameters (such as the smearing type and
broadening or the k-point integration mesh) have been fixed for all codes (see
also Box 3 in the main text). We also remind that the parameters for the two
all-electron codes are discussed instead in SI Sec. S5.

###  S13.1<EMAIL_ADDRESS>
The Abinit calculations were performed with the version 9.6.2 of the code[42,
43, 13] and v0.2a2 of the aiida-abinit plugin. All calculations were run with
tolvrs = 1e-10, Fermi–Dirac smearing (occopt = 3) of 0.0045 Rydberg (tsmear =
0.00225), a reasonable number of empty bands (fband = 2), and a minimal
k-spacing of 0.06 Å-1. Norm-conserving pseudopotentials from the
PseudoDojo-v0.5 scalar-relativistic PBE standard library have been employed;
the “high” stringency recommended energy cutoffs were used. The RMM-DIIS
diagonalization algorithm[55](rmm_diis = 1) was used for calculations with
norm-conserving scalar-relativistic PBE pseudopotentials because of its
improved computational efficiency over the default CG method. Note that
rmm_diis = 1 means that the first four SCF iterations are performed with the
CG method in order to obtain reasonably good trial states before changing to
the RMM-DIIS method. RMM-DIIS is more efficient although less stable than the
CG algorithm as there is no explicit orthogonalization while optimizing the
trial states. This led to approximately a 0.7% error rate (108 out of 15460
total from all calculations performed in the process of this study) in the
calculations, however in all cases, running the failed calculations with CG
resulted in successful convergence.

###  S13.2 BigDFT@DW|HGH-K(Valence)

The version of the BigDFT code employed for these calculation is the 1.9.2.
For the great majority of the structures presented here, the pseudopotentials
employed in the calculations are norm-conserving
Hartwigsen–Goedecker–Hutter[70] of the Krack family[72] (HGH-K). For this
verification campaign, rather than choosing the most precise pseudopotentials
for a given element, we employed the default pseudopotentials with the least
possible number of valence electrons. We have therefore used a set which
provides an overestimation of the precision error of this pseudopotential
family. For comparison we have also included some of the semicore
pseudopotentials, see SI Sec. S16.2. The BigDFT code formalism employs
Daubechies wavelets basis sets to express the Kohn-Sham (KS) orbitals, which
enable to reach precise converged results for a given pseudopotential with
moderate effort with respect to the number of degrees of freedom employed. The
wavelets grid spacing was set to a value of $0.3\leavevmode\nobreak\ a_{0}$,
with all the high-resolution degrees of freedom activated, and the k-point
integration mesh correspond to a equivalent length of $94\leavevmode\nobreak\
a_{0}$. Density mixing scheme was employed for electronic convergence, reached
for a threshold voxel accuracy of $10^{-12}$ atomic units. Also, 120 empty
Kohn–Sham states were included for each k-point. Symmetry operations were also
included to limit the calculations to the irreducible k-points. At the time of
developing the workflow, the BigDFT code was migrating its user interface into
the PyBigDFT python module, which provides a user interface to the underlying
executable. Therefore, to avoid issues in the API modifications, we fixed the
PyBigDFT version to a beta release in the aiida-bigdft plugin, which activated
limited features with respect to the stable version available nowadays. In
particular, not all the structures which were defined with a non-orthogonal
unit-cell were transformed in an orthogonal supercell, required by the code.
This compatibility problem resulted in less structures treated by this
approach. It is planned to release a stable version of the plugin compatible
with the AiiDA 2.x API.

###  S13.3 CASTEP@PW|C19MK2

CASTEP is a plane-wave pseudopotential code[45], the 20.1.1 version is used in
this work. Calculation parameters closely follow the "precision" setting in
the initial common workflow implementation[36]. The cut off energy is fixed at
800 eV for all calculations since energy comparison is needed between
different chemical systems. The reciprocal space sampling is done through
$\Gamma$ centered Monkhorst–Pack grids with a fixed spacing of 0.06
$\text{\AA}^{-1}$ (i.e., 0.00954929 $2\pi\text{\AA}^{-1}$), in line with other
codes. On-the-fly generated (OTFG) core-corrected ultrasoft pseudopotentials
from the library C19 is used for the study except for the f-block elements
(see section S16.3 for more details). The C19 library is aimed for general use
with a balance between precision and speed, and it has been the default
potential library since CASTEP version 19.1.1. The modified pseudopotential
generation strings for f-block elements are tabulated in Table S16.1. The
energy convergence threshold for electronic minimization is set to $1\times
10^{-8}$ eV per atom.

###  S13.4 CP2K/Quickstep@TZV2P|GTH

The DFT module Quickstep of the open-source simulation package CP2K is an
implementation of the Gaussian and plane wave (GPW) and the all-electron
Gaussian augmented plane wave approaches[46, 47]. Therein, the Kohn-Sham
orbitals are represented by contracted Gaussian basis functions, whereas the
electronic charge density is expanded in plane waves[109]. For the former, an
accurate molecularly optimized triple-$\zeta$ basis set with two additional
sets of polarization functions (TZV2P-MOLOPT) is employed[73], whereas for the
latter a density cutoff of 2400 Ry is utilized, which differs from a
conventional plane wave cutoff by a factor of four. Due to its GPW method,
however, CP2K/Quickstep is rather insensitive with respect to high density
cutoffs. Furthermore, four multi-grids are used to ensure an efficient mapping
of product Gaussians onto the real-space integration grids, so that wide and
smooth Gaussian functions are mapped onto a coarser grid than narrow and sharp
Gaussians. To control which product Gaussians are mapped onto which level of
the multi-grid, a relative cutoff of 80 Ry is applied that defines the plane
wave cutoff of a reference grid covered by a Gaussian with unit standard
deviation. Separable and norm-conserving Goedecker–Teter–Hutter-type
pseudopotentials including scalar relativistic effects are used to describe
the interactions between the valence electrons and the ionic cores[71, 72].

###  S13.5<EMAIL_ADDRESS>
GPAW[48, 49] is an open-source DFT code developed at the Technical University
of Denmark (DTU) and other universities and computer centers, originally
created for combining a homogeneous grid basis set with the projector
augmented wave (PAW) method. Today, the code also provides a linear
combination of atomic orbitals (LCAO) basis[110] and a plane wave mode. The
latter has been applied in this study with a plane-wave cutoff of 800 eV for
all calculations, with the exception of systems containing noble gases, where
a cutoff of 1200 eV combined with a tighter energy and density convergence was
applied. As with other codes, Fermi–Dirac smearing of 0.06122 eV was used.
Besides the mentioned parameters, default values (as per GPAW v. 21) were used
for all other keywords needed to perform the calculation in order to ensure
that the results reflect the most representative user experience. The applied
PAW potentials, included in GPAW’s PAW potential suite, were specifically
created for the PBE exchange correlation functional, by applying GPAW’s setup
creator[111]. In particular, we use the pseudopotentials included in the setup
release 0.9.20000 available at
https://wiki.fysik.dtu.dk/gpaw/setups/setups.html#atomic-paw-setups. GPAW is
tightly linked to the atomic simulation environment (ASE)[8, 9], which handles
the user interface and is developed independently.

###  S13.6 Quantum<EMAIL_ADDRESS>
All calculations have been run using version 7.0 of the Quantum ESPRESSO code
and version 3.5.1 of the AiiDA Quantum ESPRESSO plugin. For the results
presented in the main text, all pseudopotentials were selected from the SSSP
PBE Precision 1.3[76] library, with plane-wave cutoffs corresponding to the
largest recommended value from the elements in each structure. In accordance
with the verification-PBE-v1 protocol, the Brillouin zone sampling was
performed using $\Gamma$-centered meshes with a spacing of 0.06 Å and a
Fermi–Dirac smearing of 0.0045 Ry. All other inputs parameters were set via
the precision protocol as described in the SI of Ref. 37, most importantly the
energy convergence threshold was set to a very strict 0.1 $\cdot 10^{-9}$ Ry
per atom. For the comparison with Abinit and CASTEP discussed in S14, the
pseudopotentials were selected from the PseudoDojo SR PBE standard set,
version 0.4, in UPF format, after a small modification of the .upf files as
described in S14. Plane-wave cutoffs were obtained from the "high" stringency
hints provided by the PseudoDojo table, all other computational parameters
were unchanged.

###  S13.7<EMAIL_ADDRESS>
The calculations presented in this work have been carried out with Siesta
version Max-1.2.0 (https://gitlab.com/siesta-project/siesta/-/tags/MaX-1.2.0)
powered by the aiida-siesta plugin version 1.2.0. Pseudopotentials from the
PseudoDojo FR standard set, version 0.4, in PSML format, have been employed.
The real-space cutoff for the representation of charge densities and
potentials is fixed at 900 Ry. The recommended course of action regarding
basis sets in Siesta is to perform an optimization considering the key
features of the chemical environment of each system. In this project, we have
not carried out the optimization for all 960 systems. Instead, we have
attempted a partial, per element, optimization, considering only the unary
diamond crystals at their central volume. The orbitals thus optimized for each
element are then reused for all the other unary and oxide structures involving
that element. The optimization starts from a TZDP basis with the addition of
an extra $f$ orbital. For alkali metals and alkaline earth metals the addition
of a $d$ orbital shell is also necessary. The optimization is performed with
the Nelder–Mead algorithm (multidimensional optimization without derivatives),
having as variables the first-zeta radius of each orbital and the split norm
parameter that controls the ratio between the first and subsequent zetas. We
foresee using the information garnered in this verification study to develop
further heuristics and guide the development of fully automatic methods to
generate basis sets taking into account appropriate chemical environment
descriptors.

###  S13.8<EMAIL_ADDRESS>
SIRIUS is a domain-specific library, which implements pseudopotential plane
wave and full potential linearized augmented plane wave methods and is
designed for GPU acceleration[54]. As such it brings additional
functionalities to CP2K such as collinear and non-collinear magnetic systems
with or without spin–orbit coupling. It is written in C++14 with the MPI,
OpenMP, and CUDA/ROCm programming models. As shown previously, SIRIUS/CP2K
allows for energy conserving ab-initio molecular dynamics simulations with a
constant shift in the order of $\mu$Ha compared to Quantum ESPRESSO reference
calculations [47]. All SIRIUS/CP2K simulations were performed using the
pseudopotentials of SSSP PBE Precision 1.2[80] together with a plane-wave
cutoff of 55 a.u.-1 for the density and potential, as well as 10 a.u.-1 for
$|\mathbf{k}+\mathbf{G}|$, respectively.

###  S13.9 VASP@PW|GW-PAW54*

All VASP results in this work have been obtained with VASP version 6.3.0 and
AiiDA-VASP 2.1.0, using preferentially the GW VASP PBE potential set version
54.

The input parameters explicitly set for this work were (defaults in brackets):
$\mathrm{PREC}=\mathrm{Accurate}(\mathrm{Normal})$,
$\mathrm{EDIFF}=\mathrm{1E-7}(\mathrm{1E-4})$,
$\mathrm{ALGO}=\mathrm{Normal}(\mathrm{Normal})$, $\mathrm{NELM}=300(60)$, and
$\mathrm{LMAXMIX}=6(2)$. The last setting ensures that all electronic states
up to quantum numbers $l=6$ are included in the density mixer, and it is
necessary to change this in order to converge some $d$-electron systems and
most lanthanides. The plane wave cutoff was fixed at 1000 eV ($\sim 73.50$ Ry)
for all calculations. The PAW method was used[112].

We used the recommended GW potential sets, whenever possible, with exception
of oxygen, where we chose the hard O_h_GW potential instead of the recommended
O_GW since it improved results for the oxides dataset. This potential is
required if short bonds to oxygen atoms are encountered. The recommended plane
wave cutoff energy for this potential is 765 eV, however, for high precision
studies we recommend to increase the cutoff in VASP by 30% yielding the
employed 1000 eV. For the elements where no GW potential is supplied, we use
the recommended standard PBE PAW potential. (The recommendations were taken
from the VASP-wiki on April 22nd, 2022).

##  S14 Precision of plane-wave codes when using the same pseudopotential
library

In order to assess how much codes implementing the same computational
approaches agree among each other, we compare the $\varepsilon$ metric among a
subset of the plane-wave codes considered in this work, when using the same
pseudopotential library. We consider two different cases. In the first case,
we compare three different codes (Abinit, Quantum ESPRESSO, and CASTEP) using
the PseudoDojo v0.4 pseudopotential library. In the second case, we compare
Quantum ESPRESSO and SIRIUS/CP2K using the SSSP PBE precision v1.2
pseudopotential library.

For the first case (PseudoDojo v0.4 library), an initial comparison
highlighted some discrepancies among the codes for Cu, Zn, and Ne. After
investigation, we found that these discrepancies stem from differences in the
numerical treatment of the form factors, the local part of the
pseudopotential, the model core charge density or the beta projectors with
spherical Bessel functions (usage of spline or not, different integration
methods, etc.).

While we have not precisely identified the main root of the effect, which
would require additional investigation, it is clear that the treatment of the
long range part of the pseudopotential is rather sensitive to the
implementation details. For this reason, we manually truncated the radial mesh
reported in the pseudopotential files (in UPF format, normally truncated at 10
Bobr) after 6 Bohr for all elements, as it is done in the psp8 format used by
Abinit, in order to have the same radial mesh in all codes. We verified that,
except for the three cases mentioned above, this truncation had no visible
effects on the results when compared to all-electron results. In the figures
below, the tables for which the radial mesh has been truncated will be
referred to as “trim” as a suffix to the approach label reported in the
periodic-table titles.

Using this simple truncation to avoid numerical instabilities, we show in SI
Figs. S14.1, S14.2 and S14.3 the six possible pairwise comparisons between the
three codes. Remarkably, except for two noble gases, the agreement is
excellent (we have used the same color scale for $\varepsilon$ as in the main
manuscript). This agreement is on average even better than the agreement
between AE codes of this work (see Fig. 1 of the main text). We highlight
however that this is expected, since the plane-wave codes use the same basis
set and, in this SI Section, an almost identical set of numerical parameters,
while the two AE codes differ in the details of the basis set and numerical
parameters.

Figure S14.1:  Value of the comparison metric $\varepsilon$ between the CASTEP
code (with truncated radial mesh) and the Abinit code for the unaries and
oxides set using the same PseudoDojo v0.4 library.

Figure S14.2:  Value of the comparison metric $\varepsilon$ between the
Quantum ESPRESSO code (with truncated radial mesh) and the Abinit code for the
unaries and oxides set using the same PseudoDojo v0.4 library.

Figure S14.3:  Value of the comparison metric $\varepsilon$ between the
Quantum ESPRESSO and CASTEP codes (both with truncated radial mesh) for the
unaries and oxides set using the same PseudoDojo v0.4 library.

In the second case, we compared the SSSP PBE precision v1.2 pseudopotential
library for Quantum ESPRESSO and SIRIUS/CP2K. The comparison of the
$\varepsilon$ metric is shown in SI Fig. S14.4. The match is very good for all
systems except for Ba-diamond and RbO3. We note here that in the main text, we
used the SSSP PBE precision v1.3 for Quantum ESPRESSO, but we stress that it
is equivalent to v1.2 for the subset of chemical elements considered in this
section (the difference being that v1.3 also includes actinides).

Figure S14.4:  Value of the comparison metric $\varepsilon$ between the
Quantum ESPRESSO and SIRIUS/CP2K software for the unaries and oxides set using
the same SSSP PBE precision v1.2 library.

The overall excellent agreement displayed in this SI Section demonstrates that
different codes implementing the same computational approach (including, in
addition to the basis set, also the pseudopotential library and other
computational parameters) can reproduce the same results. In addition, our
investigation helps identifying the remaining numerical aspects that might
produce different outcomes and that merit further investigation, such as the
truncation of the radial mesh discussed here.

##  S15 Transferability of conclusions from the previous smaller crystal set
to the current 960 crystal set

In this section we want to show that even if two codes provide nearly
identical results for the 71 crystal set of Refs. 29, 30, they will not
necessarily provide identical results for the larger 960 crystal set of the
present paper, demonstrating the value of the larger crystal-structure set.

To perform a quantitative comparison using the same computational approach, we
use only data obtained as part of the present study. That means that we do not
use the 71 crystal set (that has several crystals that are not part of the
present set), but rather the set of 29 crystals common to both sets, as listed
in SI Sec. S17, as a proxy for the 71 crystal set. SI Fig. S15.1 shows on the
horizontal axis the $\epsilon$ averaged over these 29 crystals for any pair of
approaches or codes used in this work, and on the vertical axis the $\epsilon$
for the same pairs but now averaged over all crystals (up to 960 crystals).
The overall correlation in SI Fig. S15.1 shows that the 29 crystal set does
indeed capture part of the information. However, as expected, the data points
fall above the $y=x$ line, demonstrating that the large set probes relevant
behavior that is not probed in the small sets previously used. In particular,
there are code pairs (see shaded area in SI Fig. S15.1) that have a mutual
$\varepsilon\leq 0.2$ for the 29 unaries (good agreement), yet have a larger
$\varepsilon$ for the large set (and similarly for $\varepsilon\leq 0.06$,
that signals excellent agreement). SI Fig. S15.2 shows two examples when this
happens: panels (a) and (c) demonstrate that while values for $\varepsilon$
are small on the 29 crystals of the small set (boxed cases in SI Fig.
S15.2(a,c)), other crystals outside this set may contribute to a larger
average $\varepsilon$. This is typically the case for lanthanides, not
included in the earlier 71-crystal benchmark (that stopped at Rn). Moreover,
some elements lead to a low $\varepsilon$ for the unaries, but to a larger
$\varepsilon$ for the oxides (e.g., Cs-Ba-Fr-Ra and Te-I-Xe-Bi-Po-At-Rn for
CASTEP).

The overall conclusion of this analysis is that the 71-crystal set of unaries
– or the 29-crystal set as its proxy – gives a fair first assessment of the
comparative behavior of two DFT codes, while the complete 960-crystal set of
unaries and oxides provides a more detailed comparison, both because more
elements and more structures per element are included. Conclusions that were
based on the benchmarks of Refs. 29, 30 will therefore still hold, yet can be
refined. For testing newly developed pseudopotential libraries, a stepwise
approach can be implemented: a first quick test on one unary per element will
reveal the largest deviations; once these are fixed, the whole dataset of
unaries and oxides can be used to hunt for smaller deviations.

Figure S15.1: Vertical axis: $\varepsilon$ averaged over all crystals of the
unary and oxide set computed in this work, for all pairs of computational
approaches considered in the main text. Horizontal axis: $\varepsilon$
averaged over the 29 crystals from these 960 that appear as well in the 71
crystal set from Refs. 29, 30. The horizontal and vertical lines indicate the
threshold values for excellent (green) and good (red) agreement, as discussed
in the main text. The region above the horizontal red line and to the left of
the vertical red line, highlighted by the shaded blue area, indicates pairs
that are in good agreement according to the small crystal-structure set, but
less so when considering the full set.

Figure S15.2: Same periodic tables for the comparison metric $\varepsilon$ as
in SI S9 for codes CASTEP and Quantum ESPRESSO, where the 29 crystal
structures that appear also in the 71-crystal set of Ref. 29 are highlighted.
While the 29 structures display values with excellent agreement, other
structures result in a less good agreement. These two examples illustrate some
of the cases in SI Fig. S15.1 for which the agreement is very good on the
small crystal set, yet less good on the large crystal set.

##  S16 Additional pseudopotential datasets

This section discusses more in detail some additional datasets obtained using
the same codes of the main text, but different basis sets or pseudopotential
families, as well as the comparison with earlier versions of pseudopotentials,
before those that were optimized here.

###  S16.1 Abinit

During the process of computing the oxide verification equations of state
using the initial standard norm-conserving scalar-relativistic PBE PseudoDojo
(version 0.4) [68, 113], we observed that the results for around 11 of the
pseudopotentials were not in as good agreement with the all-electron equations
of state as other elements from the same pseudopotential family. This led to
an investigation into possible improvements to the pseudopotentials for Ba,
Bi, I, Pb, Po, Rb, Rn, S, Te, Tl, and Xe. With the exception of S, we found
that the accuracy of the pseudopotentials is significantly improved by
including a projector for the $f$ channel. In the original version, indeed,
the local part of the pseudopotential was not able to reproduce the all-
electron scattering properties of the $f$ channel in the empty region.

A particularly severe case is the one of Ba which is shown in SI Fig. S16.1(a)
where the pseudized $f$ channel (black dashed line) does not reproduce well
the all-electron reference. As a consequence, the EOS for the stable BCC phase
of Ba is slightly off (SI Fig. S16.1(c)) but the one for BaO3 is completely
wrong, as shown in SI Fig. S16.1(d). This issue is fixed by the additional
projector as shown in SI Fig. S16.1(b), which then gives excellent agreement
with the all-electron EOS.

Figure S16.1: Comparison of the scattering properties of the all-electron (AE)
atom and the pseudized (PS) Hamiltonian for two different Ba pseudopotentials
generated without (with) an $f$ projector. The subfigures (a) and (b) show the
arctangent of the $l$-dependent logarithmic derivative computed for some $R$
greater than the pseudization radius, where $\psi$ is the solution of the non-
local radial equation regular at the origin (a) without and (b) with $f$
projector. Since the transferability of a pseudopotential is directly related
to the capability of reproducing the AE logarithmic derivative over a wide
range of energies, the version with $f$ projector is expected to provide more
accurate results and we provide two equation of state example in (c) for BCC
barium and in (d) for BaO3 where the pseudopotential with the additional
projector is more accurate than the one without. This is true for all 10 cases
tested in this study.

It should be noted that the inclusion of the $f$ projector increases the
computational cost associated to the application of the non-local part of the
Hamiltonian $V_{nl}$. This is especially true if the computation of
$V_{nl}|\psi\rangle$ is obtained by projecting the wavefunction over spherical
Harmonics $Y_{lm}$. We stress, however, that in Abinit the projection is
implemented by expressing the sum over $m$ in terms of Legendre polynomials
(useylm = 0 input variable, default option when norm-conserving
pseudopotentials are employed). In this case, the computational cost of
including an additional projector for $l$ is not so high because, contrary to
the case when spherical harmonics are used, the number of floating point
operations required to apply a projector does not scale with the total number
of magnetic quantum numbers $2l+1$. This study led to the creation and
adoption of 11 new pseudopotentials for Ba, Bi, I, Pb, Po, Rb, Rn, S, Te, Tl,
and Xe in a new PseudoDojo (version 0.5). Because these pseudopotentials have
not been subject to a convergence study with respect to all-electron results,
20 Hartree were added to the “high” stringency recommended cutoffs from
ONCVPSP for safety.

In addition to norm-conserving pseudopotentials, we also investigated the JTH
PAW PBE v1.1 [63] table which has improved versions of the pseudopotentials
for H, Li, Si, Cu, Zn, Ga, Cd, Sb, Lu, Os, Ir and Bi with respect to v1.0 used
in the previous AiiDA common workflows study. Through testing, it was noted
that the recommended kinetic energy cutoff values for these PAW
pseudopotentials were not sufficient for the desired level of agreement with
all-electron codes, so twice the “high” recommended cutoffs were used. The
conjugate-gradient (CG) diagonalization algorithm was used for all
calculations with PAW pseudopotentials. Note that there do exist
pseudopotentials for lanthanide- and actinide-series elements in this PAW
table but they have not been verified against all-electron reference results.
Because of this and because the PAW potentials are included primarily as a
point of comparison against the norm-conserving PseudoDojo family (which does
not provide lanthanide nor actinide elements), these elements are not included
in the reported results and analysis.

We therefore compare in SI Fig. S16.2 the performance of the three PBE
pseudopotential libraries tested here using Abinit: the norm-conserving
PseudoDojo (version 0.4), the norm-conserving PseudoDojo (version 0.5) and the
PAW JTH v1.1 tables. We conclude that PseudoDojo (version 0.5) is the most
precise one and is therefore used in the main manuscript.

Figure S16.2: Comparison of the three tested PBE pseudopotential tables using
the Abinit software.

In summary, unary and oxide verification results were calculated for Abinit
with three sets of pseudopotentials: JTH PAW PBE v1.1, PseudoDojo norm-
conserving standard scalar-relativistic PBE v0.4, and a new PseudoDojo norm-
conserving standard scalar-relativistic PBE v0.5 based on v0.4 with improved
potentials for the 10 elements listed above. The latter set is reported in the
main text of the manuscript.

###  S16.2 BigDFT

Data production using BigDFT showed clear outliers in comparison to reference
codes (FLEUR/WIEN2k) using HGH-K Valence only pseudopotentials. We therefore
recalculated the EOS using semicore pseudopotentials for all the crystals for
which this type of pseudopotential was available and that showed a difference
in volume $\geq 0.2$ Å3 with respect to the reference AE data. The new
calculations resulted, in most cases, in a significant shift towards a closer
agreement with the reference AE, as shown in SI Fig. S16.3.

Figure S16.3: Comparative plot displaying the improvement obtained by using
HGH-K (Semicore) pseudopotentials over Valence only pseudopotentials.

###  S16.3 CASTEP

A unique feature of CASTEP is that the pseudopotentials are typically
generated on-the-fly during the calculations, although file-based potentials
are still supported. Each on-the-fly generated (OTFG) potential is defined
using a compact configuration string. An OTFG library is a collection of
predefined configurations. Many such libraries are built into the CASTEP
executable itself, with different focuses. In this study, the pseudopotentials
used are provided by the C19 library, except for the lanthanide and actinide
elements. For these elements, a new set of pseudopotentials are generated to
improve agreement with the all-electron reference data, which had not been
available during the development of the original C19 library.

The configuration strings for these elements are tabulated in Table S16.1.
Each field is separated by "|". The first field is the local angular momentum
channel. The second is the core radius in atomic units and the next three are
the recommended cut off energies (Ha) corresponding to the Coarse, Medium, and
Fine settings for CASTEP, which do not affect a calculation if a cut off
energy is specified explicitly (as in this study). This is followed by the
orbitals to be pseudized, which are separated by ":". Each orbital is defined
by a nl number and may have suffixes to indicate what kind of projectors
should be used. For example, "60U" indicates that a single ultrasoft projector
should be used for the 6s channel (the default is two ultrasoft projectors).
The "NN" suffix indicates that two norm-conserving projectors should be
included. The "U2U2" suffix indicates two ultrasoft projectors each with a
core radius of 2.0. The "L" suffix pins the local channel, and the "P" suffix
indicates that the pseudized channel is not represented by an explicit
projector or a local channel. The parameter "qc" inside the brackets controls
the smoothness of the potential. Occupations of atomic states can be further
modified by settings inside a curly bracket.

Element | Old settings | New settings
---|---|---
La | 2|2.3|5|6|7|50U:60:51:52(qc=4.5) | 2|2.3|5|6|7|50U:60:51:52:43{4f0.1}(qc=4.5)
Ce | 2|2.1|9|10|11|50U:60:51:43:52L(qc=6) | 2|2.2|8|9|10|50U:60:51:52:43{5d0.1}(qc=4.5)
Pr | 2|2.1|9|10|11|50U:60:51:43(qc=6) | 2|2.1|10|12|13|50U:60:51:52:43{5d0.1}(qc=5)
Nd | 2|2.1|9|10|11|50U:60:51:43(qc=6) | 2|2.1|10|12|13|50U:60:51:52:43{5d0.1}(qc=5)
Pm | 2|2.1|10|12|13|50U:60:51:43(qc=6) | 2|2.1|8|9|11|50U:60:51:52:43{5d0.1,4f4}(qc=5.5)
Sm | 2|2.1|10|12|13|50U:60:51:43(qc=6) | 2|2.1|9|10|12|50U:60:51:52:43{5d0.1,4f5}(qc=5.5)
Eu | 2|2.1|10|12|13|50U:60:51:43(qc=6) | 2|2.1|9|10|12|50U:60:51:52:43{5d0.1,4f6}(qc=5.5)
Gd | 2|2.1|10|12|13|50U:60:51:52L:43(qc=6) | 3|2.1|9|10|12|50U:60:51:52:43(qc=5.5)
Tb | 2|2.1|10|12|13|50U:60:51:43(qc=6) | 2|2.2|12|13|15|50U:60:51:52:43{5d0.1}(qc=5)
Dy | 2|1.9|12|14|16|50U:60:51:43(qc=6.5) | 2|2.0|12|13|15|50U:60:51:52:43{5d0.1}(qc=6.5)
Ho | 2|1.9|12|14|16|50U:60:51:43(qc=6.5) | 2|2.0|12|13|15|50U:60:51:52:43{5d0.1}(qc=6.5)
Er | 2|2.1|10|12|13|50U:60:51:43{6s0.5}(qc=6) | 2|2.1|10|12|13|50U:60:51:52:43{6s0.1,5d0.1}(qc=6)
Tm | 2|2.1|10|12|13|50U:60:51:43{4f12}(qc=6) | 2|2.1|10|12|13|50U:60:51:52:43{5d0.1,4f12}(qc=6)
Yb | 2|2.1|10|12|13|50U:60:51:434f13(qc=6) | 2|2.1|10|12|13|50U:60:51:52:43{5d0.1,4f13}(qc=6)
Ac | 2|2.5|5|6|7|60U:70NN:61:62L | 2|2.4|7|7|9|60U:70U2U2:61:62:53{6d0.1,5f0.1}(qc=5)
Th | 2|2.5|7|7|9|60U:70NN:61:62 | 2|2.2|7|7|9|60U:70U2U2:61:62:53{5f0.1}(qc=5)
Pa | 2|2.1|9|10|11|60U:70:61:53:62P(qc=6) | 2|2.2|8|9|10|60U:70U2U2:61:62:53(qc=5)
U | 2|2.1|10|12|13|60U:70:61:53:62P(qc=6) | 2|2.2|8|9|10|60U:70U2U2:61:62:53(qc=5)
Np | 2|2.1|10|12|13|60U:70:61:53:62P(qc=6) | 2|2.2|9|10|12|60U:70U2U2:61:62:53(qc=5)
Pu | 2|2.1|10|12|13|60U:70:61:53:62P{6d1,7s1}(qc=6) | 2|2.2|9|10|12|60U:70U2U2:61:62:53{6d0.1}(qc=5.5)
Am | 2|2.1|10|12|13|60U:70:61:53(qc=6) | 2|2.2|9|10|12|60U:70U2U2:61:62:53{6d0.1}(qc=5.5)
Cm | 2|2.1|10|12|13|60U:70:61:53:62L(qc=6) | 2|2.2|9|10|12|60U:70U2U2:61:62:53(qc=5.5)

Table S16.1:  Configuration strings of the on-the-fly pseudopotential
generation before and after the update for the lanthanide and actinide
elements involved in this study.

Figure S16.4: Comparing EOS curves of La and Ce in the BCC configuration. The
updated pseudopotentials agree better with the all-electron data compared with
those from the C19.

Examples EOS curves for La and Ce are displayed in SI Fig. S16.4, showing that
the updated pseudopotentials agree better with the all electron data. The
improvements may vary for other elements with $f$ electrons. One should note
that a better fit to the all electron results does not necessarily mean
smaller errors compared to the experimental results. This is due to the
inherent self-interaction errors in semi-local DFT that are present in the
description of $f$ electrons.

For some elements, the improvements are achieved through the inclusion of
additional orbital states and l channels/projectors such as La, where
previously the 4f channel was neglected. In other cases, the reference atomic
calculations include partially occupied atomic states that are otherwise empty
in the original configurations. For example {5d0.1} adds 0.1 electrons to the
5d channel. In some cases, the updated potentials contain increased core radii
and are made softer via a decreased "qc" value as far as possible. These
modifications are applied consistently for elements that neighbors to each
other in the periodic table. We also want to emphasize that the settings used
in this study are one-shot updates based on the C19 library, rather than the
outcomes of iterative optimizations. With the help of the automated test
framework and publicly available all-electron reference data, it should be
easy to adjust the strings and test them rigorously for further improvements,
as required.

###  S16.4 Quantum ESPRESSO

SSSP[75] is a library of pseudopotentials that undergoes rigorous verification
and ranking procedures. It contains two distinct families of pseudopotentials.
The first family (“efficiency”) is composed of relatively soft
pseudopotentials that are still sufficiently precise for use in high-
throughput calculations. The second family (“precision”) contains
pseudopotentials that are extremely precise with respect to all-electron
references, even if more computationally expensive. In our Quantum ESPRESSO
calculations, aiming at high precision, we have therefore used the SSSP
“precision” library. Version 1.1.2 of SSSP, available before this work, was
verified only on the unary configurations presented in Ref. 29, 30. With the
additional results of this work, we have identified that certain
pseudopotentials were not the best selection, and that more precise
pseudopotentials are available in libraries that were not included in the
previous generation of the SSSP library.

Element | SSSP precision v1.1.2 | SSSP precision v1.3
---|---|---
Te |  6|US |GBRV|v1|uspp |  6|US | PSL|v1.0.0-low|ld1
Na |  9|NC | PD|v4-std|oncvpsp3 |  9|PAW| PSL|v1.0.0-low|ld1
Cu | 19|NC | PD|v4-std|oncvpsp3 | 11|PAW| PSL|v1.0.0-low|ld1
Cs |  9|US |GBRV|v1|uspp |  9|NC | PD|v4-str|oncvpsp3
Cd | 12|US | PSL|v0.3.1|ld1 | 20|PAW| PSL|v1.0.0-high|ld1
Ba | 10|PAW| PSL|v1.0.0-high|ld1 | 10|NC | PD|v5|oncvpsp4
As |  5|US | PSL|v0.2|ld1 | 15|NC | PD|v4-std|oncvpsp3
I | 17|PAW| PSL|v1.0.0-high|ld1 | 17|NC | PD|v4-std|oncvpsp3
Hg | 20|NC |SG15|v0|oncvpsp3 | 12|US |GBRV|v1|uspp
Ne |  8|NC |SG15|v0|oncvpsp3 |  8|PAW| PSL|v1.0.0-high|ld1
Ar |  8|NC |SG15|v0|oncvpsp3 |  8|PAW| PSL|v1.0.0-high|ld1
Kr |  8|NC |SG15|v0|oncvpsp3 | 18|PAW| PSL|v1.0.0-high|ld1
Xe | 18|NC |SG15|v0|oncvpsp3 | 18|PAW| PSL|v1.0.0-high|ld1
Rn | 18|NC |SG15|v0|oncvpsp3 | 18|PAW| PSL|v1.0.0-high|ld1
Ir | 15|NC |GBRV|v1|uspp | 31|US | PSL|v1.0.0-high|ld1

Table S16.2:  List of pseudopotentials modified between SSSP v1.1.2 and SSSP v1.3. Each element in the table is composed of 5 strings separated by a | symbol, respectively: the number of electrons in the valence ($Z$), the type of pseudopotential (NC: norm-conserving, PAW: projector-augmented wave, US: ultrasoft), the source library (SG15: from Ref. 114, PSL: PSlibrary[82], GBRV: from Ref. 115, PD: PseudoDojo[68]), an internal version number that identifies the pseudopotential inside the given source library and the code used to generate them (oncvpsp3, oncvpsp4: version 3 and version 4 of the ONCVPSP code[113], ld1: the ld1.x code of Quantum ESPRESSO[51], uspp: the UltraSoft PseudoPotential (USPP) generation code[116]).

The pseudopotentials that were updated are summarized in Table S16.2. In
particular, the pseudopotentials of the noble gases (Ne, Ar, Kr, Xe, Rn) were
previously obtained using the SG15[114] library; in SSSP v1.3, we have
replaced them with “PAW-high” pseudopotentials from the PSLibrary[82] based on
our more recent verification results as they provide significantly better
agreement with the all-electron reference EOS curves. An example of the
improvement in the case of NeO2 is shown in SI Fig. S16.5.

The pseudopotentials of I, Hg were obtained from SG15[114] as well. In SSSP
v1.3, the pseudopotential of I is replaced with the one from PseudoDojo
library and the pseudopotential of Hg is replaced with the on from GBRV[115]
library. It is not only because these two pseudopotentials from SG15 library
are less precise, but also because these two pseudopotentials lead to the
electronic step convergence issue during the calculation. Using new
pseudopotentials makes all equation of state calculation finished without
issues.

The SSSP v1.1.2 libraries for Ir have been updated with the latest
pseudopotential from Pslibrary US v1.0.0 in the “high” family[82]. It should
be noted that the original Ir pseudopotential from the GBRV library contains a
ghost state at 10eV.

Figure S16.5: Left: EOS of NeO2 using the Ne pseudopotential from SSSP v1.1.2.
Right: EOS of NeO2 using the Ne pseudopotential from SSSP v1.3.

For Ba, in SSSP v1.3 we select the new pseudopotential from the PseudoDojo
v0.5 library[69] (generated in the context of this work, see SI Sec. S16.1)
that includes an $f$ projector.

For Te, Na, and Cu we consider in SSSP v1.3 the PAW or ultrasoft
pseudopotentials from the PSlibrary[82] in the “low” family. These
pseudopotentials have fewer semicore states and larger cut-off radii, making
them possibly less accurate, but optimized for lower kinetic energy cut-offs.
Tests on ten configurations, including oxides, revealed that some
pseudopotentials from the “low” family are actually even more accurate than
those from the “high” family (the latter have more semicore states and smaller
cut-off radii). We note that PSlibrary suggests to use “high” pseudopotentials
only for special applications, while “low” ones can yield sufficient precision
for regular calculations. As an example of the results with the new pseudos,
we show the results for TeO2 for the two versions of SSSP in SI Fig. S16.6.

Figure S16.6: Left: EOS of TeO2 using the Te pseudopotential from SSSP v1.1.2.
Right: EOS of TeO2 using the Te pseudopotential from SSSP v1.3.

Starting from v1.3, the SSSP library also includes pseudopotentials for
actinides (Th-Lr) developed in Ref. 81 as well as the pseudopotentials for Ac,
At, Ra, and Fr from the “high” family of PSlibrary[82].

Finally, in the case of Cs, As, while the test of the SSSP v1.1.2
pseudopotentials on unaries was resulting in good-quality EOS curves, we
obtained significant disagreements for oxides. We therefore replaced the
corresponding pseudopotentials with others that provided more precise
agreements. We illustrate the improvements in the case of Cs in SI Fig. S16.7.

Figure S16.7: Comparison (via the $\nu$ metric) between Quantum ESPRESSO
results and the all-electron reference from WIEN2k for the different crystal
structures of Cs, revealing that while the old pseudopotential was already
quite precise for unaries, only the new pseudopotential in SSSP v1.3 (from the
PseudoDojo library) provides high precision results for oxides.

###  S16.5 VASP

The lanthanides Pr, Nd, Pm, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, and Yb were
updated during the preparation of the data, since the initial ones showed
quite significant deviations from the all-electron results. For elements
containing $4f$ electrons, VASP recommends to use potentials with the $f$
electrons placed in the frozen core, to avoid the well known self-interaction
errors resulting from DFT. VASP therefore provides well-tested potentials with
frozen $f$ electrons for the lanthanides, with a valency of 2 or 3 (for Er,
Eu, and Yb, both valences are available as separate potentials). Since this
study has settled on treating $f$ electrons explicitly as valence, we used the
potentials that place all $f$ electrons in the valence. Using semi-local
functionals, these potentials lead to significant over-binding and too small
unit cells (compared to experiment as well as compared to PAW potentials that
place the $f$ electrons in the core). However, the aim of this work is code
comparison primarily. The lanthanide potentials used in this work have been
generated by G. Kresse. They are using a much smaller core radius of 2.2 a.u.
than the previous versions. Additionally, the reference electronic
configuration was altered, by placing 0.5 electrons (instead of 1 electron)
from the 4$f$ shell into the 5$d$ shell. Generally, two projectors were used
for the $f$ shells. Due to contracting $f$ shells towards the right of the
series (Tm, Er, and Yb) a third projector was found to be necessary to obtain
an accurate description of the $f$ scattering properties. The new lanthanide
potentials will be released on the VASP portal in the 6.4 PBE PAW potential
set with an _h suffix. For future reference the specific potential mapping
that was used is $\mathrm{RECOMMENDED\\_ACWF\\_LANTH}$, where $\mathrm{LANTH}$
indicates that the new lanthanide potentials have been used.

For the unaries we have compared the elevated settings optimized for this
study with the default settings of VASP (keeping only $\mathrm{LMAXMIX}=6$)
and the chosen recommended GW potential set with the recommended PBE potential
set. The resulting $\varepsilon$ metrics are plotted in SI Fig. S16.8.

(a)

(b)

(c)

Figure S16.8: Value of the comparison metric $\varepsilon$ for the unaries
using three different settings for VASP with the averaged AE results as
reference. The datasets in (a) and (b) use the recommended PBE potentials.
Computational parameters are mostly left at default values in (a), while (b)
increases the energy cutoff to 800 eV and adopts the other parameters
described in section S13.9. In (c) we show the final dataset with recommended
GW potentials, if available, and 1000 eV energy cutoff.

In the first panel of SI Fig. S16.8 we used the recommended PAW potentials for
PBE (as given on the VASP website, for the lanthanides the new PAW potentials
were used as detailed above), with the default cutoffs (dependent on the
system) and precision setting ($\mathrm{PREC}=\mathrm{Normal}$). In panel (b)
we used the same potentials, but $\mathrm{PREC}=\mathrm{Accurate}$ and a plane
wave cutoff of 800 eV alongside the other settings detailed in section S13.9,
while panel (c) corresponds to the final results presented in the paper using
the final settings (1000 eV) and the recommended GW potentials.

Clearly, changing the default parameters substantially improves the agreement
with the AE reference calculations. The noble gases in particular are
dramatically improved by using the more accurate settings. But many other
weakly bonded structures are improved as well. If then also the GW potentials
are used, some elements do improve again by quite a lot (Ba, Hg, Xe, Cd,…),
while others improve only slightly (e.g., P, S, As,…) and some get slightly
worse (e.g., Na, K, Rn,…). For the oxides, the improvement is significant for
materials with short oxygen bonds, when using the O_h_GW potential instead of
the standard PBE one (not shown).

##  S17 Consistency check with previous benchmarks for the all-electron data

All calculations in this work have been run independently and from scratch by
expert users of the respective codes, following a strict protocol. As the all-
electron calculations WIEN2k and FLEUR serve as references here, it is useful
to compare them to results obtained for the same crystals in an earlier
benchmark[29, 30]. Table S17.1 shows this comparison for the parameters of the
Birch-Murnaghan EOS: both for FLEUR and WIEN2k, there are only small relative
differences between the parameters obtained in Refs. 29, 30 and the ones
obtained in the present work. The only crystal with significant relative
deviations is Ne (FCC), which is not surprising given its very shallow
equation of state (small value of the bulk modulus, which is known[98] to lead
to a large uncertainty in the volume). This demonstrates good agreement
between the previous and current data, obtained from independent calculations.

A different view on the same data is represented in SI Fig. S17.1. It shows on
the horizontal axis the difference between the FLEUR and WIEN2k results for
these 29 crystals, expressed by the $\Delta$ metric, as obtained in Refs. 29,
30. On the vertical axis, the $\Delta$ metric for the same crystal and the
same two codes is shown, now using the data obtained in the present work.
Although the range on the horizontal axis does not extend much beyond the
threshold of good agreement of $\Delta=1$ meV/at that was used in Refs. 29,
30, the range on the vertical axis is significantly smaller.

The conclusion of this analysis is that although the agreement between the
all-electron codes FLEUR and WIEN2k was already very good in Refs. 29, 30
($\Delta\leq 1$ meV/at), and although the relative differences in the
observable properties that can be derived from the EOS are very small (Tab.
S17.1), the agreement between the two all-electron codes in the present work
is definitely even better than it was in previous works (SI Fig. S17.1).

Table S17.1: Table comparing the $V_{0}$, $B_{0}$ and $B_{1}$ parameters for the subset of 29 structures of those suggested in Ref. 29, 30 that are also present in our set of unaries, namely those that have a cubic FCC, BCC, SC or Diamond structure and that were treated without spin polarization. In particular, Fe(BCC), Mn(FCC), Ni(FCC) and Cr(BCC) have not been included because they are treated including spin polarization in Ref. 29, 30; all other structures not shown here were not considered in a cubic structure in Ref. 29, 30. The table reports the absolute percentage error on each parameter with respect to the data for the same code and the same crystal reported in Ref. 29, 30. |  | FLEUR | WIEN2k
---|---|---|---
Element | Structure |  $V_{0}$ error [%] |  $B_{0}$ error [%] |  $B_{1}$ error [%] |  $V_{0}$ error [%] |  $B_{0}$ error [%] |  $B_{1}$ error [%]
Ag | FCC | 0.01 | 1.59 | 1.98 | 0.06 | 0.86 | 6.97
Al | FCC | 0.01 | 1.22 | 5.69 | 0.10 | 0.72 | 1.15
Ar | FCC | 0.55 | 5.61 | 29.04 | 0.28 | 1.28 | 2.71
Au | FCC | 0.20 | 1.62 | 10.40 | 0.02 | 0.27 | 2.95
Ba | BCC | 0.16 | 2.18 | 9.86 | 0.26 | 0.05 | 30.67
Ca | FCC | 0.11 | 0.24 | 5.91 | 0.00 | 1.50 | 0.59
Cs | BCC | 0.19 | 0.39 | 4.82 | 0.20 | 1.45 | 38.97
Cu | FCC | 0.16 | 0.06 | 0.71 | 0.05 | 0.10 | 4.14
Ge | Diamond | 0.09 | 0.34 | 2.68 | 0.02 | 0.31 | 2.60
Ir | FCC | 0.00 | 0.82 | 0.07 | 0.03 | 0.12 | 1.38
K | BCC | 0.15 | 0.47 | 0.47 | 0.14 | 0.14 | 21.68
Kr | FCC | 0.37 | 4.43 | 38.54 | 0.62 | 5.21 | 35.59
Mo | BCC | 0.09 | 0.03 | 5.32 | 0.04 | 0.12 | 2.87
Nb | BCC | 0.12 | 0.59 | 12.32 | 0.03 | 0.97 | 3.71
Ne | FCC | 2.69 | 0.41 | 41.55 | 0.29 | 13.74 | 99.20
Pb | FCC | 0.15 | 0.45 | 1.86 | 0.09 | 0.16 | 4.52
Pd | FCC | 0.03 | 0.77 | 3.34 | 0.08 | 0.35 | 0.65
Po | SC | 0.03 | 0.60 | 9.91 | 0.04 | 0.13 | 0.39
Pt | FCC | 0.19 | 1.19 | 6.52 | 0.08 | 0.52 | 0.02
Rb | BCC | 0.07 | 0.86 | 0.83 | 0.39 | 0.55 | 53.53
Rh | FCC | 0.07 | 0.63 | 0.91 | 0.07 | 0.34 | 2.65
Rn | FCC | 0.04 | 7.79 | 14.45 | 0.48 | 4.23 | 34.33
Si | Diamond | 0.04 | 0.03 | 0.68 | 0.03 | 0.02 | 0.07
Sn | Diamond | 0.01 | 0.47 | 3.56 | 0.06 | 0.71 | 4.94
Sr | FCC | 0.71 | 3.27 | 37.98 | 0.68 | 3.43 | 5.89
Ta | BCC | 0.02 | 1.62 | 8.91 | 0.03 | 1.11 | 1.17
V | BCC | 0.07 | 1.20 | 1.87 | 0.07 | 0.16 | 2.31
W | BCC | 0.04 | 0.72 | 6.96 | 0.03 | 0.04 | 2.63
Xe | FCC | 0.11 | 4.06 | 6.18 | 0.36 | 2.00 | 12.28
Figure S17.1:  Correlation of the $\Delta$ metric on the 29 crystals listed in
Tab. S17.1. $x$ axis: $\Delta$ metric for these crystals between FLEUR and
WIEN2k, using the results from Refs. 29, 30. $y$ axis: same metric $\Delta$
for the same crystals and the same two codes, but using data from the present
work. The black solid line indicates $y=x$. We note that $\Delta\leq 1$
meV/atom was considered in Ref. 29 to indicate a good agreement.

##  S18 Discrepancies of formation energies computed from the current dataset

In SI Fig. S18.1 we report a histogram of the difference of the formation
energy obtained using the data for the two all-electron codes FLEUR and
WIEN2k, computed from the minimum-energy value of the EOS data curves. The
histogram is obtained considering the formation energy of all XnOm oxides,
using the lowest-energy unary of element X and of oxygen as the two endpoints
(in the case of oxygen, the lowest-energy non-magnetic unary in our dataset is
the simple cubic structure). The majority of the datapoints are in the visible
$x$ axis range, i.e., with an (absolute) discrepancy smaller than 50 meV/atom.
Nevertheless, several outliers are present: 52 out of the 576 materials
considered have a discrepancy larger than 50 meV/atom. The most outstanding
outliers are AtO3, PoO3, BiO3, Am2O, Pu2O, AmO and PuO. They have a
discrepancy larger than 1 eV/atom. For compounds containing Am and Pu, the
discrepancy is also due by the fact that the unary identified as having lowest
energy is different between the two codes.

As we discuss in the main text, the reason for this discrepancy is that our
workflows have been designed to guarantee consistent simulation parameters
among calculations for a given material at different volumes. However, when
considering different structures, we did not enforce any consistency between
simulation parameters, e.g., the choices of atomic radii for the AE codes
might be different in different systems. Especially, changes to the
core/valence separation from structure to structure may lead to larger
discrepancies in this comparison, because the differing relativistic
descriptions of core and valence electrons lead to different energy
contributions. Therefore, we recommend not to use our dataset to generate
plots like the one of SI Fig. S18.1 or, more generally, to avoid performing
data analysis that considers energy differences between different structures.
In this case, instead, one should design new appropriate workflows that can
ensure the consistency of simulation parameters among all relevant
calculations.

Figure S18.1: Histogram of the discrepancy between the formation energy
obtained from our reference dataset for the two codes FLEUR and WIEN2k. We
note that no correction is applied to the formation energies (as it is
typically done, for instance, for oxygen[25]). However, since we are only
considering differences in formation energies between two codes, these
corrections cancel out so they do not need to be considered explicitly. The
number of outliers outside of the visible $x$ axis range is reported in the
top right corner of the figure.