First crystals

In this exercise we will move beyond molecules to treat bulk systems, in this case, crystals.

Hint

Enter the directory ‘MgO’

Input geometry

When dealing with crystaline systems (or any kind of periodic system), it is extremely important to provide the correct structural information. This includes:

  • The atomic coordinates.

  • The size of the simulation cell.

Take a look at MgO.fdf. Within it, you will first find the following specification for the simulation cell:

LatticeConstant 4.117 Ang
%block LatticeVectors
   0.000  0.500  0.500
   0.500  0.000  0.500
   0.500  0.500  0.000
%endblock LatticeVectors

The lattice constant multiplies all of the lattice vectors, meaning that our simulation cell is defined by the cell vectors (0, 2.0585, 2.0585), (2.0585, 0, 2.0585) and (2.0585, 2.0585, 0). Note this cell is not cubic, since we are taking the primitive cell for MgO (which has 60° angles) instead of the unit cell.

Directly below the cell vectors, we find the atomic coordinates for the crystal:

AtomicCoordinatesFormat Fractional
%block AtomicCoordinatesAndAtomicSpecies
  0.000  0.000  0.000   1
  0.500  0.500  0.500   2
%endblock AtomicCoordinatesAndAtomicSpecies

By specifying “Fractional” as the format, atomic positions get defined as multiples of the cell vectors. So, while atom 1 is on (0,0,0), atom 2 is halfway through all cell vectors.

Sampling the Brillouin zone (k-points)

In bulk systems the electronic states can be characterized by a continuous quantum number k, the so-called Bloch vector or crystalline quasi-momentum of the electrons. The space of allowed Bloch vector (the Brillouin zone) is usually sampled using a finite grid. The most popular scheme to generate this integration mesh is the Monkhorst-Pack algorithm. This is also used in SIESTA. The integration grid in the Brillouin zone is specified using the block kgrid_Monkhorst_Pack:

%block kgrid_Monkhorst_Pack
   6  0  0  0.5
   0  6  0  0.5
   0  0  6  0.5
%endblock kgrid_Monkhorst_Pack

where the 3x3 matrix of integer numbers defines the mesh and the last number (an optional float) in each line defines the displacement, which might be useful to reduce the number of (symmetry irreducible) k-points used in the calculation (for this, 0.5 is a good default), or, for some applications, to exclude the Gamma point from consideration. If not specified, the displacement will be zero in each direction. The program will warn the user if the displacements chosen are not optimal.

Explore the MgO.fdf file in search for the Monkhorst Pack block.

Run the simulation, which will take a few seconds:

siesta  MgO.fdf > MgO.out

Band structure and density of states

We will now analyse the band-structure for MgO and look at the density of states.

To plot the DOS you will use the utility program ‘Eig2DOS’ (see this how-to for background and installation notes, if needed). This program reads the file MgO.EIG (always produced by Siesta), which contains all the eigenvalues for each k-point used to sample the BZ, and the file MgO.KP, which contains the k-point sampling information. Useful options to the program (type ‘Eig2DOS -h’ for a full list) are the broadening for each state in eV (a value of the order of the default (0.2 eV) is usually reasonable), the number of energy points where the DOS will be calculated (200 by default) and the Emin and Emax of the energy window where the DOS will be calculated (the default is to compute the DOS in the whole range of energies available in the EIG file).

For example:

Eig2DOS -k MgO.KP MgO.EIG > dos

will compute the DOS in the whole range, using a grid of 200 points, a broadening (“smearing”) of 0.2 eV, and using the k-point info from MgO.KP. (The k-point file information is important when different k-points have different weights).

Plot the dos using gnuplot, type gnuplot and enter the commands:

plot "dos" with lines

The result is a bit confusing, but it is just because the range of energies is very large and the number of grid points not so big. Open the MgO.out file and look for the Fermi energy. Now generate the dos only for a range of energies of about 30 eV around the Fermi energy (use the -E and -e options of Eig2DOS). The result will be more familiar.

The file MgO.fdf will also produce a file MgO.bands containing the band structure along the several high-symmetry lines in the Brillouin zone (BZ). This file is produced because in input the the block BandLines was included:

BandLinesScale       pi/a
%block BandLines
1   1.5   1.5   0.0   K             # Begin at K
38  0.0   0.0   0.0   \Gamma        # 38 points from K to Gamma
36  0.0   2.0   0.0   X             # 36 points from Gamma to X
18  1.0   2.0   0.0   W             # 18 points from X to W
26  1.0   1.0   1.0   L             # 26 points from W to L
31  0.0   0.0   0.0   \Gamma        # 31 points from L to Gamma
%endblock BandLines

Let’s understand these few lines:

  • The BandLinesScale entry specifies the scale of the k-vectors given in the BandLines block. The possible options are “pi/a” and “ReciprocalLatticeVectors”.

  • Each line of the BandLines block represents the end-point of a segment to be spanned in the reciprocal space. The comments on the code block above help to familiarize with the format of these k-space segments.

To learn more about the special high-symmetry points and lines of the BZ for this case, you can visit the Bilbao Crystallographic Server (click on ‘Optimized listing of k-vector types…’; then on ‘Brillouin Zone’).

The bands are stored in the file MgO.bands. To plot the band structure you need to use the utility program gnubands:

gnubands < MgO.bands > bands.dat

and you can plot bands.dat with gnuplot, using the bands.gplot file provided (which has extra information to place and label the ‘ticks’ for each symmetry point in the BZ):

gnuplot -persist bands.gplot

Hint

Enter the directory ‘Al’

For metals such as Al there are electronic bands that are not completely filled and therefore for an accurate description of the total energy, forces and all properties of the materials it is necessary to use a better sampling in reciprocal space (Bloch vectors) than for insulators. In the input file Al_bulk.fdf a 4x4x4 grid is used. This might be insufficient for a good description of aluminium.

Run the input as it is.

siesta < Al_bulk.fdf > Al.out

And create the DOS

Eig2DOS -k Al.KP Al.EIG > dos

If you plot the DOS with gnuplot the result does not look ok.

Create a new folder

mkdir MORE_KP
cp Al_bulk.fdf Al.psf MORE_KP
cd MORE_KP

Change in the Al_bulk.fdf file the mesh to [14 14 14] and run siesta.

siesta < Al_bulk.fdf > Al.out

Then:

Eig2DOS -k Al.KP Al.EIG > dos

Plotting the DOS, we know recognize the “free-electron-like” curve we see in textbooks. For coarse samplings, instead, the DOS was not at all like the “free-electron-like” curve since too few points were considered.

Another point to notice is that the DOS seems to have an upper limit (DOS is zero after a certain energy). This is due to the limited number of orbitals in the basis that, in turn, limits the number of produced bands. A richer basis will increase the number of bands (see tutorial on basis.)

Try to compare also the bands among the case with few kpoints and more kpoints.

A detailed tutorial on kpoints convergence is here.