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.

import matplotlib.animation as animation
import matplotlib.pyplot as plt
import numpy as np

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 read, write
from ase.mep import NEB
from ase.optimize import QuasiNewton
from ase.optimize.fire import FIRE
from ase.visualize.plot import plot_atoms

First, we create the initial state: An N₂ molecule on a Cu(111) slab

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

# Add the N2 Molecule,  oriented at 60 degrees:
d = 1.10  # N2 bond length
n2_mol = Atoms(
    'N2', positions=[[0.0, 0.0, 0.0], [0.5 * 3**0.5 * d, 0.5 * d, 0.0]]
)

Then we put the adsorbate onto the surface

add_adsorbate(slab, n2_mol, height=1.0, position='fcc')

And add a 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))

Then we relax the structure and write the trajectory into a file

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

                Step[ FC]     Time          Energy          fmax
BFGSLineSearch:    0[  0] 16:34:46       14.459066       9.8897
BFGSLineSearch:    1[  1] 16:34:46       12.244808       4.3275
BFGSLineSearch:    2[  2] 16:34:46       11.343194       2.1646
BFGSLineSearch:    3[  3] 16:34:46       11.077854       0.6379
BFGSLineSearch:    4[  4] 16:34:46       11.056067       0.2425
BFGSLineSearch:    5[  6] 16:34:46       11.049129       0.2142
BFGSLineSearch:    6[  8] 16:34:46       11.034808       0.1607
BFGSLineSearch:    7[  9] 16:34:46       11.032950       0.0397
initial state: 11.032950397754657

Now the final state. Move the second N atom to a neighboring hollow site and relax.

slab[-1].position[0] = slab[-2].position[0] + 0.25 * slab.cell[0, 0]
slab[-1].position[1] = slab[-2].position[1]

relax.run()
print('final state:  ', slab.get_potential_energy())
write('2N.traj', slab)

BFGSLineSearch:    8[ 11] 16:34:46       11.023212       0.0971
BFGSLineSearch:    9[ 13] 16:34:46       11.021738       0.1053
BFGSLineSearch:   10[ 15] 16:34:46       11.017702       0.0733
BFGSLineSearch:   11[ 16] 16:34:46       11.017285       0.0323
final state:   11.017284554450073

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.

First, we read the previous configurations

initial = read('N2.traj')
final = read('2N.traj')

Then, we make 9 images (note the use of copy) and, as before, we fix the Cu atoms

configs = [initial.copy() for i in range(8)] + [final]

constraint = FixAtoms(mask=[atom.symbol != 'N' for atom in initial])
for config in configs:
    config.calc = EMT()
    config.set_constraint(constraint)

Next, we make the NEB object, and call its interpolate() method to guess the intermediate steps

band = NEB(configs)
band.interpolate()

/home/ase/.local/lib/python3.13/site-packages/ase/mep/neb.py:329: UserWarning: The default method has changed from 'aseneb' to 'improvedtangent'. The 'aseneb' method is an unpublished, custom implementation that is not recommended as it frequently results in very poor bands. Please explicitly set method='improvedtangent' to silence this warning, or set method='aseneb' if you strictly require the old behavior (results may vary). See: https://gitlab.com/ase/ase/-/merge_requests/3952
  warnings.warn(

Then we relax the configurations in the NEB object. We can use same FIRE optimization class we might apply to optimize an Atoms object; NEB presents appropriate degrees of freedom and gradients for this to optimise multiple configurations simultaneously.

relax = FIRE(band)
relax.run()

      Step     Time          Energy          fmax
FIRE:    0 16:34:46       11.066750        0.523744
FIRE:    1 16:34:46       11.063989        0.462236
FIRE:    2 16:34:47       11.059699        0.348578
FIRE:    3 16:34:47       11.055743        0.199917
FIRE:    4 16:34:47       11.053732        0.053868
FIRE:    5 16:34:47       11.054186        0.126908
FIRE:    6 16:34:47       11.054140        0.123746
FIRE:    7 16:34:47       11.054054        0.117517
FIRE:    8 16:34:47       11.053934        0.108416
FIRE:    9 16:34:47       11.053791        0.096753
FIRE:   10 16:34:47       11.053638        0.082986
FIRE:   11 16:34:47       11.053488        0.067787
FIRE:   12 16:34:47       11.053353        0.052238
FIRE:   13 16:34:47       11.053230        0.039264

np.True_

Finally, we compare intermediate steps to the initial energy

energy_initial = initial.get_potential_energy()
n2_distances = []
energy_differences = []
for i, config in enumerate(configs):
    n2_distance = np.linalg.norm(config[-2].position - config[-1].position)
    delta_energy = config.get_potential_energy() - energy_initial
    n2_distances.append(n2_distance)
    energy_differences.append(delta_energy)
    print(
        f'{i:>3}\td = {n2_distance:>4.3f}AA '
        f'\tdelta energy = {delta_energy:>5.3f}eV'
    )

   d = 1.661AA     delta energy = 0.000eV
   d = 1.739AA     delta energy = 0.003eV
   d = 1.829AA     delta energy = 0.012eV
   d = 1.936AA     delta energy = 0.019eV
   d = 2.055AA     delta energy = 0.020eV
   d = 2.183AA     delta energy = 0.015eV
   d = 2.319AA     delta energy = 0.004eV
   d = 2.460AA     delta energy = -0.010eV
   d = 2.612AA     delta energy = -0.016eV

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.

energy_differences = [e * 1.0e03 for e in energy_differences]  # in meV

fig, axs = plt.subplots(
    nrows=2,
    gridspec_kw={'height_ratios': [0.5, 1.0]},
)

axs[0].plot(n2_distances, energy_differences)
scat = axs[0].scatter(n2_distances[0], energy_differences[0], s=10.0**2.0)
axs[0].set_ylabel(r'$E(i) - E_{\mathrm{initial}}$ (meV)')
axs[0].set_xlabel(r'$d_{\mathrm{N}-\mathrm{N}} (\AA)$')

# Plot the atomic structure, focussing on the two Nitrogen atoms
plot_atoms(config, axs[1], rotation='0x', show_unit_cell=0)
axs[1].set_axis_off()
axs[1].set_ylim(0.0, 5.642)
axs[1].set_xlim(0.0, 14.833)


def animate(i):
    scat.set_offsets((n2_distances[i], energy_differences[i]))

    # Remove the previous atomic plot
    [p.remove() for p in axs[1].patches]
    plot_atoms(configs[i], axs[1], rotation='0x', show_unit_cell=0)
    axs[1].set_xlim(0.0, 14.833)
    axs[1].set_ylim(0.0, 5.642)
    axs[1].set_axis_off()
    return (scat,)


ani = animation.FuncAnimation(
    fig, animate, repeat=True, frames=len(configs) - 1, interval=200
)

Gallery generated by Sphinx-Gallery