index | modules | gitlab | page source

# 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.