Introduction

There is a jungle of DFT tutorials all over the web, and only some of them are good. Only few of them are oriented at DMFT calculations, and even less use w2dynamics. This tutorial is written in the hope of filling this gap, combining a little bit of theory with practical how-tos and input files. It will be another DFT tutorials over the web. Hopefully among the good ones.

In this tutorial we will compute the DMFT spectral function of NiO, which is a charge transfer insulator. There is a similar tutorial on the triqs documentation which uses VASP + TRIQS. Instead, we will use Quantum ESPRESSO + w2dynamics. Compared to other DFT codes Quantum ESPRESSO has the advantage of being free and open source. Other than that, all DFT codes are essentially equivalent.

The amount of theoretical and practical knowledge needed for this tutorial is large, feel free to skip to the part that interests you through the table of contents.

Installation

Quantum ESPRESSO

can be installed in two ways, assuming you have superuser privileges. If you don't, it likely means you're using a supercluster, then you'd better bow to the cluster service and ask that they install it for you.

Easy way:

sudo apt-get install quantum-espresso

General installation:

Install the minimum needed libraries:
sudo apt-get install libopenmpi-dev -y sudo apt-get install libopenblas-dev -y sudo apt-get install libfftw3-dev -y sudo apt-get install liblapack-dev -y sudo apt-get install libscalapack-openmpi-dev -y sudo apt-get install g++ -y
Download from the download page (7.1 version)
Unpack with tar -xvf qe\-ReleasePack.tar
cd to qe-7.1
Run the ./configure script
Run sudo make all
Run sudo make install
If problems arise you have to go through the complaints of make to figure out what's wrong. Fortunately Quantum ESPRESSO has many users so you often find someone else who had the same problem.

VESTA

Can be downloaded from the jp minerals website. Installation is straightforward for all standard operative systems.

Density Functional Theory

Setting up a basic DFT calculation

General workflow of a DFT calculation

In general, a DFT calculation is well done when the total electronic energy reaches its minimum, which by the second Hohenberg and Kohn's theorem also means that the charge density is the true ground-state charge density of the system.

This is achieved by solving the Schrödinger's equation for the Kohn and Sham single-particle hamiltonian in a self consistent loop, which works as follows:

Start from an initial guess. Typically a wavefunction that is the sum of atomic wavefunctions.
Compute the charge density from the wavefunctions
Compute the potential terms of the hamiltonian from the charge density
Solve the single-particle hamiltonian with the new potential to obtain new eigenvalues and eigenfunctions
Go back to a until self-consistency is achieved (typically on total energy)

This is called the self-consistent loop. Quantum ESPRESSO is a modular code. The part of the code that does the self-consistency is called pw.x

Once self-consistency is reached one obtains the ground-state charge density, which can be re-used to compute everything else. For instance, to compute the band structure one first obtains a good ground-state charge density, and then does a new non-self-consistent calculation in which the Scrhödinger's equation for the Kohn and Sham hamiltonian is constructed using the charge density of the previous step.

Create and visualize the crystal structure using VESTA

The first thing you need is information on your crystal structure. It's a CsCl-type cubic crystal with a = 4.1674 Å. We can use this information to visualize the primitive cell.

Open VESTA, then click: File → New Structure, a prompt will open
Click Unit cell and select Cubic, Space Group No. 225 (Fm3m), and enter the lattice parameter

Now go to Structure parameters → New and add a Ni atom in (0, 0, 0) and O in (1/2, 1/2, 1/2). You should see something like this. Notice that VESTA does the work of recognizing the Wyckoff positions for us.

VESTA is a very powerful software. It can plot isosurfaces, lattice planes, polyhedral bonds.. it would need a separate tutorial.

Choose a pseudopotential

What is a pseudopotential and why you need it

The Kohn and Sham equation is a second-order Schrödinger's equation for a single electron in an external field. Any numerical solution first goes through an expansion into a basis set, which transforms the problem of solving the differential equation into the one of finding the coefficients of the expansion. The latter is a problem of linear algebra and can be done numerically, provided that the expansion is truncated at some finite value. Many DFT codes, including Quantum ESPRESSO, use plane waves as the basis for expansion.

Plane waves have mainly two advantages:

They are orthogonal to each other, and they are a complete set
Truncation into a finite set does not break Hellman-Feynman's theorem

They also have a major inconvenience: in order to accurately describe an oscillation in real space of the order of a unit length, one needs plane waves which have at least that wavelength, or smaller. All the atomic wavefunctions present fast oscillations in real space, even those of the valence shell which are not per se localized. If the DFT calculation took into account all electrons plane waves would be simply too expensive.

Pseudopotentials overcome this problem by removing the core electrons, and replacing the core region of the potential with a weakly repulsive part that mimicks the effect of the core electrons, in such a way that the electronic properties of the valence electrons are unchanged.

Where do I get good pseudopotentials and how do I choose

Choosing a pseudopotential for a DFT calculation while being fully aware what one is doing requires quite a lot of experience. I recommend to anyone without experience to simply use a well-tested norm-conserving PBE pseudopotential. Anything else risks doing more harm than good.

A well-maintained and reliable datased is found at the pseudo-dojo website. The Optimized Norm-Conserving Vanderbilt (ONCV) pseudopotentials are a little costly, but typically very reliable.

Quantum ESPRESSO uses the Unified Pseudopotential Format (UPF) format.

Set up of the input file to compute the ground-state charge density

Download the pseudopotentials used in this tutorial

Ni_ONCV_PBE_sr.upf

O_ONCV_PBE_sr.upf

Convergence criteria to watch out for: plane waves expansion and Brillouin zone integration

Plane waves expansion

In the section about pseudopotentials it was mentioned that in Quantum ESPRESSO the wavefunctions are expanded in plane waves, and this expansion is truncated at some finite value. In most DFT codes the truncation is defined in terms of an energy value such that

Where is a reciprocal lattice vector.

Using a large enough value of is necessary to get converged result. How large depends very much on which pseudopotential is used and a very little on the specific crystal structure.

As we add more plane waves the total energy will asymptotically approach the true value from above. It looks somewhat like this.

In an accurate calculation typically one wants to have within 1 meV/atom, i.e. Anything above 70 Ry is fine for NiO. You can use the following files to produce the plot above.

CONVERGE_ECUT.sh: Bash script that changes the cutoff and runs pw.x

NiO.src.in: Template input file

plot_convergence.py: Plotting script

Brillouin zone integration

The Bloch Hamiltonian is diagonal in k-space, hence the Kohn-Sham equations can be solved separately for every k point. However, the charge density (in reciprocal space) depends on the integral over the whole Brillouin zone, hence one defines a finite and discrete grid to integrate on, usually following the Monkhorst-Pack scheme. In most cases the specific scheme does not matter, what matters is if the k-mesh is dense enough so that the integrated charge is reasonable.

In insulators, integrals over the whole Brillouin zone typically converge quite quickly, because the bands over which the integrals are performed cannot change (you integrate up to the valence band).

In metals, things are more complicated because the integral is up to the Fermi energy, hence a tiny change in the electronic structure will change which bands are occupied and which are not, hence self-consistency becomes much more difficult to achieve. To avoid this circumstance, one introduces a smearing of the states, i.e. fractional occupation of the states is allowed, which is equivalent to introducing a finite temperature in the electronic states.

There are several smearing functions (Gaussian, Methfessel-Paxton, Marzari-Vanderbilt). Again, for most systems it does not matter which one is used, as long as the charge density is converged. All smearing functions have one parameter which is the width of the smearing function.

You should converge your calculation with respect to both smearing and k-mesh density.

Self-consistent calculation

Also called a scf calculation (as for self-consistent field), it is the bread and butter and the starting point of any other calculation.

It consists in starting from an initial guess of the charge density, which is used in the self-consistent loop until convergence is reached (i.e. the energy does not change up to a user-defined threshold .

The result of the self-consistent calculation provides the following useful quantities:

A converged ground-state charge density that we can re-use for subsequent calculations
The ground-state total electronic energy (or enthalpy) of the system
Thanks to the Hellman-Feynman theorem: any first derivative of the total energy. In particular: forces acting between the atoms and stresses acting on the unit cell employed
The Kohn-Sham eigenvalues for the k-points employed.

Note the naming convention: one scf calculation consists of several steps. When we refer to a scf calculation we intend that we repeat the steps above until the total energy is converged.

The calculation is run by typing:
pw.x < NiO.scf.in > NiO.scf.out

The arrows < and > redirect the input file to pw.x, and redirect the standard output of pw.x to an output file, respectively.

NiO.scf.in: input file for a scf calculation
NiO.scf.out: corresponding output file

Structural relaxation using DFT

A given input structure found in the literature or constructed by hand might not be exactly at equilibrium. The equilibrium positions also slightly change if you change exchange-correlation functional or if you apply an external pressure.

Quantum ESPRESSO can automatically relax a given crystal structure to the closest minimum, i.e. minimize forces between atoms and stresses.

The way in which this works is simply by performing one scf calculation (as in the previous paragraph), from which forces and stresses are obtained, then displace the atoms and change the lattice parameters accordingly.

This is repeated until the forces/stresses are within a desired threshold. The examples uses a very convenient flag:

cell_dofree = "ibrav"

Using this flag Quantum ESPRESSO keeps the Bravais lattice fixed, and directly outputs the lattice parameters (called celldm) in the standard output.

Note: a common mistake
After the structural relaxation the code will perform a fresh self-consistent calculation. However, you should always instead copy the result (i.e. the cell parameters and the atomic coordinates) into the input files, and perform a fresh self-consistent calculations with the relaxed settings, or you may end up with very confusing units.

The calculation is run by typing:
pw.x < NiO.relax.in > NiO.relax.out

NiO.relax.in: input file for a structural relaxation (both atoms and cell volume/shape)
NiO.relax.out: corresponding output file

DFT band structure

DFT is very often used to compute the so-called band structure, i.e. the dispersion relationship . For visualization purpose, the band structure is best plotted along high-symmetry lines in the Brillouin zone, that connect high-symmetry points. For instance, in the example provided NiO has a 225 space group. Its Brillouin zone can be seen for instance in the Bilbao crystallographic server. A good path in reciprocal space touches most of the high-symmetry lines.

In the example provided we follow the path:

What is really used to plot the DFT band structure are the Kohn-Sham eigenvalues, calculated at various wavevectors . Hence, we need to solve again the single-particle Schrödinger's equation for these points. However, these points do not form a uniform mesh over reciprocal space, hence we must use a ground-state charge density from a previous self-consistent calculation as input, and solve the Schrödinger's equation non-self consistently.

In practice, Quantum ESPRESSO does this for us. As long as we have run a scf calculation in the same folder, we can re-run the same code to do a different type of calculation (called bands) to obtain the band. Finally, a post-processing script called bands.x will read the result and combine it into a plot-friendly format.

In summary, to compute the bands we do the following steps:

A scf calculation (see above) using the pw.x code
A non-scf calculation (type bands) using the pw.x code
Post-processing using the bands.x code
Obtain the Fermi energy from the output file of the scf calculation

NiO.scf.in: scf input file (for pw.x)
NiO.bands.in: bands input file (for pw.x)
bands.in : post-processing input file (for bands.x)
NiO.bands.dat.gnu: expected output file (gnuplot-friendly)
plot_bands.ipynb: python plotting file (jupyter notebook)

The calculation is run by typing:
pw.x < NiO.scf.in > NiO.scf.out pw.x < NiO.bands.in > NiO.bands.out bands.x < bands.in > bands.out

Useful terminal commands:

Obtain the Fermi energy:grep 'Fermi' NiO.scf.out

Obtain the position of the special points in reciprocal space:
awk '{if(NR%30==1){print $1}}' NiO.bands.dat.gnu | head -6

In the end you should get something like the following (note: the Fermi energy is set as zero)

Important disclaimer: The whole world uses the Kohn-Sham eigenvalues as a good estimate of the real eigenenergies of the interacting system. In general this is not justified. The reason why people do it is that they have noticed that this often works quite well.

DFT density of states (with atom projection)

Density of states (DOS) is defined as:

Where is the band index and are wavevectors in the Brillouin zone. When you try to compute it numerically the sharpness of the function would require so sum over an impractically large number of points. Hence, you either replace it with a smoother function such as a Gaussian, or you use an alternative integration method which makes the sum finite (i.e. the tetrahedron method).

Even with a smarter integration method, the DOS requires a dense mesh to converge, much denser than what is required to converge the ground-state charge density. On the other hand, to compute the DOS we can use any converged charge density.

The common practice then is to first compute the charge density with a self-consistent calculation on a grid that is sufficient to converge it, and then re-start a non-self consistent calculation to obtain the eigenvalues on the dense grid.

An additional, very useful information, comes from the orbital projection of the DOS. In other words, we would like to break down the DOS into the contributions to single atomic orbitals. This is particularly useful when trying to understand which orbitals are relevant in a certain energy range. Different DFT codes use different schemes, each with pros and cons. Quantum ESPRESSO uses the atomic wavefunctions contained in the pseudopotential file.

In particular, suppose we have the atom with atomic wavefunction . Then we can write

Where is the Kohn-Sham eigenfunction with band index and wavevector , defined in the unit cell with the periodicity of the lattice. The code projwfc.x does this projection, which we need to do after the computation on the fine grid.

In summary, to compute the bands we do the following steps:

A scf calculation (see above) using the pw.x code on a grid such that the charge density is converged
A non-scf calculation (type nscf) using the pw.x code, on a dense grid such that the DOS is converged
Post-processing using the dos.x code
Post-processing using the projwfc.x code
Post-processing (optional) using the sumpdos.x code to sum certain projections

NiO.scf.in: scf input file (for pw.x)
NiO.nscf.in: nscf input file (for pw.x)
dos.in: dos input file (for dos.x)
projwfc.in: projection input file (for projwfc.x)
NiO.dos: DOS data file
NiO.Ni.dos.dat: DOS-Ni projection data file
NiO.O.dos.dat: DOS-O projection data file
plot_dos.ipynb: plotting script

The calculation is run by typing:
pw.x < NiO.scf.in > NiO.scf.out pw.x < NiO.nscf.in > NiO.nscf.out dos.x < dos.in > dos.out projwfc.x < projwfc.in > projwfc.outsumpdos.x NiO.pdos_*Ni* > NiO.Ni.dos.dat sumpdos.x NiO.pdos_*O* > NiO.O.dos.dat

Useful terminal commands:
grep 'electrons' NiO.scf.out

In the end you should get something like this

Wannierization

Wannier functions are a complete set of orthogonal functions. While plane waves are localized in reciprocal space, Wannier functions are localized in real space. The theory deserves a treatment in itself, for which I suggest the following resources:

As you may guess, Nicola Marzari is an important person for what concerns Wannier functions.

In short, the Wannier functions are a basis set. Going from the Bloch basis to the Wannier basis corresponds to a change of basis, and it is achieved through a unitary matrix. The unitary matrix is not uniquely defined, due to a Gauge freedom. This freedom can be exploited to find the set of Wannier functions which are most localized in real space. These are called Maximally Localized Wannier Functions (MLWF).

If we perform the transformation from the Bloch basis to the Wannier basis, we can also do the opposite, which is an excellent way of checking that the unitary matrix we are using to change basis is a good one. In addition, we can choose not to transform all the bands, but only a subset of them. This is often how Wannier functions are used to extract a minimal model from a DFT calculation.

Obtaining a good Wannierization for a small set of bands is an art which is hard to summarize. In general, do not expect to get a good Wannierization at first try. You will have to choose an energy window that includes all your relevant bands, but no more than that, then play with the initial atomic projections to start with a good guess, and then play with the so-called frozen window to help the code understand which bands it should absolutely keep. An excellent tutorial with tips and tricks can be found in this blog.

In summary, to obtain the Wannier functions we do the following steps:

A scf calculation (see above) using the pw.x code on a grid such that the charge density is converged
A non-scf calculation (type nscf) using the pw.x code, on a chosen grid. Starting with a 4 x 4 x 4 grid is a reasonable choice (but you may need to change it).
Pre-processing using wannier90.x -pp
Construction of the unitary matrices and through the pw2wannier90.x code
Running of wannier90.x to obtain the MLWF
Plotting the DFT bands and the bands obtained from the backtransform of the Wannier functions to check the quality of the Wannierization.
If necessary, change the input for Wannier90 to use a different initial projection, energy window, or frozen energy window. Then back to step 3. Sometimes you may also want to go back to step 2 and change the grid for the non-scf calculation.

At the end of the Wannierization the code also tells us how localized our Wannier functions are. From the output file NiO.wout we see:

Final State
WF centre and spread 1 ( -0.000000, 0.000000, 0.000000 ) 1.34669862 WF centre and spread 2 ( -0.000000, -0.000000, 0.000000 ) 0.79295816 WF centre and spread 3 ( 0.000000, 0.000000, 0.000000 ) 0.79295816 WF centre and spread 4 ( 0.000000, -0.000000, -0.000000 ) 1.34667241 WF centre and spread 5 ( 0.000000, 0.000000, 0.000000 ) 0.79297623

The numbers on the right indicate the spread of the Wannier functions as defined in the Marzari-Vanderbilt paper. One needs some experience to know whether a number is good or bad, which also depends on the type of orbital. As a rule of thumb a spread below 2 is rather localized, but how localized our WFs can actually be does depend on the system.

Input files:
NiO.scf.in: scf input file (for pw.x)
NiO.nscf.wan.in: nscf input file (for pw.x)
pw2wan.in: pw2wannier90 input file (for pw2wannier90.x)
NiO.win: Wannier90 input file (for wannier90.x)
NiO_band.dat: Wannier90 bands
NiO.bands.dat.gnu: Quantum ESPRESSO bands
plot_bands_wann.ipynb

The calculation is run by typing:
pw.x < NiO.scf.in > NiO.scf.out pw.x < NiO.nscf.in > NiO.nscf.out wannier90.x -pp NiO.winpw2wannier90.x < pw2wan.in > pw2wan.out wannier90 NiO.win

In the end you should get something like the picture below. In this model we we able to reproduce well our bands using 5 Wannier functions that resemble the Ni-d orbitals. The three bands around -5 eV are O\-p orbitals. We could easily get a good Wannierization without them because they are not strongly hybridized with the nickel orbital, to the point that there is a well-defined gap between the Ni and O bands. Whether we want to include the oxygen bands or not in the Wannierization depends on the type of physics that we want to describe.

Interfacing wannier90 and w2dynamics

Since Wannier functions are a basis set, once we are satisfied with our set of WFs we can write the DFT Hamiltonian in the Wannier basis. This is essentially a tight-binding Hamiltonian, since Wannier functions are extremely localized in real space. Compared to the DFT calculation this Hamiltonian now only contains the orbitals that really interest us. In the NiO example they are the five Ni-d orbitals.

Wannier90 automatically writes the Hamiltonian in real space as seedname_hr.dat, provided that in the input we set write_hr = .true.

Some versions of Wannier90 interfaced with other codes also automatically write the Fourier transform of the real-space Hamiltonian into reciprocal space. Unfortunately the Quantum ESPRESSO version does not, so before running w2dynamics one needs to transform it on its own.

The script HR_hamiltonian.ipynb provided below does exactly this (I also needs its auxiliary functions contained in DMFT_tools.ipynb). Thanks go to Martin Brass for sharing the original routine doing the Fourier transform.

DMFT_tools.ipynb
HR_hamiltonian.ipynb

The code is commented and should be straightforward to use/understand what it does. Executing it will perform a Fourier transform of the Hamiltonian from real to reciprocal space on a grid decided by n_kx, n_ky, n_kz. The result is written into a file called prefix_hk.dat, where the prefix in this case is NiO.

Dynamical Mean Field Theory

I have finished the part on which I am expert. I would like for this part to be continued by someone with more experience than me on DMFT calculations and analytic continuation. For the sake of completeness I will share here the rest of the calculations that lead to the physical quantities of interest, i.e. spectral function. I will try to keep comments to a minimum to avoid saying wrong things.

The input file for w2dynamics is below. I've found that for metal to insulator transitions keeping the chemical potential constant greatly helps with the stability and convergence of the calculation. As long as the chemical potential ends up between the gap, its precise value should not matter.

NiO_DMFT_fix_mu.in

Run the code with:
mpirun -np N DMFT.py NiO_DMFT_fix_mu.in

Where N are the cores you have available. This will produce a file in the hdf5 file format containing the local Green's function and self-energy on the imaginary axis for each iteration. In order to plot something physical we need to perform the analytic continuation onto the real axis.

Analytic continuation

Before you perform the analytic continuation, you want to have a look at the convergence of your Green's function. I've found that comparing the DOS from an actual DFT calculation and the "DFT" one produced by w2dynamics (which is actually generated from the Wannier Hamiltonian) is useful.

Then, you may perform the analytic continuation directly of the Green's function, or of the self energy. The scripts I used are provided below.

You should get something like the following.

plot_DMFT_convergence.ipynb

plot_DMFT_dos.ipynb

plot_DMFT_spectral_function.ipynb

Table of contents

Page tree

DFT + DMFT tutorial (NiO)