Calculations with LDA+U

This handson will cover the application of the DFT+U method, which can help to describe strongly correlated materials more properly. We will cover two examples here. NiO, where we improve the the description of the 3d electrons, and bulk Gd, where the description of the 4f orbitals is improved.

NiO

NiO has a very small bandgap (0.4eV) and too small magnetic moment (1.2 Bohr) in a LDA calculation. The LDA+U method cures these problems.

In the following we will calculate the antiferromagnetic structure (AFM1) of this material.

export THOME=$PWD #Save tutorial home dir
cd NiO; mkdir GGA_U; cd GGA_U

We create multiple folders for calculations here, since we want to compare results for GGA and GGA+U

First we create a inp.xml for a GGA calculation (GGA also underestimates gap and moments):

cp ../inp_NiO_AFM .; cat inp_NiO_AFM

Here we just call the inpgen

inpgen -f inp_NiO_AFM

GGA+U calculation

Now we could try to converge the GGA calculation, but this is rather unstable. Therefore, we directly turn on GGA+U.

Insert the parameters U and J (taken from literature) in the produced inp.xml for both Ni species below the energy parameters and above the line for the LO's:

<ldaU l="2" U="8.0" J="0.9" l_amf="F"/>

It also makes sense to raise the maximum number of iterations (itmax) to 25.

Open an editor and modify the inp.xml file as sketched above.

The current working directory is:

pwd

Now we can just start the fleur calculation adn converge the system

Repeat the call to fleur_MPI until the charge density is converged

fleur_MPI

In the first iteration there will be two messages that are normal. Since we did not specify an initial guess for the density matrix, in the first iteration no LDA+U term is added and we get the message:

 no density matrix found ... skipping LDA+U

The second message is:

Reset of history

This means, that fleur removes any previous mixing history for the Anderson mixing scheme, that is used in fleur by default. This is done since introducing an LDA+U term into the Hamiltonian will inevitably lead to a large distance jump which would hinder convergence if the previous history stayed around.

The DFT+U calculations will produce some additional output. For example if we look at the distances of the charge densities, we will find additional output related to the density matrix (lines starting with n_mmp)

grep "n_mmp distance" out
grep "bandgap" out | tail -1
#You might have to zoom out (ctrl + -) to have a nice looking table
grep -A8 'Magnetic moments |' out | tail -9

We should see that the bandgap is 3.3 eV/0.122 htr and the magnetic moment of the Ni atoms is 1.8 Bohr magnetons

ls

We also find a new file called n_mmp_mat_out in the directory of the calculation. This contains the current density matrix for DFT+U in the following format.

The file contains the density matrix for all LDA+U procedures in 7x7 complex matrices. Since the ouput format is 7 numbers per line, this means that every 2 rows corresponds to a row in the density matrix. The blocks are written out for each orbital and spin in the following order:

- orbital 1, spin-up
- orbital 2, spin-up
...
- orbital 1, spin-down
- orbital 2, spin-down

The density matrix can also be found in the out.xml, if you search for ldaUDensityMatrix, but the format of the n_mmp_mat_out file can be used to provide an initial guess for the density matrix in Fleur. Namely, if the file is renamed to n_mmp_mat it will be read in by Fleur in the next calculation and will initialize the density matrix

See Notebook Appendix- Creating initial density matrices.ipynb for some helpful tools for generating such a file from scratch if you are interested.

cat n_mmp_mat_out

DOS calculation

Now we want to calculate the DOS for this system and plot it. We use the same folder for this.

Set dos to T and a sigma of 0.005 works well for this system

Open an editor and modify the inp.xml file as sketched above.

fleur_MPI

There is a prepared plotting script for this system using the masci-tools library, which we will use here. It will produce the DOS with the filename dos_plot.png in the current working directory.

There should be a clearly visible bandgap of around 3.3 eV

cat $THOME/NiO/plot_dos.py
python $THOME/NiO/plot_dos.py

GGA calculation

Now we want to obtain the pure GGA values. First we create a new directory with all the required files:

We create a new input file and just take the GGA+U charge density as our starting point

mkdir $THOME/NiO/GGA; 
cp cdn.hdf inp_NiO_AFM ../GGA; 
cd $THOME/NiO/GGA

Now we invoke the input generator and start fleur. This calculation will converge much slower than the GGA+U calculation

inpgen -f inp_NiO_AFM

Repeat the call to fleur_MPI until the charge density is converged

fleur_MPI
grep "bandgap" out | tail -1
#You might have to zoom out (ctrl + -) to have a nice looking table
grep -A8 'Magnetic moments |' out | tail -9

Now we get a vanishing bandgap and a magnetic moment of 1.4 Bohr magnetons.

DOS calculation

We repeat the exact same procedure to calculate the density of states for this system

Set dos to T and a sigma of 0.005 works well for this system

Open an editor and modify the inp.xml file as sketched above.

fleur_MPI
python $THOME/NiO/plot_dos.py

We should see that the bandgap disappears and the DOS shows the characteristics of a metal

GGA+U calculation with AMF double counting

We want to calculate a third scenario for this system by calculating it with GGA+U and the around mean field limit for the double counting. For this calculation we start from scratch again.

cd $THOME/NiO;
mkdir GGA_U_AMF;
cp inp_NiO_AFM GGA_U_AMF;
cd GGA_U_AMF
inpgen -f inp_NiO_AFM

To turn on the around mean-field double counting we just set the l_amf switch to T in the ldaU tags that are added to the inp.xml

Insert the parameters U and J (taken from literature) in the produced inp.xml for both Ni species below the energy parameters and above the line for the LO's:

<ldaU l="2" U="8.0" J="0.9" l_amf="T"/>

It also makes sense to raise the maximum number of iterations (itmax) to 25.

Open an editor and modify the inp.xml file as sketched above.

The current working directory is:

pwd

Repeat the call to fleur_MPI until the charge density is converged

fleur_MPI
grep "bandgap" out | tail -1
#You might have to zoom out (ctrl + -) to have a nice looking table
grep -A8 'Magnetic moments |' out | tail -9

Now we should get a smaller bandgap of 2.3 eV and a magnetic moment of 1.8 Bohr magnetons.

DOS calculation

We repeat the exact same procedure to calculate the density of states for this system

Set dos to T and a sigma of 0.005 works well for this system

Open an editor and modify the inp.xml file as sketched above.

fleur_MPI
python $THOME/NiO/plot_dos.py