Lab 2: Quantum Espresso Input and Output

Back to Course Overview

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 to Si.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 of 2 is used to indicate an fcc lattice, while 4 is used to indicate a hexagonal lattice. Take a look at the pw.x documentation file for more details on this variable: less ~/MSE404/qe-doc/INPUT_PW.txt. You can search for ibrav by typing /ibrav and pressing enter. You can press n 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 variable C. 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:

  1. Use a denser k-point grid than you would need for a semiconductor or insulator.
  2. 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 and degauss 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 and degauss 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.