\section{Constrained random phase approximation}

The constrained random phase approximation as proposed in
\cite{aryasetiawan_frequency-dependent_2004s} has allowed for the first
principles determination of the effective local Coulomb interaction, the
``Hubbard $U$'', that is used as the bare interaction in a low-energy
Hamiltonian for a wide range of correlated materials ranging from transition
metals \cite{aryasetiawan_frequency-dependent_2004s,miyake_abinitio_2009s,
sasioglu_effective_2011s,shih_screened_2012s},
their oxides
\cite{miyake_screened_2008s,tomczak_realistic_2010s,martins_reduced_2011s,
vaugier_hubbard_2012s,sakuma_first-principles_2013s,werner_dynamical_2015s},
pnictides and chalcogenides
\cite{werner_satellites_2012s,nakamura_ab_2008s,miyake_d-_2008s,
imada_electronic_2010s,van_roekeghem_dynamical_2014s,van_roekeghem_inprep},
$f$-electron materials
\cite{karlsson_method_2010s,tomczak_rare-earth_2013s,nilsson_ab_2013s} to organic
conductors \cite{nomura_ab_2012s,nakamura_ab_2012s} and solid hydrogen
\cite{tomczak_effective_2009s}.  Probably the most important advantage of cRPA
over previous methods such as the constrained LDA
\cite{dederichs_ground_1984s,gunnarsson_density-functional_1989s,
anisimov_density-functional_1991s} or linear response
\cite{cococcioni_linear_2005s} schemes is the possibility of adapting $U$ to the
chosen low-energy model.  For reviews, see
\cite{biermann_dynamical_2014s,jiang_first-principles_2015s,imada_electronic_2010s}.

The underlying idea of the constrained random phase approximation is the
observation that the Hubbard $U$ should be obtained as the matrix element of a
{\it partially-screened} interaction in some localized basis.  This
partially-screened interaction is constructed such that screening at the RPA
level calculated for the low-energy Hamiltonian or for the full Coulomb
Hamiltonian lead to the same result for the fully-screened Coulomb interaction.

Mathematically, this can be obtained by constructing a partial polarization
function where certain screening processes are suppressed: $P_r = P - P_d$,
where $P_d$ is the polarization within the low-energy space, which at the RPA
level is a sum over particle-hole transitions that take entirely place within
the low-energy subspace. The partially-screened interaction is then obtained as
$W_r = v/(1-P_r v)$ where $v$ is the bare Coulomb interaction. Finally, the
matrix of Hubbard interactions $U_{m m^{'}m^{''} m^{'''}}$ is obtained as
matrix elements of $W_r$ within the basis of localized Wannier functions used
for the construction of the Hamiltonian.

Subtleties arise when one aims at constructing a multi-orbital Hamiltonian
which comprises more orbital degrees of freedom at the one-particle level than
in its interacting part.  The commonly applied strategy \cite{miyake_d-_2008s}
in this case is to construct Wannier functions for all orbitals, but exclude
screening transitions corresponding to the interacting degrees of freedom only.
However, ambiguities arise in this case when the correlated bands are entangled
with itinerant ones.  Several disentanglement strategies have been proposed
\cite{miyake_abinitio_2009s,nilsson_ab_2013s}, but it may appear somewhat
insatisfactory to deal with effectively an altered one-particle band structure,
in particular since by construction, the stronger the entanglement, the larger
the deviations.  Alternatively, one may think that in these cases it could be
advisable to rather reconsider the philosophy of cutting out only the few
correlated states out of a bigger number of entangled states living roughly at
the same energy.  Our work presents a well-defined solution following this
philosophy.

\section{Background on the benchmark compounds}

\subsection{Cerium}

Elemental cerium has intrigued the solid state community for decades due to its
isostructural volume-collapse transition.  This ``$\alpha$-$\gamma$''
transition is of first order below a critical temperature T$_c$ = 600K, and
ends in a second-order critical point at T$_c$.  At room temperature the
decrease in volume when going from the $\gamma$ to the $\alpha$ phase is as
large as 15$\%$. A Curie-Weiss magnetic susceptibility in the $\gamma$ phase
signals the presence of localized $4f$ electrons that delocalize in the
$\alpha$ phase where the magnetic susceptibility becomes Pauli-like.  In this
phase, the $4f$ electrons participate in the bonding (resulting in the smaller
volume) and the formation of quasiparticles, as seen in photoemission
\cite{wuilloud_electronic_1983s,wieliczka_high-resolution_1984s}.

The origin of this transition has been the subject of intense debate in the
literature. The ``Mott transition picture'' interprets it as a
delocalization-localization transition in the $4f$ shell only
\cite{johansson_a-g_1974s}.  The ``Kondo volume collapse'' (KVC) picture
\cite{allen_kondo_1982s,allen_alpha-gamma_1992s,lavagna_volume_1982s} in contrast
attributes the leading role to the hybridization between the $4f$ electrons and
($spd$-) itinerant states: the stronger hybridization in the $\alpha$-phase
leads to a high Kondo temperature and thus to screening of the $4f$ local
moment, while the high-volume $\gamma$-phase has a low Kondo temperature,
leading in practice to unscreened moments at intermediate temperatures.

Cerium has long served as a test material for methods that go beyond the LDA
\cite{sakuma_self-energy_2012s,casadei_density-functional_2012s} and in
particular for DMFT-based electronic structure calculations
\cite{held_cerium_2001s,zolfl_spectral_2001s,mcmahan_thermodynamic_2003s,
haule_cerium_2005s,amadon_cerium_2006s,streltsov_magnetic_2012s,litsarev_correlated_2012s,
bieder_thermodynamics_2014s,lanata_interplay_2014s,chakrabarti_alpha-gamma_2014s}.
Combined DFT+DMFT was found to be able to disentangle subtle
electronic-structure changes in the hybridization to elucidate the mechanism of
the transition, and give a description of spectral, optical and magnetic
properties in agreement with experiments.  It was also emphasized, however,
that at any relevant temperature the entropic contribution dominates over the
internal energy changes and pressure contribution such that the question of the
existence of a double minimum in the internal energy acquires a somewhat
academic character.  We mention very recent QMC calculations
\cite{devaux_electronic_2015s} that advocate the existence of a transition down
to zero temperature.

With the advent of more efficient Monte Carlo techniques for the DMFT equations
in the multi-orbital case in recent years, the problem has been extensively
revisited, and not only total energies, but also entropy and free energy have
been computed \cite{bieder_thermodynamics_2014s}.  Recent work has also focused
on the role of spin-orbit coupling
\cite{dong_orbital-dependent_2014s,lanata_interplay_2014s}.

While many-body calculations for cerium have thus been extremely successful,
$U$ values have generally been used as adjustable parameters, and several
authors have discussed the difficulties in determining the effective Hubbard
interactions from first principles.  This comes without surprise since this
material can be considered a ``worst-case'' scenario, precisely due to the
strong hybridization of the $4f$-states with itinerant conduction band states
present at the same energy.  Ref.~\onlinecite{amadon_screened_2014s} has
presented a detailed study of various choices of low-energy models, arriving at
the conclusion that any low-energy Hamiltonian of $f$-$spdf$ type would have
$U^{ff}$-values lower than 1 eV in cRPA. This illustrates the problem: there is
a clear consensus in the literature that a much larger value of 5 to 6 eV is
necessary to reproduce the available experimental findings.  Nilsson et
al.~\cite{nilsson_ab_2013s} proposed a disentanglement scheme which, by
essentially replacing the original band structure by an approximate one with
reduced hybridization, yields larger (though still too small) values of 4.3 eV
and 5.4 eV for $\alpha$ and $\gamma$ cerium respectively.

\subsection{NiO}

NiO has been put forward as a challenge to solid state theory and electronic
structure calculations since the very early days: Mott highlighted the
inadequacies of the band picture in describing the physics of this compound.
NiO is an antiferromagnetic insulator below a N\'eel temperature of 525 K.
Above T$_N$, the antiferromagnetic order disappears, but the material remains
insulating with a large gap of 4 eV that largely exceeds the energy scale of
the N\'eel temperature). Furthermore, its spectral properties are remarkably
insensitive to the presence or absence of magnetic order. This calls for a
theoretical description of the insulating nature without relying on magnetic
order.

The story became more involved with the seminal work of Zaanen, Sawatzky and
Allen \cite{zaanen_band_1985s}, who distinguished between Mott insulators where
the gap opens due to the Coulomb blocking between correlated states of the same
shell character, and charge-transfer insulators where the Hubbard interactions
are larger than the O $2p$ to Ni $3d$ charge-transfer energy such that the
first ionization states are in fact of O $p$ character.

In NiO, it was already known that a satellite feature at -8 eV binding energy
is mainly of $d$ character, while the low-energy states at the Fermi level have
large oxygen $2p$ character.  A series of photoemission (including resonant
photoemission) works elucidated the different transitions in detail
\cite{hufner_photoemission_1984s,sawatzky_magnitude_1984s,fender_covalency_1968s,
cheetham_magnetic_1983s,shen_aspects_1990s,shen_electronic_1991s,shen_angle_1991s,
fujimori_multielectron_1984s,tjeng_giant_1991s,fink_electron_1994s}.

NiO was treated within many-body perturbation theory within the GW
approximation by Aryasetiawan and Gunnarsson
\cite{aryasetiawan_electronic_1995s}.  The authors found the gap to open
within this approximation, but also pointed out failures concerning the gap
character which remains Mott-Hubbard like, and satellite structures related to
the strongly-atomic character of the $3d$ states.

The first DMFT works on NiO used a Hamiltonian where only the TM $3d$ states
were included \cite{ren_lda+dmft_2006s}. Soon thereafter, it was pointed out
that this was questionable due to the $d$-ligand interplay at the Fermi level.
Several subsequent works included ligand states as itinerant states in the
one-body Hamiltonian \cite{korotin_construction_2008s,kunes_nio:_2007s}.
However, none of the standard double counting corrections appeared to be
appropriate in this case \cite{karolak_double_2010s}, and NiO still appears to
be a challenge to DMFT-based or DMFT-like theories
\cite{manghi_multi-orbital_2014s}.

\subsection{Actinide oxides: UO$_2$, PuO$_2$ and NpO$_2$}

The actinide dioxides UO$_2$, PuO$_2$ and NpO$_2$ are insulating materials with
similar gaps (of
2.1 eV \cite{schoenes_optical_1978s},
2.9 eV \cite{mccleskey_optical_2013s}
and
2.8 eV \cite{mccleskey_optical_2013s}
respectively, but quite different ground states.  UO2 has been and remains the
subject of intense interest from the applied materials perspective thanks to
its application as a nuclear fuel. Electronic structure calculations thus focus
not only on questions of stability of different structures, but also attempt to
assess the physics of impurities resulting from the radioactive decay
processes.  We mention for example the study of He atom incorporation into
UO$_2$ crystals by \cite{gryaznov_ab_2009s}, where LDA+U with $U=4.6$ eV was
employed.  For an extensive review and recent work see
\cite{vathonne_etude_2014s}.  Interestingly, corrections to simple DFT-LDA
calculations turn out to be necessary not only for the description of spectral
properties but even in assessing the energetics, and UO$_2$ has become a
playground for electronic structure approaches.  While the value of $U$ has
typically been considered an adjustable parameter or guessed from spectroscopy
experiments, we note that the values employed coincide remarkably well with our
results discussed in the main text.

The series of actinide dioxides has also attracted much interest due to their
exotic low-temperature (multi-polar) ordered states.  While in PuO$_2$ the
$f^4$ configuration leads to a non-magnetic ground state, UO$_2$ has a complex
antiferromagnetic order below a N\'eel temperature of 30.8 K
\cite{wen_density_2013s}, with a transverse 3-k arrangement driven by
quadrupolar interactions. NpO$_2$ on the other hand displays a longitudinal 3-k
order below an ordering temperature of 25 K, with a rank-5 magnetic multipole
(``triakontadipole'').  Despite being difficult due to the noncollinear
character of the order and high-rank order parameters, LDA+U calculations are
successful in describing the multi-polar order with interaction values of $U=4$
eV \cite{suzuki_first-principles_2010s}.

Their complex electronic structure has always been considered a challenge also
for purely functional-based first-principles calculations, and it is not
surprising that the actinide oxides became a playground and testbed for hybrid
functionals (see e.g., \cite{kudin_hybrid_2002s,roy_dispersion_2008s}).  A
comprehensive review of various electronic structure calculations of the
actinide oxides is given in \cite{wen_density_2013s}.

\section{Detailed derivation of the shell-folded Hamiltonian}

Here we present the derivation of the shell-folded Hamiltonian by the
Peierls-Feynman-Bogoliubov variational principle in a step-by-step manner and
determine the form of our effective interactions
$\tilde{U}^{dd} = U^{dd} - U^{pd}$ and $\tilde{U}^{pp} = U^{pp} - U^{pd}$.

Consider the $dp$ Hamiltonian
\begin{align*}
 h_{\bf R} = & \frac{\tilde{U}^{dd}}{2} \sum_{\begin{subarray}{c}
                                                       (m\sigma) \ne (m'\sigma^{\prime})\\
                                                       m,m' \in {\cal C}
                                                      \end{subarray}}
                                            \hspace{-.4cm}
n_{m \sigma} n_{m' \sigma'}
+ \frac{\tilde{U}^{pp}}{2} \sum_{\begin{subarray}{c}
                                                       (m\sigma) \ne (m'\sigma^{\prime})\\
                                                       m,m' \in {\cal L}
                                         \end{subarray}}
                                            \hspace{-.4cm}
                                                      n_{m \sigma} n_{m' \sigma'}
\\  + &
  \tilde{U}^{dp}           \sum_{\begin{subarray}{c}
                                                       m \in {\cal C}, \sigma \\
                                                       m' \in {\cal L}, \sigma^{\prime}
                                         \end{subarray}}
                                            \hspace{-.4cm}
                                                      n_{m \sigma} n_{m' \sigma'},
\end{align*}

and the auxiliary Hamiltonian
\begin{align*}
 \tilde{h}_{\bf R} = & \frac{\tilde{U}^{dd}}{2} \sum_{\begin{subarray}{c}
                                                       (m\sigma) \ne (m'\sigma^{\prime})\\
                                                       m,m' \in {\cal C}
                                                      \end{subarray}}
                                            \hspace{-.4cm}
n_{m \sigma} n_{m' \sigma'}
+ \frac{\tilde{U}^{pp}}{2} \sum_{\begin{subarray}{c}
                                                       (m\sigma) \ne (m'\sigma^{\prime})\\
                                                       m,m' \in {\cal L}
                                         \end{subarray}}
                                            \hspace{-.4cm}
                                                      n_{m \sigma} n_{m' \sigma'}.
\end{align*}
where $\tilde{U}^{dd/pp}$ are the effective Coulomb interactions for the $d$
and $p$ manifolds. We find their optimal values by minimizing $\tilde{F} =
F[\tilde{h}_{\bf R}] + \ave{h_{\bf R} - \tilde{h}_{\bf R}}_0$ with respect to
$\tilde{U}^{dd}$ and $\tilde{U}^{pp}$, where $\ave{\ldots}_0$ refers to the
average taken with the auxiliary Hamiltonian.
Let us also extend the definition
\begin{eqnarray*}
D_Z^{XY} = \frac{\partial}{\partial \tilde{U}^{Z}}
         \sum_{(i\sigma) \neq (j\sigma')} \ave{ n^{X}_{i \sigma} n^{Y}_{j \sigma'}}_0 .
\end{eqnarray*}
where $X,Y,Z=d$ or $p$.

Then, from the stationarity condition
$ \partial \tilde{F} / \partial \tilde{U}^{dd} = 0 $
we obtain
\begin{align*}
\label{eq:U0_shellfolded}
\tilde{U}^{dd} = U^{dd} +
\frac{ 2 U^{dp} D_d^{dp}
       + (U^{pp} - \tilde{U}^{pp}) D_d^{pp}
      }{D_d^{dd}}.
\end{align*}
and an analogous equation for $\tilde{U}^{pp}$.

Solving the system of equations, we have for
$\Delta = \tilde{U}^{dd} - U^{dd}$,
\begin{align*}
\Delta & =
\frac{\Delta D_p^{dd} - 2 U^{dp} D_p^{dp}}
     {D_p^{pp}}
\frac{D_d^{pp}}
     {D_d^{dd}}
+
\frac{2 U^{dp} D_d^{dp} }
     {D_d^{dd}}
\end{align*}

Simplifying this expression gives
\begin{align*}
\frac{\tilde{U}^{dd} - U^{dd}}
     {U^{pd}} & =
\frac{2 D_d^{pd} D_p^{pp} - 2 D_p^{pd} D_d^{pp}}
     {  D_d^{dd} D_p^{pp} -   D_p^{dd} D_d^{pp}}.
\end{align*}

Assuming that the total number of particles is invariant with respect to
changes in the effective parameters, namely that
\begin{align*}
  \frac{\partial \ave{N(N-1)}}{\partial \tilde{U}^{dd}}
= \frac{\partial \ave{N(N-1)}}{\partial \tilde{U}^{pp}}
= 0,
\end{align*}
or equivalently that
\begin{align*}
 D_d^{dd} + 2 D_d^{dp} + D_d^{pp} = 0,\\
 D_p^{dd} + 2 D_p^{dp} + D_p^{pp} = 0,
\end{align*}
we can simply the above to
\begin{align*}
\frac{\tilde{U}^{dd} - U^{dd}} {U^{pd}} & = -1.
\end{align*}

Thus, we obtain that
\begin{align*}
\tilde{U}^{dd} = U^{dd} - U^{pd}, \\
\tilde{U}^{pp} = U^{pp} - U^{pd}.
\end{align*}

While the above derivation is presented for spin- and orbital-averaged
interactions for simplicity, it is equally applicable to the orbital-resolved
interactions $U_{m \sigma m' \sigma'}$. As long as we assume the intershell interaction
to be an orbital-independent average $U^{dp}$, the resulting effective
interactions are obtained from an effective reduced Slater integral
$\tilde{F}^0 = F^0 - U^{dp}$.  Note that quadratic terms $H_0$ are unchanged.
As we are devising an effective model for the \emph{local} $h_{\bf R}$, and as
our summations run explicitly over all orbitals contained therein, we do not
have to consider either the number of equivalent atoms within the unit cell or
its coordination.

\section{Construction of transition metal centred Wannier functions}

In the principal part of the Letter, we derived a low-energy model that we
showed to be successful in treating intershell interactions in addition to the
intrashell interactions.  Here, using NiO as an example, we describe the
construction of the ``ring'' orbitals centred on the correlated site as a
superposition of the ligand $p$ orbitals.  This ring orbital construction is
inspired by the cell-perturbation methods
\cite{jefferson_derivation_1992s,feiner_effective_1996s,raimondi_effective_1996s}
that emerged as a consequence of the Zhang-Rice construction
\cite{zhang_effective_1988s} for the Cu superconductors.

\section{NiO}

\subsection{Construction of ring orbitals}

NiO adopts the rock-salt structure with a face-centred cubic lattice of Ni atoms
intercalated with a face-centred cubic lattice of O atoms. As the conventional unit
cell contains four formula units, for simplicity, we choose to work globally in
the Cartesian coordinate system with the primitive unit cell vectors
\begin{align*}
 {\bf a} = (1/2, 1/2, 0), \\
 {\bf b} = (1/2, 0, 1/2), \\
 {\bf c} = (0, 1/2, 1/2).
\end{align*}
Our basis consists of one Ni atom at $(0,0,0)$ and one O atom at $(1/2, 1/2,
1/2)$.  We define the corresponding $k$-vectors ${\bf k_{a}}, {\bf k_{b}}, {\bf
k_{c}}$ such that ${\bf k_{i}} \cdot {\bf j} = 2\pi \delta_{ij}$.

In order to find matrix elements of $H(\mathbf{k})$, we first determine the Ni
$3d$ to O $2p$ hopping directions. Using the Slater-Koster parametrization
\cite{slater_simplified_1954s}, we write the matrix elements of the real-space
tight-binding Hamiltonian $H(\mathbf{R})$ as functions of the $p$-$p$, $p$-$d$
and $d$-$d$ two-center bond integrals that give the strength of bonding in a
given direction, reflecting the relative orbital lobe positions. We use the
standard notation of $\sigma$-bonding if the lobes point to each other and
$\pi$-bonding if the lobes are parallel. The resulting matrix is finally
Fourier-transformed to $H(\mathbf{k})$ in momentum space. We note that, besides
hopping, the cubic crystal field splitting of Ni $3d$-states into $e_{g}$ and
$t_{2g}$ (parametrized in the usual notation as $10Dq$) is also included in
$H(\mathbf{k})$. The final $H(\mathbf{k})$ in matrix form is given in Table
\ref{tab:hk}. The diagonalization of this Hamiltonian yields the tight-binding
band structure.

\begin{sidewaystable}[htdp] \scriptsize
\caption{
The Hamiltonian operator in momentum space, $H(\mathbf{k})$, in the basis
$(\text{O} \, 2p_x, \text{O} \, 2p_y, \text{O} \, 2p_z,
\text{Ni} \, 3d_{xy}, \text{Ni} \, 3d_{yz}, \text{Ni} \, 3d_{xz},
\text{Ni} \, 3d_{x^2-y^2}, \text{Ni} \, 3d_{z^2})$.
To improve readability we employ the abbreviations
$s_i = \sin(\frac{k_i}{2})$ and $c_i = \cos(\frac{k_i}{2})$.
}
\label{tab:hk}
\centering
\begin{tabular}{cccccccc}
 $ \epsilon_p+2 \text{pp$\sigma $} c_x \left(c_y+c_z\right) $ &
 $ -2 \text{pp$\sigma $} s_x s_y $ &
 $ -2 \text{pp$\sigma $} s_x s_z $ &
 $ -2 i \text{pd$\pi $} s_y $ &
 $ 0 $ &
 $ -2 i \text{pd$\pi $} s_z $ &
 $ -i \sqrt{3} \text{pd$\sigma $} s_x $ &
 $ i \text{pd$\sigma $} s_x $ \\
 $ -2 \text{pp$\sigma $} s_x s_y $ &
 $ \epsilon_p+2 \text{pp$\sigma $} (c_x c_y + c_y c_z) $ &
 $ -2 \text{pp$\sigma $} s_y s_z $ &
 $ -2 i \text{pd$\pi $} s_x $ &
 $ -2 i \text{pd$\pi $} s_z $ &
 $ 0 $ &
 $ i \sqrt{3} \text{pd$\sigma $} s_y $ &
 $ i \text{pd$\sigma $} s_y $ \\
 $ -2 \text{pp$\sigma $} s_x s_z $ &
 $ -2 \text{pp$\sigma $} s_y s_z $ &
 $ \epsilon_p+2 \text{pp$\sigma $} (c_x c_z + c_y c_z) $ &
 $ 0 $ &
 $ -2 i \text{pd$\pi $} s_y $ &
 $ -2 i \text{pd$\pi $} s_x $ &
 $ 0 $ &
 $ -2 i \text{pd$\sigma $} s_z $ \\
 $ 2 i \text{pd$\pi $} s_y $ &
 $ 2 i \text{pd$\pi $} s_x $ &
 $ 0 $ &
 $ -4 \text{Dq}+\epsilon_d+3 \text{dd$\sigma $} c_x c_y $ &
 $ 0. $ &
 $ 0. $ &
 $ 0 $ &
 $ \sqrt{3} \text{dd$\sigma $} s_x s_y $ \\
 $ 0 $ &
 $ 2 i \text{pd$\pi $} s_z $ &
 $ 2 i \text{pd$\pi $} s_y $ &
 $ 0. $ &
 $ -4 \text{Dq}+\epsilon_d+3 \text{dd$\sigma $} c_y c_z $ &
 $ 0. $ &
 $ \frac{3}{2} \text{dd$\sigma $} s_y s_z $ &
 $ -\frac{\sqrt{3}}{2} \text{dd$\sigma $} s_y s_z $ \\
 $ 2 i \text{pd$\pi $} s_z $ &
 $ 0 $ &
 $ 2 i \text{pd$\pi $} s_x $ &
 $ 0. $ &
 $ 0. $ &
 $ -4 \text{Dq}+\epsilon_d+3 \text{dd$\sigma $} c_x c_z $ &
 $ -\frac{3}{2} \text{dd$\sigma $} s_x s_z $ &
 $ -\frac{\sqrt{3}}{2} \text{dd$\sigma $} s_x s_z $ \\
 $ i \sqrt{3} \text{pd$\sigma $} s_x $ &
 $ -i \sqrt{3} \text{pd$\sigma $} s_y $ &
 $ 0 $ &
 $ 0 $ &
 $ \frac{3}{2} \text{dd$\sigma $} s_y s_z $ &
 $ -\frac{3}{2} \text{dd$\sigma $} s_x s_z $ &
 $ 6 \text{Dq}+\epsilon_d+\frac{3}{4} \text{dd$\sigma $} (c_x c_z + c_y c_z) $ &
 $ \frac{\sqrt{3}}{4} \text{dd$\sigma $} \left(c_x-c_y\right) c_z $ \\
 $ -i \text{pd$\sigma $} s_x $ &
 $ -i \text{pd$\sigma $} s_y $ &
 $ 2 i \text{pd$\sigma $} s_z $ &
 $ \sqrt{3} \text{dd$\sigma $} s_x s_y $ &
 $ -\frac{\sqrt{3}}{2} \text{dd$\sigma $} s_y s_z $ &
 $ -\frac{\sqrt{3}}{2} \text{dd$\sigma $} s_x s_z $ &
 $ \frac{\sqrt{3}}{4} \text{dd$\sigma $} \left(c_x-c_y\right) c_z $ &
 $ 6 \text{Dq}+\epsilon_d+\frac{1}{4} \text{dd$\sigma $} (c_y c_z + c_x \left(4 c_y+c_z\right)) $ \\
\end{tabular}
\end{sidewaystable}

We now construct Ni-centred ring orbitals from the three O $2p$ states. This is
similar to what is done in typical cluster calculation ligand orbitals
\cite{haverkort_multiplet_2012s}. The bonding cluster orbitals corresponding to the
NiO$_6$ octahedron are shown in Fig.~\ref{fig:nio_cluster_orb}.

\begin{figure}[htpb]
  \centering
    \includegraphics[width=1.00\columnwidth]{nio_cluster_orbitals.pdf}
    \caption{The NiO$_6$ octahedron cluster orbitals.
             \label{fig:nio_cluster_orb}}
\end{figure}

One key difference between our ring orbitals and the cluster orbitals, however,
is that the cluster orbitals are constructed entirely in real space from
combinations of a finite number of orbitals. In our case, we use a coherent
superposition of orbitals on an all (infinite) sites, giving rise to a
$k$-dependence (see the $\sin$ terms in Eq.~\ref{eq:p_to_ring}).  From the
three available oxygen states we can construct two $e_{g}$-type rings, which
locally hybridize with the Ni $e_{g}$ orbitals and a third, $a_{1g}$ (i.e.,
$s$-like) ring which is locally non-bonding with the Ni $e_{g}$ orbitals.  The
form of these rings is shown in Fig.~\ref{fig:nio_ring_orb}.

\begin{figure}[htpb]
  \centering
    \includegraphics[width=1.00\columnwidth]{nio_ring_orbitals.pdf}
    \caption{The NiO ring orbitals. The first two orbitals have $e_{g}$ symmetry
             while the last is of $a_{1g}$ symmetry.
             \label{fig:nio_ring_orb}}
\end{figure}

Mathematically, we perform a unitary transformation on the $2p$ subspace of the
original Hamiltonian $H$ to $\tilde{H}$ by $\tilde{H} = T H {T}^\dagger$. The
transformation matrix $T$ takes us from the
$(p_x,p_y,p_z,d_{xy},d_{yz},d_{xz},d_{x^2-y^2},d_{z^2})$ to
$(\tilde{p}_{x^2-y^2},\tilde{p}_{z^2},\tilde{p}_s,d_{xy},d_{yz},d_{xz},d_{x^2-y^2},d_{z^2})$
basis, and is given by

\begin{equation}
\label{eq:p_to_ring}
T = \begin{pmatrix}
  \frac{s_x}{\sqrt{2}}  &
 -\frac{s_y}{\sqrt{2}}  &
 0 & 0 & 0 & 0 & 0 & 0 \\
  \frac{s_y}{\sqrt{6}}  &
  \frac{s_x}{\sqrt{6}}  &
  \frac{-2 s_z}{\sqrt{6}}  &
 0 & 0 & 0 & 0 & 0 \\
  \frac{s_y}{\sqrt{3}}  &
  \frac{s_x}{\sqrt{3}}  &
  \frac{s_z}{\sqrt{3}}  &
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
\end{pmatrix}
\end{equation}

We have used the abbreviation $s_i = \sin(\frac{k_i}{2})$.  The transformation
is chosen such that the ring with $d_{x^2-y^2}$ symmetry is exact, while the
$d_{z^2}$ and $s$-like rings are approximate in order to orthogonalize them.

In the original O-centred $p$ basis, the hybridization elements decay with
distance as shown in Fig.~\ref{fig:p_hopping}. Note that the distances are
measured in conventional lattice units. As we choose to construct a
maximally-symmetric primitive cell as described above, the O atoms closest to
the Ni atom at (0,0,0) actually lie in the nearest-neighbour unit-cells rather
than within the (0,0,0) unit cell. This is reflected in the peak in $t_{t_{2g}
p}$ and $t_{e_{g} p}$ at (0,0,$\frac{1}{2}$). As the $d$-subspace of the
Hamiltonian is untouched by the transformation, the $t_{d d}$ hybridizations in
the original $p$ basis and in the ring basis are identical.

\begin{figure}[htpb]
  \centering
    \includegraphics[width=1.00\columnwidth]{NiO_p_t_decay.pdf}
    \caption{Matrix elements of the initial $dp$ Hamiltonian $H(\mathbf{k})$ in
             real space plotted as a function of distance (measured in conventional unit
             cell lengths).
             \label{fig:p_hopping}}
\end{figure}

In the Ni-centred ring basis, the decay of the hybridization elements is
plotted in Fig.~\ref{fig:ring_hopping}. The $s$-like ring orbital is always
non-bonding locally; however, given that the symmetry of this orbital is not
imposed exactly in order to ensure unitarity of the transformation, it
continues to hybridize weakly with the Ni $d$ orbitals. For further symmetry
reasons, the $d_{x^2-y^2}$-like ring orbital alternates between being bonding
and anti-bonding over second nearest-neighbours.
\begin{figure}[htpb]
  \centering
    \includegraphics[width=1.00\columnwidth]{NiO_ring_orb_t_decay.pdf}
    \caption{Matrix elements of the transformed $d$-ring Hamiltonian
             $\tilde{H}(\mathbf{k})$ in real space plotted as a function of distance
             (measured in conventional unit cell lengths).
      \label{fig:ring_hopping}}
\end{figure}

The transformation of the $p$ orbitals into ring orbitals also modifies the
interaction parameters $U^{pp}$. Notably, the intra-ring interaction decreases
as the rings become more delocalized, while inter-ring interactions increase as
overlap between orbitals increases. The ring-$d$ interactions are unchanged due
to symmetry reasons.

\section{UO2}

In Fig.~\ref{fig:uo2_cluster_orb}, we show the 7 bonding ligand orbitals for
the eight-fold coordinated U $5f$ states. We can construct ligand orbitals by
linear combination of O $2p$ states matching the respective symmetry of the U
$5f$ orbitals as was done for NiO above. However, unlike for a finite size
cluster, we cannot create Bloch states associated to each of the ligand
orbitals (in the unit cell we have only six oxygen states).  Hence, in order to
construct ring-orbitals one has to choose which six basis states to include, as
in NiO.

\begin{figure}[htpb]
  \centering
    \includegraphics[width=1.00\columnwidth]{uo2_cluster_orbitals.pdf}
    \caption{The 7 bonding ligand orbitals constructed as linear combinations
             of the 8 $p$ orbitals with the symmetry of the corresponding $f$ orbital.
             \label{fig:uo2_cluster_orb}}
\end{figure}

Finally, we give the computational details for the cluster model calculations
for UO$_2$ that result in the spectral functions shown in Fig.~2 of the main
text. In the calculation, both $f$- and $p$-states, separated by the
charge-transfer energy given by LDA, are included but are not allowed to
hybridize. The parameters used for the UO$_2$ cluster model are
$E_{a_u} = 0.349$,
$E_{t_{1u}} = -0.039$,
$E_{t_{2u}} = -0.077$,
$\Delta_{fp} = 5.189$,
which were obtained from LDA and $\zeta_{5f} = 0.27$, obtained from an atomic
HF calculation. The hybridization parameters $V(a_u)=V(t_{1u})=V(t_{2u})=0$.
All values are in eV.

\section{Overview of U values used in the literature}

\begin{table}[htdp]
\caption{Screened Coulomb interactions used in the literature. All values in eV.}
\label{tab:literature_u}
\centering
\begin{tabular}{ccc}
\multicolumn{2}{c}{Interaction parameters} & References \\
\hline
$U = 4.5$         &  $J = 0.51$        &   \onlinecite{dorado_dft+u_2009s,geng_structural_2007s,geng_first-principles_2011s} \\
$\bar{U} = 4.5$   &                    &   \onlinecite{dudarev_understanding_2000s} \\
$\bar{U} = 4.5$   &  $\bar{J} = 0.51$  &   \onlinecite{dudarev_electronic_1998s} \\
$U = 4.5$         &  $J = 0.54$        &   \onlinecite{dudarev_effect_1997s} \\
$U = 4.6$         &  $J = 0.5 $        &   \onlinecite{gryaznov_density_2010s} \\
$U = 7.7$         &                    &   \onlinecite{gunnarsson_electronic_1988s} \\
$U = 4.0$         &  $J = 0.51$        &   \onlinecite{hongliang_optical_2010s} \\
$U = 4.5$         &  $J = 0.5 $        &   \onlinecite{jollet_electronic_1997s} \\
$U = 3.0,6.5$     &  $J = 0.7 $        &   \onlinecite{kolorenc_electronic_2015s} \\
$U = 4.5-4.75$    &  $J = 0.54$        &   \onlinecite{kotani_systematic_1992s} \\
$U = 4.0$         &                    &   \onlinecite{petit_electronic_2010s} \\
$U = 8.0$         &                    &   \onlinecite{shim_x-ray_2009s} \\
$U = 6.0$         &  $J = 0.5-0.6$     &   \onlinecite{yin_electronic_2011s} \\
$U = 3.0$         &                    &   \onlinecite{yin_origin_2008s} \\
$U = 4.0$         &  $J = 0-0.5$       &   \onlinecite{suzuki_microscopic_2013s} \\
\end{tabular}
\end{table}

In Table \ref{tab:literature_u} we present an overview of Hubbard U values used
for calculations for UO$_2$ in the literature. It is remarkable how close most
of these values are to the result of our shell-folded theory.

Nevertheless, a recent work by Koloren\v{c} {\it et
al.}~\cite{kolorenc_electronic_2015s} presented LDA+DMFT calculations using a
larger value of $F^0=6.5$eV.  The authors show that the observed spectral
function closely resembles spectra calculated within the Hubbard I
approximation using a smaller $U$ of 3 eV. Apart from a satellite structure at
about 10 eV, the spectral features between these two theories presented an
impressive match.  They also note that these $U$ values are slightly on the
small side relative to experiment. This is consistent with the observation in
Ref.~\onlinecite{suzuki_microscopic_2013s} that the best description of the
low-temperature ordered phase is obtained within LDA+U with $U = 4$ eV.
Koloren\v{c} {\it et al.}~add an interesting discussion of the differences in
the treatment of screening in LDA+U and LDA+Hubbard I.  They argue that
screening by ligand states is absent at the LDA+U or LDA+Hubbard I levels such
that the $U$ used must already incorporate these screening processes.  For
LDA+DMFT, on the other hand, they argue that some screening comes in through
the coupling of the $f$-states to a the dynamical bath, requiring a larger $U$.
Our cRPA calculations for the reduction of $U$ by screening through the ligands
agrees well with these estimates.

For completeness we note that there is an additional effect. $U$ values will
depend on the basis used with the two extremes being a set of localized
atomic-like $f$- and $p$-orbitals and a block-diagonal basis where $f$-$p$
hybridization enters effectively. We illustrate this on a simple two-orbital
model with atomic-like $f$- and $p$-orbitals, hybridizing with strength $V$.
The eigenvalues read, in terms of the energies of the local levels $\epsilon_f$
and $\epsilon_p = \epsilon_f - \Delta$  of the $f$ and $p$ state respectively:
\begin{eqnarray}
E = \frac{1}{2}(\epsilon_f + \epsilon_p)
 \pm \frac{1}{2} \sqrt{
(\epsilon_f - \epsilon_p)^2 + 4 |V|^2
}
\end{eqnarray}
and the eigenvectors can be expressed in terms of the atomic
states $| f \rangle$, $| p \rangle$:
\begin{eqnarray}
| \pm \rangle
= \cos(\theta) |f \rangle \pm \sin(\theta) | p \rangle
\end{eqnarray}
with
$\tan(\theta) = -\frac{\Delta}{2 V} + \sqrt{(\frac{\Delta}{2 V})^2 + 1}$.
Since $U^{ff}$ transforms as $\cos(\theta)^4$, the values in an atomic-like
basis are reduced when transforming to a basis that block diagonalizes the $fp$
Hamiltonian. The precise amount of reduction depends on the orbitals since the
hybridization is strongly orbital-dependent.

\section{Relation between cluster calculations and DMFT
calculations}

The rewriting of the Hamiltonian presented in the main text of our paper also
allows for a deeper understanding of the relation between cluster model
calculations and DFT+DMFT. Indeed, by definition, the cluster model assumes the
total particle number to be fixed to specific integer values, and diagonalizes
the cluster Hamiltonian for those. In DMFT, in general, the bath introduces the
possibility of charge fluctuations and the above assumption becomes less
obvious. A cluster approximation for insulators is well justified since the
bath hybridization is small. For metallic systems, restricting the Hamiltonian
the low-energy correlated subspace avoids the problem. The Zhang-Rice-like
construction presented above provides a way to do so explicitly.  Otherwise,
when the O $p$-states are to be kept explicitly in the description, the shifts
induced by the changes in the particle number might have to be taken into
account explicitly.

\bibliographystyle{plain}\input{letter_supp.bbl}

\end{document}