Skip to main content
`Ctrl`+`K`

# Surface diffusion energy barriers using the Nudged Elastic Band (NEB) method

Note

Go to the end to download the full example code.

# Surface diffusion energy barriers using the Nudged Elastic Band (NEB) method

    import matplotlib.pyplot as plt

    from ase.build import add_adsorbate, fcc100
    from ase.calculators.emt import EMT
    from ase.constraints import FixAtoms
    from ase.io import read
    from ase.mep import NEB
    from ase.optimize import BFGS, QuasiNewton
    from ase.parallel import world
    from ase.visualize.plot import plot_atoms

First, set up the initial and final states: 2x2-Al(001) surface with 3 layers and an Au atom adsorbed in a hollow site:

    slab = fcc100('Al', size=(2, 2, 3))
    add_adsorbate(slab, 'Au', 1.7, 'hollow')
    slab.center(axis=2, vacuum=4.0)

Make sure the structure is correct:

    fig, ax = plt.subplots()
    plot_atoms(slab, ax)
    ax.set_axis_off()

Fix second and third layers:

    mask = [atom.tag > 1 for atom in slab]
    print(mask)
    slab.set_constraint(FixAtoms(mask=mask))

    [np.True_, np.True_, np.True_, np.True_, np.True_, np.True_, np.True_, np.True_, np.False_, np.False_, np.False_, np.False_, np.False_]

Use EMT potential:

    slab.calc = EMT()

Initial state:

    qn = QuasiNewton(slab, trajectory='initial.traj')
    qn.run(fmax=0.05)

                    Step[ FC]     Time          Energy          fmax
    BFGSLineSearch:    0[  0] 09:17:40        3.323870       0.2462
    BFGSLineSearch:    1[  1] 09:17:40        3.314754       0.0378

    np.True_

Final state:

    slab[-1].x += slab.get_cell()[0, 0] / 2
    qn = QuasiNewton(slab, trajectory='final.traj')
    qn.run(fmax=0.05)

                    Step[ FC]     Time          Energy          fmax
    BFGSLineSearch:    0[  0] 09:17:40        3.320051       0.1208
    BFGSLineSearch:    1[  1] 09:17:40        3.316117       0.0474

    np.True_

Note

Notice how the tags are used to select the constrained atoms

Now, do the NEB calculation:

    initial = read('initial.traj')
    final = read('final.traj')

    constraint = FixAtoms(mask=[atom.tag > 1 for atom in initial])

    images = [initial]
    for i in range(3):
        image = initial.copy()
        image.calc = EMT()
        image.set_constraint(constraint)
        images.append(image)

    images.append(final)

    neb = NEB(images, method='improvedtangent')
    neb.interpolate()
    qn = BFGS(neb, trajectory='neb.traj')
    qn.run(fmax=0.05)

          Step     Time          Energy          fmax
    BFGS:    0 09:17:40        4.219952        3.520752
    BFGS:    1 09:17:40        3.937039        2.176460
    BFGS:    2 09:17:40        3.719814        0.435144
    BFGS:    3 09:17:40        3.709652        0.230134
    BFGS:    4 09:17:40        3.708878        0.244118
    BFGS:    5 09:17:40        3.706087        0.257725
    BFGS:    6 09:17:40        3.698531        0.213294
    BFGS:    7 09:17:40        3.692120        0.246238
    BFGS:    8 09:17:40        3.692278        0.187553
    BFGS:    9 09:17:40        3.693494        0.172952
    BFGS:   10 09:17:40        3.692670        0.151720
    BFGS:   11 09:17:40        3.690813        0.073884
    BFGS:   12 09:17:40        3.690200        0.071211
    BFGS:   13 09:17:40        3.690380        0.078122
    BFGS:   14 09:17:40        3.690427        0.103668
    BFGS:   15 09:17:40        3.689897        0.100381
    BFGS:   16 09:17:40        3.689034        0.055057
    BFGS:   17 09:17:40        3.688737        0.028898

    np.True_

Visualize the results with:

    .. GENERATED FROM PYTHON SOURCE LINES 95-99

    fig, ax = plt.subplots()
    plot_atoms(final, ax)
    ax.set_axis_off()

or from the command line:

> $ ase gui neb.traj@-5:

and select Tools->NEB.

You can also create a series of plots like above, that show the progression of the NEB relaxation, directly at the command line:

    $ ase nebplot --share-x --share-y neb.traj

For more customizable analysis of the output of many NEB jobs, you can use the `ase.mep.NEBTools` class. Some examples of its use are below; the final example was used to make the figure you see above.

    import matplotlib.pyplot as plt

    from ase.io import read
    from ase.mep import NEBTools

    images = read('neb.traj@-5:')

    nebtools = NEBTools(images)

Get the calculated barrier and the energy change of the reaction.

    Ef, dE = nebtools.get_barrier()

Get the barrier without any interpolation between highest images.

    Ef, dE = nebtools.get_barrier(fit=False)

Get the actual maximum force at this point in the simulation.

    max_force = nebtools.get_fmax()

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

Create a figure like that coming from ASE-GUI.

    fig = nebtools.plot_band()

Create a figure with custom parameters.

    fig = plt.figure(figsize=(5.5, 4.0))
    ax = fig.add_axes((0.15, 0.15, 0.8, 0.75))
    nebtools.plot_band(ax)

    <Figure size 550x400 with 1 Axes>

Note

For this reaction, the reaction coordinate is very simple: The x-coordinate of the Au atom. In such cases, the NEB method is overkill, and a simple constraint method should be used like in this tutorial: Constrained Calculations - Surface diffusion energy barriers.

See also

- `ase.mep.neb`

- `ase.constraints`

- Constrained Calculations - Surface diffusion energy barriers

- `fcc100()`

## Restarting NEB

Restart NEB from the trajectory file:

read the last structures (of 5 images used in NEB)

    images = read('neb.traj@-5:')

    for i in range(1, len(images) - 1):
        images[i].calc = EMT()

    neb = NEB(images, method='improvedtangent')
    qn = BFGS(neb, trajectory='neb_restart.traj')
    qn.run(fmax=0.005)

          Step     Time          Energy          fmax
    BFGS:    0 09:17:40        3.688737        0.028898
    BFGS:    1 09:17:40        3.688734        0.026547
    BFGS:    2 09:17:41        3.688731        0.021515
    BFGS:    3 09:17:41        3.688730        0.021379
    BFGS:    4 09:17:41        3.688718        0.009577
    BFGS:    5 09:17:41        3.688718        0.008389
    BFGS:    6 09:17:41        3.688718        0.009122
    BFGS:    7 09:17:41        3.688718        0.009290
    BFGS:    8 09:17:41        3.688716        0.007568
    BFGS:    9 09:17:41        3.688714        0.005563
    BFGS:   10 09:17:41        3.688715        0.002547

    np.True_

## Parallelizing over images with MPI

Instead of having one process do the calculations for all three internal images in turn, it will be faster to have three processes do one image each. In order to be able to run python with MPI you need a special parallel python interpreter, for example gpaw python (see GPAW parallel runs) and set `parallel=True` in the NEB calculation.

The example below can then be run with `gpaw -P3 python diffusion_parallel.py`:

    initial = read('initial.traj')
    final = read('final.traj')

    constraint = FixAtoms(mask=[atom.tag > 1 for atom in initial])

    n_images = 1  # Set to number of processes you use with mpiexec
    images = [initial]
    j = world.rank * n_images // world.size  # my image number
    for i in range(n_images):
        image = initial.copy()
        if i == j:
            image.calc = EMT()
        image.set_constraint(constraint)
        images.append(image)
    images.append(final)

    neb = NEB(images, parallel=True, method='improvedtangent')
    neb.interpolate()
    qn = BFGS(neb, trajectory='neb.traj')
    qn.run(fmax=0.05)

          Step     Time          Energy          fmax
    BFGS:    0 09:17:41        4.219952        3.520752
    BFGS:    1 09:17:41        3.937039        2.176459
    BFGS:    2 09:17:41        3.730687        0.621235
    BFGS:    3 09:17:41        3.710767        0.211379
    BFGS:    4 09:17:41        3.707518        0.163942
    BFGS:    5 09:17:41        3.704890        0.187028
    BFGS:    6 09:17:41        3.699676        0.197511
    BFGS:    7 09:17:41        3.694802        0.156084
    BFGS:    8 09:17:41        3.691892        0.081850
    BFGS:    9 09:17:41        3.691165        0.050493
    BFGS:   10 09:17:41        3.690969        0.048196

    np.True_

Gallery generated by Sphinx-Gallery