Calculate on-site U parameter in a self-consistent fashion

Background

The Hubbard model is a revised tihgt-binding model suited to describe solid-state systems with prominent correlation effects (i.e. narrow band transition metal oxides - Mott Transition). Traditional tight-binding model does not consider the coulomb energy originated from the overlap between orbitals, in that sense, they treat electrons as isolated particles. However, in some systems, the overlap between atomic orbitals is large enough so that the Coulomb interaction energy cannot be ignored. To solve this, the Hubbard model adds a on-site repulsion term to represent the on-site electron-electron Coulomb repulsion (Other inter-site energy terms can also be added, they are usually called Vs and Js). The Hamiltonian can be written as:

The first term is the traditional transfer integral identical to the tight-binding model, where in the second term, the Hubbard U parameters describe the on-site Coulomb repulsion strength.

The so called “DFT+U” method can be used to give a more accurate description to the band gap and even recover the insulation nature of some systems. One thing to point out, this method is based on the a projection weight penalty term added to the Kohn-Sham Hamiltonian which means the charge density and the total energy are fully self-consistent.

Where the Hubbard term direct counts the coulomb energy with a specific orbital electron occupation. And the last term is for the double-counting correction since both exchange-correlation and the Hubbard terms have exchange contribution.

Although this method can be somehow sufficient for the description of strongly correlated materials, the U parameter are typically chosen arbitrarily or fitted to the experimental result. To obtain this value from first-principles (or “in a self-consistent fashion”), several methods like linear response, Density functioanl pertubation theory (DFPT) and constrained Random Phase Approximation (cRPA) have been proposed. Here, I’ll only give a technical description to the first one as proposed by Matteo Cococcioni[1].

Calculation procedure

Taking the BCC iron as example, the calculation will be conducted using VASP (v5.4.4) and Quantum-espresso (v6.4.1):

VASP procedure

The on-site potential shift \(\alpha\) control is a hidden feature of VASP(v5.4+) and is not documented in the official documentation. The documentations I found online was not very through: Marianetti Group VASP forum

To invoke this hidden feature:

LDAU = .TRUE.       # invoke DFT+U procedure
LDAUTYPE = 3        # potential shift mode
LDAUPRINT = 2       # print occupation and potential matrix
LDAUU = 0.0 0.0     # energy shift of spin-up electrons
LDAUJ = 0.0 0.0     # energy shift of spin-down electrons
LDAUL = 2 -1        # which orbital to perform procedure

Procedure outline:

1. Make a supercell large enough (Since the effect of onsite potential shift can also affect its neighboring atoms, under the periodic boundary conditions, we have to make a supercell large enough to eliminate interaction between changed site). Here I made a \(2\times2\times2\) supercell of BCC iron with 8 iron atoms.
2. We only want to shift the potential on one iron site, the POTCAR file should contain the two identical pseudopotential information for iron.
3. At first, we do a self consistent calculation WITHOUT potential shift, and save the charge density in CHGCAR file.

LCHARG=.TRUE.
LDAU = .TRUE.
LDAUTYPE = 3
LDAUPRINT = 2
LDAUU = 0.0 0.0
LDAUJ = 0.0 0.0
LDAUL = 2 -1
... + other scf tags

1. Then, a “bare run” non-self consistent calculation WITH potential shift should be carried out using previous saved CHGCAR.

ICHARG=11
LDAU = .TRUE.
LDAUTYPE = 3
LDAUPRINT = 2
LDAUU = 0.02 0.0
LDAUJ = 0.02 0.0
LDAUL = 2 -1
... + other scf tags

1. Similarly, a “interactive run” self consistent calculation WITH potential shift to calculate the charge change with screening.(Here, I choose to read in the charge density since it can sped up the process.)

ICHARG=1
LDAU = .TRUE.
LDAUTYPE = 3
LDAUPRINT = 2
LDAUU = 0.02 0.0
LDAUJ = 0.02 0.0
LDAUL = 2 -1
... + other scf tags

1. Repeat step 4-5 with different potential shift (typically, I use -0.08 to 0.08 eV as potential shift magnitude)

2. Check projected charge on atomic site, and plot it against the potential shift magnitude. Linear fitting will give the so-called \(\chi\) value for “bare run” and “interactive run”.

   \[\chi=\frac{\partial Q}{\partial \alpha}\]

   Where \(Q\) is the projected charge number (here d-orbital), \(\alpha\) is the potential shift used.

3. Finally, depend on accuracy requirement, either reconstruct the whole response matrix or use only on-site \(\chi\) to calculate the \(U\) value.

   \[U=\chi_0^{-1}-\chi^{-1}\]

   where \(\chi_0\) is the bare run result, and \(\chi\) is the interactive run result.

Results:

• PBE result:

• LDA result:

We can see this value is significantly higher than the published value.

• PBE with \(3\times3\times3\) supercell:

The size of supercell does not affects the final result significantly, so a supercell of \(2\times2\times2\) should be enough for this system.

Response Matrix generation:

There are 4 nearest neighbor (labeled red) and 3 second-near neighbor (white) in the \(2\times2\times2\) supercell. The diagonal terms are the onsite charge number (blue).

From the response matrix, we can see that the symmetry is broken a little especially for atom 8 and atom 7.

QE procedure

Similar procedure can also be performed by quantum espresso. To enable this feature, following tags are needed:

lda_plus_U=.true.
Hubbard_U(1)=1D-40,Hubbard_alpha(1)=-0.05

The output of this method contains 3 blocks of projected charge information.

1. initial charge (zero potential shift)
2. Bare response
3. SCF (interactive) response.

So one only need to do one calculation to extract all projected charge number needed for the calculation.

In addition, QE also allows to shift potential \(\alpha\) with U applied. Simply apply:

lda_plus_U=.true.
Hubbard_U(1)="Your U parameter",Hubbard_alpha(1)=-0.05

Results:

Final result using only onsite response:

Final result using full matrix:

Input

I’ve put all input file in a zip file for download: VASP, QE.