# ABINIT, tutorial on LDA+U :

## The projected density of states of NiO

This lesson aims at showing how to perform a LDA+U calculation using Abinit.

You will learn what is a LDA+U calculation and what are the main input variables controlling this type of calculation.

It is supposed that you already know how to do PAW calculations using ABINIT. Please follow the two lessons on PAW in ABINIT (PAW1, PAW2), if this is not the case.

This lesson should take about 1 hour to complete.

##### Copyright (C) 2008-2009 ABINIT group (SPesant,MCote,XG)
This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

Goto : ABINIT home Page | Suggested acknowledgments | List of input variables | Tutorial home page | Bibliography

Help files : New user's guide | Abinis (main) | Abinis (respfn) | Mrgddb | Anaddb | AIM (Bader) | Cut3D | Optic | Mrgscr

### Content of the LDA+U lesson

- 0. Short summary of the LDA+U method
- 1. Ground state calculation of NiO using LDA
- 2. Ground state calculation of NiO using LDA+U with the FLL double-counting method.
- 3. Initialization of the density matrix in LDA+U for the NiO crystal
- 4. Ground state calculation of NiO using LDA+U with the AMF double-counting method.
- 5. Projected density of states in LDA+U

### 0. Short summary of the LDA+U method

The standard Local Density Approximation (LDA), where the exchange and correlation energy is fit to homogenous electron gas results, is a functional that works well for a vast number of compounds. But, for some crystals, the interactions between electrons are so important that they cannot be represented by the LDA alone. Generally, these highly correlated materials contain rare-earth metals or transition metals, which have partially filled d or f bands and thus localized electrons.

The LDA tends to delocalize electrons over the crystal, and each electron feels an average of the Coulombic potential. For highly correlated materials, the large Coulombic repulsion between localized electrons might not be well represented by a functional such as the LDA. A way to avoid this problem is to add a Hubbard-like, localised term, to the LDA density functional. This approach is known as LDA+U. In the actual implementation, we separate localized d or f electrons, on which the Hubbard term will act, from the delocalized ones (s and p electrons). The latter are correctly described by the usual LDA calculation. In order to avoid the double counting of the correlation part for localized electrons (already included in the LDA, although in an average manner), another term - called the double-counting correction - is subtracted from the Hamiltonian.

In Abinit, two double-counting corrections are currently implemented :

-The Full localized limit (FLL) (see A. Lichtenstein et al PRB 52, 5467 (1995))

-The Around Mean Field (AMF) (see M. T. Czyzyk and G. A. Sawatzky PRB 49, 14211 (1994))

For some systems, the result might depend on the choice of the double-counting method. However, the two methods generally give similar results.

### 1. Ground state calculation of NiO using LDA

You might create a subdirectory of the ~abinit/tests/tutorial directory, and use it for the tutorial. In what follows, the names of files will be mentioned as if you were in this subdirectory

Copy the files ../Input/tldau_1.in and ../Input/tldau_1.files in your Work directory and run

This run should take around 7 minutes on a PowerPC 2.5 GHz G5. It calculates the LDA ground state of the NiO crystal. A low cutoff and a small number of k-points are used in order to speed up the calculation. During this time you can take a look at the input file.

The NiO crystallizes in the rocksalt structure, with one Ni and one O atom in the primitive cell (the crystallographic primitive cell). However, NiO is known to exhibit an antiferromagnetic ordering at low temperature (along the <111> direction). From the electronic point of view, the true unit cell has two Ni and two O atoms : the local magnetic moment around the first Ni atom will have a sign opposite to the one of the other Ni atom.

You should take some time to examine the values used for the input variables xred, rprim (note the last line !), typat, spinat, nsppol, and nspden, that define this antiferromagnetic ordering along the <111> direction (of a conventional cubic cell).

If you take a look at the output file (tldau_1.out), you can see the integrated total density in the PAW spheres (see the PAW1 and PAW2 tutorials on PAW formalism). This value roughly estimate the magnetic moment of NiO :

     Integrated total density in atomic spheres:
     -------------------------------------------
     Atom  Sphere radius  Integrated_up_density  Integrated_dn_density  Total(up+dn)   Diff(up-dn)
        1        2.31080             5.09436023             3.90700183    9.00136206    1.18735839
        2        2.31080             3.90700183             5.09436023    9.00136206   -1.18735839
        3        2.31080             3.44707732             3.44707732    6.89415463    0.00000000
        4        2.31080             3.44707732             3.44707732    6.89415463    0.00000000
     Note: Diff(up-dn) can be considered as a rough approximation of a local magnetic moment.

The atoms in the output file, are listed as in the typat variable (the first two are nickel atoms and the last two are oxygen atoms). The results indicate that spins are located in each nickel atom of the doubled primitive cell. Fortunately, the LDA succeeds to give an antiferromagnetic ground state for the NiO. But the result does not agree with the experimental data. The magnetic moment (the difference between up and down spin on the nickel atom) range around 1.6-1.9 according to experiments (A. K. Cheetham and D. A. O. Hope, Phys. Rev. B. 27, 6964 (1983), H. A. Alperin, J. Phys. Soc. Jpn. 17, 12 (1962), W. Neubeck et al., J. Appl. Phys. 85, 4847 (1999), G. A. Sawatzky and J. W. Allen, Phys. Rev. Lett. 53, 2339 (1984) and S. Hufner et al., Solid State Comm. 52, 793 (1984) ). Also, the band gap obtained between the last occupied (0.29691 Ha, at k point 2) and the first unoccupied band (0.32453 Ha, at kpoint 10) is approximately 0.8 eV which is lower than the measured value of 4.0-4.3 eV (This value could be modified using well-converged parameters but would still be much lower than what is expected).

Making abstraction of the effect of insufficiently convergence parameters, the reason for the discrepancy between the DFT-LDA data and the experiments is first the fact the DFT is a theory for the ground state and second, the lack of correlation of the LDA. Alone, the homogenous electron gas cannot correctly represent the interactions among d electrons of the Ni atom. That is why we want to improve our functional, and be able to manage the strong correlation in NiO.

### 2. Ground state calculation of NiO using LDA+U with the FLL double-counting

As seen previously, the LDA does not gives good results for the magnetization and band gap compared to experiments.

At this stage, we will try to improve the correspondence between calculation and experimental data. First, we will use the LDA+U with the Full localized limit (FLL) double-counting method.

FLL and AMF double-counting expressions are given in the papers listed above, and use the adequate number of electrons for each spin. For the Hubbard term, the rotationally invariant interaction is used.

It is important to notice that in order to use LDA+U in Abinit, you must employ PAW pseudopotentials.

You should run abinit on the input file tldau_2.in. This calculation takes approximately 9 minutes on a PowerPC 2.5 GHz G5.

During the calculation, you can take a look at the input file. Some variable describing the LDA+U parameters have been added to the previous file. All other parameters were kept constant from the preceeding calculation. First, you must set the variable usepawu to one (for the FLL method) and two (for the AMT method) in order to enable the LDA+U calculation. Then, with lpawu, you give for each atomic species (znucl) the values of angular momentum (l) for which the LDA+U correction will be applied. The choices are 2 for d-orbitals and 3 for f-orbitals. You cannot treat s and p orbitals with LDA+U in the present version of ABINIT. Also, if you do not want to apply LDA+U correction on a species, you can set the variable to -1.For the case of NiO, we put lpawu to 2 for Ni and -1 for O.

Finally, as described in the article cited above for FLL and AMF, we must define the screened Coulomb interaction between electrons that are treated in LDA+U, with the help of the variable upawu, and the screened exchange interaction, with jpawu. Note that you can choose the energy unit by indicating at the end of the line the unit abbreviation (e.g. eV or Ha). For NiO, we will use variables that are generally accepted for this type of compound:

upawu = 8.0 eV

jpawu = 0.8 eV (10 % of U)

You can take a look at the result of the calculation. The magnetic moment is now :

     Integrated total density in atomic spheres:
     -------------------------------------------
     Atom  Sphere radius  Integrated_up_density  Integrated_dn_density  Total(up+dn)   Diff(up-dn)
        1        2.31080             4.53843823             4.53872511    9.07716334   -0.00028688
        2        2.31080             4.53872511             4.53843823    9.07716334    0.00028688
        3        2.31080             3.42606450             3.42606450    6.85212900    0.00000000
        4        2.31080             3.42606450             3.42606450    6.85212900    0.00000000
     Note: Diff(up-dn) can be considered as a rough approximation of a local magnetic moment.

And the system is an low band gap insulator (the 16th band at k point 6 has an eigenenergy of 0.30415 Ha, that is a little bit lower than the eigenenergy of the 17th band at k point 10, namely 0.30603 ).

As you can see, the magnetic moment has not been improved compared to experimental results. In fact, it is even worse. This result highlights the problem with convergency which can occur in LDA+U, and in magnetic system in general. In this particular case, this problem comes from the fact that poor convergency parameters have been used (ecut, and kpoints). With a higher energy cut-off, and a finer kpoint grid, it is possible to obtain a band gap of 3.3 eV and a reasonable magnetic moment (this calcution takes approximately 6 hours) :

     Integrated total density in atomic spheres:
     -------------------------------------------
     Atom  Sphere radius  Integrated_up_density  Integrated_dn_density  Total(up+dn)   Diff(up-dn)
        1        2.31080             5.30568572             3.59317402    8.89885974    1.71251170
        2        2.31080             3.59317402             5.30568572    8.89885974   -1.71251170
        3        2.31080             3.48210428             3.48210428    6.96420857    0.00000000
        4        2.31080             3.48210428             3.48210428    6.96420857    0.00000000
     Note: Diff(up-dn) can be considered as a rough approximation of a local magnetic moment.

As we can see, this value of magnetization is in good agreement with experimental values. Thiscalculation also reveals another difficulty in LDA+U : it can be delicate to find the good ground state. For NiO, ABINIT does not succeeded to find the lowest energy state, which is antiferromagnetic, for the insufficient convergence parameters that we used. (For this tutorial, better parameters would have taken too much computational time, which is the reason for using unconverged parameters).

### 3. Initialization of the density matrix in LDA+U for the NiO crystal

You should begin by running the tldau_3.in file before continuing.

In order to help the LDA+U find the ground state, you can define the initial density matrix for correlated orbitals with dmatpawu. To enable this feature, usedmatu must be set to a non-zero value (default is 0). When positive, the density matrix is kept to the dmatpawu value for the usedmatpu value steps. For our calculation(tldau_3.in) , usedmatpu is 5, thus the density matrix is kept constant for 5 SCF steps.

In the log file (not the usual output file), you might find for each step, the calculated density matrix, followed by the imposed density matrix. After the first 5 SCF steps, the initial density matrix is no longer imposed.

    ========== LDA+U DATA ===================================================

    ====== For Atom    2, occupations for correlated orbitals. lpawu =   2

    Atom   2. Occ. for lpawu and for spin  1 =   3.91002
    Atom   2. Occ. for lpawu and for spin  2 =   4.77891
    => On atom   2,  local Mag. for lpawu is      0.86889325E+00

    == Calculated occupation matrix for correlated orbitals:

    Calculated occupation matrix for spin  1
         0.78200    0.00000    0.00000    0.00000    0.00000
         0.00000    0.78200    0.00000    0.00000    0.00000
         0.00000    0.00000    0.78200    0.00000    0.00000
         0.00000    0.00000    0.00000    0.78200    0.00000
         0.00000    0.00000    0.00000    0.00000    0.78200

    Calculated occupation matrix for spin  2
         0.95578    0.00000    0.00000    0.00000    0.00000
         0.00000    0.95578    0.00000    0.00000    0.00000
         0.00000    0.00000    0.95578    0.00000    0.00000
         0.00000    0.00000    0.00000    0.95578    0.00000
         0.00000    0.00000    0.00000    0.00000    0.95578

    ======= LDA+U Energy terms (in Hartree) ====

     For Atom   2
       Contributions to the direct expression of energy:
         Double counting  correction   =    0.9387965491E+01
         Interaction energy            =    0.9528687788E+01
         Total LDA+U Contribution      =    0.1407222972E+00

       For the "Double-counting" decomposition:
         LDA+U Contribution            =    0.1008802376E+01

    -------------------------------------------------------------------------

    Occupation matrix for correlated orbitals is kept constant
    and equal to initial one !
    ----------------------------------------------------------

    == Atom   1 == Imposed occupation matrix for spin 1 ==
         0.90036    0.00000   -0.00003    0.00000    0.00000
         0.00000    0.90036   -0.00001    0.00000    0.00002
        -0.00003   -0.00001    0.91309   -0.00001    0.00000
         0.00000    0.00000   -0.00001    0.90036   -0.00002
         0.00000    0.00002    0.00000   -0.00002    0.91309

    == Atom   1 == Imposed occupation matrix for spin 2 ==
         0.89677   -0.00001    0.00011   -0.00001    0.00000
        -0.00001    0.89677    0.00006    0.00001   -0.00010
         0.00011    0.00006    0.11580    0.00006    0.00000
        -0.00001    0.00001    0.00006    0.89677    0.00010
         0.00000   -0.00010    0.00000    0.00010    0.11580

    == Atom   2 == Imposed occupation matrix for spin 1 ==
         0.89677   -0.00001    0.00011   -0.00001    0.00000
        -0.00001    0.89677    0.00006    0.00001   -0.00010
         0.00011    0.00006    0.11580    0.00006    0.00000
        -0.00001    0.00001    0.00006    0.89677    0.00010
         0.00000   -0.00010    0.00000    0.00010    0.11580

    == Atom   2 == Imposed occupation matrix for spin 2 ==
         0.90036    0.00000   -0.00003    0.00000    0.00000
         0.00000    0.90036   -0.00001    0.00000    0.00002
        -0.00003   -0.00001    0.91309   -0.00001    0.00000
         0.00000    0.00000   -0.00001    0.90036   -0.00002
         0.00000    0.00002    0.00000   -0.00002    0.91309

This time, by imposing an initial density matrix to the correlated orbitals, the results are in better agreement with experiments. In fact, the magnetization is 1.713 and the band gap is 4.7 eV. Also, the total energy is lower than the previous calculation by approximately 0.16 Ha. This confirms that the calculation was not converged. Moreover, one could show that the antiferromagnetic calculation is the most stable one in LDA+U.

Generally, the LDA+U functional meets the problem of multiple local minima, much more than the usual LDA or GGA functionals. One often gets trapped in a local minimum. Trying different starting points might be important...

### 4. Ground state calculation of NiO using LDA+U with the AMF double-counting method.

Now we will use the other implementation for the double-counting term in LDA+U (in Abinit), known as AMF. As the FLL method, this method uses the number of electrons for each spin independently and the complete interactions U(m1,m2,m3,m4) and J(m1,m2,m3,m4).

As in the preceding run, we will start with a fixed density matrix for d orbitals. You might now Start your calculation, with the tldau_4.in and tldau_4.files, or skip the calculation, and rely on the reference file provided in the ~abinit/tests/tutorial/Refs directory. Examine the tldau_4.in file. The only difference in the input file compared to tldau_3.in is the value of usepawu = 2. We obtain a band gap of 4.1 eV. The value of the band gap with AMF and FLL is different. However, we have to remember that these results are not well converged. By contrast, the magnetization,

     Integrated total density in atomic spheres:
     -------------------------------------------
     Atom  Sphere radius  Integrated_up_density  Integrated_dn_density  Total(up+dn)   Diff(up-dn)
        1        2.31080             5.28943941             3.60071406    8.89015347    1.68872535
        2        2.31080             3.60071406             5.28943941    8.89015347   -1.68872535
        3        2.31080             3.49862376             3.49862376    6.99724752    0.00000000
        4        2.31080             3.49862376             3.49862376    6.99724752    0.00000000
     Note: Diff(up-dn) can be considered as a rough approximation of a local magnetic moment.

is very similar to the LDA+U FLL. In fact, this system is not very complex. But for other systems, the difference can be more important. In fact, FLL is designed to work well for crystal with diagonal occupation matrix with 0 or 1 for each spin. The AMF should be used when orbital occupations are near the average occupancies.

### 5. Projected density of states in LDA+U

Using prtdos 3, you can now compute the projected d and f density of states. For more information about projected density of states, for more details see the PAW1 tutorial.

Goto : ABINIT home Page | Suggested acknowledgments | List of input variables | Tutorial home page | Bibliography

Help files : New user's guide | Abinis (main) | Abinis (respfn) | Mrgddb | Anaddb | AIM (Bader) | Cut3D | Optic | Mrgscr