Electronic Structure Calculations of Bulk Crystals Using Quantum ESPRESSO

Tutorial Created with Software Release: 2025-1
Topics: Catalysis & Reactivity, Energy Capture & Storage, Metals, Alloys & Ceramics, Thin Film Processing
Methodology: Periodic Quantum Mechanics
Products Used: MS Maestro, Quantum Espresso

Tutorial files

78 MB
This tutorial is written for use with a 3-button mouse with a scroll wheel.
Words found in the Glossary of Terms are shown like this: Workspacethe 3D display area in the center of the main window, where molecular structures are displayed

 

Tip: You can hover over a glossary term to display its definition. You can click on an image to expand it in the page.
Abstract:

 

In this tutorial, we will learn the basics of the Quantum ESPRESSO (QE) interface for periodic density functional theory (DFT) calculations of bulk solids, including geometry optimization, electronic band structures, density of states (DOS), projected density of states (PDOS), and convergence testing.

 

Tutorial Content

  1. Introduction to Quantum ESPRESSO

  1. Creating Projects and Importing Structures

  1. Optimizing Atomic Positions and the Crystal Cell

  1. Calculating Electronic Properties of Relaxed Structures

  1. Analyzing Electronic Properties

  1. Running Convergence Tests for Energy Cutoffs and K-point Grids 

  1. Analyzing Convergence Calculations    

  1. Conclusion and References

  1. Glossary of Terms                

1. Introduction to Quantum ESPRESSO

Quantum ESPRESSO (QE) is a plane-wave pseudopotential DFT code that uses fully periodic boundary conditions. The plane-wave basis set and periodic boundary conditions make QE especially valuable for quantum-level simulations of 3D ordered solid materials and other periodic materials, including 2D or 1D periodic materials such as surfaces or nanomaterials. Simulations of molecules and 0D materials are sometimes appropriate as well, typically for exploring interactions of these materials with periodic systems.

Figure 1: Example of a wavefunction in a 1D crystal, illustrating the extended nature of wavefunctions in crystals. The wavefunction in fact extends infinitely in the ±x directions. The green circles represent atom sites on a periodic lattice.

In this tutorial, we will cover the most straightforward case: a bulk crystalline material. First, we will optimize the crystal structure for the rutile and anatase phases of TiO2 at a quantum mechanical level. Then, we will analyze the electronic properties of the relaxed structures, including electronic band structures, density of states (DOS), and projected density of states (PDOS). At the end, we will begin to learn about convergence testing to find a balance between accuracy and computational expense for our simulations.

For users new to periodic DFT calculations we recommend reading the references in Section 8 to understand the most important concepts in solid state electronic structure theory.

For users who are new to using Materials Science Maestro with Quantum ESPRESSO, see the Installing and Configuring Quantum ESPRESSO guide and the Quantum ESPRESSO Calculation Panel help documentation. It is important to note that Quantum ESPRESSO jobs are not supported on Windows or on Mac operating systems. To run the calculations in this tutorial, a Linux host is required. For those who are unable to or choose not to run the calculations, the input and output files are provided (in Section 2) and referenced throughout the tutorial.

For a follow-up tutorial involving surface and slab calculations with QE, see Modeling Surfaces.

2. Creating Projects and Importing Structures

At the start of the session, change the file path to your chosen Working Directorythe location where files are saved in MS Maestro to make file navigation easier. Each session in MS Maestro begins with a default Scratch Projecta temporary project in which work is not saved, closing a scratch project removes all current work and begins a new scratch project, which is not saved. A MS Maestro project stores all your data and has a .prj extension. A project may contain numerous entries corresponding to imported structures, as well as the output of modeling-related tasks. Once a project is saved, the project is automatically saved each time a change is made.

Structures can be built in MS Maestro or can be imported using File > Import Structures (or drag-and-dropped), and are added to the Entry Lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion and Project Tabledisplays the contents of a project and is also an interface for performing operations on selected entries, viewing properties, and organizing structures and data. The Entry Lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion is located to the left of the Workspacethe 3D display area in the center of the main window, where molecular structures are displayed. The Project Tabledisplays the contents of a project and is also an interface for performing operations on selected entries, viewing properties, and organizing structures and data can be accessed by Ctrl+T (Cmd+T) or Window > Project Table if you would like to see an expanded view of your project data.

  1. Double-click the Materials Science icon

Figure 2-1. Change Working Directory option.

  1. Go to File > Change Working Directory
  2. Find your directory, and click Choose
  3. Pre-generated input and results files are included for running jobs or examining output. Download the zip file here: schrodinger.com/sites/default/files/s3/release/current/Tutorials/zip/qe_basic_rutile.zip
  4. After downloading the zip file, unzip the contents in your Working Directorythe location where files are saved for ease of access throughout the tutorial

Figure 2-2. Save Project panel.

  1. Go to File > Save Project As
  2. Change the File name to QE_bulk_crystal_tutorial and click Save
    • The project is now named QE_bulk_crystal_tutorial.prj

Figure 2-3. Importing structure file.

  1. Go to  File > Import Structures
  2. Navigate to where you saved the tutorial files and Open TiO2_input.mae
    • Two new entries titled rutile and anatase are selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries in the Workspacethe 3D display area in the center of the main window, where molecular structures are displayed and the first is includedthe entry is represented in the Workspace, the circle in the In column is blue

Note: Please refer to the Glossary of Terms for the difference between includedthe entry is represented in the Workspace, the circle in the In column is blue and selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries.

Figure 2-4. The two entries  in the entry list.

  1. The imported entries are now available

3. Optimizing Atomic Positions and the Crystal Cell

In this tutorial, we will be calculating and analyzing the electronic structure of bulk rutile and anatase. The crystal files imported in Section 2 are experimental X-ray structures. In order to calculate the electronic structure of these materials, we need to first optimize their crystal structures at the quantum mechanical level. In this section, we will relax the crystal structure using DFT, ensuring that it is prepared to undergo further calculations. The optimization will adjust the atomic positions until the interatomic forces and stresses on the unit cell are zero (within specified tolerances). Here, we will use parameters for the wavefunction cutoff energy and k-point grid that are chosen on the basis of convergence testing. More about convergence testing and how these parameters were chosen can be found in Section 6.

Figure 3-1. Selecting and including rutile and anatase crystal structures.

  1. Select(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries and includethe entry is represented in the Workspace, the circle in the In column is blue both rutile and anatase from the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion (Shift + click)
  2. Open the Workspace Configuration Panel 
  3. Click on Tile ()
    • The two crystal structures are tiled in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed

Figure 3-2. The asymmetric unit cell of the rutile TiO2  in the Workspace as imported from the cif file. Grey: Ti; Red:O.

  1. Change the representation to ball-and-stick by clicking on the Style menu () and choosing Apply ball-and-stick representation
  2. In Color Atoms, choose Element from the dropdown menu

Note: The rutile and anatase structures includedthe entry is represented in the Workspace, the circle in the In column is blue in the Workspacethe 3D display area in the center of the main window, where molecular structures are displayed depict the asymmetric unit cell (as described in the .mae file). Using the information about the space group, which is also provided in the given .mae file, the full unit cell can be  generated as explained in the following step.

Figure 3-3. Above: Show periodic structure tool settings. Below: The full unit cell of the rutile and anatase. TiO2 : Grey: Ti; Red: O.

  1. Click Show periodic structure tool window () in the footer bar
  2. Hover over Build Cell in the pop-up menu
  3. Select Translate to First Unit Cell, Recalculate Connectivity, and Recalculate Bond Orders
  4. Click Apply
    • The full unit cell is generated

Note: By choosing Extents and selecting values greater than 1, one can visualize the extended unit cell in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed

For a complete overview of crystals and periodic boundary conditions, the Building and Manipulating Crystal Structures tutorial is recommended

Figure 3-4. Selecting rutile and anatase in the entry list.

  1. Open the Workspace Configuration Panel 
  2. Click on Tile ()
    • The two crystal structures are no longer tiled in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed
  3. Ensure that rutile and anatase are both selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries in the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion. Includethe entry is represented in the Workspace, the circle in the In column is blue rutile in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed. The geometry optimization we will perform in the next few steps will take the selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries structures as inputs for the calculation. This way, we can run calculations on rutile and anatase simultaneously.

Figure 3-5. Setting up the geometry optimization.

We are now ready to optimize the geometry of both anatase and rutile using Quantum ESPRESSO

  1. Go to Tasks > Materials > Quantum Mechanics > Quantum ESPRESSO > Quantum ESPRESSO Calculations
  2. Ensure that Use structures from shows Project Table (2 selected entries)
  3. Ensure that Convergence tests is unchecked
  4. In Geometry select the Optimize atomic positions and cell radio button
    • This option allows relaxation of the positions of the atoms and the unit cell parameters to a stable structure with no residual forces. In this case, no constraints are applied
  5. Uncheck all properties in the Property section of the panel
  6. Click Pseudopotentials

Figure 3-6. Pseudopotentials panel settings.

  1. In Path click Browse
  2. Navigate to where you saved the tutorial file and choose the all_pbe_UPF_v1.5 directory
    • The Pseudopotentials panel is automatically filled with the paths to the directory and to the Ti and O specific pseudopotential files
  3. Click OK to close the Pseudopotentials panel
    • A message appears that the energy cutoffs were not defined. Click OK to close the message as we will define the cutoffs in the next step

Note: Pseudopotentials are not included with the Quantum ESPRESSO installation and must be downloaded separately. There are several libraries available (see this link for more information). The pseudopotentials used in this tutorial can be downloaded from GBRV pseudopotential site

Note: Pseudopotentials are usually generated for a specific type of density functional. It is therefore recommended to employ the same type of density functional, if possible, in all QE calculations that was used in generation of the pseudopotentials. The functional used for the generation of the pseudopotential is stated in the Pseudopotentials panel (see the Figure)

Figure 3-7. Advanced Options, Theory tab settings.

We will now select specifics for the level of theory for our geometry optimization

  1. In the Quantum ESPRESSO Calculations panel, click Advanced Options
  2. In the Theory tab, ensure the following settings are used:
    • Choose GGA for Density functional type
    • Choose PBE for Density functional
    • Check Use symmetry
    • Uncheck Use primitive cell

Note: Regarding the primitive unit cell option, QE always calculates the geometry shown in the selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries or includedthe entry is represented in the Workspace, the circle in the In column is blue entries except in one special case, when the “Use Primitive” option is checked. With this option QE first uses symmetry analysis to establish whether the smaller (primitive) unit cell exists for the selected entry. If such a cell exists the calculations will be performed for the primitive cell

  • In Brillouin zone partition options, select Grid plane distance and set it to 0.05
  • Check Include Γ-point

Note: Rutile has a tetragonal unit cell, i.e. the lattice constant c in z-direction is smaller than the lattice constants a/b in x/y direction (c=2.9587 Å versus a/b= 4.5937 Å). Therefore, we need more k-points in the z-direction to obtain the same k-point density in all 3 directions of the reciprocal cell. By specifying the grid-plane distance, such an evenly-distributed k-point grid is automatically constructed. Therefore, this option is very useful for non-cubic cells

Note: Include Γ-point is selected to also include the Γ-point ((0,0,0) in the reciprocal cell) in the k-point mesh

  1. Confirm that the rest of the Theory settings match those that are shown in the Figure, and then click on the SCF tab

Note: By clicking Update K-point mesh with a single entry includedthe entry is represented in the Workspace, the circle in the In column is blue in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed, you can see the Monkhorst-Pack k-point grid corresponding to the chosen grid plane distance settings for that structure

Figure 3-8. Advanced Options, SCF tab settings.

  1. Set the Custom energy cutoff for wavefunction to 60 Ry
    • We select the value of 60 Ry because it is shown to converge the energy of the system in Section 6
  2. Set the Custom energy cutoff for charge density to 300 Ry

Note: The ultrasoft pseudopotentials used in this tutorial have been optimized for a 5 times higher charge density cutoff than wavefunction cutoff. However, depending on the type of pseudopotentials the multiplier recommended to determine the charge density cutoff could be different

  1. Click Save

Figure 3-9. Running the geometry optimization.

  1. For the Job name, input cell_optimization
  2. Adjust the job settings () as needed
    • This job requires a CPU host. The job can be completed in less than 30 minutes on a CPU host
  3. If you would like to run the job, click Run. Otherwise, we will proceed with pre-generated results in the next step
  4. Once the job is successfully finished, a new cell_optimization1 (2) group, with 2 entries, is selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries in the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion and the first structure is includedthe entry is represented in the Workspace, the circle in the In column is blue in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed

Figure 3-10. DFT-relaxed rutile in the workspace.

  1. If you did not run the optimization yourself, from the main menu, choose File > Import Structures
  2. Navigate to where you downloaded the tutorial file and choose the Section_03 > cell_optimization > cell_optimization.maegz
  3. Click Open
    • A new entry group is added to the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion. The first entry is automatically selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries and includedthe entry is represented in the Workspace, the circle in the In column is blue in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed

After the job is finished, you should compare the output with the initial structure (which was an experimental X-ray structure). The DFT-relaxed structure has slightly longer cell vectors than the experimental structure for both anatase and rutile. Different computational methods can lead to slightly different optimized geometries. For example, GGA functionals like the PBE used in this tutorial tend to slightly overestimate bond lengths and cell vectors of inorganic materials, whereas LDA functionals are known to underestimate them. However, DFT-relaxed structures are typically within 1% of experimental structures, even considering the caveats discussed here. The small difference between the experimental and computed structures here indicates that the DFT-relaxed structures are good models to use for further calculations.

 

Note also that for materials with significant Van der Waals (VdW) interactions (such as molecular crystals) the use of dispersion corrections or Van der Waals functionals is strongly suggested.

4. Calculating Electronic Properties of Relaxed Structures

We are now ready to explore the electronic properties of the DFT-relaxed anatase and rutile structures. In this section, we will calculate the density of states (DOS), projected density of states (PDOS), and band structure for the two crystal phases of TiO2. The properties calculated here can help us better understand the electronic properties of the materials and the band gap can be compared to experiment. For information about computing other properties such as phonons, elastic constants, and NMR shifts, see the help documentation.

Figure 4-1. QE Calculations panel settings.

  1. Ensure that the DFT-relaxed anatase and rutile structures are both selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries, then open the Quantum ESPRESSO Calculations panel again using the search function of the Task menu
    • Alternatively, go to Tasks > Materials > Quantum Mechanics > Quantum ESPRESSO > Quantum ESPRESSO Calculations
  2. For Geometry, select Frozen
    • For the property calculation, we are using the structures that we optimized in Section 3. These optimized geometries should be frozen during the property calculations
  3. For Property, check Density of states, Projected density of states, and Band structure
  4. Click Advanced Options
    • The Advanced Options panel opens

Figure 4-2. Advanced Options, (P)DOS tab settings.

  1. In the Theory tab, ensure the following settings are used:
    • Check Use symmetry
    • Check Use primitive cell
    • Retain all other selections from  Section 3
  2. Switch to the (P)DOS tab
  3. In Brillouin zone partition options, select Grid Plane distance and set the distance to 0.025
    • This is half of the grid plane distance we used for the geometry optimization. This corresponds to a k-point grid that is twice as dense as the one used during the relaxation
  4. Check Include Γ-point
  5. Click Save

Note: Here we used a k-point grid for the DOS that is two times denser than the original one. In order to capture the rapid variations in the conduction band, a tight k-point sampling is necessary for DOS and band structure calculations. Since the conduction band does not affect the total energy of the system, structural relaxations are often done at a comparably looser grid, whereas for the calculation of electronic properties, like the (P)DOS, a tighter grid is chosen (see Kratzer P and Neugebauer J (2019) The Basics of Electronic Structure Theory for Periodic Systems. Front. Chem. 7:106. DOI:10.3389/fchem.2019.00106).

Figure 4-3. Running a property calculation.

  1. For the Job name, input properties
  2. Adjust the job settings () as needed
    • This job requires a CPU host. The job can be completed in less than 30 minutes on a CPU host
  3. If you would like to run the job, click Run from the Quantum ESPRESSO Calculations panel. Otherwise, we will proceed with pre-generated results in the next section
  4. Close the Quantum ESPRESSO Calculations panel

Having computed the electronic properties of two bulk crystals, tools for analyzing them are presented in the next section.

5. Analyzing Electronic Properties

This section presents graphical interface tools for analyzing the electronic structure of periodic systems.  In this section, we will use pre-generated results.  Alternatively, you can use the electronic properties that you calculated in Section 4 above.

Figure 5-1. Viewing the output of the property calculation.

  1. From the main menu, choose File > Import Structures
  2. Navigate to where you downloaded the tutorial file and choose the properties.maegz
  3. Click Open
    • A new entry group with 10 structures is added to the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion. There are five entries corresponding to SCF, DOS, PDOS, DOS Plotting, and Band calculations respectively for rutile and anatase. The group is automatically selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries and the first entry includedthe entry is represented in the Workspace, the circle in the In column is blue in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed

5.1 Using the Density of States Viewer

The DOS plot shows the energies at which the electronic states are located. States at energies below the Fermi level are occupied (or bound) states and those above the Fermi level are unoccupied states.  It is usual to refer to the highest-energy band of occupied states as the valence band and the lowest-energy band of unoccupied states as the conduction band.

Figure 5-2. Launching the Density of States Viewer panel from the WAM.

Now, let's evaluate the DOS for rutile and then compare it to that of anatase

The most convenient method to access the Density of States Viewer panel is by using the Workflow Action Menu (WAM) button which appears next to the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion

 

  1. With the properties (10) entry group selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries, click on the WAM button and select Density of States Viewer
    • The Density of States Viewer panel opens
    • Alternatively, go to Tasks > Materials > Quantum Mechanics > Quantum ESPRESSO > Density of States Viewer

Figure 5-3. Density of States Viewer panel.

  1. Check Zero to Fermi energy
    • This will set the zero of energy on the horizontal axis to the Fermi energy. The actual value of the Fermi energy is reported in the status bar at the foot of the panel
  2. Check Plot integrated DOS

Figure 5-4. Opening multiple instances of the Density of States Viewer panel.

We can also take a side-by-side look at the DOS of anatase and rutile

  1. To view and compare the DOS of anatase to rutile, Ctrl+Click (Cmd+Click) the entry titled properties..anatase...dos_plot_3 in the properties (10) entry group
  2. In the Density of States Viewer Panel click Open Viewers for Selected Entries
    • Two viewer panels open, one with the DOS of rutile and the other with the DOS of anatase
  3. Repeat Steps 5 and 6 for each viewer panel

 

You can continue on to examine the DOS plots, change the axes limits, and save the figure. See the help documentation for more.

  1. Close the Density of States Viewer panels

The zero density of states at the Fermi level tells us that both of our materials are insulators.

The band gap is the energy difference between the highest occupied state of the valence band and the lowest unoccupied state of the conduction band. The Fermi level and the band gap are printed at the bottom of the panel. We will discuss the band gap further in Section 5.3.

5.2 Using the Projected Density of States Viewer

While the DOS plot gives us information about the total density of states, analyzing the projected/partial DOS helps us to understand how the individual atoms, or even orbitals, contribute to the total DOS.

Figure 5-5. Launching the Projected Density of States Viewer panel from the WAM.

The most convenient method to access the Projected Density of States Viewer panel is by using the Workflow Action Menu (WAM) button which appears next to the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion

  1. With the properties (10) entry group selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries, click on the WAM button and select Projected Density of States Viewer
    • The Projected Density of States Viewer panel opens
    • Alternatively, go to Tasks > Materials > Quantum Mechanics > Quantum ESPRESSO > Projected Density of States Viewer
    • The entry associated with the projected density of states is selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries by default in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed when the panel is opened using the WAM button

Figure 5-7. Initial settings of the Projected Density of States Viewer.

We will now analyze how individual atoms and orbitals contribute to the DOS. To do this, we will view the contributions of an O-2p orbital and Ti-3d orbital for rutile. We will first create sets of atoms for which we want to view the PDOS, Ti and O. Then we will plot the contributions of the 3d and 2p orbitals of Ti and O respectively.

  1. Change Selected atoms to Atom set Ti
  2. In the workspacethe 3D display area in the center of the main window, where molecular structures are displayed, in the atom selection menu under Quick Select, click on the three dots and choose Metal Atoms. All of the Ti atoms in the selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries structure will be highlighted
    • You can also manually select all Ti atoms by shift-clicking the Ti atoms in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed
  3. Back in the Projected Density of States Viewer, you will see Ti1 and Ti4, the selected atoms, listed under Selected atoms. Click Add to List

Figure 5-8. Projected Density of States Viewer panel final settings.

  1. Repeat steps 12 to 14 for O atoms (name the Selected Atoms Atom Set O)
    • Use the Invert button in the atom selection menu to choose all non-Ti atoms quickly
  2. Check Zero to Fermi energy
  3. In Atomic PDOS, select Atom set Ti
  4. In Atomic Orbitals PDOS, select Atom set Ti (n = 3, l = 2) and Atom set O (n = 2, l = 1) simultaneously using Ctrl+Click (Cmd+Click)
    • You have now plotted the contributions of an O-2p orbital and Ti-3d orbital for rutile

The PDOS can tell us about the relative contribution of a particular atom or atomic orbital to the total DOS. Looking at the PDOS of rutile, we see that the valence band is dominated by an O-2p orbital (Atom set O, n=2, l=1,  depicted in green), while the conduction band is dominated by the Ti-3d orbital (Atom set Ti, n=3, l=2, depicted in orange). Repeat the steps taken in this section for anatase if interested. The PDOS for anatase is similar to that of rutile, the major contributions to the valence and conduction bands are from the O-2p and Ti-3d orbitals respectively.

5.3 Using the Band Structure Viewer

The electronic band structure of a solid shows the allowed bands of energy levels in reciprocal (k) space, including special high symmetry k-points. Band structures reveal information about band gaps, carrier mobility, and other electronic properties. We can obtain the magnitude and type of band gap (i.e. direct or indirect) from the band structure diagram. Here, we will look at the band structures for anatase and rutile.

Figure 5-9. Launching Band Structure Viewer panel from the WAM.

  1. Close the Projected Density of States Viewer panel
  2. With the properties (10) entry group selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries, click on the WAM button and select Band Structure Viewer
    • The Band Structure Viewer panel opens
    • Alternatively, go to Tasks > Materials > Quantum Mechanics > Quantum ESPRESSO > Band Structure Viewer

Figure 5-10. Band Structure Viewer panel settings.

  1. In Band Structure Viewer panel, check Zero to Fermi energy
    • The panel records that the Fermi level is at 9.634 eV, which we have shifted to 0 eV on our plot, and you can see that this locates the Fermi level in the gap between bands, meaning that rutile is an insulator
    • The panel also records that the band gap is 1.853 eV, again indicating that rutile is an insulator, and that it has a direct band gap at the Gamma point

Note: The Fermi level and the band gap are printed at the bottom of the panel as shown in the Figure

Figure 5-11. Band Structure Viewer panel.

  1. Repeat Step 20 and 21 with the properties (10) entry group and select the band structure of anatase
    • The panel records that the Fermi level is at 8.290 eV, which we have indicated as 0 eV on our plot, and that the indirect band gap is 2.110 eV. This indicates that anatase is an insulator as well
  2. Close the Band Structure Viewer panels

The experimentally measured band gaps for rutile and anatase are 3.03 eV and 3.20 eV, respectively (Scanlon et al. Nature Materials volume 12, pages 798–801(2013)), so our calculated values of 1.853 eV and 2.109 eV are lower than experiment by almost 40%. This underprediction of the band gaps is expected: GGA functionals, and local functionals in general, usually severely underestimate band gaps because of the approximations made at this level of theory. Accurate band gaps can be obtained using corrections such as GGA+U or more computationally expensive hybrid functionals. However, these more advanced techniques are often unnecessary, since GGA functionals are accurate for other properties of materials such as optimized geometries and reaction energies.

 

Rather than focus on the absolute values of the band gaps, it is often more useful to observe trends in the band gap between systems and other qualitative features of the band structures. In these calculations, we have correctly determined that anatase has a larger band gap than rutile. Furthermore, features of the band structure are predicted correctly: rutile has a direct band gap and anatase has an indirect band gap. For anatase, our calculations show the conduction band minimum  at the Gamma point and the valence band maximum near the X point, in agreement with angle-resolved photoemission spectroscopy data (Baldini, E., Chiodo, L., Dominguez, A. et al. Strongly bound excitons in anatase TiO2 single crystals and nanoparticles. Nat Commun 8, 13 (2017)). The nature of the band gaps for these materials has implications for their performance in applications such as photocatalysis.

 

Note: we have defined the band gap in our calculations as the difference between energies of the highest occupied and lowest unoccupied Kohn-Sham states, consistent with convention in the literature.

6. Running Convergence Tests for the Energy Cutoffs and K-point Grids

In this section, we will run convergence tests to investigate the accuracy of the total energy calculation with respect to integration grids in reciprocal (k) space and with respect to size of plane wave basis for wavefunctions and charge density. The convergence of total energy will allow us to make an appropriate choice for the k-point grid and wavefunction/density energy cutoff parameters.

We now provide some background on the significance of converging the calculations with respect to these parameters. If you are new to periodic DFT calculations, please see additional references in Section 8 to learn more.

According to Bloch’s theorem, in an infinite periodic crystal, the wavefunction can be written as the product of a plane wave and a factor, , which is periodic with regard to the Bravais lattice.

Eq. 1:

can be expressed as an infinite Fourier sum over the reciprocal lattice vectors, .

Eq. 2:

Therefore, the problem of finding the wavefunction in the plane wave formalism is reduced to finding the expansion coefficients and of the plane wave basis expansion (Eq. 2). Note, that the sum in Eq. 2 is infinite. In any practical calculation it has to be truncated, which means truncating the plane wave basis via defining an energy cutoff. The representation of the charge density is truncated in a similar way.

Once these coefficients are defined, any observable (such as energy) can be expressed as the expectation value of the corresponding operator, , as follows, where indicates 3D integration over the Brillouin Zone.

Eq. 3:

Similarly, in order to evaluate the integral over the Brillouin Zone in Eq. 3, the integral has to be approximated by a sum over a finite number of appropriately-chosen points within the Brillouin Zone that we call the k-point grid.

Thus, the accuracy of the calculations critically depends on the wavefunction energy cutoff and the k-point grid. For a more detailed description  please refer to the standard solid state textbooks such as Kratzer P and Neugebauer J (2019) The Basics of Electronic Structure Theory for Periodic Systems. Front. Chem. 7:106. doi: 10.3389/fchem.2019.00106 and Density Functional Theory: A Practical Introduction David Sholl, Janice A Steckel, Wiley, 2009, Chapter 2-4.

Recall in Section 3, when we optimized the geometry of the structure, we selected a value for the wavefunction energy cutoff and the k-point grid. These values were determined by running the following convergence tests for the unoptimized rutile cell. Now, let's run the convergence tests so we can select reasonable values for the wavefunction energy cutoff and the k-point grid parameters in Section 7. This tutorial will go through the process of running and analyzing the convergence tests for rutile. If you are interested, feel free to perform a similar analysis for anatase.

Figure 6-1. QE Calculations panel settings.

For the convergence tests, we will use the experimental crystal structure of rutile imported in Section 2. You may carry out convergence tests on any other reasonable structure, but in that case the exact results will differ from what is presented here

  1. With rutile includedthe entry is represented in the Workspace, the circle in the In column is blue and selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries in the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion, go to Tasks > Materials > Quantum Mechanics > Quantum ESPRESSO > Quantum ESPRESSO Calculations
    • The Quantum ESPRESSO Calculations panel opens
    • Ensure that the initial rutile structure that was imported into the project is selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries and not the DFT-relaxed output
  1. Ensure that Use structures from shows Project Table (1 selected entry)
  2. Check the Convergence tests box
  3. Click Advanced Options

Figure 6-2. Advanced Options, Theory tab settings.

  1. In the Theory tab, ensure the following are selected:
    • GGA for Density functional type
    • PBE for Density functional
    • Uncheck Use primitive cell

Figure 6-3. Advanced Options, SCF tab settings.

  1. Switch to the SCF tab
  2. For Custom energy cutoff for wavefunctions, input 40
  3. For Custom energy cutoff for charge density, input 200

Note: The values of the wavefunction and charge density cutoffs will be overwritten by the values provided in the Convergence tests tab

Figure 6-4. Advanced Options, Convergence tests tab settings.

  1. Switch to the Convergence tests tab
  2. Make the following selections:
    • Minimum energy cutoff to 40 Ry for wavefunctions
    • Maximum energy cutoff to 80 Ry for wavefunctions
    • Step size to 10 Ry
    • Charge density multiplier to 5
    • Check Grid plane distance, with the following settings:
    • Max: 0.12
    • Decrement:  0.01
    • Steps: 10
    • Check Include Γ-point
  3. Click Save to return to the Quantum ESPRESSO Calculations panel

Figure 6-5. Convergence job settings.

  1. Change the job Name to convergence_rutile
  2. Adjust the job settings () as needed
    • This job requires a CPU host. The job can be completed in approximately 30 minutes on a CPU host

Note: with the current settings the workflow will perform a number of fixed cell and frozen atom QE calculations with the energy cutoff and k-points value ranges specified in the Convergence Test subpanel. Each calculation will be performed on 1 thread with a maximum of 5 calculations run simultaneously. So, the workflow will request 5 CPU cores 

Note: The QE parallel options -npools must be a divisor of Threads to evaluate several k-points at once. Other systems might require different parallelization settings, for more information check this help topic

  1. If you would like to run the job, click Run. Otherwise, we will proceed with pre-generated results in the next section
  2. Close the Quantum ESPRESSO Calculations panel

7. Analyzing Convergence Calculations

By analyzing the convergence calculations, we can select reasonable values for the wavefunction energy cutoff and the k-point grid parameters to use for further calculations. In this section, we will use pre-generated results.  Alternatively, you may skip importing the results and use those that you have computed in Section 6.

Figure 7-1. Importing Structures.

  1. From the main menu, choose File > Import Structures
  2. Navigate to where you downloaded the tutorial file and choose convergence_rutile.maegz file
  3. Click Open
    • A new group with 50 structures is added to the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion. The entire group is automatically selected(1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries and the first entry includedthe entry is represented in the Workspace, the circle in the In column is blue in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed
    • Maintain the entire entry group selection

Figure 7-2. Launching Convergence Tests Viewer panel.

  1. From Tasks in the toolbar, search for convergence
    • Alternatively, go to Tasks > Materials > Quantum Mechanics > Quantum ESPRESSO > Convergence Tests Viewer
  2. Choose Convergence Tests Viewer from the list

Figure 7-3. Loading data from Workspace.

The Convergence Tests Viewer panel will be used to analyze the results of the convergence calculations

  1. Click Load data from Workspace
    • This loads in the entry group from the previous calculation
  2. Check Convert energy to eV
  3. Check Show relative energy
    • Energy values are now displayed relative to the lowest energy value in the data set

Figure 7-4. Convergence with the wavefunction cutoffs.

  1. In Wavefunction cutoffs (Ry) column, shift-click all entries
    • The graph shows convergence with k-points for the different wavefunction cutoffs. We observe a fast convergence with the number of k-points already after 18 k-points. The convergence is slightly improved by going up to 24 k-points which we will use for our example. This means that we can use the 5 x 5 x 7 k-point grid (or 4 x 4 x 6), or the corresponding grid plane distance of 0.05 /Å, without significant loss of accuracy for the calculation of the total energy

Figure 7-5. Convergence with the wavefunction cutoff for a given k-point grid

  1. In # of k-points column, click 24 (5,5,7)
    • The graph shows convergence with wavefunction cutoff for our chosen k-point grid (24 in total, 5 x 5 x 7 k-point grid)  

 

 

Let us examine the results of the convergence test. Figure 7-4 depicts the total energy as a function of k-point grid density. Using the Zoom to rectangle tool to get a closer look at the convergence, we observe that with k-point sampling 5 x 5 x 7 the energy per unit cell is converged to below 5 meV. However, the converged energy value strongly depends on wavefunction energy cutoff. Choosing this k-point sampling and examining energy dependence on wavefunction cutoff (Figure 7-5), we observe that at the cutoff value 60 Ry the total energy is converged to below 20 meV/cell or ≤ 3 meV/atom. This accuracy is sufficient for our calculations. However, depending on the type of calculation/property of interest, you may require a tighter energy convergence.

One might naively expect a k-point grid of 5 x 5 x 7 to give a total number of 175 k-points, which is more than the 24 k-points actually reported in this calculation. The reduction in k-points is achieved by exploiting the lattice symmetry and sampling the k-points only inside the irreducible part of the Brillouin Zone using the Monkhorst-Pack approach. For more information, see Kratzer P and Neugebauer J (2019) The Basics of Electronic Structure Theory for Periodic Systems. Front. Chem. 7:106. doi: 10.3389/fchem.2019.00106.

Based on the convergence test, we choose a wavefunction cutoff of 60 Ry, a charge density cutoff of 300 Ry, which is 5 times the wavefunction cutoff and a k-point grid plane distance of 0.05 /Å which corresponds to a k-point mesh of 5 x 5 x 7 in the rutile unit cell. This is why we use these parameters in Sections 3 and 4 of the tutorial.

If you are interested, feel free to perform a similar analysis for anatase. You will find that the parameters we choose for rutile are sufficient for use with anatase as well.

8. Conclusion and References

In this tutorial, we learned the basics of the Quantum ESPRESSO (QE) interface for periodic density functional theory (DFT) calculations of bulk solids.

For further learning:

For introductory content, focused on navigating the Schrödinger Materials Science interface, an Introduction to Materials Science Maestro tutorial is available. Please visit the materials science training website for access to 70+ tutorials. For scientific inquiries or technical troubleshooting, submit a ticket to our Technical Support Scientists at help@schrodinger.com. 

For self-paced, asynchronous, online courses in Materials Science modeling, including access to Schrödinger software, please visit the Schrödinger Online Learning portal on our website.

For some related practice, proceed to explore other relevant tutorials:

For further reading:

For the interested user without much experience in periodic DFT calculations, we recommend taking a look into the following references to complement the tutorial. Please be aware that those references are suggested in the context of this tutorial. The list is far from being complete and there are many more review articles and textbooks available that discuss solid state DFT.

9. Glossary of Terms

Entry List - a simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion

Included - the entry is represented in the Workspace, the circle in the In column is blue

Project Table - displays the contents of a project and is also an interface for performing operations on selected entries, viewing properties, and organizing structures and data

Recent actions - This is a list of your recent actions, which you can use to reopen a panel, displayed below the Browse row. (Right-click to delete.)

Scratch Project - a temporary project in which work is not saved, closing a scratch project removes all current work and begins a new scratch project

Selected - (1) the atoms are chosen in the Workspace. These atoms are referred to as "the selection" or "the atom selection". Workspace operations are performed on the selected atoms. (2) The entry is chosen in the Entry List (and Project Table) and the row for the entry is highlighted. Project operations are performed on all selected entries

Working Directory - the location where files are saved

Workspace - the 3D display area in the center of the main window, where molecular structures are displayed