Solution LAB1 kpt convergence

• Prev: Structural and electronic properties of semiconductors and metals#Exercises
• Next: LabQSM#Link to Selected Laboratories

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).