Introduction to Nudged Elastic Band (NEB) calculations

This tutorial describes how to use the NEB method to calculate the diffusion barrier for an Au atom on Al(001). If you are not familiar with the NEB method, some relevant references are listed on ASE’s documentation page on the ase.mep.neb module.

The tutorial uses the EMT potential instead of DFT, as this is a lot faster. It is based on a tutorial found on the ASE webpage.

import matplotlib.pyplot as plt
from ase.build import fcc100, add_adsorbate
from ase.constraints import FixAtoms
from ase.calculators.emt import EMT
from ase.optimize import BFGS
from ase.visualize import view
from ase.mep import NEB
from ase.io import read

First we set up the initial state and check that it looks ok:

# 2x2-Al(001) surface with 3 layers and an
# Au atom adsorbed in a hollow site:
slab = fcc100('Al', size=(2, 2, 3))
add_adsorbate(slab, 'Au', 1.7, 'hollow')
slab.center(axis=2, vacuum=4.0)
view(slab)

Then we optimise the structure and save it.

# Fix second and third layers:
mask = [atom.tag > 1 for atom in slab]
slab.set_constraint(FixAtoms(mask=mask))

# Use EMT potential:
slab.calc = EMT()

# Optimise initial state:
qn = BFGS(slab, trajectory='initial.traj')
qn.run(fmax=0.05)

We make the final state by moving the Au atom by one lattice constant and optimise again.

# Optimise final state:
slab[-1].x += slab.get_cell()[0, 0] / 2
qn = BFGS(slab, trajectory='final.traj')
qn.run(fmax=0.05)

Now we make a NEB calculation with 3 images between the inital and final states. The images are initially made as copies of the initial state and the call interpolate() performs a linear interpolation between the initial and final states. As always, we check that everything looks ok before we run the calculation.

NOTE: The linear interpolation works well in this case but not for e.g. rotations. In this case an improved starting guess can be made with the IDPP method.

initial = read('initial.traj')
final = read('final.traj')

constraint = FixAtoms(mask=[atom.tag > 1 for atom in initial])

n_im = 3  # number of images
images = [initial]
for i in range(n_im):
    image = initial.copy()
    image.calc = EMT()
    image.set_constraint(constraint)
    images.append(image)

images.append(final)

neb = NEB(images, method='improvedtangent')
neb.interpolate()
view(images)

After verifying the structure, we can run the optimization:

qn = BFGS(neb, trajectory='neb.traj')
qn.run(fmax=0.05)

We visualize the final path with:

view(images)

You can find the barrier by selecting Tools->NEB in the gui, or you can make a script using ase.mep.neb.NEBTools, e.g.:

from ase.mep import NEBTools

images = read('neb.traj@-5:')

nebtools = NEBTools(images)

# Get the calculated barrier and the energy change of the reaction.
Ef, dE = nebtools.get_barrier()
print('Barrier:', Ef)
print('Reaction energy:', dE)

# Generate new matplotlib axis - otherwise nebtools plot double.
fig, ax = plt.subplots()
nebtools.plot_band(ax)

Exercise

Now you should set up your own NEB using the configuration with N2 lying down as the initial state and the configuration with two N-atoms adsorbed on the surface as the final state.

Parallelisation

The NEB should be parallelised over images. An example can be found in this GPAW tutorial. The script enumerates the CPUs and uses this number (rank) along with the total number of CPUs (size) to distribute the tasks. Note that the NEB class now needs to be initialized by using NEB(images, parallel=True).

# This code is just for illustration
from gpaw.mpi import world
n_im = 4              # Number of images
n = world.size // n_im      # number of cpu's per image
j = 1 + world.rank // n     # image number on this cpu
assert n_im * n == world.size

For each image, we assign a set of CPUs identified by their rank. The rank numbers are given to the calculator associated with this image.

# This code is just for illustration
from gpaw import GPAW, PW
images = [initial]
for i in range(n_im):
    ranks = range(i * n, (i + 1) * n)
    image = initial.copy()
    image.set_constraint(constraint)
    if world.rank in ranks:
        calc = GPAW(mode=PW(350),
                    nbands='130%',
                    xc='PBE',  # student: ...,
                    txt=f'{i}.txt',
                    communicator=world.new_communicator(ranks))
        image.calc = calc
    images.append(image)
images.append(final)

When running the parallel NEB, you should choose the number of CPU cores properly. Let Ncore = N_im * Nk where N_im is the number of images, and Nk is a divisor of the number of k-points; i.e., if there are 6 irreducible k-points, Nk should be 1, 2, 3 or 6. Keep the total number of cores to 24 or less, or your job will wait too long in the queue.

Input parameters

Some suitable parameters for the NEB are given below:

  • Use the same calculator and constraints as for the initial and final images, but remember to set the communicator as described above

  • Use 6 images. This gives a reasonable description of the energy landscape and can be run e.g. on 12 cores.

  • Use a spring constant of 1.0 between the images. A lower value will slow the convergence

  • Relax the initial NEB until fmax = 0.1 eV/Å, then switch on the climbing image and relax until fmax = 0.05 eV/Å.

Once the calculation is done, you should check that the final path looks reasonable. What is the N—N distance at the saddle point? Use NEBTools to calculate the barrier. Is N2 likely to dissociate on the surface at room temperature?