Dissociation of a molecule using the NEB method

In this tutorial we provide an illustrative example of a nudged-elastic band (NEB) calculation. For more information on the NEB technique, see ase.mep.neb. We consider the dissociation of a nitrogen molecule on the Cu (111) surface.

The first step is to find the relaxed structures of the initial and final states.

from ase import Atoms
from ase.build import add_adsorbate, fcc111
from ase.calculators.emt import EMT
from ase.constraints import FixAtoms
from ase.io import write
from ase.optimize import QuasiNewton

# Find the initial and final states for the reaction.

# Set up a (4 x 4) two layer slab of Cu:
slab = fcc111('Cu', size=(4, 4, 2))
slab.set_pbc((1, 1, 0))

# Initial state.
# Add the N2 molecule oriented at 60 degrees:
d = 1.10  # N2 bond length
N2mol = Atoms('N2', positions=[[0.0, 0.0, 0.0],
                               [0.5 * 3**0.5 * d, 0.5 * d, 0.0]])
add_adsorbate(slab, N2mol, height=1.0, position='fcc')

# Use the EMT calculator for the forces and energies:
slab.calc = EMT()

# We don't want to worry about the Cu degrees of freedom,
# so fix these atoms:

mask = [atom.symbol == 'Cu' for atom in slab]
slab.set_constraint(FixAtoms(mask=mask))

# Relax the structure
relax = QuasiNewton(slab)
relax.run(fmax=0.05)
print('initial state:', slab.get_potential_energy())
write('N2.traj', slab)

# Now the final state.
# Move the second N atom to a neighboring hollow site:
slab[-1].position[0] = slab[-2].position[0] + 0.25 * slab.cell[0, 0]
slab[-1].position[1] = slab[-2].position[1]
# and relax.
relax.run()
print('final state:  ', slab.get_potential_energy())
write('2N.traj', slab)

Having obtained these structures we set up an NEB calculation with 9 images. Using interpolate() provides a guess for the path between the initial and final states. We perform the relaxation of the images and obtain the intermediate steps.

import numpy as np

from ase.calculators.emt import EMT
from ase.constraints import FixAtoms
from ase.io import read
from ase.mep import NEB
from ase.optimize.fire import FIRE as QuasiNewton

# Read the previous configurations
initial = read('N2.traj')
final = read('2N.traj')

#  Make 9 images (note the use of copy)
configs = [initial.copy() for i in range(8)] + [final]

# As before, fix the Cu atoms
constraint = FixAtoms(mask=[atom.symbol != 'N' for atom in initial])
for config in configs:
    config.calc = EMT()
    config.set_constraint(constraint)

# Make the NEB object, interpolate to guess the intermediate steps
band = NEB(configs)
band.interpolate()

relax = QuasiNewton(band)

# Do the calculation
relax.run()

# Compare intermediate steps to initial energy
e0 = initial.get_potential_energy()
for config in configs:
    d = config[-2].position - config[-1].position
    print(np.linalg.norm(d), config.get_potential_energy() - e0)

After the calculation is complete, the energy difference with respect to the initial state is given for each image, as well as the distance between the N atoms.