My Community

Is it possible in VASP to calculate a value for the Hubbard U using linear response theory, as described in PRB 71 035105?

I'd like more information on this too.

This is implemented but not released and not documented.
The developers provide following instructions:
The calculation of U from constrained DFT in the spirit of the work of sec. 3.4 of the thesis work of Coccocioni.

(1) The first step is a well-converged DFT calculation without any additional Hubbard potentials.
Probably it is wise to save the charge density and wave function files for later restarts.
Use something like the following INCAR:

LDAU = .TRUE.
LDAUTYPE = 3
LDAUL = 2 -1 -1
LDAUU = 0.00 0.00 0.00
LDAUJ = 0.00 0.00 0.00
LDAUPRINT= 2
LMAXMIX = 4 # for d-electrons, set to 6 when you're dealing with f-electrons
LORBIT = 11

LDAUTYPE=3 is the one you're going to use to add a shifting potential on the atomic sites.
The present INCAR was used to compute U for NiO.
You see 3 entries for LDAUL, the first pertains to the Ni site were I want to put a shifting potential (on the d-electrons: L=2),
and the second and third entries pertain to the other Ni sites and the oxygen that do not get a shifting potential (L=-1).

So the site were you put your shifting potential must be a species of its own in your POSCAR and POTCAR files.

(2) Then one switches on a Hubbard potential and converges the wave functions for the "bare perturbation" i.e. you do not account
for changes in the potential due to charge redistribution [beware you need to restart from the CHGCAR file of the calculation (1) !!!]:

ICHARG = 11
LDAU = .TRUE.
LDAUTYPE = 3
LDAUL = 2 -1 -1
LDAUU = 0.10 0.00 0.00
LDAUJ = 0.10 0.00 0.00
LDAUPRINT= 2
LMAXMIX = 4 # for d-electrons, set to 6 when you're dealing with f-electrons
LORBIT = 11

In the example above a shifting potential of 0.1 eV is put in both the spin-up as well as the spin-down channel of the d-electrons
(for LDAUTYPE=3, LDAUU and LDAUJ denote the shifts acting on the up and down electrons, respectively).

(3) Do a fully self-consistent calculation with the shifting potential turned on:

LDAU = .TRUE.
LDAUTYPE = 3
LDAUL = 2 -1 -1
LDAUU = 0.10 0.00 0.00
LDAUJ = 0.10 0.00 0.00
LDAUPRINT= 2
LMAXMIX = 4 # for d-electrons, set to 6 when you're dealing with f-electrons
LORBIT = 11

Now it's basically a matter of computing

\chi_IJ = d n^I_d / d \alpha_J

from the calculations (3), and the corresponding \chi^0_IJ from the calculations (2), and U from eq. 3.25 in the thesis of Coccocioni.

(LORBIT=11 will make VASP write the site and L-decomposed charge density, n^I_d, in the OUTCAR file).

The change in n^I_d due to the shifting potential you might get maybe best compute using central differences, i.e., from the
differences between a shift by +V and -V.

Just as a comment, by "LDAUTYPE = 3" you probably mean "LDAUTYPE = 1 or 2", because no "LDAUTYPE = 3" is available as I know just "LDAUTYPE = 1 or 2 or 4". The "LDAUTYPE = 1 and 2" are able to do LSAD+U (spin polarized) but "LDAUTYPE = 4" can just do LDA+U.

LDAUTYPE = 3
is a hidden feature of vasp.5

Thanks for your help, admin!