Large difference between ΔN_NSCF and ΔN_SCF in Cococcioni linear-response U calculation

[tasdiquearman]

Dear VASP developers and users,

I am performing a Cococcioni linear-response calculation of the Hubbard U parameter for an Fe site in a polyanionic Na–Fe–Mn–phosphate material. I follow the procedure in the VASP wiki (“Calculate U for LSDA+U”) and obtain well-converged SCF and NSCF runs for different perturbation amplitudes. My tabulated occupancies for the Fe 3d orbitals are:

Code: Select all

alpha (eV)    n_nscf     n_scf
-0.20         5.449      6.011
-0.15         5.607      6.025
-0.10         5.719      6.026
-0.05         5.880      6.039
+0.05         6.220      6.051
+0.10         6.381      6.060
+0.15         6.523      6.068
+0.20         6.640      6.073

The ground-state (unperturbed) d-occupancy is 6.049. From these data, I observe that, The NSCF response is large, e.g. ΔN_NSCF(0.10 eV) ≈ 6.381 − 6.049 ≈ 0.33, and the SCF (screened) response is almost an order of magnitude smaller, e.g. ΔN_SCF(0.10 eV) ≈ 6.060 − 6.049 ≈ 0.011. This difference is consistent across all perturbation strengths. When I fit N vs α to obtain χ₀ and χ, I find an effective U ≈ 6.2 eV, which is somewhat higher than reported values (~4–5 eV) for similar Fe-based compounds.
My questions are:
(a) Is it normal in VASP’s linear-response workflow for ΔN_NSCF to be much larger than ΔN_SCF (sometimes by a factor of 20–30)?
(b) Can a larger U (≈6 eV) be physically reasonable for Fe in a polyanionic framework, or does this point to a convergence/projection/mixing issue in the SCF+NSCF steps?
(c) Are there recommended convergence settings or projection parameters to ensure robust linear-response occupancies (e.g., LORBIT, LMAXMIX, ALGO, EDIFF)?
(d) Could such a large difference indicate that my perturbation amplitudes (±0.2 eV) are too large, i.e., the system is outside the linear regime?

Any guidance on whether my results indicate correct behavior or a setup issue would be greatly appreciated.
Thank you very much.