Dear Krystian,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?
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?
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:
SPIN 1 eigenvalues: 0.471 0.511 0.928 0.928 0.957 SPIN 2 eigenvalues: 0.039 0.407 0.459 0.905 0.905As 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.
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).
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 *71*, 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.
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.
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).
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.
All the best for you and your work, Eric -- 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 http://www.hmi.uni-bremen.de/
smime.p7s
Description: S/MIME Cryptographic Signature
_______________________________________________________________________________ The Quantum ESPRESSO Foundation stands in solidarity with all civilians worldwide who are victims of terrorism, military aggression, and indiscriminate warfare. -------------------------------------------------------------------------------- Quantum ESPRESSO is supported by MaX (www.max-centre.eu) users mailing list [email protected] https://lists.quantum-espresso.org/mailman/listinfo/users
