Metals and the Density of States
This week we will study the electronic structure of metals and the electronic density of states. Solving the Kohn-Sham equation for metals requires some additional considerations compared to the case of semiconductors and insulators which we studied in the last lab.
As before, all the input files and scripts you need can be found in
/opt/MSE404-MM/docs/labs/lab05
and you should copy this folder to
your home directory.
Density of States
The electronic density of states (DOS) contains information about the number of electronic states with certain energies. Mathematically, it is defined as
where \(\epsilon_{n\mathbf{k}}\) are the Kohn-Sham eigenvalues for band \(n\) and k-point \(\mathbf{k}\), \(w_{\mathbf{k}}\) is the k-point weight to keep the total integrated DOS to be the number of bands. I.e.,
Note that the sum over k-points can be replaced by an integral if the grid is sufficiently fine.
For a molecular system, the DOS consists of a series of discrete peaks at the energies of the molecular Kohn-Sham orbitals, since we only use one k-point (the \(\Gamma\) point) for DFT calculations of isolated molecules. In contrast, for DFT calculations of crystals, we use many k-points to sample the first Brillouin zone and in this case the discrete peaks merge to form a continous curve. For example, here are graphs of the DOS of a water molecule (isolated molecule) and carbon diamond (periodic crystal):
The shape of the DOS is intimately connected to the band structure of a crystal: bands with a strong dispersion (i.e. the Kohn-Sham energies change rapidly along a path in k-space) result in low DOS spread over a large interval, whereas less dispersive (i.e. flatter) bands result in high DOS in a smaller energy interval.
Van Hove Singularity
While the DOS is a continuous function, there can be sharp peaks in the DOS that breaks the smoothness of the curve. These are called Van Hove singularities and they occur when the band structure has a saddle point or an inflection point. They are important because they can lead to interesting properties in optical absorption spectra of materials.
For insulators and semiconductors, the DOS is zero inside the band gap, as there are no Kohn-Sham states in that energy range. Hence, the DOS allows us to determine the band gap of a crystal (unlike a band structure plot which only shows the KS energies along a specific path in the Brillouin zone, the DOS contains information about the KS energies at all k-points in the first Brillouin zone).
Smearing
To evaluate the equation for the DOS of a crystal numerically, we have to replace the delta-functions by functions with a finite width, such as a Gaussian (otherwise the DOS will always look like a set of sharp peaks). This is called "smearing".
The width of the Gaussian is a parameter that needs to be tuned to ensure that a physically meaningful result is obtained:
- If the broadening is too large, you may smear out important features of the DOS.
- If the broadening is too small, you will see some unphysical peaks in the DOS and it will look very spikey.
- In principle you want the smearing to be comparable to the typical difference of Kohn-Sham energies at neighboring k-points. In practice, however, it is often easier to just try different values for the broadening until the DOS looks physically meaningful.
Tetrahedron Method
Another way to obtained an accurate DOS is to linearly interpolate the calculated Kohn-Sham energies using the so-called tetrahedron method. Such a calculation is noticeably slower than using a broadening but there is no need to to worry about using the correct smearing. Your only convergence parameter is the fineness of the k-point grid you use which determines how accurate the interpolation is.
It's important to note that in a real measurement of the density of states of material (for example, using photoemission spectroscopy), there is always some broadening coming from
-
Electron-phonon coupling: since the atoms in a real material vibrate around their equilibrium positions, the energies of the electronic states will be smeared out.
-
Any measurement apparatus will have a finite energy resolution which will further broaden the measured DOS.
Calculating the DOS
In a similar way to the electronic band structure, we produce the density of states plot in three steps.
Step 1 - SCF Calculation
Perform a self-consistent calculation as before to produce a converged charge density.
Task 1.1 - SCF Calculation
Run pw.x
using the input file
01_C_diamond_scf.in
for diamond.
1 |
|
Step 2 - NSCF Calculation
Take the density calculated in the previous step and use it to perform a non-self-consistent calculation on a denser k-point grid. We want an accurate description of how the KS energies vary as we move around the Brillouin zone so we use a much denser grid here than we need to obtain the converged density in the previous step.
The difference between this and the band structure calculation is that here we use a uniform sampling of the Brillouin zone, rather than a path between k-points. The input file for this calculation can be found at 02_C_diamond_nscf.in:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
|
calculation = nscf
specifies that we are calculating the non-self-consistent calculation.nbnd = 8
specifies that we want to calculate 8 bands.K_POINTS automatic
specifies that we are using an automatically generated k-point grid. We've increased the k-point sampling to a 20x20x20 grid, and we have removed the shift (without a shift, Quantum Espresso uses the convention that the gide always includes the origin, for more information see this). Many systems have a valence band maximum or conduction band minimum at the gamma point, so it is good to ensure it's explicitly included in the grid.
Task 1.2 - NSCF Calculation
Run pw.x
using the input file
02_C_diamond_nscf.in
for diamond.
1 |
|
Step 3 - Density of States Calculation
From the Kohn-Sham energies calculated on the dense k-point grid we then
calculate the density of states using dos.x
.
03_C_diamond_dos.in is the input
file for dos.x
and contains just a DOS
section. What this input file
instruct dos.x
to do is to discritise an energy range
Emin
and
Emax
using an
interval of
DeltaE
and at
each energy calculate the DOS using the formula given
above. Each delta function is replaced using a gaussian
function (which can be replace by other functions using
ngauss
tag) with
a width of
degauss
.
1 2 3 4 |
|
degauss
specifies the Gaussian broadening (in Rydberg) to use in the density of states calculation.DeltaE
specifies the spacing between energies at which the DOS is calculated, in eV.
Note
We have picked similar values for degauss
and DeltaE
(after converting
them to the same units). In fact if degauss
is not specified, and no
broadening scheme is used in the DFT calculation, degauss
will take the
value of DeltaE
by default. You can check the documentation
INPUT_DOS for more
details.
Task 1.3 - Plotting Density of States
Run dos.x
using the input file
03_C_diamond_dos.in
for diamond.
1 |
|
The final step produces a file named pwscf.dos
by default. This is a simple
text file you can plot. It has three columns: Energy (eV), Density of States
(states/eV), and Integrated Density of States (states) up until that energy. And
in the first line, the estimated Fermi energy is also given.
1 2 3 4 5 6 |
|
It is customary to shift the x-axis in the plot such that the Fermi energy
or valence band maximum is at 0. A value for the Fermi level is given in
the file header of pwscf.dos
. This value is determined in a simple way from
the integrated density of states. Sometimes, it may be worth calculating the
Fermi level in a separate calculation using a relatively small broadening if
you're studying a metallic system (as we shall see later). For semiconductors
and insulators you can determine the valence band maximum energy from the output
file.
Fermi Energy From The Integrated Density of States
The integrated density of states is the number of states up to a certain energy. I.e., \(\int_{-\infty}^{E_1} DOS(E) d(E)\). One can easily found the Fermi energy by looking at the energy where the integrated density of states is equal to the number of electrons in the system. However, it is worth noting that since we are using a smearing scheme and a denser k-point grid in the NSCF step (which is not the case for the SCF calculation), the Fermi energy calculated in this way may not be accurate. For accurate Fermi energy, one should always refer to the output of the SCF calculation.
The directory 03_densityofstates
contains a python script that can be used to
plot the shifted DOS together with the integrated DOS.
Task 1.4 - Density of States Calculation
Plot the density of states using the script provided.
1 |
|
DOS Plot
Now try re-running Tasks 1.3 and 1.4 with different degauss
values and see
how the DOS plots change. Does the change match your expectations?
Metals
Metals have a Fermi surface in k-space which separates the occupied from the unoccupied Kohn-Sham states. The shape of the Fermi surface can be quite complicated and is not known a priori. This means that in contrast to insulators or semiconductors where every k-point has the same number of occupied states, in a metal the number of occupied states can vary from k-point to k-point. This makes it more difficult to calculate the electron density which involves a summation of the squared magnitudes of the wavefunctions of all occupied states.
Tackling Discontinuities
Generally, there are two things that we typically do for metals to help with the convergence of the SCF calculation:
-
Use a denser k-point grid than you would need for a semiconductor or insulator. This allows us to have a better description of the complex Fermi surface that metals have.
-
Use a smearing scheme for the calculation of the occupation number of Kohn-Sham states at each k-point. Instead of having a integer occupation number of each band (0 or 1), we now allow a fractional occupation number so that the occupation number varies smoothly across k-points.
Charge Sloshing Effect
Without occupation smearing, each band is either fully occupied or unoccupied. At each step of the SCF loop, bands may switch between being occupied and occupied, leading to a large change in the electron density which can slow down the convergence of the SCF calculation.
One way to obtain smeared occupation numbers is to calculate them using the Fermi-Dirac distribution \(f_T(\epsilon_{n\mathbf{k}})\) at a finite temperature \(T\)
\[ f_T(E) = \frac{2}{e^\frac{E-E_F}{k_B T} + 1}, \]where \(E_F\) denotes the Fermi energy which we set to the energy of the highest occupied Kohn-Sham state and the factor of 2 results from the Pauli principle. We can calculate the Fermi level using the condition that the sum over all occupation numbers must be equal to the total number of electrons \(N_e\) in the cell(\(w_\mathbf{k}\) is the weight of the k-point):
\[ N_e = \sum_n \sum_\mathbf{k} w_\mathbf{k} f_T(\epsilon_{n\mathbf{k}}). \]At low temperatures, the Fermi-Dirac distribution becomes a discontinous step function. So we must choose the temperature to be sufficiently high to ensure smooth occupation numbers near the Fermi surface (using the tag
degauss
in the input file). It is worth noting that other smearing methods, such as Gaussian smearing.Using a smearing scheme for the occupation numbers helps significantly in achieving a smooth SCF convergence for metals, as otherwise a small change in the energy of a KS state from one self-consistent iteration to the next one could lead to a very large change in its occupation number and to the electron density. We set the smearing scheme (for both DOS and occupation function) and width with the
occupations
anddegauss
variables in the input file.Additional read
For a thorough review of the most typical smearing method used in DFT, have a read at this paper.
Example: Aluminium
Aluminium has an FCC crystal structure with one atom per unit cell, which we know how to deal with at this point. The thing about Aluminium that makes it more complicated within DFT is that it is a metal.
Here is an example input file for a SCF calculation of Aluminium:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
|
- The
verbosity
variable is set tohigh
to give more information about the SCF calculation, including occupation of the bands at each k-point. - The
occupations
variable is set tosmearing
to tell Quantum Espresso to use a smearing scheme input description. - The
smearing
variable is set tofermi-dirac
to tell Quantum Espresso to use a Fermi-Dirac smearing scheme. input description. - The
degauss
variable is set to 0.1d0 to set the width of the smearing. see input description.
Task 2.1 - Smearing for Metals
First, run the pw.x
calculation with the supplied input file in
02_aluminium/Al.in.
Then, take a look at the output file. Can you find the Kohn-Sham eigenvalues at each k-point and the corrresponding occupations number.
Example
``` k =-0.0625 0.0625 0.0625 ( 59 PWs) bands (ev):
3.5809 40.9721 43.7477 43.7477 49.9157 53.1364
occupation numbers 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000
1 |
|
4.7127 35.0485 40.5915 44.7949 50.9368 50.9459
occupation numbers 1.0000 0.0001 0.0000 0.0000 0.0000 0.0000 ... ```
For a metal the default number of bands that are used in the SCF calculation
is at least four more than are needed to accommodate all electrons (even
without us setting the nband
tag). The pseudopotential we have used for
aluminium has 3 valence electrons, which could be accommodated in two bands.
Adding the four extra bands, this gives a total of 6 bands which are used in
the calculation.
Now, try to play around with the degauss
value and see how the
occupations change. Does it match your expectations?
Now you know how to use smearing to help with the convergence of a metal, you can obtain the density of states and band structure of a metal just as easy.
Task 2.2 - DOS of Aluminium
Try calculating the density of states of Aluminium by following Tasks 1.1 to 1.4. Remember to change the Fermi energy in the python script to the correct value. Can you find an energy gap? Where does the Fermi level lie?
DOS Plot
Task 2.3 - Band Structure of Aluminium
Try calculate the band structure of Aluminium following the path of:
1 2 3 4 5 6 7 8 9 10 |
|
Band Structure Plot
Summary
In this tutorial, we have learned:
- How to use the
dos.x
code from the Quantum Espresso package. - How to treat a metallic system.