
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples_generated/tutorials/neb_dissociation.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_generated_tutorials_neb_dissociation.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_generated_tutorials_neb_dissociation.py:

.. _dissociation:

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 :mod:`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.

.. GENERATED FROM PYTHON SOURCE LINES 15-30

.. code-block:: Python


    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








.. GENERATED FROM PYTHON SOURCE LINES 31-33

First, we create the initial state: An N\ :sub:`2`\  molecule on a
Cu(111) slab

.. GENERATED FROM PYTHON SOURCE LINES 33-44

.. code-block:: Python


    # 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]]
    )








.. GENERATED FROM PYTHON SOURCE LINES 45-47

Then we put the adsorbate onto the surface


.. GENERATED FROM PYTHON SOURCE LINES 47-50

.. code-block:: Python


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








.. GENERATED FROM PYTHON SOURCE LINES 51-52

And add a calculator for the forces and energies:

.. GENERATED FROM PYTHON SOURCE LINES 52-55

.. code-block:: Python


    slab.calc = EMT()








.. GENERATED FROM PYTHON SOURCE LINES 56-58

We don't want to worry about the Cu degrees of freedom,
so fix these atoms:

.. GENERATED FROM PYTHON SOURCE LINES 58-62

.. code-block:: Python


    mask = [atom.symbol == 'Cu' for atom in slab]
    slab.set_constraint(FixAtoms(mask=mask))








.. GENERATED FROM PYTHON SOURCE LINES 63-64

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

.. GENERATED FROM PYTHON SOURCE LINES 64-70

.. code-block:: Python


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





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

                    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




.. GENERATED FROM PYTHON SOURCE LINES 71-73

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

.. GENERATED FROM PYTHON SOURCE LINES 73-81

.. code-block:: Python


    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)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    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




.. GENERATED FROM PYTHON SOURCE LINES 82-89

Having obtained these structures we set up an NEB
calculation with 9 images.  Using :func:`~neb.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

.. GENERATED FROM PYTHON SOURCE LINES 89-93

.. code-block:: Python


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








.. GENERATED FROM PYTHON SOURCE LINES 94-96

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

.. GENERATED FROM PYTHON SOURCE LINES 96-104

.. code-block:: Python


    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)








.. GENERATED FROM PYTHON SOURCE LINES 105-107

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

.. GENERATED FROM PYTHON SOURCE LINES 107-111

.. code-block:: Python


    band = NEB(configs)
    band.interpolate()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

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




.. GENERATED FROM PYTHON SOURCE LINES 112-116

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.

.. GENERATED FROM PYTHON SOURCE LINES 116-120

.. code-block:: Python


    relax = FIRE(band)
    relax.run()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

          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_



.. GENERATED FROM PYTHON SOURCE LINES 121-122

Finally, we compare intermediate steps to the initial energy

.. GENERATED FROM PYTHON SOURCE LINES 122-136

.. code-block:: Python


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





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

      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




.. GENERATED FROM PYTHON SOURCE LINES 137-141

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.


.. GENERATED FROM PYTHON SOURCE LINES 141-176

.. code-block:: Python


    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
    )



.. video:: /examples_generated/tutorials/images/sphx_glr_neb_dissociation_001.mp4
   :class: sphx-glr-single-img
   :height: 480
   :width: 640
   :autoplay:
   :loop:







.. _sphx_glr_download_examples_generated_tutorials_neb_dissociation.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: neb_dissociation.ipynb <neb_dissociation.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: neb_dissociation.py <neb_dissociation.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: neb_dissociation.zip <neb_dissociation.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
