Activation Energies for Reactivity in Solids and on Surfaces

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

Tutorial files

25 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 to locate the transition state for the reaction of a small molecule on a surface via the nudged elastic band method. Periodic quantum mechanics (specifically Density Functional Theory, DFT) is used to compute the activation energy at the transition state, which dictates the reaction rate in chemical kinetics.

 

Tutorial Content

  1. Introduction to Nudged Elastic Band Calculations

  1. Catalyzing Ammonia Oxidation with Platinum 

  1. Creating Projects and Importing Structures 

  1. Running Nudged Elastic Band Calculations

  1. Z-matrix Interpolation

  1. Conclusion and References

  1. Glossary of Terms

1. Introduction to Nudged Elastic Band Calculations

The transition stateor activated complex, corresponds to a first order saddle point on the potential energy surface linking the initial reactant and final product (TS), or activated complex, corresponds to a first order saddle point on the potential energy surface linking the initial reactant and final product. Locating a transition state is essential for computing the activation energythe energy difference between the transition state of a reaction and the ground state of the reactants of a reaction, and thereby the reaction rate. To learn more about transition states and their importance in chemical kinetics, please visit the Locating Transition States: Part 1 tutorial for gas-phase, molecular systems. The Nudged Elastic Band Calculations panel is used to locate the transition state for an elementary reactiona reaction for which no reaction intermediates have been detected or need to be postulated in order to describe the chemical reaction on a molecular scale. An elementary reaction is assumed to occur in a single step and to pass through a single transition state in a periodic system, such as a surface, given an initial and final configuration. In Schrödinger’s Materials Science (MS) Maestro with Quantum ESPRESSO, we can locate the transition state and evaluate the activation energy for the reaction of a molecule on a surface following a straightforward workflow as summarized here:

The Nudged Elastic Band (NEB) approach is commonly used to calculate activation energies for elementary chemical reactions including but not limited to: adsorption/desorption of small molecules on surfaces, gas diffusion at metal surfaces, and diffusion of atoms between vacancy sites in a bulk crystal. First, the user must provide structures of the initial reactant and final product connected via an elementary reaction step - there should only be one transition state between the reactant and product. Hidden intermediates are a common reason for problems in locating transition states so it is best to use your chemical knowledge to choose the initial and final structures wisely. It is best practice to provide geometry optimized structures, however they can be optimized during the NEB calculation if preferred. Then, candidate structures (‘images’) that are approximately on the reaction path are generated through interpolation between the provided initial and final structures. A sufficient number of images must be generated on this path to span the path and reliably locate the transition state. However, there is a balance to be struck between accuracy and computational expense, as more images require more computational time.

By performing simultaneous DFT optimizations of the images, with their geometries constrained through fictional forces (the ‘elastic band’) to a continuous sequence from reactant to product, one of the images should converge to the transition state along that path. Each image also feels a force orthogonal to the path, the true force, which is minimized. In this method, images are not optimized independently; instead, the force on each image depends on its neighbors. At every step, the force components parallel to the reaction path are replaced by a "spring force" that tries to keep each image between its neighbors and does not allow the images to optimize to the nearest minima i.e. the initial or final reaction state. Thus, the total force is the sum of the parallel component of the spring force and the perpendicular component of the true force, where the latter is mandatory for an accurate prediction of the minimum energy path.

NEB calculations allow us to locate transition states and predict activation energies, which can be useful for many applications, including catalysis, thin film processing, and energy storage.

In this tutorial, we will study one elementary step in the dissociation of ammonia, *NH3 → *NH2 + *H (where asterisks indicate surface-bound states), while adsorbed on the Pt(111) surface using the Nudged Elastic Band Calculations panel. The dissociation reaction of interest is described in more detail in Section 2. We will begin with optimized structures of the reactant and product of the dissociation reaction on the Pt(111) surface. If you are unfamiliar with modeling the adsorption of small molecules on a surface, reviewing the Modeling Surfaces tutorial is available and recommended. For a more in-depth explanation of transition states and other tools available to locate them within MS Maestro, see the Locating Transition States: Part 1 and 2 tutorials.

2. Catalyzing Ammonia Oxidation with Platinum

Ammonia (NH3) oxidation reactions play an essential role in removing NH3 from waste systems and in the production of nitric acid, which has a wide range of applications in the chemical industry. For these reasons, NH3 oxidation on platinum (Pt) is a reaction of interest to scientists across chemical industries.

In this tutorial, we will consider a Pt(111) surface as the heterogeneous catalyst for NH3 decomposition (*NH3 → *NH2 + *H). As a surface model we choose (√7 × √7) Pt(111), which denotes a supercell with an extent of √7 in the a- and b-directions of the (111) crystallographic plane of platinum. The (√7 × √7) system is chosen because it is neither too small (which would cause periodic interactions between the adsorbates) nor too large (which would significantly increase the computational costs of the NEB calculations). The image below of the product structure (NH2 + H) for a  (2 × 2) Pt(111) and (√7 × √7) Pt(111) surface model demonstrates our first point. These structures are available in the provided tutorial files under the optional directory. The dissociated H is colored orange in order to help visualize its position. In the smaller (2 × 2) cell, we see that the distance wise, the dissociated H atom is very close to NH2 and its periodic image, which will lead to anomalies in the computed energies as H and NH2 move apart from each other. Increasing the cell size to (√7 × √7) gives the ammonia molecule space to dissociate fully which allows us to avoid such artificial correlations between periodic images during the NEB calculation.  

Supercell of Product Structure for (2 × 2) Pt(111)

Supercell of Product Structure for (√7 × √7) Pt(111)

The NEB calculation requires two minima as input (reactant and product), between which we will search for the single transition state of the elementary reaction that links these two minima. The two input geometries must consist of the same atoms, arranged in different geometries, so that a one-to-one correspondence of atoms can be established for the reaction. 

Though not explicitly demonstrated in this tutorial, determining whether two minimum-energy structures are linked is non-trivial. In particular, to understand how the reactant (NH3) and products (NH2 + H) bind to the Pt(111) surface, we need to compare the energetics of the structures associated with the different possible adsorption sites for each adsorbate on the surface (Ammonia Dissociation on Pt{100}, Pt{111}, and Pt{211}:  A Comparative Density Functional Theory Study. DOI:10.1021/jp073083f). The possible adsorption sites are top, fcc (face centered cubic), hcp (hexagonal close-packed) and bridge. Previous DFT calculations revealed an energetic preference for the NH3 molecule to adsorb on the top site. On the other hand, the NH2 molecule prefers the bridge site over the top site by 0.47 eV. The H atom dissociated from the NH3 → NH2 + H reaction is 0.03 eV more stable on the bridge site than on any other site. Therefore, we use the NH3 (top) and NH2 (bridge) + H (bridge) as our initial and final structures, respectively. For a review of the methods used to construct and optimize the provided minima, revisit the Modeling Surfaces tutorial.

Having optimized the minima, we can compare the DFT energies of the product and the reactant and we find that adsorbed NH2 + H is 0.66 eV higher in energy than adsorbed NH3. This means that the dissociation reaction is thermodynamically unfavorable (endothermic in terms of enthalpy at zero Kelvin).

3. 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 3-1. Change Working Directory option.

  1. Go to File > Change Working Directory
  2. Find your directory, and click Choose
  3. Pre-generated files are included for running jobs or examining output. Download the zip file here: schrodinger.com/sites/default/files/s3/release/current/Tutorials/zip/neb.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 3-2. Save Project panel.

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

Figure 3-3. Entry list and workspace after importing.

We will import geometry optimized initial and final structures for our reaction of interest, the dissociation of ammonia on Pt(111). Let’s import the structures now:

  1. Go to File > Import Structures
  2. Select input_structures.mae
  3. Click Open
    • A new entry group containing three entries is added to the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion

4. Running Nudged Elastic Band Calculations

In this section, we will run an NEB calculation for the dissociation of NH3 on Pt(111) using the reactant and product structures imported in Section 3. We will generate the structures that are expected to be near the reaction path (‘images’) through linear interpolation and then run the NEB calculation. For an example of Z-matrix interpolation, see Section 5. Finally, we will analyze the outputs of the calculation: the activation energy at the transition state and the transition state geometry.

Figure 4-1. Opening the Nudged Elastic Band Calculations panel. 

We will start by inputting structures and generating images for the ammonia dissociation reaction pathway using the Nudged Elastic Band Calculations panel:

 

  1. Go to Tasks > Materials > Quantum Mechanics > Quantum ESPRESSO > Nudged Elastic Band Calculations

 

Note: Many features found at the bottom of this panel are similar to those of the Quantum ESPRESSO Calculations panel. To learn more about these features, please visit the Electronic Structure Calculations of Bulk Crystals Using Quantum ESPRESSO tutorial or the help documentation.

Figure 4-2. Loading the initial and final structure.

First, let's load our geometry-optimized initial and final structures into the panel.

  1. In the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion, ensure initial_structure_NH3 and final_structure_NH2_H_bridge-reordered 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
  2. Back in the Nudged Elastic Band Calculations panel, click the Load from Workspace button
    • Both structures are loaded into the images table

Figure 4-3. Setting the parameters for interpolation.

Next, let’s generate structures of images that are linearly interpolated between the initial and final state.

Note: Consistent atom ordering between the initial and final structure is necessary for a successful calculation. The Reorder Atoms button allows you to reorder the atoms in the final structure to match the ordering in the initial structure. In this section, the provided initial and final structures are already ordered appropriately. In Section 5, we will practice reordering atoms for a different system.

  1. Select the two rows in the images table corresponding to the initial and final structure using Shift + Click
  2. For Number of structures, enter 5
    • We will generate five images between the reactant and product structures; note that more images may be needed for NEB convergence in the case of more complex reactions
    • Use of the Z-Matrix option will be covered in Section 5 below

Figure 4-4. Interpolating images and exporting them to the workspace.

  1. Click Interpolate Images
    • 5 images are added in between the initial and final structures in the images table.
    • The order of the images in the table corresponds to the path from reactant (first image) to product (last image).
  2. We can visualize the images generated by exporting to the workspace. Click Export Images to Workspace
    • A new entry group titled periodic_dft_1-images (7) 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 Nudged Elastic Band Calculations panel is not synced with the workspace. Structures from the entry list can be loaded into the panel using the Load Selected Entries button while structures generated in the panel, through interpolation, can be copied to the workspace using the Export Images to Workspace button. Any changes made to structures in the workspace will not be reflected in the panel unless the edited structures are loaded into the panel.

Best practice is to check the images to ensure that a reasonable pathway for the reaction of interest has been generated. For instance, sometimes interpolation leads to unphysically close atoms in some images. Note that image structures may be edited in the workspace manually and re-loaded into the panel. Images can also be added or deleted in the panel.

 

In the Figure below, we see that the seven total structures (reactant  + 5 images + product) show a reasonable initial pathway for the NH3 molecule to dissociate into NH2 and H. The dissociating H atom is close enough to the Pt(111) surface by the sixth structure for a Pt-H bond to be visible. Feel free to further visualize the interpolated structures.

In general, images can be interpolated linearly, using the Z-matrix, or using IDPP. In this section, we linearly interpolated images between the initial and final structure in Step 6. Use of the Z-Matrix option will be covered in Section 5. Images can also be generated between any two contiguous structures. For example, if you have a guess for the transition state, the initial, transition state guess, and final structure can be loaded into the panel. Images can be interpolated between the initial structure and transition state as well as the transition state and final structure.

 

It is important to note that no quantum mechanical calculations have been performed at this stage. The images will be used as input for the subsequent calculation.

Figure 4-5. Opening the Atomic Constraints Dialog Box.

Prior to NEB, the bare Pt(111) slab was generated from bulk Pt with optimization of atom positions. It is important to use consistent constraints to those used when optimizing the reactant and product structures. We will now use Atomic Constraints to practice freezing the bottom two layers of Pt atoms in the slab during NEB. We leave the first two layers unfixed in these optimizations to allow these Pt atoms to relax in response to the dissociation of ammonia on the surface. As is often the case, the bottom layers are already constrained from previous calculations; here we will clear these restraints to practice using the panel. If you are familiar with setting atomic constraints, feel free to skip ahead to Step 15.

 

  1. In the Nudged Elastic Band Calculations panel, click on Atomic Constraints
  2. In the Atomic Constraints Dialog Box, click on Delete All
    • We are deleting any pre-existing constraints on the structure

Figure 4-6. Freezing the bottom two layers of the slab.

  1. For Type, choose Cartesian - XYZ from the drop-down
    • The atom will be frozen in the x, y and z directions
  2. In the workspacethe 3D display area in the center of the main window, where molecular structures are displayed, 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 all of the Pt atoms in the bottom two layers of the slab (Pt19-Pt32), 14 atoms in total. This can be achieved in different ways including manual selection of each atom or highlighting the area in the workspace and holding Shift
  3. Click Add Workspace Selection
    • The panel populates with the selected atoms
  4. Click Save

Figure 4-7. Opening the Pseudopotential panel.

  1. Back in the Nudged Elastic Band Calculations panel, click on Pseudopotentials
  2. For the top Path, click Browse

Figure 4-8. Choosing the pseudopotential directory.

  1. Navigate to where you saved the tutorial file and choose the uspp_gbrv_all_pbe_UPF_v1.5 directory

 

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 or the Schrödinger Github repository

Figure 4-9. Closing the Pseudopotential panel.

The Pseudopotentials panel is automatically filled with the paths to the directory and atom specific pseudopotential files

 

  1. Click OK to close the Pseudopotentials panel

Figure 4-10. Opening the Advanced Options panel.

  1. Click on Advanced Options

 

If you would like to learn more about the Advanced Options, visit the Electronic Structure Calculations of Bulk Crystals Using Quantum ESPRESSO tutorial or the help documentation

Figure 4-11. Setting parameters in the Theory tab.

  1. For Density functional type, choose GGA
  2. For Density functional, choose PBE
  3. For Dispersion correction, choose None
  4. Uncheck Use symmetry
  5. For the Brillouin zone partition options, select Monkhorst-Pack and set it to 5 x 5 x 1
    • A high density of k-points is required for accurate calculation of metals. A higher density of k-points will increase the computational expense
    • However, slabs should only have one k-point in the z-direction
  6. Go to the SCF tab

Figure 4-12. Setting parameters in the SCF tab.

  1. Set the Custom energy cutoff for wavefunctions to 40 Ry
  2. Set the Custom energy cutoff for charge density to 320 Ry
  3. Set the Max steps in SCF to 100
  4. Set the Energy convergence threshold to 1e-06
  5. For Occupations, select Smearing and choose Marzari-Vanderbilt from the drop-down
    • We use Marzari-Vanderbilt as it is a preferred method for smearing in metallic systems

 

Note: The parameters in the Theory and SCF tabs are chosen to match those used to optimize the bulk structure, bare slab, and reactant and product structures. The level of theory used to compute absolute energies must be the same for the reactant, images, and product in order to generate a single reaction energy profile for the transformation.

Figure 4-13. Setting parameters in the NEB tab.

  1. Go to the Custom NEB tab
  2. Enter 0.1000 for the Force norm threshold
    • The convergence threshold is defined as the norm of the force orthogonal to the path. The simulation stops when the norm is smaller than the threshold value entered here
  3. Increase the Number of steps to 200
    • We increase the maximum number of steps to lower the chance that the calculation will fail due to non-convergence
  4. Ensure Climbing image box is checked
  5. Click Save
    • All of our Advanced Options are saved in the Nudged Elastic Band Calculations panel

Figure 4-14. Running the NEB calculation.

  1. Change the Job name to NEB_NH3_Pt
  2. Due to the long calculation time for this job, we will proceed with the provided files. Go to File > Import Structures, navigate to the provided tutorial files and Open  Section_04 > NEB_NH3_Pt > NEB_NH3_Pt.maegz
    • This job requires a CPU host and must be run on a Linux machine. The job can be completed in 3 days on a CPU host with 8 threads and 8 subjobs
  3. Close the Nudged Elastic Band Calculations panel

Figure 4-15. Visualizing the output.

  1. After importing, 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 the NEB_NH3_Pt1 (7) entry group and includethe entry is represented in the Workspace, the circle in the In column is blue the first entry, NEB_NH3_Pt_1
    • Feel free to visualize the output system in the workspace, such as changing the representation style to ball-and-stick
    • It is possible for the unit cell to become misaligned with the periodic structure; this is purely visual, but can be easily adjusted using the Periodic Structure Tool Window if needed

Figure 4-16. NEB convergence monitor.

We can visualize the energy convergence for our calculation.

  1. Go to Tasks > Materials > Quantum Mechanics > Quantum ESPRESSO > Nudged Elastic Band Convergence Monitor
  2. Click Select and navigate to the provided tutorial files and Open  Section_04 > NEB_NH3_Pt > NEB_NH3_Pt.csv

Figure 4-17. Viewing the error convergence for the NEB calculation.

The error convergence plot shows that the forces were converged for all five intermediate images. The dashed line is the force convergence criteria.

  1. Click Energy profile

Figure 4-18. Viewing the energy profile for the NEB calculation.

We see that we have a well defined reaction pathway with one transition state (image 2)

In this example, our calculation successfully converged. However, if the calculation runs into any convergence issues the convergence monitor is very useful. See the help documentation for more information on the NEB Monitor panel.

Figure 4-19. Viewing the Project Table.

We can also view properties of the image structures using 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.

 

  1. Use the Property Tree to add the NEB Energy (eV) (under All > Materials Science > Primary > …)

The Nudged Elastic Band Monitor panel monitors the progress of jobs that have been submitted from the Nudged Elastic Band Calculations panel as they are running or once they are completed. This panel helps one decide how to proceed if the jobs appear to be problematic. For example, if the calculation reaches the maximum number of optimization steps but convergence is not achieved, you may want to restart the job with a larger number of steps. Alternatively, if the energy profile shows two transition states, you may want to separate the calculation into two NEB jobs.

 

Another useful tool to view NEB results is the Reaction Profile Viewer panel. In this panel, you can plot an energy diagram of a multistep reaction that shows how the energy changes through the course of the reaction. Let’s use our pre-optimized reactant and product molecules as well as our located transition state to plot the overall reaction energy profile using the Reaction Profile Viewer panel.

Figure 4-20. Importing the reactant in the Reaction Profile Viewer.

  1. Go to Tasks > Materials > Quantum Mechanics > Tools > Reaction Profile Viewer
  2. Select Quantum ESPRESSO at the top of the panel and keep the dropdown as Total energy
  3. Click on Reactant and then Delete
    • The table is cleared of all rows. We will import all entries at once in the next few steps
  4. In the entry list, 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 the NEB_NH3_Pt1 (7) entry group

Figure 4-21. Importing steps in the Reaction Profile Viewer.

The table on the left takes the various components of each step as input. We will have seven steps: reactant, images, transition state and product. We will import all the steps at once here. Note that multistep reactions can be plotted here as well

  1. Ensure the NEB_NH3_Pt1 (7) entry group is selected
  2. Click Import to New Steps
    • All 7 structures are pulled into the table. The first structure is labeled Reactant and the remaining 6 are labeled Product
  3. Double-click Product for the 5 optimized images and change the Title to IMG (an abbreviation for Image)
    • Your table in the panel should match that which is shown in the Figure

Figure 4-22. Viewing the reaction energy profile.

  1. Change the Reaction Energy units to eV
  2. Click Plot
    • The energy diagram is updated to display the reaction energy profile

Figure 4-23. Labeling the transition state.

From this reaction profile, we can see that the geometry with the highest energy (1.29 eV) is the saddle point or transition state (TS) along the reaction pathway.  We can indicate this finding by labeling the reaction profile as follows:

  1. Change the Title of the second image to TS for Transition State
  2. Click Plot
    • The energy diagram is updated to display the reaction energy profile

Feel free to adjust the plot as you please. A plot with only the reactant, TS, and product can also be generated

After the NEB calculation ends, we can visualize the optimized images for the NH3 dissociation.

 

The transition state (structure shown above) can be thought of as the moment when the N-H bond breaks, which leads to the NH2 still near the initial adsorbate position and the free H from the N-H bond near its final position. The DFT energy of the transition state (TS) relative to the reactant is the activation energy of the forward reaction, +1.29 eV. The energy of the TS relative to the product is the activation energy of the backward reaction and is also provided as a result of the calculation as 0.66 eV.

5. Z-matrix Interpolation

The Nudged Elastic Band method is based on structures called ‘images’ that approximately span the reaction path between reactant and product. The closer the generated images are to the actual reaction pathway, the faster the calculation will converge. Linear interpolation of Cartesian coordinates (x, y, z) is the simplest way to generate images and is a good approximation to reactions where atoms move in straight lines. Periodic boundary conditions are taken into account in the linear interpolation scheme.

However reactions sometimes involve the rotation of entire groups of atoms, which is not treated realistically by the linear algorithm. Instead, interpolation can be performed relative to the bond lengths, bond angles and torsional degrees of freedom of the molecule, which together are referred to as ‘internal coordinates’ and are represented by the Z-matrix. Using Z-matrix interpolation for the reactive atoms of a system can therefore produce more chemically-reasonable images that converge faster to the TS. The remaining atoms, including any that cross the periodic boundaries, should be interpolated linearly (the Z-matrix does not recognize periodic boundary conditions).

In either case, you can inspect and edit the interpolated structures to ensure that they are reasonable before starting the NEB calculation.

In this section, we will demonstrate how Z-matrix interpolation can provide more reasonable structures than linear interpolation for the rearrangement of a formate molecule, HC(=O)O─, from monodentate to bidentate configurations on a copper surface. Before doing so, we will need to use the Reorder Final Structure Atoms dialog box to order the atoms in the initial and final structure.

Figure 5-1. The entry list and workspace after importing.

We will import geometry optimized initial and final structures for our reaction of interest, the rearrangement of HC(=O)O- on copper. Let’s import the structures now:

  1. Go to File > Import Structures
  2. Select Section_05 > HCOO_Cu.mae
  3. Click Open
    • A new entry group containing three entries is added to the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion

Figure 5-2. Opening the Nudged Elastic Band Calculations panel.

Once again, we will proceed to generate images for the reaction pathway using the Nudged Elastic Band Calculations panel:

 

  1. Go to Tasks > Materials > Quantum Mechanics > Quantum ESPRESSO > Nudged Elastic Band Calculations

Figure 5-3. Loading the initial structure.

First, let's load our geometry-optimized initial and final structures into the panel

 

  1. In the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion, ensure HCOO_side_initial and HCOO_top-final 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
  2. Back in the Nudged Elastic Band Calculations panel, click the Load Selected Entries button
    • Both structures are added to the images table

Figure 5-4. Opening the Reorder Final Structure Atoms Dialog Box.

The atoms in the initial and final structures loaded here are not in the same order. We must reorder the atoms before we can generate images. It is good practice to check if the atoms are ordered to ensure no errors arise in setting up or running the calculation

 

  1. Select the final structure in the images table
  2. Click on Reorder Atoms

Figure 5-5. Visualizing the reference and comparison atoms in the workspace.

  1. Uncheck Automatically map hydrogens and monatomic elements
  2. Check Pick reference and comparison atoms in the workspace
  3. Arrange the panel and MS Maestro window in a way that you can see the structures in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed and the panel simultaneously

Figure 5-6. Selecting atoms in the Reference section.

Let’s begin with the formate molecule as the corresponding atoms are easy to spot on both structures. The crucial mapping here is that of the two oxygen atoms since those are binding to the copper surface. You can get two different reactions (with different activation energies) if you map the oxygen atoms in different ways.

 

  1. Scroll to the bottom of the Reference and Comparison sections
  2. Click on O 50 in the Reference section
    • The corresponding atom in the Comparison section is guessed and shown in gray
    • The guess in the Comparison section and the workspacethe 3D display area in the center of the main window, where molecular structures are displayed is incorrect. The H 52 atom is selected instead of the corresponding O. Let’s remedy that in the next step

Figure 5-7. Selecting atoms in the Comparison section.

  1. Click on O 50 in the Comparison section
    • The atoms that are mapped are now highlighted in green
  2. Click on O 51 in the Reference section

Figure 5-8. Mapping formate.

  1. Click on O 51 in the Comparison section

 

Continue to map the atoms comprising the formate molecule in the Reference and Comparison sections by iterating between steps 13 and 14

Figure 5-9. Mapping Cu atoms.

Now, let's move onto mapping the copper atoms. Continue to map the atoms by iterating between steps 16 and 17 for Cu as seen in the Figure.

 

Use the Hide Mapped Atoms button as needed to hide the atoms that have already been mapped. This will allow for a clearer view of the remaining atoms to be mapped.

 

Figure 5-10. Guessing how to map atoms.

  1. After mapping a sufficient number of atoms, such as ~20 for this example, click on Guess
    • The panel makes a guess on how additional atoms are mapped and highlights them in blue

Figure 5-11. Accepting a guess.

  1. Click Accept Guess if the guess is reasonable
    • The rows highlighted in blue are now mapped and highlighted in green

Figure 5-12. Finishing the atom mapping process.

  1. Continue to map the atoms manually and by using the Guess button until all 52 atoms are mapped
  2. Click OK

Figure 5-13. Viewing the reordered final structure.

The reordered 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 and also replaces the previous entry as the Final structure

Figure 5-14. Generating linearly interpolated structures.

First, for comparison, let’s generate images that are linearly interpolated between the initial and final state:

 

  1. Select the two rows in the images table corresponding to the initial and final structure using Shift + Click
  2. For Number of structures, enter 5
    • We will generate five structures along the reaction pathway
  3. Change the Job name to linear_interpolation
  4. Click Interpolate Images
  5. Click Export Images to Workspace
    • A new entry group titled linear_interpolation-images (7) is added to the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion. Feel free to visualize the images

Figure 5-15. Deleting the linearly interpolated images.

Next, let’s generate structures using Z-matrix interpolation:

 

  1. Before being able to interpolate new images, we must delete the old ones. Select the 5 linearly interpolated images
  2. Click Delete Selected Images

 

Figure 5-16. Opening the Atom Selection Dialog Box for Z-matrix interpolation.

  1. Check Use Z-matrix
  2. Next to with ASL, click Select

Figure 5-17. Selecting atoms.

We can now define the atoms that will be moving in the reaction during our NEB calculation

 

  1. At the bottom of the panel, click Clear
  2. Back in the workspacethe 3D display area in the center of the main window, where molecular structures are displayed, select all of the atoms that make up the formate molecule (Ctrl or Cmd + Click) as shown in the Figure

Figure 5-18. Adding the atoms to ASL.

  1. Click Add
  2. Click OK

 

Figure 5-19. Generating structures by Z-matrix interpolation.

  1. Back in the NEB Calculation panel, change the Job name to z-matrix_interpolation
  2. Select the two rows in the images table corresponding to the initial and final structure using Shift + Click
  3. Click Interpolate Images
  4. Click Export Images to Workspace
    • A new entry group titled z-matrix_interpolation-images (7) is added to the entry lista simplified view of the Project Table that allows you to perform basic operations such as selection and inclusion

Let’s look at the generated structures in more detail. Below are the initial and final optimized structures along with the five unoptimized intermediate structures as generated by linear interpolation. 

We can see that linear interpolation does not properly capture the rotation of the formate molecule, instead causing atoms to move very close to each other as they swap locations. If we were to run an NEB calculation using these initial structures, the DFT steps would likely fail to converge. We can check if the Z-matrix interpolated structures are any better in the Figure below.

Here we see that the interpolated structures are better suited to find the transition state. By comprehending rotation as one of the internal coordinates of the HC(=O)O- molecule, the Z-matrix approach handled this reaction in a much more reasonable manner than the linear interpolation. The images generated by Z-matrix interpolation would give us a better starting point for an NEB transition state search. Having generated images in this way, the procedure for running and analyzing the NEB calculation is the same as in Section 3.

6. Conclusion and References

In this tutorial, we learned how to generate images for, run and analyze a Nudged Elastic Band Calculation for locating a transition state and computing the associated activation energy. The sample systems were reactions of small molecules on a surface.

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 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.

7. Glossary of Terms

Activation Energy - the energy difference between the transition state of a reaction and the ground state of the reactants

Elementary Reaction - a reaction for which no reaction intermediates have been detected or need to be postulated in order to describe the chemical reaction on a molecular scale. An elementary reaction is assumed to occur in a single step and to pass through a single transition state

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

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

Transition State - or activated complex, corresponds to a first order saddle point on the potential energy surface linking the initial reactant and final product

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