by Pablo Garcia Risueno & Pasquale Pavone for exciting neon
(Jupyter notebook by Mara Voiculescu & Martin Kuban)
Purpose: In this tutorial, you will investigate how the choice of computational parameters can affect the result of an exciting
calculation. The procedure which is used in this case is often called convergence study and allows to obtain the best value of such parameters. Explicit examples are presented for the convergence of the total energy of silver and diamond with respect to the choice for the k-points sampling (related to the parameter ngridk
) and the dimension of the basis set (related to the parameter rgkmax
).
Table of Contents
1. Convergence of Total-Energy Calculations
Read the following paragraphs before starting with the rest of this tutorial!
Before running any Jupyter tutorials, please refer to the 00_before_starting.md
document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting
's root directory:
source tools/excitingjupyter/venv/excitingvenv/bin/activate
Here is a list of the scripts found within the excitingscripts
module which are relevant for this tutorial with a short description.
setup.convergence_test.py
: Python script for generating input files with different values of the main computational parameters.
execute.convergence_test.py
: Python script for running a series of exciting calculations.
plot.convergence.py
: Python script for visualization of convergence results.
As a first step, you may set the exciting_root directory and create a running directory for the notebook.
%%bash
mkdir -p run_simple_convergence_tests
A ground-state calculation using a DFT code like exciting
fundamentally depends on 2 parameters
groundstate
attribute: ngridk
);groundstate
attribute: rgkmax
).In order to be able to rely on your calculation, you need to understand what these two parameters mean and what is the effect of a change of their values on the physical quantities which are relevant for you. In this tutorial, we focus on the effect of such changes on the total energy of a crystal only for pedagogical reasons. In the general case, this kind of analysis must be made for all the properties of interest (e.g., lattice parameters, bulk modulus, equation of state, bandgap energies, etc.).
For more informations on these parameters, please expand the links below.
The mesh of k-points
From basic considerations of solid-state theory, it can be shown that the Schrödinger (or, in our case, Kohn-Sham) equation for a periodic system has a bunch of solutions, where each solution is characterized by a vector k in reciprocal space. This k-vector is essentially related to the periodicity of the corresponding solution of the Schrödinger equation: Roughly speaking, the solution will have the same value at coordinates in real space that, along the direction of k, have a distance of 2$\pi$/|k|.
Many properties of a solid, including the total energy, are represented as integrals, performed over all possible k-vectors. Obviously, the direct calculation of these integrals is very demanding (one should calculate the solution of the Schrödinger equation for a huge number of k-points). Therefore, the integrals are approximated by sums performed on a set of k-points distributed on a finite grid. The spacing of the points on the grid is a measure of the accuracy of the calculation of the integrals which also depends on how fast the integrand (the solution of the Schrödinger equation) varies by changing k. The challenge is to find a good number of k-points: Large enough to capture the physical properties of your system well, but small enough to keep calculations as fast as possible.
Some qualitative considerations:
In exciting
the parameter which is associated to the choiche of the mesh of k-points is ngridk
, see Input Reference for further details.
The basis-set size
Most of the DFT codes solve the Kohn-Sham equation in reciprocal space, expanding the wavefunctions in terms of suitable periodic basis functions. In our case, the basis functions are the linearized augmented plane waves (LAPW). The larger the size of the basis set, the more accurate is the calculation. The corresponding groundstate
attribute in the exciting code is rgkmax
(see Input Reference). The default value of rgkmax
is 7.0. In some cases, (e.g., for a good description of forces), larger values (8.0 — 9.0) may be necessary to get reliable results. In some other cases, (e.g., for describing very light materials like hydrogen), already smaller values of rgkmax
may be sufficient.
Scaling of the calculation time
The calculation time scales
rgkmax
to the power of 9.In this first example, we investigate the dependence of the total energy on the parameters ngridk
and rgkmax
for the Ag crystal.
We start by creating an exciting
(xml) input file called input.xml which should appear as the one below.
<input>
<title>Ag: Convergence test</title>
<structure speciespath="$EXCITINGROOT/species">
<crystal scale="7.7201">
<basevect>0.5 0.5 0.0</basevect>
<basevect>0.5 0.0 0.5</basevect>
<basevect>0.0 0.5 0.5</basevect>
</crystal>
<species speciesfile="Ag.xml">
<atom coord="0.0 0.0 0.0" />
</species>
</structure>
<groundstate
xctype="GGA_PBE_SOL"
ngridk="4 4 4"
rgkmax="5.0">
</groundstate>
</input>
The next step is writing the complete input as a string and saving it in your working directory as input.xml.
N.B.: Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT using the command
%%bash
cd run_simple_convergence_tests
python3 -m excitingscripts.setup.excitingroot
cd ..
In order to run exciting
from the terminal, you simply need to execute the exciting_smp binary in the running directory. After a few seconds, the calculation should be finished.
Here we used the time
command before exciting_smp in order to get, at the end of the run, the elapsed time explicitly written on the screen.
%%bash
cd run_simple_convergence_tests
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..
After a few seconds, the calculation will be completed. Indeed, the calculation you have just performed uses very rough values for the parameters ngridk
and rgkmax
. The corresponding total energy (in Hartree) can be found as the last non empty line of the file TOTENERGY.OUT. You can save the value of parameters and total energy at a text file called convergence-test as indicated below (here, only the first value of the parameter ngridk
is considered).
4 5.0 -5314.29682116 # ngridk rgkmax total-energy
In order to improve the calculation, we should change the value of ngridk
and rgkmax
and repeat the calculation in a systematic way. By changing the values in the input file of one (or both) parameters, e.g.
...
<groundstate
xctype="GGA_PBE_SOL"
ngridk="6 6 6"
rgkmax="5.0">
</groundstate>
...
and performing a new calculation using the command exciting_smp, we can extract the value of the calculated total energy and add it to the file convergence-test, which now should appear as:
4 5.0 -5314.29682116 # ngridk rgkmax total-energy
6 5.0 -5314.31785599
Then, the same procedure can be repeated for other values of the parameters.
The procedure described above can be easily performed with the help of some useful script. For instance, the script excitingscripts.setup.convergence_test
can be used for preparing inputs for a sequence of calculations. You can execute the script in your running directory in the following way.
%%bash
cd run_simple_convergence_tests
python3 -m excitingscripts.setup.convergence_test 4 20 5 5
cd ..
where the first 2 arguments are the minimum and maximum values (with a step of 2) for ngridk
, respectively. The last 2 (integer!) entries are the minimum and maximum values (with a step of 1) for rgkmax
, respectively. This means that, in this example, we are considering the convergence of the total energy as a function of ngridk
at the constant value of rgkmax
= "5".
The complete run is performed by the script excitingscripts.execute.convergence_test
.
%%bash
cd run_simple_convergence_tests
python3 -m excitingscripts.execute.convergence_test 4 20 5 5
cd ..
The result of this sequence of calculation is summarized automatically in the file convergence-test, which now should appear as:
!cat run_simple_convergence_tests/convergence-test
To visualize these results you can use the following command (the line entry k indicates the running on different values of ngridk
).
%%bash
cd run_simple_convergence_tests
python3 -m excitingscripts.plot.convergence k
cd ..
Then, you can use this procedure to plot the convergence for the following examples. In order to visualize the different results, we create a subdirectory for each of the examples.
ngridk
at fixed rgkmax
= "9"%%bash
cd run_simple_convergence_tests
mkdir -p silver_4_20_9_9 && cd silver_4_20_9_9
cp ../input.xml .
python3 -m excitingscripts.setup.convergence_test 4 20 9 9
python3 -m excitingscripts.execute.convergence_test 4 20 9 9
python3 -m excitingscripts.plot.convergence k
cd ../..
rgkmax
at fixed ngridk
= "12 12 12"%%bash
cd run_simple_convergence_tests
mkdir -p silver_12_12_5_9 && cd silver_12_12_5_9
cp ../input.xml .
python3 -m excitingscripts.setup.convergence_test 12 12 5 9
python3 -m excitingscripts.execute.convergence_test 12 12 5 9
python3 -m excitingscripts.plot.convergence r
cd ../..
rgkmax
and ngridk
%%bash
cd run_simple_convergence_tests
mkdir -p silver_4_20_5_9 && cd silver_4_20_5_9
cp ../input.xml .
python3 -m excitingscripts.setup.convergence_test 4 20 5 9
python3 -m excitingscripts.execute.convergence_test 4 20 5 9
python3 -m excitingscripts.plot.convergence rk
cd ../..
In this subsection we refer to the results of Example 1. In this case the convergence-test file is:
4 9.00 -5314.79400821
6 9.00 -5314.79814418
8 9.00 -5314.79964620
10 9.00 -5314.80025331
12 9.00 -5314.80060003
14 9.00 -5314.80081026
16 9.00 -5314.80084371
18 9.00 -5314.80090872
20 9.00 -5314.80092356
Due to the fact that the difference between energies in contiguous rows becomes (increasingly) smaller going to large values of ngridk
, this difference can be interpreted as a measure of the accuracy in the value of the calculated energy.
Repeat what you have done in the previous section for diamond. In contrast to silver, which is a metal, diamond an insulator. The corresponding input for exciting
(input.xml) could look like the following.
<input>
<title>Diamond: Convergence test</title>
<structure speciespath="$EXCITINGROOT/species">
<crystal scale="6.719">
<basevect>0.5 0.5 0.0</basevect>
<basevect>0.5 0.0 0.5</basevect>
<basevect>0.0 0.5 0.5</basevect>
</crystal>
<species speciesfile="C.xml">
<atom coord="0.00 0.00 0.00" />
<atom coord="0.25 0.25 0.25" />
</species>
</structure>
<groundstate
xctype="GGA_PBE_SOL"
ngridk="4 4 4"
rgkmax="5.0"
gmaxvr="14">
</groundstate>
</input>
Finally, on the basis of what you learned above, try to reproduce the results for Diamond displayed in the following.
%%bash
cd $EXCITINGROOT/tools/excitingjupyter/excitingjupyter/01_getting_started/reference_data_simple_convergence_tests/diamond_4_20_5_5
python3 -m excitingscripts.plot.convergence k
cd ../..
%%bash
cd $EXCITINGROOT/tools/excitingjupyter/excitingjupyter/01_getting_started/reference_data_simple_convergence_tests/diamond_4_20_9_9
python3 -m excitingscripts.plot.convergence k
cd ../..
%%bash
cd $EXCITINGROOT/tools/excitingjupyter/excitingjupyter/01_getting_started/reference_data_simple_convergence_tests/diamond_12_12_5_9
python3 -m excitingscripts.plot.convergence r
cd ../..
%%bash
cd $EXCITINGROOT/tools/excitingjupyter/excitingjupyter/01_getting_started/reference_data_simple_convergence_tests/diamond_4_20_5_9
python3 -m excitingscripts.plot.convergence rk
cd ../..
ngridk
faster for silver or diamond? Why?rgkmax
?