Quantum Espresso
Quantum Espresso is a freely available package of open-source codes for electronic-structure calculations and materials modelling at the nanoscale. It is based on density-functional theory, plane waves, and pseudopotentials, which you will be learning about in the theoretical part of the course.
Quantum espresso is used via the command line. There is no graphical interface by default. This is typical of most electronic structure codes, where you are often interacting with a remote HPC system solely via ssh, and submitting calculations and text-file scripts to a queueing system which takes your calculation and executes it on some set of compute server nodes when they become available.
To run a calculation you first need to make an input file, describing the various calculation parameters along with giving the location of any other input files that will be used. Then you run the code giving it your input file (redirected from stdin in the case of quantum espresso - see lab 1), and it will create one more files with the result of your calculation and potentially intermediate data also.
First you’ll need to copy the material for this Lab to your home directory.
You should already have a directory named MSE404
in your home directory; you
can confirm this by going to your home directory and checking the contents
with the following two commands:
cd ~
ls
If you don’t see an MSE404
folder listed, you can create it with
mkdir MSE404
The above command requests the creation of a directory named MSE404
within
whatever directory you are currently. If you weren’t in your home directory
you could still create a directory there, by giving the absolute path: mkdir
~/MSE404
. Once you’ve confirmed you have the directory as expected, you can
make a copy of the lab material there so that you can run the calculations for
this lab, with cp -r /opt/Courses/MSE404/lab02 ~/MSE404
.
A set of basic input files for a variety of systems have been set up in the
different directories here. In the material that follows you’ll need to run
the calculations in your own copy of this directory. You’ll get an error if
you try to run calculations anywhere under the /opt
directory as you don’t
have the permissions needed to create files here.
You should also ensure you have a copy of the directory with the pseudopotentials used for the labs:
cd ~/MSE404
cp -r /opt/Courses/MSE404/pseudo .
A basic input file
A simple example that will use the pw.x
code from quantum espresso to
calculate the total energy of carbon in the diamond structure is given in the
directory 01_carbon_diamond
. Let’s first open the input file
C_diamond.in
and read through it.
In this case we are relying on many default values that would normally be
specified, but this lets us focus on the most important things to begin. There
is an annotated version of this file called
C_diamond_detailed.in
which will let us go through the input in more detail. Let’s look at this now.
Running the calculation
The quantum espresso package has been compiled as a module on the mt-student server. As discussed in the previous lab, modules are often used on HPC systems to make different versions of various packages as compiled with different compilers available to users. To add quantum espresso to your environment, along with its dependencies type the following in a terminal:
module load gcc mkl espresso
Now to run the code: make sure you are in the 01_carbon_diamond
directory
then do:
pw.x < C_diamond.in
Here we redirected our input file to the stdin of pw.x
, as described in the
IO Redirection section of lab 1. (We could also
have used the -i
flag to specify an input file as pw.x -i C_diamond.in
.)
You’ll see a lot of output has been generated in your terminal, but a simple
calculation like this will finish quite quickly. It’s better to explicitly
save the output (and any error information) to a file. To do this, we can
instead run the calculation a redirect the output to a file with
pw.x < C_diamond.in &> C_diamond.log
The output files
Take a look in the directory with ls
. You’ll see some additional files have
been generated. We have pwscf.xml
, which has all the details of the system
and results of the calculation in an xml format. If you skip to the end of
this file (if you have opened it with less
, then pressing G
will go
straight to the end, while g
will go back to the start) you will see the
eigenvalues and other details at each k-point. Note, that while the band
energies listed in the output file are in eV, those in the xml file are in
Hartree units. And we have a folder called pwscf.save
. This has some other
files such as the charge density, another copy of the xml file we saw above, a
copy of the psuedoptential that was used, and files with the calculated
wavefunction at each k-point stored in a binary format that can be used as
input to other calculations.
Now let’s take a quick look through the output that we generated. (e.g.
less C_diamond.log
).
- First we have details of the code, version and date of the calculation.
- There’s some info about the calculation details, including some numbers where we had just accepted the defaults.
- We have our list of automatically generated k-points. We have 10 here since crystal symmetries have been taken into account after generating the 64 points on the 4x4x4 grid we requested.
- Then we go into our self-consistent calculation, starting from randomized atomic wavefunctions. At the end of each step it outputs the total energy.
- After four steps we have converged to the default level (ethr < 1.0E-6).
- The eigenvalues at each k-point are output in eV. Note although we have 8 electrons, we have treated them as 4 doubly occupied states, so only four numbers are output.
- The total energy is listed, as well as its decomposition into several terms.
- And finally some timing information. If calculations are taking a long time this can be helpful in seeing where it is spent.
Calculating other systems
To start with let’s focus on the input parameters that are important for defining our system. There are input files for a number of other types of systems here that we’ll go through.
Silicon
Silicon crystallises in the same diamond structure that we have already seen in the carbon example. So in this case, all we need to do is change the lattice length to the silicon value, and change the atoms from carbon to silicon, which also means we need to specify a silicon pseudopotential rather than use the carbon one.
If you want to see what we’ve changed in the silicon input
file relative to the carbon diamond input
file a useful tool is
diff
. If you’re in the lab02
directory
you can type diff 01_carbon_diamond/C_diamond.in 02_silicon/Si.in
. You’ll be
able to see that we’ve changed four lines in our input file and everything
else is the same. So switching to a system with the same structure but
different atoms can be quite simple. Note we’ve just used the measured
lattice constant for both systems; in later labs we’ll see how to find the
lattice constant predicted by DFT.
Task
- Run
pw.x
using this input and save the output toSi.out
.
Pseudopotentials
You’ll be learning more about pseudopotentials later your lectures. For now
you can think of these as files describing the approximation you want the DFT
code to use for each atomic species. As you’ve seen we set the input files to
look in the current directory for the required pseudopotential file. This is
by setting pseudo_dir = '.'
in the CONTROL
section of the input, where .
represents the current directory as we saw in the last lab. We then use a link
to this file from some central directory rather than making multiple copies of
the pseudopotential for different calculations.
An alternative way to do this would be to specify the central pseudopotential
directly in the input file. We’ve set this in the input file in
02a_silicon/Si.in
.
- Take a look at the directory contents and you’ll see there’s only an input file there.
- Use the
diff
command again to see what’s changed compared to the previous silicon input file. - Try to run the calculation again and confirm that you get the same output as previously.
- Take a look at the pseudopotential file we’ve used for silicon. The header has some useful information regarding how the pseudopotential was generated, such as what states are included, and what approximations are used for exchange and correlation which you will learn more about in your lectures.
Exchange & Correlation Functional
You’ll learn more about the exchange and correlation term in DFT in your
lectures. This is a key part of DFT. The functional we use determines how we
approximate all the many body interactions. By default, within quantum
espresso, the exchange and correlation functionals that are used are taken
from the header of the pseudopotential file as mentioned above. It’s possible
to override this using the input_dft
variable in the system section, but
it’s best to use the same approximation as was used in the pseudopotential
generation.
- In the directory
02b_silion_pbe
, an input file is present that is identical to the previous calculation, except we specify a different pseudopotential here:Si.pbe-rrkj.UPF
. - Look at this file, and compare the header section to the pseudopotential you
used previously. You’ll notice a line mentioning “PBE Exchange-Correlation
functional” in the pseudopotential we’re using here, where it said “PZ
Exchange-Correlation functional previously”.
- Here PZ refers to a particular parametrisation of the Local Density Approximation (LDA) - the simplest approximation for the exchange and correlation, while PBE is more advanced approximation (though not necessarily better performing), which is classed as a Generalized Gradient Approximation (GGA). The letters PZ and PBE are the initials of the authors of the papers in which the particular approximations were published.
Task
- Run
pw.x
with the provided input file. How does the bandwidth of the calculated valence states compare with that of the previous calculation? (You hopefully calculated this at the end of the last lab - check back if you’re not sure what’s looked for).
Carbon Graphite
An input file for carbon graphite is given in
03_carbon_graphite/C_graphite.in
.
Look at the diff
of this input file with the carbon diamond input file.
Now we have more differences than in the silicon case.
- The input
ibrav
has changed. The value of2
is used to indicate an fcc lattice, while4
is used to indicate a hexagonal lattice. Take a look at thepw.x
documentation file for more details on this variable:less ~/MSE404/qe-doc/INPUT_PW.txt
. You can search foribrav
by typing/ibrav
and pressing enter. You can pressn
to go to the next match if you don’t go straight to the main entry. (The input description is also viewable online.) A
has changed, and we have the additional variableC
. For a hexagonal lattice these variables are used to define the three lattice vectors as follows:v1 = A(1,0,0), v2 = A(-1/2,sqrt(3)/2,0), v3 = C(0,0,1)
.- The basis is increased from 2 atoms to 4.
- The basis has changed to a different set of fractional atomic positions.
Carbon Graphene
As you may have noticed, Quantum Espresso is set up to run periodic 3D crystals. So what do we do for systems that aren’t periodic in three dimensions? In this case, we just make the non-periodic dimension large enough that the atoms are too far apart to interact in this dimension, or more accurately that the interaction is sufficiently small relative to the other interactions.
So to move from graphite to graphene, we want to take a single layer of
graphite and extend the interlayer spacing so that the inter-layer interaction
is small. The input for this is in
04_carbon_graphene/C_graphene.in
. diff
the graphite and graphene input files and you’ll see these changes. The other
change that you’ll notice is to the k-point sampling. Since the system we are
trying to calculate is no longer periodic in the z-direction, we only have a
single k-point in this direction.
Task
- Compare the total energy per atom that you have calculated for graphene and graphite. Which is lower? Why do you think that is?
- How does this compare to diamond? Does this seem correct?
Methane
So now we have seen a two-dimensional system, the next question is what do we need to do for molecules or systems with no periodicity. In this case, we just put the entire system in a large box, that we hope (we check) is big enough that periodic images don’t interact, or interact weakly compared to the other interactions we’re interested in.
The input file for a single methane molecule is in
05_methane/CH4.in
. Here by setting ibrav = 1
we have
chosen to use a cubic box, and A = 15.0
makes the box 15 Angstrom along an
edge. And then we define the 5 atomic positions. In this case, it doesn’t
make sense to define the atomic positions in terms of an arbitrarily chosen
box size, so we use the setting ATOMIC_POSITIONS angstrom
to indicate that
the we are giving the atomic positions in Cartesian coordinates in Angstrom.
Note there are a number of equivalent ways of defining the unit cell vectors
within Quantum Espresso. Instead of relying on choosing the correct ibrav
value, we can set ibrav = 0
and give the unit cell lattice vectors
explicitly in a CELL_PARAMETERS
section. This is shown in
05a_methane/CH4.in
.
Task
- Try running
pw.x
with both methane input files and confirm you get the same total energy in both cases.
Aluminium
Aluminium forms in a standard fcc structure with one atom per cell, which we know how to deal with at this point. The thing about Aluminium that makes it more complicated within DFT is that it is a metal. Metals have a Fermi surface that can be quite complex in k-space. This means that in contrast to an insulator or semiconductor where every k-point has the same number of occupied states, in a metal the number of occupied states can vary from k-point to k-point. This makes them a little more difficult to converge than other systems. In short, there are generally two things you need to do:
- Use a denser k-point grid than you would need for a semiconductor or insulator.
- Use some smearing scheme. This means that a state is no longer simply
occupied or empty, but has instead some probability of occupation. In a
semiconductor, if we have a system with 4 electrons, then at every k-point
the first 2 doubly occupied states will be populated. But in a metal,
whether or not a state is occupied depends on its energy and the value of
the Fermi level (which in turn will depend on the energies of all the
electronic states we calculate). Adding a smearing helps significantly in
achieving a smooth SCF convergence, as otherwise a small change in a state
energy from once cycle to the next could lead to a very large change in its
occupation and to the total energy in turn. You will learn more about this
in the lectures. We set the smearing scheme and width with the
occupations
anddegauss
variables in the input file.
Run the pw.x
calculation with the supplied input file in
06_aluminium/Al.in
and take a look through the output.
You’ll need to look in the pwscf.xml
file and find the various ks_energies
entries towards the end. These give the various k-points used in the
calculation and the energies and occupations of each state for this k-point.
Note, for a metal the default number of bands is at least four more than are
needed for the number of electrons per cell. The pseudopotential we have used
has 3 valence electrons, which could be represented with two potentially
doubly occupied bands, so we have four more bands in the calculation for a
total of 6.
- Try removing the
occupations
anddegauss
variables from the input file and see what happens when you try to run the calculation.
Summary
- In this lab we’ve looked at how to create an input file for many different
types of system:
- Carbon in the diamond structure
- Silicon, which has the same diamond structure, but the atomic species are different and the bond length is different.
- Carbon in the graphite structure, where the unit cell is different and the number of atoms per unit cell is different with respect to carbon in the diamond structure.
- Carbon in the graphene structure, where we now have a 2D system rather than a 3D system.
- Methane, which is a molecule in contrast to the previous periodic examples.
- Aluminium, which is a metal and needs some additional inputs to take care of the complications associated with this.
- We’ve looked at how we specify the pseudopotentials for our atomic species.
- We’ve seen an example of using different approximations for exchange and correlation within our DFT calculations.
Extra - Visualising Structures
There are many different tools that can be used to visualize atomic
structures. xcrysden
is installed as a module on the server you’re using for
this course, and conveniently can read Quantum Espresso input files. Try
loading the module with module load xcrysden
, running the command xcrysden
and opening the input files for the various structures we’ve looked at in this
lab. There are many options to control how the structure looks, and you can
grab and rotate the structure with your mouse.
Task
- See if you can figure out how to save an image of each structure as a png file.