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 settingdisk_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 theconv_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
orawk
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 called03_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?