Linear response calculation of U?

[xamnidar]

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

[nsolomat]

I'd like more information on this too.

[admin]

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.

[szvasp]

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.

[admin]

LDAUTYPE = 3
is a hidden feature of vasp.5

[xamnidar]

Thanks for your help, admin!