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
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 fcc111, add_adsorbate from ase.calculators.emt import EMT from ase.constraints import FixAtoms from ase.optimize import QuasiNewton from ase.io import write # 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 = slab[-2].position + 0.25 * slab.cell[0,0] slab[-1].position = slab[-2].position # 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
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.constraints import FixAtoms from ase.calculators.emt import EMT from ase.neb import NEB from ase.optimize.fire import FIRE as QuasiNewton from ase.io import read # 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.