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
Messages in the output related to LDA+U
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