
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples_generated/tutorials/diffusion_neb.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_diffusion_neb.py>`
        to download the full example code.

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

.. _sphx_glr_examples_generated_tutorials_diffusion_neb.py:

.. _diffusion tutorial:

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

.. GENERATED FROM PYTHON SOURCE LINES 10-21

.. code-block:: Python

    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








.. GENERATED FROM PYTHON SOURCE LINES 22-25

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

.. GENERATED FROM PYTHON SOURCE LINES 27-31

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 32-33

Make sure the structure is correct:

.. GENERATED FROM PYTHON SOURCE LINES 35-39

.. code-block:: Python

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




.. image-sg:: /examples_generated/tutorials/images/sphx_glr_diffusion_neb_001.png
   :alt: diffusion neb
   :srcset: /examples_generated/tutorials/images/sphx_glr_diffusion_neb_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 40-41

Fix second and third layers:

.. GENERATED FROM PYTHON SOURCE LINES 43-47

.. code-block:: Python

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





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

 .. code-block:: none

    [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_]




.. GENERATED FROM PYTHON SOURCE LINES 48-49

Use EMT potential:

.. GENERATED FROM PYTHON SOURCE LINES 51-53

.. code-block:: Python

    slab.calc = EMT()








.. GENERATED FROM PYTHON SOURCE LINES 54-55

Initial state:

.. GENERATED FROM PYTHON SOURCE LINES 57-60

.. code-block:: Python

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





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

 .. code-block:: none

                    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_



.. GENERATED FROM PYTHON SOURCE LINES 61-62

Final state:

.. GENERATED FROM PYTHON SOURCE LINES 64-68

.. code-block:: Python

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





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

 .. code-block:: none

                    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_



.. GENERATED FROM PYTHON SOURCE LINES 69-72

.. note::  Notice how the tags are used to select the constrained atoms

Now, do the NEB calculation:

.. GENERATED FROM PYTHON SOURCE LINES 74-93

.. code-block:: Python

    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)





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

 .. code-block:: none

          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_



.. GENERATED FROM PYTHON SOURCE LINES 94-95

Visualize the results with::

.. GENERATED FROM PYTHON SOURCE LINES 95-99

.. code-block:: Python

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




.. image-sg:: /examples_generated/tutorials/images/sphx_glr_diffusion_neb_002.png
   :alt: diffusion neb
   :srcset: /examples_generated/tutorials/images/sphx_glr_diffusion_neb_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 100-105

or from the command line:

   $ ase gui neb.traj@-5:

and select Tools->NEB.

.. GENERATED FROM PYTHON SOURCE LINES 108-116

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 :class:`ase.mep.NEBTools` class. Some examples of its use are below; the
final example was used to make the figure you see above.

.. GENERATED FROM PYTHON SOURCE LINES 118-127

.. code-block:: Python

    import matplotlib.pyplot as plt

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

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

    nebtools = NEBTools(images)








.. GENERATED FROM PYTHON SOURCE LINES 128-129

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

.. GENERATED FROM PYTHON SOURCE LINES 131-133

.. code-block:: Python

    Ef, dE = nebtools.get_barrier()








.. GENERATED FROM PYTHON SOURCE LINES 134-135

Get the barrier without any interpolation between highest images.

.. GENERATED FROM PYTHON SOURCE LINES 137-139

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 140-141

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

.. GENERATED FROM PYTHON SOURCE LINES 143-145

.. code-block:: Python

    max_force = nebtools.get_fmax()





.. 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 146-147

Create a figure like that coming from ASE-GUI.

.. GENERATED FROM PYTHON SOURCE LINES 149-151

.. code-block:: Python

    fig = nebtools.plot_band()




.. image-sg:: /examples_generated/tutorials/images/sphx_glr_diffusion_neb_003.png
   :alt: $E_\mathrm{f} \approx$ 0.374 eV; $E_\mathrm{r} \approx$ 0.373 eV; $\Delta E$ = 0.001 eV
   :srcset: /examples_generated/tutorials/images/sphx_glr_diffusion_neb_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 152-153

Create a figure with custom parameters.

.. GENERATED FROM PYTHON SOURCE LINES 155-159

.. code-block:: Python

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




.. image-sg:: /examples_generated/tutorials/images/sphx_glr_diffusion_neb_004.png
   :alt: $E_\mathrm{f} \approx$ 0.374 eV; $E_\mathrm{r} \approx$ 0.373 eV; $\Delta E$ = 0.001 eV
   :srcset: /examples_generated/tutorials/images/sphx_glr_diffusion_neb_004.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    <Figure size 550x400 with 1 Axes>



.. GENERATED FROM PYTHON SOURCE LINES 160-174

.. 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: :ref:`constraints diffusion tutorial`.

.. seealso::

  * :mod:`ase.mep.neb`
  * :mod:`ase.constraints`
  * :ref:`constraints diffusion tutorial`
  * :func:`~ase.build.fcc100`


.. GENERATED FROM PYTHON SOURCE LINES 176-180

Restarting NEB
==============

Restart NEB from the trajectory file:

.. GENERATED FROM PYTHON SOURCE LINES 182-183

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

.. GENERATED FROM PYTHON SOURCE LINES 185-194

.. code-block:: Python

    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)





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

 .. code-block:: none

          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_



.. GENERATED FROM PYTHON SOURCE LINES 195-207

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 <https://gpaw.readthedocs.io/documentation/parallel_runs/parallel_runs.html>`_)
and set ``parallel=True`` in the NEB calculation.

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

.. GENERATED FROM PYTHON SOURCE LINES 209-229

.. code-block:: Python

    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)




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

 .. code-block:: none

          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_




.. _sphx_glr_download_examples_generated_tutorials_diffusion_neb.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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