Nudged Elastic Band (NEB)

In this tutorial we will use the climbing-image Nudged Elastic Band method (NEB) to identify the minimum energy path of:

Transition state of a chemical reaction

In this first example, we will find the minimum energy reaction path and transition state of the following reaction:

/scm-uploads/doc.2023/Tutorials/_images/NEB_reaction.png

We will be using the DFTB engine with the GFN1-xTB model. This is a computationally fast method (which is convenient for the purposes of a tutorial) but it is not very accurate in predicting reaction energies. For better accuracy, using the DFT engine ADF (or BAND) is recommended.

Setting up and running the calculation

You can set up and run the calculations for this tutorial using either the Graphical User Interface (GUI), the python library PLAMS or a bash script.

Start AMSjobs
1. in the menu bar select SCM → New input

This will open a new AMSinput window.

In AMSinput:
1. Switch to the DFTB panel: ADFPanel → DFTBPanel
2. Select Task → NEB

We will use the GFN1-xTB model, which is the default model as of the AMS2020 release. Your AMSinput window should look like this:

/scm-uploads/doc.2023/Tutorials/_images/NEB_Main_DFTB_panel.png

When setting up an NEB calculation we need to specify two systems: the initial and final states of the reaction. The NEB algorithm will then generate a set of images by interpolating between the initial and final systems. This will be the initial approximation of the reaction path, which will be optimized during the NEB calculation.

1. Download the following xyz files: file NEB_initial.xyz and NEB_final.xyz
2. Import the coordinates of the initial system: in the menu bar, select File → Import Coordinates (System)… and select the file NEB_initial.xyz
3. Create a new molecule-tab: in the menu bar, select Edit → New Molecule…
4. Import the coordinates of the final system in the newly created molecule tab: in the menu bar, select File → Import Coordinates (System)… and select the file NEB_final.xyz

You can switch between the two molecule tabs by clicking arrows at the bottom of the molecule drawing area.

Important

In NEB calculations, the order of the atoms in the initial and final system should be the same (if you provide an intermediate system, you should use a consistent atom-ordering for that too). The order of the atoms should be consistent because the images-interpolation algorithm maps the n-th atom of the initial system to the n-th atom of the final system.

You can see the indices of the atoms by clicking in the menu bar on View → Atom Info → Name → Show. It is possible to change the order of the atoms in the Coordinates panel (in the panel bar: Model → Coordinates) using the Move atom(s) option.

Now, go to the NEB details panel where we will set up the NEB calculation:

1. Click on MoreBtn next to Task → NEB to go to the NEB details panel
2. From the drop-down menu next to initial system, select Mol-1
3. From the drop-down menu next to final system, select Mol-2
4. Check the Characterize PES point checkbox

Your AMSinput window should look like this:

/scm-uploads/doc.2023/Tutorials/_images/NEB_ready_to_run.png

Tip

From most AMSinput panels you can jump to the relevant section of the user manual by clicking on InfoBtn, which is located in the top-right corner of the panel.

We are now ready to run the calculation:

In AMSinput:
1. In the menu bar, click on File → Save and give it the name “NEB_tutorial”
2. In the menu bar, click on File → Run . This will bring the AMSjobs window to the front
3. Wait for the calculation to finish. It should take just a few seconds

In the logfile you can monitor the progress of your NEB calculation:

In AMSjobs:
1. Select the job “NEB_tutorial” and in the menu bar click on SCM → Logfile. This will open the logfile

A NEB calculation consists of several steps, which are automatically executed one after the other:

  • first, the two end points (the initial and final molecules) are optimized (in the logfile, look for Optimizing initial state and Optimizing final state)

  • then the NEB reaction path will be optimized (in the logfile, look for NEB Path Optimization). During the reaction path optimization, the highest-energy image on the path will climb to the transition state

  • finally, a single point calculation for the TS is performed (in the logfile: Final calculation for highest-energy image). If the option Characterize PES point is on, then the lowest-lying normal modes will be calculated in order to validate the TS (the TS should have exactly one imaginary frequency). Some information on the reaction path is printed at the end of the logfile:

    NEB found a transition state!
TS barrier height from the left           0.02576078 Hartree
                                         16.165 kcal/mol
                                         67.635 kJ/mol
TS barrier height from the right          0.08632064 Hartree
                                         54.167 kcal/mol
                                        226.635 kJ/mol

To improve the initial approximation of the reaction path in an NEB calculation, you can (optionally) provide an intermediate system.

Another important NEB option is the number of images. In case of problematic NEB path optimization convergence, using more images might help (note that the computation time increases with the number of images used).

You can read more about the various NEB options in the AMS manual.

Results of the calculation

Now, let’s visualize the reaction path computed by NEB:

In AMSjobs:
1. Select the job “NEB_tutorial” and, in the menu bar, click on SCM → Movie. This will open the AMSmovie module
/scm-uploads/doc.2023/Tutorials/_images/NEB_AMSmovie.png

In AMSmovie, you can click on play (or drag the slide-bar) so see the reaction happening. On right-hand side, the energy and gradients of the images in the NEB reaction path are plotted.

The transition state is characterized by having one imaginary frequency. Let’s visualize the normal modes of the transition state geometry with AMSspectra:

In AMSmovie:
1. click on SCM → Spectra. This will open the AMSspectra
/scm-uploads/doc.2023/Tutorials/_images/NEB_AMSspectra.png

Here you will see the computed normal modes for the TS geometry. Notice that there is one negative frequency (imaginary frequency are shown as negative numbers).

In AMSspectra:
1. In the table, click on the line with the negative frequency

The corresponding normal mode will be displayed in the molecule-visualization area. This normal mode gives you an idea of how the atoms are moving as they cross the transition state.

Ionic diffusion in a solid

In this second example, we will determine the minimum energy path and the activation energy of Li diffusion in a typical cathode materials LiTiS2.

We will use the ML potential engine with the M3GNet universal machine learning potential. Although M3GNet was trained to reproduce ground state structures of the Materials Project database, it can also be accurate to describe transition states, as we will demonstrate in this tutorial. This tutorial can be reproduced with the periodic DFT engine BAND.

First, make sure to install the M3GNet package inside AMS as described in this tutorial.

Layer LiTiS2 cathode material

In some Li-based cathode materials, the diffusion of Li ions is tightly connected to the content of Li in the materials. We would like to demonstrate this by computing the diffusion barrier of Li in LiTiS2 with 1 and 2 Li vacancies. In layer LiTiS2, the lithium ions are octahedrally coordinated between the S−Ti−S sandwich layers.

Li diffusion in LiTiS2

To perform the NEB calculation we will use the graphical user interface and perform a 3-step calculation: initial state, final state, and NEB.

Step 1: Initial state

Download the crystal structure of layer LiTiS2.

Start AMSjobs
1. in the menu bar select SCM → New input

This will open a new AMSinput window. Let’s perform a geometry and cell optimization of the initial structure.

In AMSinput
1. Switch to the ML Potential panel: ADFPanel → MLPotentialPanel
2. Select M3GNet potential ANI-2x → M3GNet-UP-2022
3. Click File → Import System.. and navigate to the crystal file you just downloaded named LiTiS2_pristine.cif
4. Click the MoreBtn next to Task → Geometry Optimization
5. Check the Optimize lattice option
6. Click Details → Expert AMS and scroll down to ExpertAddons D3Dispersion and check the box next to Enabled

Note

When using M3GNet it is highly recommended to include dispersion via the addon option. In the following, we will build the final state, and the NEB calculation from the initial state calculation and therefore the D3 dispersion will be included throughout.

Save the job (File → Save As..), give it a name (LiTiS2_relax), and run the job Job → Run. This should take a few minutes. When asked to update coordinates with bonds and data say yes.

/scm-uploads/doc.2023/Tutorials/_images/NEB_relax.png

We can now delete the final position of the Li atom. Orient the crystal along the x-axis (ctrl+1) and select the middle Li atom of the top layer (index 25) and tap the delete key to remove it. Let’s relax only the ions.

1. Click the MoreBtn next to Task → Geometry Optimization
2. Uncheck the Optimize lattice option

Save the job (File → Save As..), give it a name (LiTiS2_init), and run the job Job → Run. When asked to update coordinates with bonds and data say yes. Save the new job (File → Save As..) as LiTiS2_final and jump to the next section.

/scm-uploads/doc.2023/Tutorials/_images/NEB_delete.png

Step 2: Final state

First, re-orient the crystal along the y-axis (ctrl+2) and select the middle Li atom of the top layer (index 1). With the atom selected, right click on the atom and drag the atom to its final position (the position of the atom we previously deleted).

/scm-uploads/doc.2023/Tutorials/_images/NEB_final.png

You can do Bonds → Guess Bonds afterward in order to re-evaluate the bonds.

Note

When you perform NEB in a crystal it is recommended to keep the cell dimension between the initial and final states.

Save and run the calculation without changing the settings to perform a geometry optimization of the final state. Check that the energy of the initial and final state are almost identical.

Step 3: NEB

We now have the initial and finial state geometries and we can proceed to the NEB calculation. We will modify the initial state.

1. From AMSjobs, click the ML button left of the LiTiS2_init job
2. Say yes when asked to update bonds and data
3. Save the calculation as LiTiS2_NEB

We can now define the settings for the NEB calculation.

1. From the Main tab, select the Task → NEB
2. Click the MoreBtn next to Task → NEB
3. Set the Number of images to 17
4. Select Mol-1 from the drop-down menu Initial system
5. Select New molecule from the drop-down menu Final system

Notice that the molecule area has been cleared. What happened is that AMS created a new empty molecule, named Mol-2 by defaults. You can switch between the initial state Mol-1 and the final state Mol-2 by clicking the small black triangles at the bottom of the molecule area. There are several ways to define the final system. We will copy and paste the final system we have previously prepared. Keep the LiTiS2_NEB input window opened, with the empty molecule area Mol-2.

1. From AMSjobs, click the ML button left of the LiTiS2_final job
2. Say yes when asked to update bonds and data
3. Hold the shift key and drag the mouse to select all the atoms of the final system
4. Click Edit → Copy (or ctrl+c)
5. Select the window of the LiTiS2_NEB input and paste the final system to the empty molecule area Mol-2 (Edit → Paste or ctrl+v)

Note

Alternatively, you can save the coordinates of the final state and load it to Mol-2.

The NEB is now ready. Save the calculation. This will automatically rename the molecules as initial and final.

/scm-uploads/doc.2023/Tutorials/_images/NEB_settings.png

Run the calculation. This might take a few more minutes as AMS will perform a simultaneous optimization of all the 17 images by minimizing the forces parallel to the reaction path. Once the calculation is finished, select LiTiS2_NEB in AMSjobs and click SCM → Movie. You can re-orient the crystal and use the slidebar under the molecule area to appreciate the diffusion along the minimum energy path. If you double click the y-axis label of the graph, you can switch to your prefered units. We found the activation energy EA = 0.59 eV.

/scm-uploads/doc.2023/Tutorials/_images/NEB_movie.png

Li diffusion in LixTiS2

You can now reproduce the 3-step calculation by removing one Li neighbor to appreciate the effect of Li concentration on its diffusion. We found that with a neighbor Li vacancy, the activation energy is decreased by approximately a factor of 2. An interesting note is that while the trajectory of the Li migration is linear when all Li neighbor sites are occupupied, the trajectory path becomes curved, passing throught a tetrahedral site (resulting in a shallow minima), when a neighbor Li is vacant. You can also further reproduce these calculations with BAND to validate the accuracy of M3GNet. In the end, we found that M3GNet reproduces well the diffusion pathway of Li in LixTiS2. This can be explained by the large variety of Li-Ti-S crystals in the Materials Project database spaning numerous Li environments.

/scm-uploads/doc.2023/Tutorials/_images/NEB_plot.png