Searching for a shell
Calculate onsite U parameter in a selfconsistent fashion
Background
The Hubbard model is a revised tihgtbinding model suited to describe solidstate systems with prominent correlation effects (i.e. narrow band transition metal oxides  Mott Transition). Traditional tightbinding 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 onsite repulsion term to represent the onsite electronelectron Coulomb repulsion (Other intersite 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 tightbinding model, where in the second term, the Hubbard U parameters describe the onsite 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 KohnSham Hamiltonian which means the charge density and the total energy are fully selfconsistent.
\[E_{LDA+U}[n(\textbf{r})] = E_{LDA}[n(\textbf{r})] + E_{Hub}[{n^{I \sigma}_{m}}]  E_{DC}[{n^{I \sigma}_{m}}]\]Where the Hubbard term direct counts the coulomb energy with a specific orbital electron occupation. And the last term is for the doublecounting correction since both exchangecorrelation 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 firstprinciples (or “in a selfconsistent 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 Quantumespresso (v6.4.1):
VASP procedure
The onsite 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 spinup electrons
LDAUJ = 0.0 0.0 # energy shift of spindown electrons
LDAUL = 2 1 # which orbital to perform procedure
Procedure outline:
 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.
 We only want to shift the potential on one iron site, the
POTCAR
file should contain the two identical pseudopotential information for iron.  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
 Then, a “bare run” nonself 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
 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
 Repeat step
45
with different potential shift (typically, I use0.08 to 0.08 eV
as potential shift magnitude) 
Check projected charge on atomic site, and plot it against the potential shift magnitude. Linear fitting will give the socalled \(\chi\) value for “bare run” and “interactive run”.
\[\chi=\frac{\partial Q}{\partial \alpha}\]Where \(Q\) is the projected charge number (here dorbital), \(\alpha\) is the potential shift used.

Finally, depend on accuracy requirement, either reconstruct the whole response matrix or use only onsite \(\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 secondnear 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)=1D40,Hubbard_alpha(1)=0.05
The output of this method contains 3 blocks of projected charge information.
 initial charge (zero potential shift)
 Bare response
 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:
\[U=3.29 eV\]Final result using full matrix:
\[U=2.94 eV\]Input
I’ve put all input file in a zip file for download: VASP, QE.