Skip to main content
`Ctrl`+`K`

# NEB: Dissociation of a molecule using the NEB method

Note

Go to the end to download the full example code.

# NEB: 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 N2 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] 09:41:36       14.459066       9.8897
    BFGSLineSearch:    1[  1] 09:41:36       12.244808       4.3275
    BFGSLineSearch:    2[  2] 09:41:36       11.343194       2.1646
    BFGSLineSearch:    3[  3] 09:41:36       11.077854       0.6379
    BFGSLineSearch:    4[  4] 09:41:36       11.056067       0.2425
    BFGSLineSearch:    5[  6] 09:41:36       11.049129       0.2142
    BFGSLineSearch:    6[  8] 09:41:36       11.034808       0.1607
    BFGSLineSearch:    7[  9] 09:41:36       11.032950       0.0397
    initial state: 11.03295039775445

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] 09:41:36       11.023212       0.0971
    BFGSLineSearch:    9[ 13] 09:41:36       11.021738       0.1053
    BFGSLineSearch:   10[ 15] 09:41:36       11.017702       0.0733
    BFGSLineSearch:   11[ 16] 09:41:36       11.017285       0.0323
    final state:   11.01728455444998

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()

    /build/python-ase-iXAd1o/python-ase-3.28.0/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 09:41:36       11.066750        0.523744
    FIRE:    1 09:41:36       11.063989        0.462236
    FIRE:    2 09:41:36       11.059699        0.348578
    FIRE:    3 09:41:36       11.055743        0.199917
    FIRE:    4 09:41:36       11.053732        0.053868
    FIRE:    5 09:41:36       11.054186        0.126908
    FIRE:    6 09:41:36       11.054140        0.123746
    FIRE:    7 09:41:36       11.054054        0.117517
    FIRE:    8 09:41:36       11.053934        0.108416
    FIRE:    9 09:41:36       11.053791        0.096753
    FIRE:   10 09:41:36       11.053638        0.082986
    FIRE:   11 09:41:37       11.053488        0.067787
    FIRE:   12 09:41:37       11.053353        0.052238
    FIRE:   13 09:41:37       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'
        )

    0     d = 1.661AA     delta energy = 0.000eV
    1     d = 1.739AA     delta energy = 0.003eV
    2     d = 1.829AA     delta energy = 0.012eV
    3     d = 1.936AA     delta energy = 0.019eV
    4     d = 2.055AA     delta energy = 0.020eV
    5     d = 2.183AA     delta energy = 0.015eV
    6     d = 2.319AA     delta energy = 0.004eV
    7     d = 2.460AA     delta energy = -0.010eV
    8     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