DFT+U

The objective of this tutorial is to give an idea of how DFT+U in Octopus works. As a prototypical example, we will consider bulk NiO in its anti-ferromagnetic configuration.

Details on the implementation of the DFT+U can be found in Ref.¹ .

Bulk NiO with PBE

We will start by setting up an input file for the system we want to simulate, but without DFT+U, as there are a few important things to consider in order to do this correctly.

Input

The input file we will use is the following one:

CalculationMode = gs
PeriodicDimensions = 3
BoxShape = parallelepiped

a = 7.8809
%LatticeParameters
  a | a | a
%
%LatticeVectors
 0.0 | 1/2 | 1/2
 1/2 | 0.0 | 1/2
 1.0 | 1.0 | 0.0
%

%Species
"Ni1" | species_pseudo | set | hscv_pbe
"Ni2" | species_pseudo | set | hscv_pbe
%

%ReducedCoordinates
 "Ni1" | 0.0 | 0.0 | 0.0
 "Ni2" | 0.0 | 0.0 | 0.5
 "O"  | 0.5 | 0.5 | 0.25
 "O"  | 0.5 | 0.5 | 0.75
%
Spacing = 0.4
%KPointsGrid
2 | 2 | 2
%
KPointsUseSymmetries = yes
KPointsUseTimeReversal = no

ExtraStates = 3

SpinComponents = polarized
GuessMagnetDensity = user_defined
%AtomsMagnetDirection
 8.0
-8.0
 0.0
 0.0
%

Let’s now look more in detail at some of the other input options:

• LatticeVectors and LatticeParameters: Here we have decided to neglect the small lattice distortion of bulk NiO in its anti-ferromagnetic configuration and consider only its cubic cell. Note that, as we are interested in the antiferromagnetic order, the primitive cell is doubled along the last lattice vector.

• Species and ReducedCoordinates: Octopus uses the spglib library for detecting crystal symmetries, but spglib cannot detect the magnetic space group of a system at the present time. Therefore, in order to obtain the correct symmetries, we need to explicitly make the atoms different. To do this, we defined two species, Ni1 and Ni2, which are both Ni atoms, but having different magnetic moments.

• KPointsUseTimeReversal = no: Although we obviously want to use the ‘‘k’’-point symmetries in our calculation, we need to deactivate the use of time-reversal symmetry, as this is currently not implemented when doing a spin-polarized calculation.

• GuessMagnetDensity and AtomsMagnetDirection: An initial guess is added to break the system symmetries and help the convergence.

We also added a few extra state to the calculation. While this is not needed to obtain the ground state, it is interesting to do so in order to get the value of the electronic bandgap, which we might want to compare to the experiment. Note that depending on the system degeneracy, one might need to add more than one extra state to get an estimate of the bandgap.

Output

After running Octopus, we can have a look at the output and notice a few things. The information about the symmetries of the system is printed at the start of the calculation:

***************************** Symmetries *****************************
Space group No.  166
International: R-3m
Schoenflies: D3d^5
  Index                Rotation matrix                      Fractional translations
    1 :     1   0   0     0   1   0     0   0   1      0.000000    0.000000    0.000000
    2 :    -1   0   0     0  -1   0     0   0  -1      0.000000    0.000000    0.000000
    3 :    -1   1   0    -1   0   0    -2   0   1      0.000000    0.000000    0.000000
    4 :     1  -1   0     1   0   0     2   0  -1      0.000000    0.000000    0.000000
    5 :     0  -1   0     1  -1   0     0  -2   1      0.000000    0.000000    0.000000
    6 :     0   1   0    -1   1   0     0   2  -1      0.000000    0.000000    0.000000
    7 :     1   0   0     1  -1   0     2   0  -1      0.000000    0.000000    0.000000
    8 :    -1   0   0    -1   1   0    -2   0   1      0.000000    0.000000    0.000000
    9 :     0  -1   0    -1   0   0     0   0  -1      0.000000    0.000000    0.000000
   10 :     0   1   0     1   0   0     0   0   1      0.000000    0.000000    0.000000
   11 :    -1   1   0     0   1   0     0   2  -1      0.000000    0.000000    0.000000
   12 :     1  -1   0     0  -1   0     0  -2   1      0.000000    0.000000    0.000000
Info: The system has    12 symmetries that can be used.
**********************************************************************

Here we see that the symmetry finder found a different spacegroup for our cubic supercell than the one of cubic NiO where all Ni atoms are considered equivalent (i.e. spacegroup 225).

As the calculation is a spin-polarized calculation, during the SCF cycle Octopus outputs the total magnetic moment as well as the local magnetic moments around the atoms, obtained by intergrating the magnetization density on a sphere centered around each atom. For the last SCF iteration, it should look like this:

Total Magnetic Moment:
 mz =  -0.000000
Local Magnetic Moments (sphere radius [b] =   1.970):
 Ion                    mz
   1       Ni1       1.501271
   2       Ni2      -1.501272
   3         O       0.000001
   4         O       0.000001

As we can see, we indeed obtained an antiferromagnetic state, where the total magnetization is zero and the local magnetic moments are non-zero and of opposite signs.

Finally, let’s have a look at the values of the direct and indirect bandgaps, which can be found in the static/info file:

Direct gap at ik=    2 of  0.0304 H
Indirect gap between ik=    8 and ik=    2 of  0.0177 H

It will be instructive to compare these values with the ones obtained with DFT+U.

Bulk NiO with DFT+U

Input

We will now run the same system, but with a Hubbard U correction. For this we modify our previous input file so that it looks like this:

 CalculationMode = gs
 PeriodicDimensions = 3
 BoxShape = parallelepiped
 FromScratch = yes

 a = 7.8809
 %LatticeParameters
   a | a | a
 %
 %LatticeVectors
  0.0 | 1/2 | 1/2
  1/2 | 0.0 | 1/2
  1.0 | 1.0 | 0.0
 %

 %Species
 "Ni1" | species_pseudo | set | hscv_pbe | hubbard_l | 2 | hubbard_u | 5.0*eV
 "Ni2" | species_pseudo | set | hscv_pbe | hubbard_l | 2 | hubbard_u | 5.0*eV
 %

 DFTULevel = dft_u_empirical

 %ReducedCoordinates
  "Ni1" | 0.0 | 0.0 | 0.0
  "Ni2" | 0.0 | 0.0 | 0.5
  "O"  | 0.5 | 0.5 | 0.25
  "O"  | 0.5 | 0.5 | 0.75
 %
 Spacing = 0.4
 %KPointsGrid
 2 | 2 | 2
 %
 KPointsUseSymmetries = yes
 KPointsUseTimeReversal = no

 ExtraStates = 3

 SpinComponents = polarized
 GuessMagnetDensity = user_defined
 %AtomsMagnetDirection
  8.0
 -8.0
  0.0
  0.0
 %

 %Output
   occ_matrices
 %

The main differences compared to the previous input file are the following ones:

• DFTULevel = dft_u_empirical: This variable specifies the level of DFT+U used. Here we have chosen to use an empirical correction. The other two available options are dft_u_none, which is the default and corresponds to no +U correction, and dft_u_acbn0, which corresponds to the ab initio U correction based on the ACBN0 functional².

• Species: As we have chosen to use an empirical correction, we need to specify the value of the effective Hubbard U and the orbitals this will be aplied to. This is done by adding two options to the Species block: hubbard_l specifies the orbitals (l=0 for s orbitals, l=1 for p orbitals, …) and hubbard_u is used to set the value of the effective Hubbard U. Here we add a Hubbard U of 5eV on the 3d orbitals (corresponding to the quantum number l=2).

• Output: We ask the code to output the density matrix of the selected localized subspaces, usually called occupation matrix. This will create a file named static/occ_matrices .

Output

After running Octopus using the above input, we can look at the output. Some information specific to DFT+U is printed at the start of the calculation:

******************************* DFT+U ********************************
Input: [DFTUDoubleCounting = dft_u_fll]
Input: [DFTUPoissonSolver = dft_u_poisson_fft]

Method:
  [1] Dudarev et al., Phys. Rev. B 57, 1505 (1998)
Implementation:
  [1] Tancogne-Dejean, Oliveira, and Rubio, Phys. Rev. B 69, 245133 (2017)

Input: [AOTruncation = ao_full]
Input: [AOThreshold = 0.1000E-01]
Input: [AONormalize = yes]
Input: [AOSubmesh = no]
Input: [AOLoewdin = no]
Building the LDA+U localized orbital basis.
Found   2 orbital sets.
Orbital set  1 has a value of U of  0.18375 Ha.
It contains  5 orbitals.
The radius is  4.40000 Bohr,  with   8007 grid points.
Orbital set  2 has a value of U of  0.18375 Ha.
It contains  5 orbitals.
The radius is  4.40000 Bohr,  with   8007 grid points.
**********************************************************************

Among other details, this output confirms that we are indeed using 2 sets of d orbitals, each corresponding to the orbitals centered around a Ni atom.

Lets now look at the magnetic moments:

Total Magnetic Moment:
 mz =   0.000000
Local Magnetic Moments (sphere radius [b] =   1.970):
 Ion                    mz
   1       Ni1       1.698444
   2       Ni2      -1.698445
   3         O       0.000000
   4         O       0.000000

and bandgaps:

Direct gap at ik=    2 of  0.0828 H
Indirect gap between ik=    5 and ik=    2 of  0.0621 H

We see that the magnetic moments on the Ni atoms are larger than at the PBE level and that the bandgap is also found to be larger. Clearly, by adding an Hubbard U correction, we have increased the electron localization in the d orbitals, which leads to an increase of the local magnetic moment and to the opening of the electronic bandgap.

Double counting term

The DFT+U method, like DFT+DMFT, needs to specify an approximate treatment of the so called double-counting term, i.e., the part of the electron-electron interaction that is already included in the DFT part and would be double counted if this would not be considered. By default, Octopus uses the fully-localized limit (FLL) double counting. It is also possible to use the around mean-field (AMF) double counting, as well as a “mixed” double counting term proposed in Ref.³. This is controlled by the variable DFTUDoubleCounting.

Choice of basis

In this tutorial, we have seen how to use DFT+U with atomic-center orbitals obtained from the pseudopotential. It is also possible to construct a localization subspace from states defined on the mesh, obtained from a different calculation. For instance, one can perform an DFT calculation to obtain a localized state (core state, flat band, …). Then is it possible to perform a DFT+U calculation based on these state, using the variable called DFTUBasisFromStates, together with DFTUBasisStates. It is also possible to use oct-wannier90 to obtain the Wannier state and then to use it for the localized subspace. This method will be presented in a different tutorial.

References