In this section, we provide two examples for self-consistently calculating the U parameter. One example is with the Ruddlesden Popper perovskite, Sr $_3$ Ti $_2$ O $_7$ , where the Hubbard sites are symmetrically equivalent, and therefore the $U^i$ is equivalent for each. In the other example, a doped Ruddlesden Popper perovskite is considered with inequivalent Hubbard sites.

Example 1

First an initial Self-Consistent Field (SCF) calculation must be done with regular DFT or experimental lattice parameters to obtain the KS states.The example input for Sr $_3$ Ti $_2$ O $_7$ . This should, of course, be done after performing convergence tests with respect to k-points and plane wave cutoffs. Please note that the SSSP efficiency V. 1.0 pseudopotential library is used for all calculations within this example and are available here.

The inputs within this section constitute one iteration in the self-consistent procedure.

The input file for the SCF calculation:

&control
calculation='scf'
restart_mode='from_scratch',
prefix='pwscf'
pseudo_dir = '/path/to/pseudopotentials'
outdir='./temp'
verbosity='high'
/
&system
ibrav = 0,
nat = 24,
ntyp = 3,
ecutwfc = 50.0
ecutrho = 400.0
lda_plus_u = .true.,
Hubbard_U(2) = 1.d-8
/
&electrons
conv_thr = 1.d-15
mixing_beta = 0.7
/
ATOMIC_SPECIES
Sr 87.62 Sr_pbe_v1.uspp.F.UPF
Ti 47.867 ti_pbe_v1.4.uspp.F.UPF
O 16.0 O.pbe-n-kjpaw_psl.0.1.UPF

ATOMIC_POSITIONS {crystal}
Ti 0.0 0.0 0.097768
Ti 0.5 0.5 0.402232
Ti 0.5 0.5 0.597768
Ti 0.0 0.0 0.902232
Sr 0.0 0.0 0.315737
Sr 0.5 0.5 0.0
Sr 0.5 0.5 0.184263
Sr 0.5 0.5 0.815737
Sr 0.0 0.0 0.5
Sr 0.0 0.0 0.684263
O 0.5 0.0 0.403661
O 0.5 0.5 0.306192
O 0.0 0.0 0.0
O 0.5 0.0 0.096339
O 0.0 0.5 0.096339
O 0.0 0.5 0.403661
O 0.0 0.0 0.193808
O 0.0 0.5 0.903661
O 0.0 0.0 0.806192
O 0.5 0.5 0.5
O 0.0 0.5 0.596339
O 0.5 0.0 0.596339
O 0.5 0.0 0.903661
O 0.5 0.5 0.693808
CELL_PARAMETERS {angstrom}
3.940908 0.0 0.0
0.0 3.940908 0.0
0.0 0.0 20.511886
K_POINTS {automatic}
7 7 1 0 0 0

Note the stringent electronic convergence threshold, ‘conv_thr’ is required to ensure well-converged Kohn Sham states; the U parameter for all Ti atoms (symmetrically equivalent) is set to a small number to initialize the DFT+U calculation.This calculation is run with the pw.x code in quantum espresso. The results of this calculation are needed to perform the DFPT calculation.

The input for the DFPT calculation is:

&inputhp
prefix = 'pwscf',
outdir = './temp',
nq1 = 2, nq2 = 2, nq3 = 2,
conv_thr_chi = 1.0d-8,
iverbosity = 2
/

The ‘prefix’ and the ‘outdir’ must be the same as in the SCF input. The value ‘conv_thr_chi’ is the convergence threshold for the trace of the response matrix, inversely related to the second derivatives of the energy with respect to occupation. As with other DFPT calculations such as phonon calculations, the response should be converged with respect to supercell size or in Quantum Espresso, the q-point grid. A q-point grid of 2x2x2 is used in this example for the sake of brevity. This calculation is ran with the hp.x code in Quantum Espresso after the SCF calculation. The results, including the parameters and response matrices, are printed to the ‘<prefix>.Hubbard_parameters.dat’ file.

=-------------------------------------------------------------------=

Hubbard U parameters:

site n. type label spin new_type new_label Hubbard U (eV)
1 2 Ti 1 2 Ti 4.1634
2 2 Ti 1 2 Ti 4.1634
3 2 Ti 1 2 Ti 4.1634
4 2 Ti 1 2 Ti 4.1634

=-------------------------------------------------------------------=
...

The effective U parameters are printed per site. With this initial estimate of the U parameter the convergence of the U is not met ( $U_n = 4.1634$ – $U_{n-1} = 1.e-8$ > $1.e-5$ ) and a cell optimization should be performed. This is the second step in the self-consistent procedure. To do this the SCF calculation input can be adjusted and ran with the pw.x code in Quantum Espresso. The calculation flag is changed to a variable-cell relaxation. Additional flags are added to ensure convergence during the calculation and change the calculation type, but the key difference is the updated U parameter.

&control
calculation='vc-relax'
restart_mode='from_scratch',
prefix='pwscf'
pseudo_dir = '/path/to/pseudopotentials'
outdir='./temp'
verbosity='high'
/
&system
ibrav = 0,
nat = 24,
ntyp = 4,
occupations='smearing'
smearing = 'mv'
degauss = 0.01
ecutwfc = 50.0
ecutrho = 400.0
lda_plus_u = .true.,
Hubbard_U(2) = 4.1634
/
&electrons
conv_thr = 1.d-10
mixing_beta = 0.7
/
&ions
ion_dynamics='bfgs'
/
&cell
cell_dynamics='bfgs'
/

ATOMIC_SPECIES
Sr 87.62 Sr_pbe_v1.uspp.F.UPF
Ti 47.867 ti_pbe_v1.4.uspp.F.UPF
La 138.9 La.GGA-PBE-paw-v1.0.UPF
O 16.0 O.pbe-n-kjpaw_psl.0.1.UPF
ATOMIC_POSITIONS {crystal}
Ti 0.0 0.0 0.097768
Ti 0.5 0.5 0.402232
Ti 0.5 0.5 0.597768
Ti 0.0 0.0 0.902232
Sr 0.0 0.0 0.315737
Sr 0.5 0.5 0.0
Sr 0.5 0.5 0.184263
Sr 0.5 0.5 0.815737
Sr 0.0 0.0 0.5
Sr 0.0 0.0 0.684263
O 0.5 0.0 0.403661
O 0.5 0.5 0.306192
O 0.0 0.0 0.0
O 0.5 0.0 0.096339
O 0.0 0.5 0.096339
O 0.0 0.5 0.403661
O 0.0 0.0 0.193808
O 0.0 0.5 0.903661
O 0.0 0.0 0.806192
O 0.5 0.5 0.5
O 0.0 0.5 0.596339
O 0.5 0.0 0.596339
O 0.5 0.0 0.903661
O 0.5 0.5 0.693808
CELL_PARAMETERS {angstrom}
3.940908 0.0 0.0
0.0 3.940908 0.0
0.0 0.0 20.511886

K_POINTS {automatic}
7 7 1 0 0 0

This marks the end of the iteration in the iterative procedure.

After structural optimization, the primary difference are the lattice parameters: $a = 3.954898889 \AA$ and $c = 20.535525478 \AA$ . We return the first step in the self-consistent procedure. We perform another SCF calculation with the updated lattice parameters and the $U = 4.1634$ parameter. From here, the procedure is the same is in the block above marked by the indented text, the input is provided for completion though.

&control
calculation='scf'
restart_mode='from_scratch',
prefix='pwscf'
pseudo_dir = '/path/to/pseudopotentials'
outdir='./temp'
verbosity='high'
/
&system
ibrav = 0,
nat = 24,
ntyp = 3,
ecutwfc = 50.0
ecutrho = 400.0
lda_plus_u = .true.,
Hubbard_U(2) = 4.1634
/
&electrons
conv_thr = 1.d-15
mixing_beta = 0.7
/
ATOMIC_SPECIES
Sr 87.62 Sr_pbe_v1.uspp.F.UPF
Ti 47.867 ti_pbe_v1.4.uspp.F.UPF
O 16.0 O.pbe-n-kjpaw_psl.0.1.UPF
CELL_PARAMETERS (angstrom)
3.954898889 0.000000000 0.000000000
0.000000000 3.954898889 0.000000000
0.000000000 0.000000000 20.535525478

ATOMIC_POSITIONS (crystal)
Ti 0.000000000 0.000000000 0.097986454
Ti 0.500000000 0.500000000 0.402013539
Ti 0.500000000 0.500000000 0.597986461
Ti 0.000000000 0.000000000 0.902013546
Sr 0.500000000 0.500000000 0.184879041
Sr 0.000000000 0.000000000 0.315120943
Sr 0.500000000 0.500000000 0.000000000
Sr 0.500000000 0.500000000 0.815120959
Sr 0.000000000 0.000000000 0.500000000
Sr 0.000000000 0.000000000 0.684879057
O 0.500000000 0.000000000 0.403281390
O 0.500000000 0.500000000 0.305549045
O 0.000000000 0.000000000 0.000000000
O 0.500000000 0.000000000 0.096718599
O 0.000000000 0.500000000 0.096718599
O 0.000000000 0.500000000 0.403281390
O 0.000000000 0.000000000 0.194450945
O 0.000000000 0.500000000 0.903281401
O 0.000000000 0.000000000 0.805549055
O 0.500000000 0.500000000 0.500000000
O 0.000000000 0.500000000 0.596718610
O 0.500000000 0.000000000 0.596718610
O 0.500000000 0.000000000 0.903281401
O 0.500000000 0.500000000 0.694450955

K_POINTS {automatic}
7 7 1 0 0 0

This process is repeated until the desired convergence is reached. In this example, the following data are obtained with U parameters given in eV:

Iteration	0	1	2	3
Effective U	0.00000001	4.1634	3.9907	3.9983

$U_3$ – $U_2$ = 8 meV. From here, the effective U parameter $U_3 = 3.9983$ and corresponding lattice constants should be used in calculations to obtain physical quantities. In this example, we highlight the effect on the charge density. By taking the difference in the charge densities with and without the U parameter, we can highlight the relocalization of electrons near transition metal atom centers.

Fig. 1) (Left) View perpendicular to the (110) plane in the Sr $_3$ Ti $_2$ O $_7$ structure after relaxation with the self-consistent effective U correction. (Right) Electronic charge density difference with and without the effective U parameter. The blue isosurface levels indicate regions of reduced electronic density, and red isosurface levels indicate increased electronic density. The green isosurface levels show small levels of charge re-localization along the metal-oxygen bonds.

Example 2

Since Example 1 was in a structure with symmetrically equivalent Hubbard sites, the effective U parameter is the same for all sites. There may be systems where the sites are inequivalent such as in alloy, defect, and low symmetry cells. In these, the effective U corrections can vary per site; the response in the electronic occupations to perturbations varies with local environment. This example is for such a system. The perovskite in Example 1, is often doped with La to obtain certain electronic properties. We take a simple La doped Sr $_3$ Ti $_2$ O $_7$ cell and compute the effective U parameters for each site.

The first SCF calculation requires the initialization of each Hubbard site, to do this a different ‘atom type’ is assigned to each:

&control
calculation='scf'
restart_mode='from_scratch',
prefix='pwscf'
pseudo_dir = '/path/to/pseudopotentials'
outdir='./temp'
verbosity='high'
/
&system
ibrav = 0,
nat = 24,
ntyp = 7,
occupations='smearing'
smearing = 'mv'
degauss = 0.01
ecutwfc = 50.0
ecutrho = 400.0
lda_plus_u = .true.,
Hubbard_U(2) = 1.d-8
Hubbard_U(3) = 1.d-8
Hubbard_U(4) = 1.d-8
Hubbard_U(5) = 1.d-8
Hubbard_U(6) = 1.d-8

/
&electrons
conv_thr = 1.d-15
mixing_beta = 0.7
/
ATOMIC_SPECIES
Sr 87.62 Sr_pbe_v1.uspp.F.UPF
Ti1 47.867 ti_pbe_v1.4.uspp.F.UPF
Ti2 47.867 ti_pbe_v1.4.uspp.F.UPF
Ti3 47.867 ti_pbe_v1.4.uspp.F.UPF
Ti4 47.867 ti_pbe_v1.4.uspp.F.UPF
La 138.9 La.GGA-PBE-paw-v1.0.UPF
O 16.0 O.pbe-n-kjpaw_psl.0.1.UPF

ATOMIC_POSITIONS {crystal}
Ti1 0.0 0.0 0.097768
Ti2 0.5 0.5 0.402232
Ti3 0.5 0.5 0.597768
Ti4 0.0 0.0 0.902232
La 0.5 0.5 0.184263
Sr 0.0 0.0 0.315737
Sr 0.5 0.5 0.0
Sr 0.5 0.5 0.815737
Sr 0.0 0.0 0.5
Sr 0.0 0.0 0.684263
O 0.5 0.0 0.403661
O 0.5 0.5 0.306192
O 0.0 0.0 0.0
O 0.5 0.0 0.096339
O 0.0 0.5 0.096339
O 0.0 0.5 0.403661
O 0.0 0.0 0.193808
O 0.0 0.5 0.903661
O 0.0 0.0 0.806192
O 0.5 0.5 0.5
O 0.0 0.5 0.596339
O 0.5 0.0 0.596339
O 0.5 0.0 0.903661
O 0.5 0.5 0.693808
CELL_PARAMETERS {angstrom}
3.940908 0.0 0.0
0.0 3.940908 0.0
0.0 0.0 20.511886
K_POINTS {automatic}
7 7 1 0 0 0

The Hubbard site closest to the La defect is Ti1, and the other sites go in order of increasing z-position. Following the procedure in the first example, the following U parameters are obtained (in eV).

Iteration	0	1	2	3
Ti1	0.00000001	4.9778	4.5800	4.5815
Ti2	0.00000001	4.5982	4.6048	4.606
Ti3	0.00000001	4.3059	4.5378	4.5277
Ti4	0.00000001	4.3440	4.6418	4.6474
La	0.00000001	0.9294	0.7769	0.7794

Note that the effective U is different for each site, and intuitively it should be. The response for a Hubbard site closer to a defect should not be the same as in a pristine perovskite crystal. We highlight the charge density differences again, but this time showing that the extent of charge relocalization changes per Hubbard site.

Fig. 2) (Left) View perpendicular to the (110) plane in the doped Sr $_3$ Ti $_2$ O $_7$ structure after relaxation with the self-consistent effective U corrections. (Right) Electronic charge density difference with and without the effective U parameters.

General notes and remarks

In Quantum Espresso, this self-consistent Hubbard U calculation can be done for most d-shell transition metals.The self-consistent Hubbard parameters can be calculated for transition metal systems with arbitrary doping concentration and alloy configurations. By construction, the calculated U parameter values will change per configuration. A rigorous approach to calculating corrected properties in doped or alloy systems is to recalculate the U parameter for every symmetrically inequivalent structure used. A more common and less computationally intensive approach is to calculate the self-consistent U parameter for a single concentration/configuration, and using this for all others. Similarly, when used for slab systems, it is common to calculate the U parameter at a given adsorbate coverage and applying this parameter for all other adsorbate coverages considered.

Self-consistent Hubbard-U corrections for transition metal oxides

Using Quantum Espresso to calculate U parameters

Example 1

Example 2

Fig. 2) (Left) View perpendicular to the (110) plane in the doped Sr $_3$ Ti $_2$ O $_7$ structure after relaxation with the self-consistent effective U corrections. (Right) Electronic charge density difference with and without the effective U parameters.

General notes and remarks

Example 1

Example 2

Fig. 2) (Left) View perpendicular to the (110) plane in the doped SrTiO structure after relaxation with the self-consistent effective U corrections. (Right) Electronic charge density difference with and without the effective U parameters.

General notes and remarks

Fig. 2) (Left) View perpendicular to the (110) plane in the doped Sr $_3$ Ti $_2$ O $_7$ structure after relaxation with the self-consistent effective U corrections. (Right) Electronic charge density difference with and without the effective U parameters.