Calculation of Hubbard U parameters for Fe in fayalite

[Heseri]

Dear VASP developers,

I am trying to calculate the Hubbard U parameters for Fe in fayalite using the linear response approach implemented in VASP. Fayalite is the iron-rich end member of (Mg,Fe)2SiO4 olivine, and was used as an example in the original paper by Cococcioni and Gironcoli (2005, Phys. Rev. B 71, 035105). I have followed the VASP tutorial “Calculate U for LSDA+U” (https://vasp.at/wiki/Calculate_U_for_LSDA%2BU), but I have two questions.

First, as described by Cococcioni and Gironcoli (2005), there are two distinct Fe sites in fayalite, and I need to calculate two site-specific U values. In the VASP tutorial on Ni in NiO, only one Ni atom/cation was specified for the linear response calculations. In this way, the POSCAR file will be like:

Code: Select all

Fe Fe Fe Si O
1 31 32 128

For my case, I am wondering whether, for each Fe site, it would be better to include all Fe atom/cation at that site in the linear response calculations, i.e., using the POSCAR file like

Code: Select all

Fe Fe Si O
32 32 32 128

then perform linear fits of the number of d-electrons as a function of the additional potential for each Fe atom at that site, and finally average the derived U values over the 32 Fe atoms?

Second, to be rigorous, if I use only one Fe atom at a given site in the calculations, which LDAUL value (2 or -1) should I assign to the other Fe atoms at the same site? For example, should I use:

Code: Select all

LDAU = .TRUE.
LDAUTYPE = 3
LDAUL = 2 -1 -1 -1 -1
LDAUU = $v 0.0 0.0 0.0 0.0
LDAUJ = $v 0.0 0.0 0.0 0.0

or

Code: Select all

LDAU = .TRUE.
LDAUTYPE = 3
LDAUL = 2 2 -1 -1 -1
LDAUU = $v 0.0 0.0 0.0 0.0
LDAUJ = $v 0.0 0.0 0.0 0.0

Or would these two settings produce in fact the same results?

Thank you very much for your time and consideration.
Sincerely,
Tian Hua

[martin.schlipf]

We checked both questions numerically on the NiO example from the tutorial: AFM-II, LSDA, α = ±0.1 and ±0.2 eV.

1. Perturb only one Fe atom per inequivalent site. If you apply the potential to all 32 equivalent Fe atoms at once, each atom's occupation also responds to the potential on its perturbed neighbours. You then measure a response of the whole sublattice rather than the on-site response to an isolated perturbation. That calculation is in fact identical to perturbing the corresponding atoms in the primitive cell. For NiO:

   Code: Select all

   | Setup                                          | χ₀    | χ     | U (eV) |
   |------------------------------------------------|-------|-------|--------|
   | primitive cell (4 atoms)                       | 0.432 | 0.120 | 6.02   |
   | 2×2×2 supercell, whole Ni sublattice perturbed | 0.432 | 0.120 | 6.02   |
   | 2×2×2 supercell, single Ni perturbed           | 0.532 | 0.130 | 5.81   |
   | 3×3×3 supercell, single Ni perturbed           | 0.537 | 0.130 | 5.83   |
   

   Perturbing the whole sublattice in the supercell exactly reproduces the primitive cell, so the supercell gains you nothing. Perturbing a single atom converges with supercell size. Averaging over the 32 atoms adds no information either, because they are symmetry-equivalent and respond identically. For fayalite, run two separate sets of calculations: one with a single M1 Fe perturbed and one with a single M2 Fe perturbed. Check convergence with respect to the supercell size.

2. Both LDAUL settings give the same result. With LDAUTYPE = 3, LDAUU/LDAUJ set the strength of the perturbing potential. Species with LDAUU = LDAUJ = 0 feel no additional potential, whether LDAUL = 2 or -1. In our test, the two settings gave total energies that agree to within 3·10⁻⁸ eV and identical d occupations on all atoms, in every run. With LORBIT = 11, the occupations of all atoms appear in the "total charge" table of the OUTCAR in either case. Use that table for the fits, because it also contains the response in the non-self-consistent (ICHARG = 11) runs.

Sign convention. A positive LDAUU with LDAUTYPE = 3 increases the d occupation in VASP. With χ₀ and χ taken as the fitted slopes, U = 1/χ − 1/χ₀, as in the wiki tutorial.

[Heseri]

Dear Dr. Schlipf,

Thank you very much for your detailed tests and helpful reply! Following your guidance, I continued calculating the U parameters for other Fe-bearing minerals, and encountered another related issue. Could you please offer some further guidance on the following question?

The issue arose when I calculated the U value for Fe2+ in Fe-bearing MgSiO3 perovskite. I first optimized the structures (ISIF=2, with fixed lattice parameters) of 1FeSiO3+15MgSiO3 and 1FeSiO3+31MgSiO3, and then calculated the U values using the linear response approach following the NiO tutorial. I obtained U values of 2.66 eV and 2.76 eV for the smaller and larger systems, respectively, suggesting that the finite-size effect is negligible.

Then I doubled the 1FeSiO3+15MgSiO3 structure to obtain a system consisting of 2FeSiO3+30MgSiO3, in which the two Fe2+ occupy equivalent sites. After optimizing this structure, I calculated the U values for Fe2+ in two ways: (1) perturbing only one Fe2+ and (2) perturbing both two Fe2+ simultaneously. However, these two calculations yielded U values of 5.88 eV and 2.62 eV, respectively: the value obtained by perturbing both two Fe2+ simultaneously is much closer to that calculated from the 1FeSiO3+15MgSiO3 system.

I am quite puzzled by this significant difference between the two values. Could you please help me understand why these two approaches yield such different results? Furthermore, which U value should be employed for Fe2+ (i.e., ~2.6 eV or ~5.9 eV) when performing simulations of systems such as 2FeSiO3+30MgSiO3 and even 4FeSiO3+28MgSiO3?

Thank you very much for your time and consideration. I would greatly appreciate any insights or suggestions you could share.

Sincerely,
Tian Hua

P.S. I used ALGO=All and reduced mixing parameters to obtain reasonable magnetic configurations: with these settings and an initial local magnetic moment of 4 μB for Fe2+, the Fe2+ ions remained in the high-spin state in all simulations. For your references, I have uploaded the POSCAR files for the four linear response calculations and the key files for the linear response calculations for 2FeSiO3+30MgSiO3, and would be happy to provide any additional details or files if needed.

[martin.schlipf]

We expect the calculation with 2 Fe sites and both being perturbed to reproduce the U you obtain with a single Fe site in half the cell. This is the same argument about using a supercell and perturbing every site. However, we still believe that a U of ~2.7eV is more realistic in your case. The likely issue is that the d charges move to the second unperturbed Fe site. You can look at this in the output of VASP. If you see that the second site has similar changes to its charge occupancies, the approximation that only the perturbed site responds is no longer valid. You would need to go beyond that approximation by
\(U = \chi^{-1}-\chi_0^{-1} \approx \left(\frac{\partial N^{\rm SCF}_{I}}{\partial V_{I}}\right)^{-1} - \left(\frac{\partial N^{\rm NSCF}_{I}}{\partial V_{I}}\right)^{-1}\)
with a 2x2 matrix using both sites instead. This should be already possible with the data you computed already. You did not provide the output for the calculation with one Fe perturbed, so I cannot verify.

[Heseri]

Dear Dr. Schlipf,

Thank you very much for your prompt reply. For your reference, I have also uploaded the key files for the linear response calculations in which only one Fe2+ was perturbed. As shown below, the effect of the spherical potential on the d-charge of the unperturbed Fe2+ is negligible in the self-consistent calculations, but significant in the non-self-consistent calculations:

For the unperturbed Fe2+:

Code: Select all

LDAUU/LDAUJ // d charges (self-consistent) // d charges (non-self-consistent)
-0.20  6.221  6.353
-0.15  6.220  6.319
-0.10  6.219  6.285
-0.05  6.218  6.250
0.0  6.217  6.217
0.05  6.216  6.181
0.10  6.215  6.146
0.15  6.214  6.112
0.20  6.213  6.078

For the perturbed Fe2+:

Code: Select all

LDAUU/LDAUJ // d charges (self-consistent) // d charges (non-self-consistent)
-0.20  6.189  6.046
-0.15  6.196  6.089
-0.10  6.203  6.132
-0.05  6.210  6.175
0.0  6.217  6.217
0.05  6.224  6.263
0.10  6.231  6.306
0.15  6.239  6.350
0.20  6.246  6.393

Could you please provide some more detailed guidance on how to calculate the U value using the matrix approach you mentioned?

Additionally, I apologize if I have misunderstood something or am being overly persistent. If two Fe2+ separated by more than 6 angstrom can still influence each other and lead to unreasonable U values, would it be reliable to perturb only one cation when calculating U for systems in which Fe or Ni is major elements, such as NiO, Fe2SiO4, or Fe3O4?

Thank you very much for your time and further guidance. I truly appreciate your patience and help.

Sincerely,
Tian Hua

[martin.schlipf]

Thank you for the additional data. With this we can compute the 2x2 matrix of the response. Assume that χ is a matrix [[a, b], [b, a]] and likewise for χ0. Then using your data, we get non-selfconsistent a0 = 0.868 and b0 = -0.688 whereas self consistent yields a = 0.143, b = -0.020. You get this by computing the differences and dividing by the spread e.g. (6.246 - 6.189) / 0.40 = 0.143.

The method in the tutorial corresponds to neglecting b (U = 1/a - 1/a0 = 5.8). You can go beyond that and invert χ and χ0 the diagonal of these is 7.1 and 3.1, respectively. Therefore you get U = 4.0 which is a bit larger than what I would have expected for this system. You would need to converge this with respect to the supercell size.

However, I also noticed that you used a very large smearing. Is there a reason to use such large values? SIGMA = 0.34 corresponds to a temperature of 4000K. You should probably reduce this and/or switch to a narrower distribution like ISMEAR = 0. At the very least you should study how strongly your results depend on this value.