<!DOCTYPE html>
<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
</head>
<body text="#000000" bgcolor="#deddda">
<p>Dear Krystian,</p>
<p>Calculating PBAs with "standard" DFT+U formulations is a quite
tricky thing to do for fundamental reasons. While the method works
well for many other battery materials, it tends to fail for PBAs,
and LR-cDFT/DFPT indeed often yield inexplicably large Hubbard
parameters. So what is happening?</p>
<p>The LS sites (the transition metals sitting at the
C-terminations) of PBAs are subject to an extremely strong
(~octahedral) ligand field, and carbon itself is an element with a
much lower electronegativity than oxygen, which most of the DFT+U
community is used to seeing as a ligand. This results in a very
appreciable covalent character of the metal-ligand bond (Fe--C in
your case); so far so good. But what does this mean for DFT+U and
for the calculation of U parameters from perturbative methods?</p>
<p>1. While you didn't mention it, I strongly assume that you used
atomic or ortho-atomic Hubbard projectors for your calculations.
Now lets look at the eigenvalue spectrum of the d-orbital
occupation matrix in Mn[Fe(CN)6], as outputted by QE:</p>
<p> SPIN 1<br>
eigenvalues:<br>
0.471 0.511 0.928 0.928 0.957<br>
SPIN 2<br>
eigenvalues:<br>
0.039 0.407 0.459 0.905 0.905<br>
</p>
<p>As expected for a LS-Fe(III) species, we see that there are 5
almost full d-orbitals, namely those with occupation eigenvalues
>0.9. We also see that in the spin-down manifold, there is an
eigenvalue that is very close to 0, and this is the unoccupied
"t2g-like" orbital (careful that the point symmetry of Fe in
Mn[Fe(CN)6] is not really O_h, but we'll approximate it as such).
But then what are the other 4 orbitals with eigenvalues ~0.5?
Well, these should be the e_g orbitals, but as you can see, they
are half-occupied even though fundamental chemistry tells us they
should be empty. What's going on here is that the atomic-like
projectors we're using are simply very very bad approximations for
a state that is so hybridized. The e_g orbitals of Fe are forming
sigma bonds with the CN ligand MOs, and this hybridization leads
to an MO with a shape vastly different from that of a d-orbital in
a free atom calculation (for which the projector was generated!).
So it is the mismatch between and projector and the state we seek
to represent that explains these deceptively "fractional"
occupation values.</p>
<p>2. This has direct implications for SCF calculations with a
finite U value. If you apply, say, U(Fe-3d)=6eV to these states
using the Dudarev functional, then every single one of these
orbitals will give rise to a punitive Hubbard energy of
U/2*(0.5*(1-0.5) = 0.75eV, which is a huge amount of energy.
Consequently, the system will try to avoid these fractional
occupations at all cost, either by unphysically switching to the
HS state or by unphysically extending to bond lengths far beyond
their experimental values if you allow for relaxations (see Mariano
et al., J. Chem. Theory Comput. 2020, 16, 11, 6755–6762 for a
detailed analysis performed on Fe(II) hexacomplexes including
Fe[CNH]6).</p>
<p>3. The problem is aggravated when you start computing the U
values from first principles using DFPT or the equivalent LR-cDFT
approach (Cococcioni & de Gironcoli, Phys. Rev. B <b>71</b>,
035105). When you calculate the response of a manifold (such as
Fe-3d) by perturbing it to obtain U, you essentially calculate how
much the non-perturbed states are able to screen this
perturbation. If your perturbation is well-screened, you'll get
low U parameters, and vice versa. Perturbing the Fe-3d shell in
Mn[Fe(CN)6], you'd expect that the screening should be mainly due
to the C and N 2p orbitals. However, and here is the problem, if
you 3d projector already includes a large chunk of the C-2p (and
eventually even N-2p) orbitals, then these orbitals are INSIDE the
perturbed space, and therefore cannot screen the perturbation.
Thus, you obtain huge Hubbard U parameters that are only partially
screened when they should be fully screened.</p>
<p>So what can you do? One option is to change your Hubbard
projector basis to one that offers a better (more localized!)
representation of your hybridized states, such as Wannier
functions. This way you'd get rid of these deceptive fractional
occupations and also find lower U values. However, the drawback is
that the Wannierization can be quite cumbersome and then you
cannot compute forces and stresses anymore.</p>
<p>An easier ad-hoc fix is orbital-resolved DFT+U, which has been
published by us in J. Chem. Theory Comput. 2024, 20, 11,
4824–4843. Here the idea is to circumvent the problem of non-ideal
projector orbitals by simply excluding the hardly projectable
eigenstates from the Hubbard manifold. This leads to
orbital-resolved Hubbard parameters, such as U(Fe-t2g)=3.5eV
instead of U(Fe-3d)=8.0eV. The difficulty here is to find a
suitable Hubbard manifold, for which you need to invoke group
theory and also analyze the electronic structure (see Warda et
al., J. Chem. Theory Comput. 2026, 22, 2, 1016–1029 for an
example).</p>
<p>In any case, I would advise against using "standard" DFT+U
calculating PBAs, especially not in combination with U parameters
obtained from shell-averaged DFPT. For simple tests, it might be
best to simply use an empirical value of 2eV>U>5eV for most
d-shells, for anything more serious I'd use OR-DFT+U or DFT+U with
Wannier projectors. </p>
<p>All the best for you and your work,</p>
<p>Eric</p>
<pre class="moz-signature" cols="72">--
Eric Macke
PhD Student
Hybrid Materials Interfaces Group
University of Bremen
Faculty of Production Engineering
TAB-Building, Room 3.29
Am Fallturm 1
28359 Bremen, DE
<a class="moz-txt-link-freetext" href="http://www.hmi.uni-bremen.de/">http://www.hmi.uni-bremen.de/</a></pre>
</body>
</html>