2.3. DFT+U

In DFT, the interactions between electrons are discribed by the direct Coulomb intereaction term (the Hartree term) and the exchange-correlation (\(xc\)) term. In LDA and GGA, the exchange-correlation functional only depends on electron denisty and gradients of the density, any orbital dependence is neglected. The absence of orbital dependence may lead to accuracy issues in the predicted results, since LDA and GGA tend to over-delocalize valence electrons [KCS].

The DFT+U method takes into account orbital dependence of the Coulomb and exchange interactions. In particular, the Hubbard (U) model is formulated to more propertly describe the behavior of highly localized d and f electrons, for example. The U parameter in typical DFT+U can be determined from ab-initio first principles calculations, but is usually obtained semi-empirically. [CG]

In DFT+U, there is an additional energy item added to the exchange correlation energy:

\[E_{U}=\sum\left(\frac{1}{2}\right)U_{\mu}(n_{\mu}-n^{2}_{\mu}),\]

where, \(n_{\mu}\) is the projection on the atomic shell, \(U_{\mu}\) is the U value corresponding to the atomic shell, and the \(E_{U}\) energy term is zero for fully occupied or unoccupied shells, and positive for slightly occupied shells. Therefore, if the state is fully occupied, the energy is reduced. This can occur if the energy level is moved away from the Fermi level, i.e. the energy band gap is increased, or if the broadening of the states is reduced, i.e. the electrons are localized. In this way, Hubbard U improves the accuracy of LDA and GGA functionals.

In this tutorial, we will calculate the electronic density of states (DOS) using the Hubbard U correction for transition metal oxide NiO as shown in Fig. 2.3.1. Specifically, the band gap of NiO crystal is too low in LDA or GGA predictions, and DFT+U improves it.

../../_images/NiO-struct.svg

Fig. 2.3.1 Schematic view of Ni crystal.

The calculation steps are:

  1. Build:

    • NiO unit cell;

    • Expand the primitive NiO unit cell;

  2. Self-consistent calculation:

    • Generate the input files;

    • SCF calculation;

    • SCF calculation including Hubbard correction;

2.3.1. Importing the crystal

For purposes of this tutorial, we will import a NiO crystal from Device Studio data base.

Begin by starting the Device Studio software and creating a new project.

  1. Open Device Studio software ➟ Create a new Project ➟ File name NiO ➟ press Save;

  2. In the menu bar click File ➟ Import ➟ Import Local to pop up the search box:

    • Go to DeviceStudio folder ➟ material ➟ 3Dmaterials ➟ metal_oxides ➟ select NiO.hzw ➟ click open.

  3. In the menu bar click Build ➟ Redefine Crystal;

    • In the cell vector set u = (2a, 0b, 0c), v = (0a, 1b, 0c), and w = (0a, 0b, 1c), click Preview or Build to expand the cell 2 times in the a direction, click on Preview and Build the structure.

Tip

  • It is imperative to expand the unit cell in order to modify the magnetic moment and achieve an antiferromagnetic arrangement, as the original cell configuration results in a ferromagnetic structure.

  • Only to to make the calculations quicker we’ll not optimize the atom postitions and cell vectors. For the propouse of your project be sure about the equilibrium geometry.

The atomic coordenates for NiO primitive cell can be download here: NiO.hzw.

NiO.hzw

% The number of probes
0
% Uni-cell vector
0.00000000  2.00000000      2.00000000
2.00000000  0.00000000      2.00000000
2.00000000  2.00000000      0.00000000
%Total number of device_structure
2
%Atom site
Ni  1  1  1
O  3  3  3

2.3.2. SCF for NiO

Generate the input file

The user migth be able to set parameter and perform the standard scf calculation. The steps for this calculation can be found in our previous tutorial SCF calculation. Now we will set the parameters to perform the DFT+U with NanoDCAL module simulator, as shown in Fig. 2.3.2.

  1. In the menu bar of Device Studio, click Simulator ➟ NanoDCAL ➟ SCF Calculation to pop up the interface.

    • On Basic settings change the parameter k-point Sampling to n1 = 15, n2 = 15, and n3=15.

    • Go to the Basis set check the LDA+U option ➟ set for Ni Atom s = 0, p = 0, d = 4.6, f = 0

    • Go to Spin Type ➟ select CollinearSpin. On Spin polarization set Ni = 1 and O = 0.

  2. Click Generate files

../../_images/DFT-U-scf.svg

Fig. 2.3.2 Graphical user interface of Device Studio to set the calculation parameters.

Scf input file

%%What quantities should be calculated
calculation.name = scf

%Basic setting
calculation.occupationFunction.temperature = 300
calculation.realspacegrids.E_cutoff = 80 Hartree
calculation.xcFunctional.Type = LDA_PZ81
calculation.k_spacegrids.number = [ 15 15 15 ]'

system.centralCellVectors = [[0 2 2]' [2 0 2]' [4 4 0]']
system.spinType = CollinearSpin

%Hubbard U setting
calculation.Hubbard.isIncluded = true
calculation.Hubbard.parameterBlock = 2
        Ni  LDA-DZP    0.00    0.00    4.60    0.00
        O  LDA-DZP    0.00    0.00    0.00    0.00
end

%Iteration control
calculation.SCF.monitoredVariableName = {'rhoMatrix','hMatrix','totalEnergy','bandEnergy','gridCharge','orbitalCharge','spinPolar'}
calculation.SCF.convergenceCriteria = {1e-04,1e-04,[],[],[],[],[]}
calculation.SCF.maximumSteps = 200
calculation.SCF.mixMethod = Pulay
calculation.SCF.mixRate = 0.1
calculation.SCF.mixingMode = H
calculation.SCF.startingMode = H
%calculation.SCF.donatorObject = NanodcalObject.mat

%Basic set
system.neutralAtomDataDirectory = '../'
system.atomBlock = 4
AtomType OrbitalType X Y Z SpinPolarization
Ni  LDA-DZP 1.00000000      1.00000000      1.00000000      1
Ni  LDA-DZP 3.00000000      3.00000000      1.00000000      1
O   LDA-DZP 3.00000000      3.00000000      3.00000000      0
O   LDA-DZP 5.00000000      5.00000000      3.00000000      0
end

During DFT+U calculations on a crystalline system, the symmetry used is often lower than the actual crystallographic symmetry. This mismatch can reduce computational efficiency, requiring more k-points to accurately represent the system in reciprocal space.

The SpinPolarization column is useful to break the symmetry between spin up and spin down, allowing magnetic solutions. In this case the ferromagnetic (FM) ordering is the default value. Therefore, in order to obtain the antiferromagnetic (AFM) ordering, we set a different SpinPolarization pattern in the scf input files, as follows:

  1. Returning to the Project, navigate to the folder ➟ right-click on the scf.input file to pop up the file;

    • Edit the SpinPolarization column, Ni = 1 and Ni= -1 ➟ Press save.

  2. Close the file.

AtomType OrbitalType X Y Z SpinPolarization
Ni  LDA-DZP 1.00000000      1.00000000      1.00000000      1
Ni  LDA-DZP 3.00000000      3.00000000      1.00000000      -1
O   LDA-DZP 3.00000000      3.00000000      3.00000000      0
O   LDA-DZP 5.00000000      5.00000000      3.00000000      0

The keywords for U correction are specify the following:

  • calculation.Hubbard.isIncluded If it is true or 1, a Hubbard-like term, the so called LDA+U, is added to the potential.

  • calculation.Hubbard.parameterBlock This block is used to overwrite the default Hubbard (U) values of specified atomic basis. There must be exactly N lines in the block, and followed by a line of end used as a end mark of the block. N is the number of atomic species whose Hubbard values needs to be re-defined. Each line has two strings, which specify the AtomType and OrbitalType of the atomic basis respectively, and enough U values for s, p, d, and f orbitals, respectively. Note that the unit of the U values is defined by the parameter calculation.control.energyUnit.

The generated scf input file for NiO for DFT and DFT+U couldbe downloaded here:

Important

  • It is extremally important note that there is no default values for Hubbard projectors. Moreover, the pseudopotentials used for the DFT+U has to be carefully evaluated.

  • The value specified for Hubbard projectors can be extracted by fitting various experimental properties. Or by atomic calculations used to generate the respective pseudopotentials.

  • In magnetic systems can be found fifferents local minima. Thus, the inicial occupation has be carefully selected.

To perform the calculation the users can also use Device Studio, as following:

On the Device Studio navigate to the Project panel:

  1. Go to the NiO folder ➟ right click over scf.input file ➟ run to pop up the Run dialog window;

  2. In the Run window, choose the Gateway location, number of cores and press run;

  3. Once the calculation has completed, the results will be returned to the NiO folder.

Optionally, several ground state properties of the system could also be calculated from the obtained SCF results. In the next step, we’ll show how to calculate the electronic density of states.

2.3.3. DOS for NiO

  1. In the Device Studio navigate to the menu bar, click on Simulator ➟ NanoDCAL ➟ Analysis to set physical quantity, as shown in Fig. 2.3.3.

  2. Navigate to the Analysis panel (left side) ➟ select DensityOfStates ➟ click on the Right arrow button to transfer DensityOfStates calculator to the Calculation Selected panel.

    • Navigate to the Analysis panel ➟ Calculation Selected click on select DensityOfStates ➟ in k-pint sampling set n1 = 15, n2 = 15, n3 = 15

  3. click on Generate files.

../../_images/DOS-u.svg

Fig. 2.3.3 Physical properties interface of the Device Studio.

The Device Studio will create new files in the NiO_Rede directory that will be the DensityOfStates input file.

DOS input file

system.object = NanodcalObject.mat
calculation.name = densityOfStates
calculation.densityOfStates.kSpaceGridNumber = [ 15 15 15 ]'
calculation.densityOfStates.numberOfEnergyPoints = 401
calculation.densityOfStates.energyRange = [-10,10]
calculation.densityOfStates.whatProjected = 'Atom'
%calculation.densityOfStates.plot = true
calculation.control.xml = true

The keywords specifies the following:

  • calculation.densityOfStates.kSpaceGridNumber The number of small k-space grids in each direction which, together with kSpaceGridShift, are used to produce the parameter kSpacePoints.

  • calculation.densityOfStates.numberOfEnergyPoints When isIntegrated is true, this parameter is used to perform the integration. When isIntegrated is false, this parameter is used together with energyRange, to define the parameter energyPoints.

  • calculation.densityOfStates.energyRange When isIntegrated is true, this parameter defines the energy range over which dos is integrated. When isIntegrated is false, this parameter is used together with numberOfEnergyPoints, to define the parameter energyPoints. Note that for bulk system its values are measured from the Fermi energy, and for system with leads its values are measured from the chemical potential of a lead that has zero applied voltage.

  • calculation.densityOfStates.whatProjected If ‘Atom’ or ‘Orbital’, the dos is projected on the given orbitals or atoms. If ‘Ell’, the dos is projected on orbitals with the given angular momentum. If ‘None’, not projected.

The generated input file could be downloaded or users can use Device Studio to peforme the calculations, as we shown above (run SCF).

Run the simulation with the DensityOfStates input file. This calculation will produced DensityOfStates dataset (.xml and .mat).

DOS analysis

Next, you should analyses the results from the produced density of states dataset (.xml and .mat) that can then be plotted using the MATLAB or Device Studio by:

  • In menu bar of Device Studio, click on Simulator ➟ NanoDCAL ➟ Analysis Plot ➟ import DensityOfStates.xml file.

For the NiO system, the Hubbard correction has a significant impact on the density of states (DOS) near the Fermi level. Without the Hubbard correction, standard DFT calculations tend to underestimate the bandgap of NiO, as shown in Fig. 2.3.4.

../../_images/DOS-spec.svg

Fig. 2.3.4 The calculated DOS for NiO crystal by DFT and DFT+U.

By introducing the Hubbard U term into the calculations, the on-site electron-electron interactions are properly accounted for, leading to a more accurate description of the electronic structure of NiO. Specifically, the Hubbard correction tends to shift the 3d states of Ni towards higher energies, effectively increasing the bandgap of NiO and promoting its insulating behavior. This results in a better agreement between theoretical calculations and experimental observations of the electronic properties of NiO, such as its insulating nature and magnetic ordering. In summary, the Hubbard correction plays a crucial role in accurately describing the electronic structure of transition metal oxides like NiO, particularly in capturing the insulating behavior and magnetic properties that arise from strong electron-electron interactions in these systems.

Tip

As an additional task the users should be able to:

  • Calculate the ferromagnetic (FM) ordoring of NiO structure. Without the Hubbard correction, standard DFT calculations the FM ordoring tend predict it to be metallic.

  • The user can repeat the procedure and import different structures from our data base to verify the influence of the Hubbard (U) correction.

[KCS]

    1. Kulik, M. Cococcioni, D. A. Scherlis, and N. Marzari. Density Functional Theory in Transition-Metal Chemistry: A Self-Consistent Hubbard U Approach Phys. Rev. Lett. 97 (2006), 103001.