# Density analysis for estimating the degree of on-site correlation on
transition-metal atoms in extended systems

Leila Kalantari  Fabien Tran  Peter Blaha  Institute of Materials Chemistry,
Vienna University of Technology, Getreidemarkt 9/165-TC, A-1060 Vienna,
Austria

###### Abstract

In the context of the modified Becke-Johnson (mBJ) potential, we recently
underlined that $\bar{g}$, the average of $\left|\nabla\rho\right|/\rho$ in
the unit cell, has markedly different values in transition-metal oxides and
pure transition metals [Tran et al., J. Appl. Phys. 126, 110902 (2019)].
However, since $\bar{g}$ is a constant it is not able to provide local
information about a particular atom in the system. Furthermore, while
$\overline{g}$ can be used only for periodic bulk solids, a local (i.e.,
position-dependent) version would allow us to consider also low-dimensional
systems and interfaces. Such a local function has been proposed by Rauch et
al. [J. Chem. Theory Comput. 16, 2654 (2020)] for the local mBJ potential.
Actually, a local version of $\overline{g}$, or of another similar quantity
like the reduced density gradient $\overline{s}$, could also be used in the
framework of other methods. Here, we explored the idea to use such a local
function $\tilde{g}$ (or $\tilde{s}$), defined as the average of $g$ (or $s$)
over a certain region around a transition-metal atom, to estimate the degree
of on-site correlation on this atom. We found a large difference in our
correlation estimators between non-correlated and correlated materials,
proving its usefulness and reliability. Our estimators can subsequently be
used to determine whether or not a Hubbard $U$ on-site correction in the
DFT+$U$ method should be applied to a particular atom. This is particularly
interesting in cases where the degree of correlation of the transition-metal
atoms is not clear, like interfaces between correlated and non-correlated
materials or oxygen-covered metal surfaces. In such cases, our estimators
could also be used for an interpolation of $U$ between correlated and non-
correlated atoms.

##  I Introduction

Density functional theory (DFT) [1, 2] is the main computational tool to
perform electronic structure calculations on systems of realistic complexity,
however, one has to be careful when dealing with strongly correlated systems
like Mott insulators. Prototypical such examples are the antiferromagnetic
(AFM) $3d$-transition-metal oxides (TMO), for which DFT predictions using
standard semilocal approximations for the exchange-correlation energy are very
inaccurate (typically, a too small or even absent band gap, and a too small
magnetic moment [3]). More accurate DFT approximations, like hybrid
functionals [4, 5, 6, 7], or beyond-DFT methods, like DFT plus dynamical mean
field theory [8, 9, 10] or reduced density matrix functional theory [11, 12],
can improve significantly the description of correlated systems, but they are
prohibitively expensive for large systems. However, at a practical level there
exist DFT methods that are barely more expensive than standard semilocal
functionals, and therefore suitable for calculations on large systems, but
also more reliable for systems with strongly correlated electrons. This is the
case for instance of the modified Becke-Johnson (mBJ) potential [13, 14] and
DFT+$U$ [15, 16, 17], both being as accurate as the hybrid functionals for
strongly correlated systems.

In DFT+$U$, the Hubbard-model parameters to represent the on-site screened
Coulomb (Hubbard $U$) and exchange (Hund $J$) interactions are used. $U$,
which is usually larger than $J$, is a measure of the on-site interaction
between two electrons in the same electronic shell and is defined as

$U=E(d^{n+1})+E(d^{n-1})-2E(d^{n}).$ (1)

It is the Coulomb energy cost to transfer one electron from one site to the
other such that the number of on-site interactions between two electrons is
increased by one (from $n^{2}-n$ to $n^{2}-n+1$). In Eq. (1), $E(d^{n})$ is
the energy of an atom with $n$ localized electrons in the $d$ (or $f$) shell.
In order to achieve reasonable results with DFT+$U$, two important points
should be considered.

First, although a qualitative improvement of the results for correlated
systems can be obtained by DFT+$U$, the results depend of course on the value
of the parameters $U$ and $J$. Their values are often empirically calibrated
such that the result for a property (e.g., band gap or oxidation energy)
matches experiment. However, they can also be obtained ab initio by some
methods like constrained-LDA [18, 19, 20], constrained-RPA [21, 22, 23, 24,
25], or from linear response [26, 27]. Nevertheless, there are still some
ambiguity as well as freedom in the numerical implementation of these methods.
For instance, localized $d$ or $f$ orbitals in solids usually hybridize with
other valence $sp$ orbitals, causing an entangled band structure with
considerable band widths, and it is difficult to uniquely define the localized
states in solids.

Second, it is not even always clear, or known in advance, if $U$ should be
applied or not. In contrast to AFM TMO where a large $U$ value of 6$-$8 eV is
required, in a pure transition metal (TM) the $d$ electrons are itinerant,
i.e. only weakly correlated, and no $U$ needs to be used in principle.
However, there are of course intermediate cases, and furthermore in a given
(complicated) system the degree of correlation may vary from one atom to the
other even of the same type. This is the case, for instance, when a surface
system consists of a TMO layer adsorbed on top of a pure metal; the degree of
correlation on a TM atom is expected to decrease when going from the surface
(TMO-like) deep into the bulk (pure TM-like).

In the present work we will show that quantities depending on the electron
density $\rho$ (that we will call correlation estimators) can be used to
distinguish between correlated and non-correlated $3d$ TM atoms in any kind of
systems. Bulk solids, interfaces, and surfaces will be considered. Such a
quantity could in principle be used to determine, at least qualitatively,
whether or not a Hubbard $U$ correction should be applied on a TM atom.

We mention that Wang and Jiang proposed in Ref. [28] to use $\rho$ to
parameterize $U$ and $J$. More specifically, a simple average of $\rho$,
calculated either in the whole unit cell or only in a sphere surrounding the
atom, was used as screening parameter in the Slater integrals that enter into
the expression of $U$ and $J$. They applied the method for calculating the
energy difference between the AFM and ferromagnetic (FM) states of strongly
correlated bulk solids.

The paper is organized as follows. In Sec. II, the theory and the
computational details are given. In Sec. III, the results are presented and
discussed, and in Sec. IV a summary is given.

##  II Theory and computational details

###  II.1 Correlation estimators

####  II.1.1 $\tilde{g}(\mathbf{r})$

Our first correlation estimator is based on

$g(\mathbf{r})=\frac{|\nabla\rho(\mathbf{r})|}{\rho(\mathbf{r})},$ (2)

where $\rho$ is the electron density and $\nabla\rho$ the first derivative. In
Ref. [29], it was proposed to use $g$ to calculate the screening parameter
$\omega$ ($\omega\propto g$), which determines the separation of short- and
longe-range exchange in the screened hybrid functional HSE [7]. Then, the
average of $g$ in the whole unit cell of volume $V_{\text{cell}}$,

$\bar{g}=\frac{1}{V_{\text{cell}}}\int\limits_{\text{cell}}g(\mathbf{r}^{\prime})d^{3}r^{\prime},$
(3)

was used to define the parameter $c$ that specifies the relative weights of
the two terms in the mBJ potential [13]. $c$ was parameterized as
$c=\alpha+\beta\sqrt{\bar{g}}$, where $\alpha$ and $\beta$ are constants. In
Ref. [30], we discussed about the ability of $\bar{g}$ to distinguish between
strongly correlated TMO and itinerant elemental TM by having clearly different
values. However, since $\bar{g}$ is a constant for a particular system it can
not be used as a local probe and distinguish between different atoms in the
same system. Furthermore, on the technical side it is not applicable to non-
periodic solids, interfaces, and systems with vacuum (low-dimensional
systems), since in such systems averaging a quantity in the unit cell has no
meaning.

In the present work, we are searching for a quantity that is local (i.e.,
position dependent) and can provide an indication about the strength of
correlation on a particular atom. Using simply Eq. (2) would not really work
since some (local) average (as in $\bar{g}$) is still necessary in order to
have a function that is able as $\bar{g}$ to distinguish between TMO and pure
TM (see discussion in Sec. III.1). To this end we will use the smeared local
estimator first suggested in Ref. [31] and then implemented by Rauch et al.
[32, 33]. This local correlation estimator, $\tilde{g}$, is a local average of
$g$:

$\tilde{g}(\mathbf{r})=\frac{1}{(2\pi\sigma^{2})^{3/2}}\int
g(\mathbf{r}^{\prime})e^{-\frac{|\mathbf{r}-\mathbf{r}^{\prime}|^{2}}{2\sigma^{2}}}d^{3}r^{\prime},$
(4)

where the smearing parameter $\sigma$ determines the size of the region
(centered around $\mathbf{r}$) over which $g$ is averaged. The details of how
$\sigma$ can be determined and set as a parameter will be discussed in Sec.
III. The expression for $\tilde{g}$ is very advantageous since it can be very
easily calculated by using the convolution theorem if $g$ is expanded in plane
waves. Note that in periodic bulk systems, $\tilde{g}$ becomes a constant and
recovers the value of Eq. (3) when $\sigma$ is large enough.

However, for surfaces and other systems with vacuum numerical issues with
$|\nabla\rho|/\rho$ that becomes very large close to the surface region need
to be solved. We followed the prescription proposed in Ref. [32], which
consists of modifying Eq. (2) as follows:

$g(\mathbf{r})=\frac{1-\alpha}{\beta}\bigg{[}1-\text{erf}\bigg{(}\frac{\rho(\mathbf{r})}{\rho_{\text{th}}}\bigg{)}\bigg{]}+\frac{|\nabla\rho(\mathbf{r})|}{\rho(\mathbf{r})}\text{erf}\bigg{(}\frac{\rho(\mathbf{r})}{\rho_{\text{th}}}\bigg{)},$
(5)

where $\rho_{\text{th}}$ is a threshold for very low densities. For
$\rho\gg\rho_{\text{th}}$, $|\nabla\rho|/\rho$ is obtained, while for
$\rho\ll\rho_{\text{th}}$, $g$ becomes
$\left(1-\alpha\right)/\beta=\left(1-(-0.012)\right)/1.023=0.989$. Equation
(5) was proposed in the framework of the local mBJ potential to cope with the
aforementioned problem, but also to have the correct asymptotic behavior of
the local mBJ potential in the vacuum region. Here, the goal of the damping
with the function $\text{erf}(\rho/\rho_{\text{th}})$ is more to have a faster
convergence of the plane-wave expansion of $g$. Although the first term in Eq.
(5) is not necessary and could be discarded for the purpose of the present
work, we decided to use the full original expression from Ref. [32]. How
$\rho_{\text{th}}$ is chosen in the present work will be explained in Sec.
III.

####  II.1.2 $\tilde{s}(\mathbf{r})$

The second correlation estimator that we will consider is based on the reduced
density gradient $s$, which is used in the enhancement factor of exchange
functionals of the generalized gradient approximation (GGA) [34]. $s$ reads

$s(\mathbf{r})=\frac{|\nabla\rho(\mathbf{r})|}{2(3\pi^{2})^{1/3}\rho^{4/3}(\mathbf{r})}.$
(6)

Besides a constant factor, $s$ differs from Eq. (2) by the power in the
denominator that makes $s$ dimensionless. Far from the nuclei, $s$ goes to
infinity (while $g$ goes to a constant), but does not show the large values
close to surface regions like $g$. Similarly to the second term in Eq. (5) for
$g$, we will damp $s$ in the vacuum:

$s(\mathbf{r})=\frac{|\nabla\rho(\mathbf{r})|}{2(3\pi^{2})^{1/3}\rho^{4/3}(\mathbf{r})}\text{erf}\bigg{(}\frac{\rho(\mathbf{r})}{\rho_{\text{th}}}\bigg{)}$
(7)

and then use it in

$\tilde{s}(\mathbf{r})=\frac{1}{(2\pi\sigma^{2})^{3/2}}\int
s(\mathbf{r}^{\prime})e^{-\frac{|\mathbf{r}-\mathbf{r}^{\prime}|^{2}}{2\sigma^{2}}}d^{3}r^{\prime}$
(8)

to get our second correlation estimator.

###  II.2 Computational details

All calculations presented in this work were carried out using the all-
electron WIEN2k code [35, 36], which is based on the full potential linearized
augmented plane-wave and local orbitals [FP-(L)APW+lo] method [37, 38]. In the
FP-(L)APW+lo method, the wave functions, electron density, and potential are
expanded in spherical harmonics inside the atomic spheres and in plane waves
in the interstitial region. The GGA functional of Perdew, Burke, and Ernzerhof
(PBE) [34] has been used to treat the exchange-correlation effects in the
self-consistent calculations. The parameters of the calculations, like the
number of $\mathbf{k}$-points or the size of the basis set, have been chosen
such that the results are well converged. Typically, we used a basis-set size
corresponding to $R_{\text{MT}}^{\text{min}}K_{\text{max}}=9$, where
$R_{\text{MT}}^{\text{min}}$ is the smallest atomic sphere radius in the
system and $K_{\text{max}}$ is the magnitude of the largest reciprocal lattice
vector.

An important point concerning the calculation of the correlation estimators
$\tilde{g}$ and $\tilde{s}$ is the following. WIEN2k is an all-electron code,
therefore a plane-wave expansion in the entire unit cell of the wave functions
and all derived quantities like $\rho$ is in principle practically impossible.
However, $g$ and $s$, and more particularly the damped versions Eqs. (5) and
(7), are smooth enough such that it is possible to expand them in plane waves
in the entire unit cell with a Fourier series that is not prohibitively large.
Therefore, the calculation of the correlation estimators $\tilde{g}$ and
$\tilde{s}$ can be performed efficiently using the convolution theorem. If the
convolution theorem could not be used, the calculation of $\tilde{g}$ and
$\tilde{s}$ would be very cumbersome (see discussion in Ref. [39] for the
weighted-density approximation).

##  III Applications

###  III.1 Bulk solids

Figure 1: Box plot of the values of $\tilde{g}$ (upper panel) and $\tilde{s}$ (lower panel) on the $3d$ TM atom in correlated and non-correlated systems obtained with different values of the smearing parameter $\sigma$ (in bohr). The inner horizontal line in each box shows the median and the outer horizontal lines are the minimum and maximum values, while the outliers are shown with plus signs. Table 1: Values of $\tilde{g}$ and $\tilde{s}$ on the $3d$ TM atom in correlated and non-correlated bulk solids calculated with $\sigma=1.78$ bohr. The magnetic state is indicated in parenthesis. Non-correlated | $\tilde{g}$ | $\tilde{s}$ | Correlated | $\tilde{g}$ | $\tilde{s}$
---|---|---|---|---|---
Ti (NM) | 1.23 | 0.41 | TiO2-anatase (NM) | 2.02 | 0.81
V (NM) | 1.27 | 0.39 | TiO2-rutile (NM) | 1.98 | 0.76
Cr (AFM) | 1.35 | 0.40 | Ti2O3 (NM) | 1.90 | 0.71
Mn (NM) | 1.40 | 0.47 | V2O3 (AFM) | 1.93 | 0.70
Fe (FM) | 1.43 | 0.48 | SrVO3 (NM) | 1.90 | 0.68
FeAl (FM) | 1.23 | 0.38 | Cr2O3 (AFM) | 1.97 | 0.75
FeNi (FM) | 1.50 | 0.49 | CrO2 (FM) | 2.00 | 0.72
Fe3Ni (FM) | 1.48 | 0.49 | MnO (AFM) | 1.85 | 0.66
Fe2P (FM) | 1.42 | 0.46 | MnO2 (AFM) | 2.02 | 0.75
FeSb2(FM) | 1.50 | 0.52 | Mn2O3 (AFM) | 1.99 | 0.72
Co (FM) | 1.51 | 0.49 | FeO (AFM) | 1.83 | 0.65
Ni (FM) | 1.57 | 0.49 | Fe2O3 (AFM) | 1.96 | 0.73
Cu (NM) | 1.58 | 0.49 | Fe3O4 (FM) | 1.92 | 0.70
Cu2Sb (NM) | 1.52 | 0.49 | FeF2 (AFM) | 2.18 | 0.91
Cu3P (NM) | 1.56 | 0.54 | CoO (AFM) | 1.86 | 0.66
CuAu (NM) | 1.61 | 0.52 | NiO (AFM) | 1.90 | 0.66
Cu3Au (NM) | 1.62 | 0.54 | CuO (AFM) | 1.91 | 0.77
Zn (NM) | 1.47 | 0.49 | Cu2O (NM) | 1.79 | 0.68
|  |  | CuI (NM) | 1.79 | 0.79
|  |  | ZnO (NM) | 1.97 | 0.80
|  |  | YBa2Cu3O6 (FM) - Cu1,Cu2 | 1.81,1.89 | 0.78,0.70
|  |  | YBa2Cu3O7 (NM) - Cu1,Cu2 | 1.81,1.89 | 0.77,0.70

In order to find a relation between the degree of on-site correlation and
$\tilde{g}$ or $\tilde{s}$ on a $3d$ TM atom, we calculated these two
correlation estimators for a series of bulk solids with a wide range of
correlation strength: from pure TM (itinerant non- or weakly correlated) to
AFM oxides (strongly correlated). All considered systems are listed in Table 1
along with their magnetic state: non-magnetic (NM), FM, or AFM. For the
non-/weakly correlated systems, not only pure TM are considered, but also
other systems such as FeAl or Cu3P. For the more correlated systems, AFM and
NM oxides are considered, including YBa2Cu3O6 and its parent compound
YBa2Cu3O7, which is a high-$T_{c}$ superconductor.

Figure 1 shows the calculated $\tilde{g}$ and $\tilde{s}$ on the TM atoms for
different smearing parameters $\sigma$ in Eqs. (4) and (8), respectively. The
largest considered value of $\sigma$ is 3.78 bohr (2 Å), as originally chosen
in Ref. [32] for the local mBJ potential, so that $\tilde{g}$ and $\tilde{s}$
are averaged over a region that covers typical interatomic distances. Then we
gradually reduced $\sigma$ until 0.78 bohr, which is more representative of
the size of an atom. Since the results for $\sigma=3.78$ and 2.78 bohr lead to
the same results for the bulk systems, we do not show the values for
$\sigma=3.78$ bohr. The main observations that can be made are the following.
The most interesting one is that except for the smallest values of $\sigma$
(0.78 and 1.03 bohr), the values of the correlation estimators in correlated
and non-correlated systems are clearly different. For a value of $\sigma$ that
is too small, the average is done only around the atomic core region and the
resulting value becomes system-independent and no distinction between
correlated and non-correlated systems can be made. We can also observe that
the range of values of $\tilde{g}$ is for the correlated systems quite narrow.
Also, $\tilde{g}$ shows a pronounced variation with respect to $\sigma$ and
increases when $\sigma$ decreases. On the other hand, for the correlated
systems $\tilde{s}$ has a much larger spread than $\tilde{g}$ and is quite
independent of $\sigma$. In any case, both $\tilde{g}$ and $\tilde{s}$ seem to
be good candidates to be used as correlation estimators. In the rest of this
section, the discussion will be mainly based on the results obtained with
$\sigma=1.78$ bohr.

The numerical values of $\tilde{g}$ and $\tilde{s}$ obtained with
$\sigma=1.78$ bohr are shown in Table 1. For the non-correlated systems the
values of $\tilde{g}$ range from 1.2 to 1.6, while they are clearly larger for
the correlated systems, between 1.8 and 2.2. For $\tilde{s}$, the range is
0.4$-$0.55 for the non-correlated systems and 0.65$-$0.9 for the correlated
solids. Thus, for instance $\tilde{g}=1.7$ and $\tilde{s}=0.6$ for
$\sigma=1.78$ could be used as boundary values to distinguish between
correlated and non-correlated systems, and to decide if, at least
qualitatively, a particular TM atom would need or not a Hubbard on-site
correction $U$ in a DFT+$U$ calculation. However, it seems to be difficult to
estimate a value of $U$ from the estimators, because we do not really see any
systematic trend neither across the $3d$ series nor within a certain class of
compounds.

Among the correlated solids, we note that FeF2 and CuI are not oxides. While
the former is a typical correlated AFM solid, the latter, which has the zinc
blend structure, is a highly mobile $p$-type wide band gap nonmagnetic
semiconductor (similar as ZnO or TiO2). Considering it as correlated is
somehow consistent with the results from Refs. [40, 41], where it is shown
that the mBJ potential alone is not accurate enough to yield the experimental
band gap, so that adding an effective Hubbard term $U$ is necessary.

Table 2: Values of $\tilde{g}$ and $\tilde{s}$ on the Ni atom in FM fcc Ni and AFM NiO calculated with different smearing parameters $\sigma$ (in bohr). $\sigma$ | 2.78 | 1.78 | 1.53 | 1.28 | 1.03 | 0.78
---|---|---|---|---|---|---
| $\tilde{g}$
Ni | 1.52 | 1.57 | 1.67 | 1.95 | 2.50 | 3.37
NiO | 1.81 | 1.90 | 2.03 | 2.29 | 2.77 | 3.51
| $\tilde{s}$
Ni | 0.48 | 0.49 | 0.50 | 0.54 | 0.60 | 0.65
NiO | 0.65 | 0.66 | 0.67 | 0.70 | 0.72 | 0.71
Figure 2: Two-dimensional plots of $g$ in FM fcc Ni (left) and AFM NiO (right)
in a (001) plane. A Ni atom is always at the center. The scale of the ruler is
in bohr and on top Gaussian functions for $\sigma=0.78$ (blue) and 1.78 (red)
bohr are shown. Figure 3: Two-dimensional plots of $s$ in FM fcc Ni (left) and
AFM NiO (right) in a (001) plane. A Ni atom is always at the center. The scale
of the ruler is in bohr and on top Gaussian functions for $\sigma=0.78$ (blue)
and 1.78 (red) bohr are shown.

In order to discuss about $\tilde{g}$ and $\tilde{s}$ in more detail, and to
show the effect of $\sigma$, we take a closer look at two examples: FM face-
centered cubic (fcc) Ni and AFM NiO (in rocksalt structure) as non-correlated
and correlated systems, respectively. Figures 2 and 3 show two-dimensional
plots of $g$ and $s$ for Ni and NiO within a (001) plane with a Ni atom at the
center of the plot. Gaussian functions with $\sigma=0.78$ and 1.78 bohr are
added to the figures to show which area is covered by the integration in Eqs.
(4) and (8). With $\sigma=0.78$ bohr the integration is done only in the
region close to the Ni atom, thus giving values of $\tilde{g}$ and $\tilde{s}$
in Ni and NiO that are rather close (see Table 2). With $\sigma=1.78$ bohr the
low $s$ and $g$ area is also covered, which leads to more different values for
$\tilde{g}$ and $\tilde{s}$ in non-correlated and correlated systems, as
discussed above. Table 2 shows the values of $\tilde{g}$ and $\tilde{s}$ in Ni
and NiO for other values of $\sigma$. In order to sufficiently take into
account the environment of the Ni atom (i.e., for NiO the effect due to the
oxygen atoms) for calculating the average, $\sigma$ has to be large enough,
let us say at least 1.2 bohr.

###  III.2 Interfaces and surfaces

Moving to more complex systems, the values of $\tilde{g}$ and $\tilde{s}$ on
TM atoms in interface and surface systems are discussed below.

####  III.2.1 Interfaces

We calculated $\tilde{g}$ and $\tilde{s}$ on Mn atoms in an interface of AFM
body-centered cubic (bcc) Mn with MnO2 in rutile structure. The (001)
interface model consists of 5 Mn and 7 MnO2 layers and has only a small
$3.5\%$ lattice mismatch. As shown in Fig. 4, the interface model consists of
Mn atoms in four different situations; one in the middle of the Mn slab
without oxygen neighbors, one in the middle of the MnO2 slab with six oxygen
neighbors, and two in the interface region with two and four oxygen neighbors,
respectively. Table 3 shows $\tilde{g}$ and $\tilde{s}$ on the Mn atoms
obtained with different values of the smearing parameter $\sigma$. There is a
clear relation between the correlation estimators and the oxygen coordination,
except for small values of $\sigma$ (0.78 bohr for $\tilde{g}$ and up to 1.28
bohr for $\tilde{s}$). For comparison, $\tilde{g}$ and $\tilde{s}$ obtained
from bulk calculations of Mn and MnO2 (both AFM) are also shown in Table 3. In
most cases, the bulk values agree very well with the values on the
corresponding side of the interface system [Mn(no O neighbor) and Mn(6 O
neighbors), respectively]. However, some differences can also be noted like in
the case of $\tilde{g}$ with $\sigma=3.78$ bohr for Mn or $\tilde{s}$ with the
small values of $\sigma$ for MnO2. Figure 5 shows the perfectly linear
relation between our correlation estimators and the oxygen coordination number
for $\sigma=1.78$ bohr. Thus with our correlation thresholds of 1.7 for
$\tilde{g}$ or 0.6 for $\tilde{s}$ we can classify the Mn atoms with two
oxygen neighbors still as non-correlated, while those with four oxygen
neighbors should already be considered as correlated. As an alternative
interpretation, one could also use $\tilde{g}$ or $\tilde{s}$ to interpolate
the value of $U$ between 0 and 5 eV (assuming a $U$ of 5 eV for MnO2). In
Figure 5 we also show the magnetic moment of the Mn atoms. Metallic Mn has a
larger moment than MnO2, but the non-monotonic behavior indicates that the
magnetic moment cannot be used as a correlation estimator.

Figure 4: Side view of the Mn/MnO2 interface system. The color coding of the atoms is as follows: red (O), dark blue (Mn with no O neighbors), light green (Mn with two O neighbors), dark green (Mn with four O neighbors), light blue (Mn with six O neighbors). Table 3: Values of $\tilde{g}$ and $\tilde{s}$ on Mn atoms in the Mn/MnO2 interface system calculated with different values of $\sigma$ (in bohr). The values in bold for $\sigma$=1.78 indicate the correlated atoms. $\tilde{g}$ and $\tilde{s}$ in bulk Mn and MnO2 (both AFM) are also shown for comparison. $\sigma$ | 2.78 | 1.78 | 1.53 | 1.28 | 1.03 | 0.78
---|---|---|---|---|---|---
| $\tilde{g}$
Mn(bulk) | 1.30 | 1.40 | 1.57 | 1.91 | 2.49 | 3.28
Mn(no O neighbor) | 1.36 | 1.43 | 1.58 | 1.89 | 2.46 | 3.25
Mn(2 O neighbors) | 1.53 | 1.61 | 1.76 | 2.05 | 2.56 | 3.29
Mn(4 O neighbors) | 1.72 | 1.84 | 1.96 | 2.19 | 2.60 | 3.25
Mn(6 O neighbors) | 1.88 | 2.01 | 2.12 | 2.32 | 2.67 | 3.27
MnO2(bulk) | 1.89 | 2.02 | 2.13 | 2.33 | 2.69 | 3.29
| $\tilde{s}$
Mn(bulk) | 0.44 | 0.47 | 0.50 | 0.57 | 0.68 | 0.75
Mn(no O neighbor) | 0.45 | 0.46 | 0.50 | 0.56 | 0.66 | 0.73
Mn(2 O neighbors) | 0.55 | 0.56 | 0.59 | 0.64 | 0.72 | 0.76
Mn(4 O neighbors) | 0.66 | 0.64 | 0.63 | 0.63 | 0.63 | 0.60
Mn(6 O neighbors) | 0.75 | 0.73 | 0.71 | 0.69 | 0.66 | 0.62
MnO2(bulk) | 0.77 | 0.75 | 0.75 | 0.74 | 0.74 | 0.72
Figure 5: $\tilde{g}$, $\tilde{s}$, and magnetic moment $M$ ($\mu_{B}$) on Mn
atoms in the Mn/MnO2 interface system calculated with $\sigma=1.78$ bohr
plotted against the number of oxygen atoms surrounding a Mn atom.

####  III.2.2 Surfaces

As mentioned in Sec. II.1 and in Ref. [32], $g$ as given by Eq. (2) is quite
large in the vacuum region, while $s$ [Eq. (6)] even goes to infinity far from
the nuclei. Therefore, they were damped by multiplying them by
$\text{erf}(\rho/\rho_{\text{th}})$ [32], see Eqs. (5) and (7). The threshold
density $\rho_{\text{th}}$ is a parameter that has to be chosen suitably,
i.e., not too small so that one has damping of $g$ and $s$ in the vacuum, but
also not too large so that they are not affected inside the bulk. In fact for
our purpose, $\rho_{\text{th}}$ should be chosen such that the values of
$\tilde{g}$ and $\tilde{s}$ on the surface atoms are the same/similar as the
values on atoms deeper in the bulk. Taking the Fe-(001) surface as example,
the effect of the damping on $g$ and $s$ is displayed in Figs. 6 and 7,
respectively. The figures show the electron density $\rho$ of sub-subsurface
and surface Fe atoms and the corresponding $g$ or $s$, undamped and damped
with $\rho_{\text{th}}$ = 0.002, 0.01, or 0.015 e/bohr3.

Figure 6: Electron density $\rho$, original $g$ [Eq. (2)], and $g$ damped with
different threshold densities $\rho_{\text{th}}$ [Eq. (5)] plotted for the
Fe(001) surface along the [001] direction into the vacuum. Fesss and Fes
indicate the sub-subsurface and surface Fe atoms, respectively. Figure 7:
Electron density $\rho$, original $s$ [Eq. (6)], and $s$ damped with different
threshold densities $\rho_{\text{th}}$ [Eq. (7)] plotted for the Fe(001)
surface along the [001] direction into the vacuum. Fesss and Fes indicate the
sub-subsurface and surface Fe atoms, respectively.

We calculated $\tilde{g}$ and $\tilde{s}$ on the Ni and Fe atoms in surface
systems with and without oxygen coverage. Different values of
$\rho_{\text{th}}$ and $\sigma$ were considered, but we show in Table 4 only
the results for $\tilde{g}$ and $\tilde{s}$ obtained with
$\rho_{\text{th}}=0.015$ e/bohr3 and $\sigma=2.78$ bohr. With this choice the
correlation estimators for the plain surfaces are only very little enhanced as
compared to the bulk, and the method is optimally sensitive to distinguish
correlated and non-correlated atoms.

We first discuss Ni(111) surfaces. In the case of (full) oxygen coverage, the
oxygen atom is located at the fcc hollow site. From Table 4 we can see that
for the plain Ni(111) surface the correlation estimators are very close to the
bulk values. However, with full oxygen coverage the Ni atoms at the surface
have values of $\tilde{g}$ and $\tilde{s}$ that are clearly larger than for
the Ni atoms in the subsurface and deeper into the bulk. By comparing
$\tilde{g}$ and $\tilde{s}$ with the corresponding values of bulk Ni and NiO
we clearly can classify the surface Ni atom in the fully O covered surface as
correlated.

For Fe(001) surfaces, we studied plain Fe(001) and Fe(001) with different
oxygen coverages, namely full, full with one additional oxygen atom in the
subsurface (in a $2\times 2$ supercell), half, and 1/9 ($3\times 3$
supercell). In these cases the surface Fe atoms have 4, 4 or 5, 2 and 1 oxygen
neighbors, respectively. The results in Table 4 show that for both full oxygen
coverages all atoms on the surface have the largest values of $\tilde{g}$ and
$\tilde{s}$. In the case of half coverage $\tilde{g}$ and $\tilde{s}$ on the
Fe atoms at the surface are smaller, whereas a further reduction is obtained
for the Fe atoms in the subsurfaces (except in the case of full coverage(2))
and for all Fe atoms in the case of 1/9 coverage. In the subsurface of full
coverage(2) each Fe atom has two oxygen nearest neighbors and is therefore
expected to be more correlated. For comparison we also show the values for
bulk Fe and FeO. For convenience we also show the results for the surface and
subsurface Fe atoms graphically in Fig. 8 and identify the surface atoms in
the full coverage case as correlated. Whether the Fe atoms on the surface with
half oxygen coverage or in the subsurface with full oxygen coverage should be
considered as correlated or not depends on the chosen boundary value (e.g.,
see dashed line in Fig. 8 as a possible choice) or one could use our
estimators to interpolate the value of $U$ for the different Fe atoms.

Table 4: $\tilde{g}$ and $\tilde{s}$ on TM atoms in plain Ni(111) and Fe(001) surfaces, as well as in Ni(111) with full oxygen coverage and Fe(001) with different oxygen coverage. The values in bold indicate the correlated atoms. $\sigma=2.78$ bohr and $\rho_{\text{th}}=0.015$ e/bohr3 were used in the calculations. The results for bulk Ni (FM), NiO (AFM), Fe (FM), and FeO (AFM) are also shown for comparison. Fe-full coverage(2) is a $2\times 2$ Fe surface with full oxygen coverage and one extra oxygen in the subsurface and (*) indicates the Fe atom on the surface which also has one oxygen neighbor in the subsurface. System | $\tilde{g}$ | $\tilde{s}$
---|---|---
Ni (Ni@surface) | 1.52 | 0.51
Ni (Ni@subsurface) | 1.54 | 0.49
Ni (Ni@middle) | 1.52 | 0.48
Ni-full coverage (Ni@surface) | 1.71 | 0.56
Ni-full coverage (Ni@subsurface) | 1.56 | 0.49
Ni-full coverage (Ni@middle) | 1.54 | 0.49
Ni(bulk) | 1.52 | 0.48
NiO(bulk) | 1.81 | 0.65
Fe (Fe@surface) | 1.40 | 0.47
Fe (Fe@subsurface) | 1.42 | 0.49
Fe (Fe@middle) | 1.38 | 0.47
Fe-full coverage (Fe@surface) | 1.60 | 0.55
Fe-full coverage (Fe@subsurface) | 1.48 | 051
Fe-full coverage (Fe@middle) | 1.39 | 0.47
Fe-full coverage(2) (Fe@surface*) | 1.75 | 0.56
Fe-full coverage(2) (Fe@surface) | 1.61 | 0.56
Fe-full coverage(2) (Fe@subsurface) | 1.52 | 0.53
Fe-full coverage(2) (Fe@middle) | 1.41 | 0.49
Fe-half coverage (Fe@surface) | 1.55 | 0.52
Fe-half coverage (Fe@subsurface(near O)) | 1.47 | 0.50
Fe-half coverage (Fe@subsurface(far O)) | 1.46 | 0.50
Fe-half coverage (Fe@middle) | 1.39 | 0.47
Fe-1/9 coverage (Fe@surface(near O)) | 1.45 | 0.49
Fe-1/9 coverage (Fe@surface (far O)) | 1.40 | 0.47
Fe-1/9 coverage (Fe@subsurface(near O)) | 1.46 | 0.50
Fe-1/9 coverage (Fe@subsurface(far O)) | 1.42 | 0.49
Fe-1/9 coverage (Fe@middle) | 1.38 | 0.47
Fe(bulk) | 1.37 | 0.47
FeO(bulk) | 1.73 | 0.64
Figure 8: $\tilde{g}$ and $\tilde{s}$ for different Fe atoms on the surface
and subsurface of Fe(001) surfaces with different O coverage calculated with
$\sigma=2.78$ bohr. The down triangle are the results for the Fe(001) surface
with full oxygen coverage and one extra oxygen in the subsurface. The dashed
line shows our selected value to separate correlated and non-correlated atoms.
The green triangles are $\tilde{g}$ and $\tilde{s}$ for bulk Fe and bulk FeO.

We also calculated $\tilde{g}$ and $\tilde{s}$ on the Cu atoms in the more
complex system shown in Fig. 9. It consists of a nonmagnetic Cu5O cluster
adsorbed on the anatase TiO2(101) ($3\times 2$) surface. This and similar
systems studied in Ref. [42] constitute very irregular surfaces with Cu atoms
ranging from neutral to +2-charged. The results obtained with $\sigma=2.78$
bohr are shown in Table 5. As discussed in Sec. III.1, the correlation
estimators are mainly determined by the environment of the TM atom and the
results show that the correlated Cu atoms (Cu1$-$Cu4) are those which are
bonded to O atoms of the surface or cluster and have values above 1.7 for
$\tilde{g}$ and at $\sim 0.6$ for $\tilde{s}$. They are thus larger than for
the non-correlated Cu5 atom ($\tilde{g}=1.59$ and $\tilde{s}=0.48$) that is
bonded only to another Cu atom. The results for the bulk solids Cu (NM), CuO
(AFM), and Cu2O (NM) are also shown in Table 5, and we can see that the values
for bulk Cu are similar to the values for the Cu5 atom. The values of
$\tilde{g}$ for CuO and Cu2O are larger and smaller than for Cu1$-$Cu4,
respectively. For $\tilde{s}$, the values are larger in both CuO and Cu2O
compared to Cu1$-$Cu4.

Obviously, also the oxidation state depends on the environment and it is
interesting to calculate the Bader charges [43] of the Cu atoms as an
additional quantity to compare with. The calculation of the Bader charge is
based on the gradient vector field of $\rho$, and the surfaces of zero flux
define the atoms and therefore their charge (nuclear charge minus integrated
$\rho$). As shown in Table 5 the Bader charge is almost zero for the non-
correlated Cu5, while it amounts to 0.60 for Cu1 and Cu2 and 0.41 for Cu3 and
Cu4, corresponding to Cu+ ions, which were estimated as correlated atoms.

Table 5: $\tilde{g}$, $\tilde{s}$, and Bader charge of the five Cu atoms of the Cu5O cluster adsorbed on the anatase TiO2(101) surface. The values in bold indicate the correlated atoms. $\sigma$ = 2.78 bohr and $\rho_{\text{th}}=0.015$ e/bohr3 were used in the calculations. The results for Cu in bulk Cu (NM), CuO (AFM), and Cu2O (NM) are also shown for comparison. System | $\tilde{g}$ | $\tilde{s}$ | Bader charge
---|---|---|---
TiO2-Cu5O(Cu1) | 1.74 | 0.58 | 0.60
TiO2-Cu5O(Cu2) | 1.74 | 0.58 | 0.60
TiO2-Cu5O(Cu3) | 1.72 | 0.59 | 0.41
TiO2-Cu5O(Cu4) | 1.72 | 0.59 | 0.41
TiO2-Cu5O(Cu5) | 1.59 | 0.48 | 0.07
Cu(bulk) | 1.51 | 0.48 | 0
CuO(bulk) | 1.80 | 0.77 | 1.07
Cu2O(bulk) | 1.64 | 0.67 | 0.55
Figure 9: Side view of a Cu5O cluster adsorbed on the anatase TiO2(101)
surface. The color coding of the atoms is as follows: red (surface O), yellow
(O bonded to Cu), blue (Ti), brown (correlated Cu), orange (non-correlated
Cu).

Thus, overall there are very clear trends in the values of the correlation
estimators on the TM atoms also in complicated systems. Thus, this
demonstrates that $\tilde{g}$ or $\tilde{s}$ can be efficiently used to
estimate the correlation on TM atoms in complex systems like interfaces or
surfaces.

##  IV Summary

In this work, we have shown that the value of the correlation estimators
$\tilde{g}$ and $\tilde{s}$ at the nucleus of a TM atom, which are local
averages of density-dependent quantities around the corresponding atom, can be
used to estimate the strength of correlation of the TM atom. In bulk solids,
where we usually know from experience in which systems the TM atoms are
correlated, there is a very clear difference in the values of the correlation
estimators between correlated (e.g., in oxides) and non-correlated (e.g., in
pure metals) TM atoms. In more complicated systems, like at interfaces or
surfaces, it may be unclear whether a certain TM atom is correlated or not,
however we showed that our correlation estimators are very reliable in
providing a very good hint on the correlation strength. We have demonstrated
the power of the estimators for oxygen-covered TM surfaces, for Cu5O clusters
adsorbed on the TiO2-anatase surface, and a Mn/MnO2 interface. Thus,
$\tilde{g}$ or $\tilde{s}$ could be used to determine for which atoms a
Hubbard $U$ correction should be applied in a DFT+$U$ calculation. According
to the results shown in this work, we would favour $\tilde{g}$ as a more
reliable estimator. However, it does not seem possible to go to a more
quantitative level and to find a relation between the estimators and a
specific value of $U$ in general, although it might be possible to use
$\tilde{g}$ or $\tilde{s}$ in systems having several TM atoms of the same type
but in different environments as demonstrated for the Mn/MnO2 interface or the
Fe and Ni surfaces.

###### Acknowledgements.

L.K. and P.B. acknowledge support by the TU-D doctoral college (TU Wien).

## References

  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Terakura  _et al._ [1984] K. Terakura, T. Oguchi, A. R. Williams, and J. Kübler, Phys. Rev. B 30, 4734 (1984).
  * Bredow and Gerson [2000] T. Bredow and A. R. Gerson, Phys. Rev. B 61, 5194 (2000).
  * Perry  _et al._ [2001] J. K. Perry, J. Tahir-Kheli, and W. A. Goddard, III, Phys. Rev. B 63, 144510 (2001).
  * Muscat  _et al._ [2001] J. Muscat, A. Wander, and N. M. Harrison, Chem. Phys. Lett. 342, 397 (2001).
  * Heyd  _et al._ [2003] J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118, 8207 (2003), 124, 219906 (2006).
  * Metzner and Vollhardt [1989] W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989), 62, 1066 (1989).
  * Kotliar  _et al._ [2006] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
  * Held [2007] K. Held, Adv. Phys. 56, 829 (2007).
  * Sharma  _et al._ [2013] S. Sharma, J. K. Dewhurst, S. Shallcross, and E. K. U. Gross, Phys. Rev. Lett. 110, 116403 (2013).
  * Di Sabatino  _et al._ [2015] S. Di Sabatino, J. A. Berger, L. Reining, and P. Romaniello, J. Chem. Pyhs. 143, 024108 (2015).
  * Tran and Blaha [2009] F. Tran and P. Blaha, Phys. Rev. Lett. 102, 226401 (2009).
  * Tran and Blaha [2017] F. Tran and P. Blaha, J. Phys. Chem. A 121, 3318 (2017).
  * Anisimov  _et al._ [1993] V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyżyk, and G. A. Sawatzky, Phys. Rev. B 48, 16929 (1993).
  * Ylvisaker  _et al._ [2009] E. R. Ylvisaker, W. E. Pickett, and K. Koepernik, Phys. Rev. B 79, 035103 (2009).
  * Himmetoglu  _et al._ [2014] B. Himmetoglu, A. Floris, S. de Gironcoli, and M. Cococcioni, Int. J. Quantum Chem. 114, 14 (2014).
  * Dederichs  _et al._ [1984] P. H. Dederichs, S. Blügel, R. Zeller, and H. Akai, Phys. Rev. Lett. 53, 2512 (1984).
  * Hybertsen  _et al._ [1989] M. S. Hybertsen, M. Schlüter, and N. E. Christensen, Phys. Rev. B 39, 9028 (1989).
  * Madsen and Novák [2005] G. K. H. Madsen and P. Novák, Europhys. Lett. 69, 777 (2005).
  * Vaguier  _et al._ [2012] L. Vaguier, H. Jiang, and S. Biermann, Phys. Rev. B. 86, 165105 (2012).
  * Springer and Aryasetiawan [1998] M. Springer and F. Aryasetiawan, Phys. Rev. B 57, 4364 (1998).
  * Aryasetiawan  _et al._ [2004] F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann, and A. I. Lichtenstein, Phys. Rev. B 70, 195104 (2004).
  * Aryasetiawan  _et al._ [2006] F. Aryasetiawan, K. Karlsson, O. Jepsen, and U. Schönberger, Phys. Rev. B 74, 125106 (2006).
  * Şaşıoğlu  _et al._ [2012] E. Şaşıoğlu, C. Friedrich, and S. Blügel, Phys. Rev. Lett. 109, 146401 (2012).
  * Cococcioni and de Gironcoli [2005] M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71, 035105 (2005).
  * Pickett  _et al._ [1998] W. E. Pickett, S. C. Erwin, and E. C. Ethridge, Phys. Rev. B 58, 1201 (1998).
  * Wang and Jiang [2019] Y.-C. Wang and H. Jiang, J. Chem. Phys. 150, 154116 (2019).
  * Krukau  _et al._ [2006] A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, J. Chem. Phys. 125, 224106 (2006).
  * Tran  _et al._ [2019] F. Tran, J. Doumont, L. Kalantari, A. W. Huran, M. A. L. Marques, and P. Blaha, J. Appl. Phys. 126, 110902 (2019).
  * Marques  _et al._ [2011] M. A. L. Marques, J. Vidal, M. J. T. Oliveira, L. Reining, and S. Botti, Phys. Rev. B 83, 035119 (2011).
  * Rauch  _et al._ [2020a] T. Rauch, M. A. L. Marques, and S. Botti, J. Chem. Theory Comput. 16, 2654 (2020a).
  * Rauch  _et al._ [2020b] T. Rauch, M. A. L. Marques, and S. Botti, Phys. Rev. B 101, 245163 (2020b), 102, 119902(E) (2020).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996), 78, 1396(E) (1997).
  * Blaha  _et al._ [2018] P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, J. Luitz, R. Laskowski, F. Tran, and L. D. Marks,  _WIEN2k: An Augmented Plane Wave plus Local Orbitals Program for Calculating Crystal Properties_ (Vienna University of Technology, Austria, 2018).
  * Blaha  _et al._ [2020] P. Blaha, K. Schwarz, F. Tran, R. Laskowski, G. K. H. Madsen, and L. D. Marks, J. Chem. Phys. 152, 074101 (2020).
  * Singh and Nordström [2006] D. J. Singh and L. Nordström,  _Planewaves, Pseudopotentials and the LAPW Method, 2nd ed._ (Springer, Berlin, 2006).
  * Karsai  _et al._ [2017] F. Karsai, F. Tran, and P. Blaha, Comput. Phys. Commun. 220, 230 (2017).
  * Singh [1993] D. J. Singh, Phys. Rev. B 48, 14099 (1993).
  * Jiang [2013] H. Jiang, J. Chem. Phys. 138, 134115 (2013).
  * Rubel  _et al._ [2021] O. Rubel, F. Tran, X. Rocquefelte, and P. Blaha, Comput. Phys. Commun. 261, 107648 (2021).
  * [42] J. S. Schubert, L. Kalantari, A. Lechner, A. Giesriegl, P. S. Nandan, P. A. Leiva, S. Kashiwaya, M. Sauer, A. Foelske, J. Rosán, P. Blaha, A. Cherevan, and D. Eder, Unpublished.
  * Bader [1985] R. F. W. Bader, Acc. Chem. Res. 18, 9 (1985).