Lab 3: Converging your Calculations

Back to Course Overview

In this lab we’ll be going through the various things you need to check to understand how well converged your calculation.

Plane-wave energy cut-off

Regardless of the type of system you’re looking at, you’ll need to check how well converged your result is (whatever it is your calculating) with respect to the plane-wave energy cut-off. This governs how many planewaves are used in the calculation.

  • In Quantum Espresso this is controlled with the parameter ecutwfc.
  • Different systems will converge differently - you shouldn’t expect carbon and silicon to be converged to the same accuracy with the same energy cut-off despite having the same structure and same number of valence electrons.
  • Different pseudopotentials for the same atomic species will also converge differently. Often pseudopotential files will suggest an energy cut-off.
  • Different calculated parameters will converge differently.
    • You should be particularly careful when calculating parameters that depend on volume, as the number of planewaves for a given energy cut-off is directly proportional to the volume so this can introduce an additional variation. We’ll see more about this later.

Silicon total energy

An example showing the total energy convergence with respect to energy cut-off is in the 01_ecut/01_silicon directory. We have already set up a series of input files which are all identical except we systematically increase the value of ecutwfc.

Examine one of these input files. You’ll notice we’ve specified some additional variables:

  • In the CONTROL section we’ve added the setting disk_io = 'none'. This suppresses the generation of the wavefunction file and the folder with the charge density etc. If we only care about the total energy we don’t need to generate these files, and the calculation will run a little faster if it’s not trying to write unnecessary files as disk io can often be a bottleneck in calculations.
  • In the ELECTRONS section we have added the conv_thr setting, though it is set to its default value. This variable controls when the self consistency cycle finishes. You should be aware of this variable, as there is little point in trying to converge to greater accuracy than we are converging self-consistently.

So now we want to run pw.x for each of these input files. Don’t forget you’ll need to load the espresso module and its dependencies before the pw.x command will work.

It’s a little tedious to type out the pw.x command for each input file in turn manually. Instead we can automate this with a small script.

Shell scripting

A shell script is a collection of Linux shell commands in a file, and when that file is executed (in the same way as any Linux command is executed) then those commands are executed. We could make a script that explicitly ran input file with pw.x in turn as follows:

#!/bin/bash

pw.x < Si_10.in &> Si_10.out
pw.x < Si_15.in &> Si_15.out
pw.x < Si_20.in &> Si_20.out
pw.x < Si_25.in &> Si_25.out
pw.x < Si_30.in &> Si_30.out

Here the first line says what command is used to interpret the script, which in this case is bash. Bash is the default shell on the system we’re using, so the same commands you type in a terminal can be used in a bash script. So we can then follow this with a series of commands as we’ve been typing directly. After we save this in a file called say run_set.sh, (e.g. by copying the above to whichever text editor you prefer such as gedit), we also need to set the file as executable using the chmod command as chmod +x run_set.sh. This command sets a flag associated with the file that tells the system this file is an executable, and can be run like any other program or command (such as pw.x or gedit or ls). To run the script, assuming you are in the same directory as the script, you’ll need to write ./run_set.sh. The preceding ./ tells the shell the path to find the executable. Recall in the first lab we discuss the $PATH variable, which has a list of directories the shell looks in for executables. If the executable is in one of these directories you don’t need to specify the path to it. Our script is not in one of these directories, so we need to explicitly tell it to run the script called run_set.sh that is contained in the current directory by typing ./run_set.sh.

There are a number of features available in bash to make scripts more general than an explicit set of commands. For example, let’s say we want to instead find whatever input files are in the current directory, and run pw.x with them, saving the output in an appropriate file, we could make our script as follows:

#!/bin/bash

for input_file in $(ls *.in)
do
    pw.x < "$input_file" &> "${input_file%.*}.out"
done

In this script, we save the output of an ls command listing all the input files as a variable. Then we loop over those input files, running pw.x with each, and saving the output in a file that has the same name as the input file but with the extension .out instead of .in. The part ${input_file%.*} returns the value stored in the $input_file variable but with the extension stripped away. This lets us make our script much more general: if we add additional input files, they’ll automatically be picked up and run without us needing to modify our script.

Save this script as run_all.sh, make it executable as described earlier, and run it within the directory with the silicon input files. At this point you should check one of your output files to make sure quantum espresso ran as expected. If your file contains only an error message, you likely forgot to load the modules required to use quantum espresso on our server.

Once the calculations have run, we want to see how the total energy changes with the input plane-wave energy cut-off. We could go through each output file and find the resulting total energy and gather them in a file, but again that would be tedious so we can write a script to do it for us (or extend our earlier script to also do this).

There are two different commands we could use to extract the resulting total energy from file. The first is often simpler to use and is called grep. For example we could use the following to print the line containing the final total energy from each output file using the fact that pw.x helpfully starts this line with a !: grep '^!.*total energy' *out. Try running this now in the directory containing your output files.

The other command we could use is awk, which is a little more powerful, and lets us pick out both the total energy value and the energy cut-off that was used with a single command. We could use this in a simple script as follows:

#!/bin/bash

awk '/kinetic-energy/{ecut=$4} /^!.*total/{print ecut, $5}' *out

Here we use the awk command to find the line with kinetic-energy and save the fourth word ($4) to a variable called ecut. Then when it finds a line that starts with a ! and has the word total, it will output the value stored in the ecut variable, followed by the fifth word ($5) on that line. Save this as a script called etot_v_ecut.sh. Make it executable, run it in the directory with the output files, saving the output as etot_v_ecut.dat. From the output we can see to how many significant figures the total energy is converged for a given energy cut-off.

If you’re interested in reading more about scripting, you can take a look at the shellscripting section. We’ll keep things relatively simple in the course, but you may find some of the more advanced functionality useful in your own work.

Plotting with Gnuplot

It’s often useful to be able to generate a quick plot when testing the relation between variables. gnuplot is a useful tool for generating plots, particularly as it is also scriptable, so we could extend our earlier script to automatically generate a plot from from the extracted data. There is a more detailed overview in the gnuplot section.

We can launch gnuplot by typing gnuplot in a terminal. Once it opens, we can, for example, plot a data file by typing plot "etot_v_ecut.dat" (provided we are in the directory containing that file). This will only plot the points by default. If we want to join these up with lines we can type plot "etot_v_ecut.dat" with linespoints. We’ll be doing more advanced things in later labs.

Carbon diamond total energy

Now let’s repeat this process to see how the carbon diamond total energy converges with the plane-wave energy cut-off. Rather than manually generate all the input files, we can modify the script that we generated earlier to run pw.x for all our input files, to take a base input file and modify it for each calculation, run it, and finally parse the output to a data file.

This is in the 01_ecut/02_carbon_diamond directory.

The script is as follows:

#!/bin/bash

template="C_diamond_base.in"
repstr="xxxx"

for val in {10..60..5}
do
  inp="C_diamond_${val}.in"
  sed "s/$repstr/$val/" $template > $inp
  pw.x < $inp &> ${inp%.*}.out
done

awk '/kinetic-energy/{ecut=$4}
     /^!.*total/{print ecut, $5}' *out > etot_v_ecut.dat

Here we’ve combined the two scripts we created for silicon above, and also automated the generation of input files using the sed command. This can be used to search for and replace some text in a file. We have set up a template input file C_diamond_base.in where we have used the placeholder text xxxx as the text we’ll search for and replace with energy cut-off we want for our input file. The bash construction for val in {10..60..5} will create a loop where the value stored in the variable $val runs from 10 to 60 in steps of 5.

  • Run this script and see what input and output files are generated.
  • See how the total energy converges with respect to ecutwfc.
  • How different is it to silicon?

Task - Graphene

  • Create a directory called 01_ecut/03_carbon_graphene and within this find the convergence of total energy with respect to the plane-wave energy cut-off.
  • Generate a plot of the difference between the total energy and it’s most converged value versus plane wave energy cut-off for both graphene and carbon diamond.
  • What do you notice about all the total energy versus energy cut-off behaviours we have seen?

K-Point convergence

The other core variable to test for convergence in a crystal system is the k-point sampling. The folder 02_kpt/01_silicon contains a template file and script to check the convergence of the total energy with respect to the k-point sampling in silicon.

Take a look at the auto_run.sh script and see how we have changed it to with respect to the previous script.

Task - Diamond

  • Create a directory called 02_kpt/02_carbon_diamond and within this find the convergence of total energy with respect to the number of k-points.
  • How does the behaviour compare to the energy cut-off convergence.

Summary

  • In this lab we looked at how to check the convergence of the total energy with respect to the plane-wave energy cut-off and k-point sampling.
    • This is done by systematically increasing the corresponding calculation parameter and looking at how the result changes.
  • We saw how we can use bash scripts to automate this process.
    • This means we don’t need to manually create a number of almost identical input files, and manually go through each one to find the values we want.
    • We can use a bash for loop to perform a calculation for a number of input files.
    • We can use grep or awk to parse results or parameters from our output files.
    • We can use sed to replace values in a template input file.
  • We can quickly generate a plot of a data file with gnuplot.

Extra - InterLayer Spacing and Box-size

For systems which are not periodic in three dimensions, such as graphene which is a 2D system, or a molecule which we are calculating within a box we need to test the convergence with respect to the inter-layer spacing or box size. This indicates that the spurious interaction between periodic images is sufficiently small for the accuracy we desire.

Task - Graphene and Methane

  • Create a folder called 03_spacing/01_graphene and using the values of energy cut-off and kpt sampling that give us a total energy accurate to within 0.001 Ry calculate the interlayer spacing that gives us a total energy converged to within the same accuracy.
  • Create a folder called 03_spacing/02_methane_20 and test the convergence of total energy versus box dimension for an energy cut-off of 20 Ry. And then create a folder called 03_spacing/03_methane_60 and test the convergence of total energy versus box dimension for an energy cut-off of 60 Ry.
    • How do these compare?