Solution LAB1 kpt convergence
- Back to the previous page: Structural and electronic properties of semiconductors and metals #Exercises
Convergence wrt k-points
A possible script to run multiple scf calculations using pw.x is
LabQSM/LAB_1/test_diamond/run_diamond.sh
As a first step, we inspect the script. The most relevant parts are highlighted here:
[...] # # set vars # ecut=40.0 alat=6.741 # nk_list="2 4 6 8 12 16"
#
# main loop
#
for nk in $nk_list ; do
#
label="nk${nk}"
#
filein=scf_${label}.in
fileout=scf_${label}.out
#
# generate the input file
#
cat > $filein <<EOF
&CONTROL
[...] # note here you can use shell variables
/
[...]
EOF
#
# running
echo "Running $label"
#
$para_prefix $bindir/pw.x $para_postfix < $filein > $fileout
done
Running the calculations.
Once we have understood how to run multiple calculations at once, we can produce the actual data. In order to run the calculations, create e dedicated folder, cd into it and run the script:
$> mkdir RUNS_diamond $> cd RUNS_diamond $> cp $HOME/LabQSM/LAB_1/test_diamond/run_diamond.sh ./ $> ./run_diamond.sh # actual running
You should get a message similar to the following:
$> ./run_diamond.sh Running nk2 Running nk4 Running nk6 Running nk8 Running nk12 Running nk16
The next step is to check that the calculations have been run successfully, without crashes, and properly converging.
In order to do this inspect one or more output files.
$> ls scf_nk* # check whether the files exist scf_nk12.in scf_nk16.in scf_nk2.in scf_nk4.in scf_nk6.in scf_nk8.in scf_nk12.out scf_nk16.out scf_nk2.out scf_nk4.out scf_nk6.out scf_nk8.out
Extracting data.
Assuming all is good, we now want to extract the computed total energies as a function of the k-point grid used.
Converged total energies are marked by ! in the output of pw.x
$> grep ! scf_nk*out scf_nk12.out:! total energy = -22.78078456 Ry scf_nk16.out:! total energy = -22.78074173 Ry scf_nk2.out:! total energy = -22.54574059 Ry scf_nk4.out:! total energy = -22.77088457 Ry scf_nk6.out:! total energy = -22.78012390 Ry scf_nk8.out:! total energy = -22.78070087 Ry
You can either redirect this information to file and then edit the file, or use some unix commands to directly extract the data you want
$> grep ! scf_nk*out | sed 's/scf_nk//' $> grep ! scf_nk*out | sed 's/scf_nk//' | sed 's/.out:!.*=//' $> grep ! scf_nk*out | sed 's/scf_nk//' | sed 's/.out:!.*=//' | sort -n > conv_kpt.dat
At this point we obtain
$> cat ./conv_kpt.dat 2 -22.54574059 Ry 4 -22.77088457 Ry 6 -22.78012390 Ry 8 -22.78070087 Ry 12 -22.78078456 Ry 16 -22.78074173 Ry
We may want to further edit this file also adding the total number of reducible k-points
(i.e. nk^3) and the number of irreducible k-points.
This last information can be directly extracted from the input file, searching for the
string number of k points= , either manually, or using grep.
After a bit of editing of the conv_kpt.dat file, ve obtain:
# grid kpt_BZ kpt_IBZ Etot [Ry] 2 8 3 -22.54574059 4 64 8 -22.77088457 6 216 16 -22.78012390 8 512 29 -22.78070087 12 1728 72 -22.78078456 16 4096 145 -22.78074173
Plotting.
We are now in the condition to plot the data.
In particular we can plot the total energy (col 4) as a function of the grid generator parameter (col 1)
$ gnuplot > plot "conv_kpt.dat" u 1:4 w lp Or $ xmgrace → data → import → ascii...
Some documentation concerning plotting can be found:
- ./LabQSM/tools (GitHub / training material) or
- Plotting and visualization tools (this wiki)
The plot should look like:
Comments:
- The use of symmetries, allowing one to use kpts in the irreducible BZ (IBZ) instead of in the whole BZ largely reduce the number of kpts
- The calculation time is linearly proportional to the number of actual kpts used
- Using 8 kpts (automatic grid 8 8 8 1 1 1) seems to be a good compromise giving an error bar on the total energy of the order of 2*10^-5 Ry.
How to choose such threshold ? It depends on what we want to compute. K-point grids are a numerical parameter, and we want the convergence wrt such parameter not to impact the description or accuracy of computed physical quantities. E.g. If we are interested in the binding energy of a molecule adsorbed on a surface, we want the numerical error bars related to numerical parameters to be smaller (one to 2 orders of magnitude, at least) than the physical binding energies (the physical quantities of interest).
Only in this way we can say to have computed the physical quantity of interest within a given level of theory (e.g. KS-DFT using LDA, implemented using plane0waves and pseudopotentials).