# Third tutorial¶

## Crystalline silicon.¶

This tutorial aims at showing you how to get the following physical properties, for an insulator:

• the total energy,
• the lattice parameter,
• the band structure (actually, the Kohn-Sham band structure).

You will learn about the use of k-points, as well as the smearing of the plane-wave kinetic energy cut-off.

Visualisation tools are NOT covered in this tutorial. Powerful visualisation procedures have been developed in the Abipy context, relying on matplotlib. See the README of Abipy and the Abipy tutorials.

This tutorial should take about 1 hour.

Note

Supposing you made your own installation of ABINIT, the input files to run the examples are in the ~abinit/tests/ directory where ~abinit is the absolute path of the abinit top-level directory. If you have NOT made your own install, ask your system administrator where to find the package, especially the executable and test files.

In case you work on your own PC or workstation, to make things easier, we suggest you define some handy environment variables by executing the following lines in the terminal:

export ABI_HOME=Replace_with_absolute_path_to_abinit_top_level_dir # Change this line
export PATH=$ABI_HOME/src/98_main/:$PATH      # Do not change this line: path to executable
export ABI_TESTS=$ABI_HOME/tests/ # Do not change this line: path to tests dir export ABI_PSPDIR=$ABI_TESTS/Psps_for_tests/  # Do not change this line: path to pseudos dir


Examples in this tutorial use these shell variables: copy and paste the code snippets into the terminal (remember to set ABI_HOME first!) or, alternatively, source the set_abienv.sh script located in the ~abinit directory:

source ~abinit/set_abienv.sh


The ‘export PATH’ line adds the directory containing the executables to your PATH so that you can invoke the code by simply typing abinit in the terminal instead of providing the absolute path.

To execute the tutorials, create a working directory (Work*) and copy there the input files of the lesson.

Most of the tutorials do not rely on parallelism (except specific tutorials on parallelism). However you can run most of the tutorial examples in parallel with MPI, see the topic on parallelism.

## Computing the total energy of silicon at a fixed number of k-points¶

Before beginning, you might consider working in a different subdirectory, as for tutorial 1 or 2. Why not Work3?

The following commands will move you to your working directory, create the Work3 directory, and move you into that directory as you did in the first and second tutorials. Then, we copy the file tbase3_1.abi inside the Work3 directory. The commands are:

while $ABI_TESTS/tutorial/Refs/tbase3_5.abo is a reference output file. You should find the band structure starting at (second dataset): Eigenvalues ( eV ) for nkpt= 39 k points: kpt# 1, nband= 8, wtk= 1.00000, kpt= 0.5000 0.0000 0.0000 (reduced coord) -4.83930 -2.21100 3.66138 3.66138 6.36920 8.18203 8.18203 12.44046 kpt# 2, nband= 8, wtk= 1.00000, kpt= 0.4500 0.0000 0.0000 (reduced coord) -4.97880 -2.00874 3.67946 3.67946 6.39165 8.20580 8.20580 12.47444 kpt# 3, nband= 8, wtk= 1.00000, kpt= 0.4000 0.0000 0.0000 (reduced coord) -5.30638 -1.49394 3.73328 3.73328 6.45364 8.26444 8.26444 12.56455 kpt# 4, nband= 8, wtk= 1.00000, kpt= 0.3500 0.0000 0.0000 (reduced coord) -5.69306 -0.79729 3.82286 3.82286 6.55602 8.33970 8.33970 12.65080 ....  One needs a graphical tool to represent all these data. In a separate file (_EIG), you will find the list of k-points and the eigenenergies (the input variable prteig is set by default to 1). Even without a graphical tool we will have a quick look at the values at L, $\Gamma$, X and $\Gamma$ again: kpt# 1, nband= 8, wtk= 1.00000, kpt= 0.5000 0.0000 0.0000 (reduced coord) -4.83930 -2.21100 3.66138 3.66138 6.36920 8.18203 8.18203 12.44046 kpt# 11, nband= 8, wtk= 1.00000, kpt= 0.0000 0.0000 0.0000 (reduced coord) -7.22396 4.87519 4.87519 4.87519 7.42159 7.42159 7.42159 8.26902 kpt# 23, nband= 8, wtk= 1.00000, kpt= 0.0000 0.5000 0.5000 (reduced coord) -3.01262 -3.01262 1.97054 1.97054 5.46033 5.46033 15.02324 15.02382 kpt# 39, nband= 8, wtk= 1.00000, kpt= 1.0000 1.0000 1.0000 (reduced coord) -7.22396 4.87519 4.87519 4.87519 7.42159 7.42159 7.42159 8.26902  The last $\Gamma$ is exactly equivalent to the first $\Gamma$. It can be checked that the top of the valence band is obtained at $\Gamma$ (=4.87519 eV). The width of the valence band is 12.1 eV, the lowest unoccupied state at X is 0.585 eV higher than the top of the valence band, at $\Gamma$. The Si is described as an indirect band gap material (this is correct), with a band-gap of about 0.585 eV (this is quantitatively quite wrong: the experimental value 1.17 eV is at 25 degree Celsius). The minimum of the conduction band is even slightly displaced with respect to X, see kpt # 21. This underestimation of the band gap is well-known (the famous DFT band-gap problem). In order to obtain correct band gaps, you need to go beyond the Kohn-Sham Density Functional Theory: use the GW approximation. This is described in the first tutorial on GW. For experimental data and band structure representation, see the book by M.L. Cohen and J.R. Chelikowski [Cohen1988]. Important There is a subtlety that is worth to comment about. In non-self-consistent calculations, like those performed in the present band structure calculation, with iscf = -2, not all bands are converged within the tolerance tolwfr. Indeed, the two upper bands (by default) have not been taken into account to apply this convergence criterion: they constitute a buffer. The number of such buffer bands is governed by the input variable nbdbuf. It can happen that the highest (or two highest) band(s), if not separated by a gap from non-treated bands, can exhibit a very slow convergence rate. This buffer allows achieving convergence of important, non-buffer bands. In the present case, 6 bands have been converged with a residual better than tolwfr, while the two upper bands are less converged (still sufficiently for graphical representation of the band structure). In order to achieve the same convergence for all 8 bands, it is advised to use nband=10 (that is, 8 + 2). ## Using AbiPy to automate the most boring steps¶ The AbiPy package provides several tools to facilitate the preparation of band structure calculations and the analysis of the output results. First of all, one can use the abistruct script with the kpath command to determine a high-symmetry k-path from any file containing structural information (abinit input file, netcdf output files etc.). The high-symmetry k-path follows the conventions described in [Setyawan2010]. Let’s try with: abistruct.py kpath tbase3_5.abi # Abinit Structure natom 2 ntypat 1 typat 1 1 znucl 14 xred 0.0000000000 0.0000000000 0.0000000000 0.2500000000 0.2500000000 0.2500000000 acell 1.0 1.0 1.0 rprim 0.0000000000 5.1080000000 5.1080000000 5.1080000000 0.0000000000 5.1080000000 5.1080000000 5.1080000000 0.0000000000 # K-path in reduced coordinates: # tolwfr 1e-20 iscf -2 getden ?? ndivsm 10 kptopt -11 kptbounds +0.00000 +0.00000 +0.00000 #$\Gamma$+0.50000 +0.00000 +0.50000 # X +0.50000 +0.25000 +0.75000 # W +0.37500 +0.37500 +0.75000 # K +0.00000 +0.00000 +0.00000 #$\Gamma\$
+0.50000  +0.50000  +0.50000 # L
+0.62500  +0.25000  +0.62500 # U
+0.50000  +0.25000  +0.75000 # W
+0.50000  +0.50000  +0.50000 # L
+0.37500  +0.37500  +0.75000 # K
+0.62500  +0.25000  +0.62500 # U
+0.50000  +0.00000  +0.50000 # X


To visualize the band structure stored in the GSR.nc file, use the abiopen script and the command line:

abiopen.py tbase3_5o_DS2_GSR.nc --expose -sns=talk


It is also possible to compare multiple GSR files with the abicomp script and the syntax

abicomp.py gsr tbase3_5o_DS1_GSR.nc tbase3_5o_DS2_GSR.nc -e -sns=talk


to produce the following figures:

For further details about the AbiPy API and the GSR file, please consult the GsrFile notebook .