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

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

.. _sphx_glr_examples_generated_tutorials_neb_selfdiffusion.py:

.. _selfdiffusion_example:

NEB and Dimer method for Self-diffusion on the Al(110) surface
==============================================================

.. GENERATED FROM PYTHON SOURCE LINES 10-17

In this exercise, we will find minimum-energy paths and transition states
using the :mod:`Nudged Elastic Band <ase.mep.neb>` method. We will illustrate
how NEB can be used in ASE to compute and compare three different
diffusion pathways for an Al atom
on a Al(110) surface. Finally, another method
for finding the transition state (i.e. the highest-energy state), the Dimer
method, will also be explored.

.. GENERATED FROM PYTHON SOURCE LINES 19-24

Initialize the system
---------------------

Al(110) surface can be generated with ASE code


.. GENERATED FROM PYTHON SOURCE LINES 25-35

.. code-block:: Python

    from math import sqrt

    import numpy as np

    from ase import Atom, Atoms

    a = 4.0614
    b = a / sqrt(2)
    h = b / 2








.. GENERATED FROM PYTHON SOURCE LINES 36-37

Create  :class:`~ase.Atoms` object and

.. GENERATED FROM PYTHON SOURCE LINES 38-46

.. code-block:: Python


    initial = Atoms(
        'Al2',
        positions=[(0, 0, 0), (a / 2, b / 2, -h)],
        cell=(a, b, 2 * h),
        pbc=(1, 1, 0),
    )








.. GENERATED FROM PYTHON SOURCE LINES 47-48

Multiply the unit cell to make it larger in x,y,z

.. GENERATED FROM PYTHON SOURCE LINES 48-50

.. code-block:: Python

    initial *= (4, 4, 2)








.. GENERATED FROM PYTHON SOURCE LINES 51-52

You can visualize the surface using

.. GENERATED FROM PYTHON SOURCE LINES 52-60

.. code-block:: Python

    import matplotlib.pyplot as plt

    from ase.visualize.plot import plot_atoms

    fig, ax = plt.subplots()
    plot_atoms(initial, ax, rotation=('-60x, 10y,0z'))
    ax.set_axis_off()




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





.. GENERATED FROM PYTHON SOURCE LINES 61-62

Then, lets add the Al addatom which will be the one moving on the surface.

.. GENERATED FROM PYTHON SOURCE LINES 62-64

.. code-block:: Python

    initial.append(Atom('Al', (a / 2, b / 2, 3 * h)))








.. GENERATED FROM PYTHON SOURCE LINES 65-66

Center the cell in vacuum along the z axis

.. GENERATED FROM PYTHON SOURCE LINES 66-68

.. code-block:: Python

    initial.center(vacuum=4.0, axis=2)








.. GENERATED FROM PYTHON SOURCE LINES 69-70

Visualize the new atom in the cell

.. GENERATED FROM PYTHON SOURCE LINES 70-75

.. code-block:: Python

    fig, ax = plt.subplots()
    plot_atoms(initial, ax, rotation=('-60x, 10y,0z'))
    ax.set_axis_off()





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





.. GENERATED FROM PYTHON SOURCE LINES 76-86

Perform a NEB calculation
-------------------------

The adatom can jump along the rows (into the picture) or across
the rows (to the right inthe picture).
We are going to compute this motion to find out which of the two jump
will have the largest energy barrier.
To do this, you need to create an image with the atoms at
their final position.
First copy the initial :class:`~ase.Atoms` object

.. GENERATED FROM PYTHON SOURCE LINES 86-88

.. code-block:: Python

    final = initial.copy()








.. GENERATED FROM PYTHON SOURCE LINES 89-91

Then move the last atom of the :class:`~ase.Atoms` object "final"
(the one atom we just added before) of +b along the second positional array

.. GENERATED FROM PYTHON SOURCE LINES 91-93

.. code-block:: Python

    final.positions[-1, 1] += b








.. GENERATED FROM PYTHON SOURCE LINES 94-95

Visualize the new atom in the cell

.. GENERATED FROM PYTHON SOURCE LINES 95-100

.. code-block:: Python

    fig, ax = plt.subplots()
    plot_atoms(final, ax, rotation=('-60x, 10y,0z'))
    ax.set_axis_off()





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





.. GENERATED FROM PYTHON SOURCE LINES 101-105

Let us fix the atoms that are not moving by creating a constraint
and setting this constraint to the images.
To do this, we create a mask of boolean array that select
fixed atoms (the two bottom layers):

.. GENERATED FROM PYTHON SOURCE LINES 105-113

.. code-block:: Python


    from ase.calculators.emt import EMT
    from ase.constraints import FixAtoms

    mask = initial.positions[:, 2] - min(initial.positions[:, 2]) < 1.5 * h
    constraint = FixAtoms(mask=mask)
    print(mask)





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

 .. code-block:: none

    [ True  True False False  True  True False False  True  True False False
      True  True False False  True  True False False  True  True False False
      True  True False False  True  True False False  True  True False False
      True  True False False  True  True False False  True  True False False
      True  True False False  True  True False False  True  True False False
      True  True False False False]




.. GENERATED FROM PYTHON SOURCE LINES 114-117

Set the :class:`~ase.constraints.FixAtoms` to the :class:`~ase.Atoms` objects,
and in the same loop, set the calculator (in this example we use EMT,
but you can use any calculator supported by ASE)

.. GENERATED FROM PYTHON SOURCE LINES 117-123

.. code-block:: Python

    initial.calc = EMT()
    initial.set_constraint(constraint)
    final.calc = EMT()
    final.set_constraint(constraint)









.. GENERATED FROM PYTHON SOURCE LINES 124-127

Use :class:`~ase.calculators.emt` calculator and
:class:`~ase.optimize.QuasiNewton` Algorithm to optimize
the geometry of the initial and final states

.. GENERATED FROM PYTHON SOURCE LINES 127-133

.. code-block:: Python

    from ase.optimize import QuasiNewton

    QuasiNewton(initial).run(fmax=0.05)
    QuasiNewton(final).run(fmax=0.05)






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

 .. code-block:: none

                    Step[ FC]     Time          Energy          fmax
    BFGSLineSearch:    0[  0] 09:41:45       17.939914       0.2540
    BFGSLineSearch:    1[  2] 09:41:45       17.897422       0.1631
    BFGSLineSearch:    2[  4] 09:41:45       17.884347       0.0522
    BFGSLineSearch:    3[  5] 09:41:45       17.881710       0.0269
                    Step[ FC]     Time          Energy          fmax
    BFGSLineSearch:    0[  0] 09:41:45       17.939914       0.2540
    BFGSLineSearch:    1[  2] 09:41:45       17.897422       0.1631
    BFGSLineSearch:    2[  4] 09:41:45       17.884347       0.0522
    BFGSLineSearch:    3[  5] 09:41:45       17.881710       0.0269

    np.True_



.. GENERATED FROM PYTHON SOURCE LINES 134-136

Then, construct a list of images by copying the first image several time in
an array and append to this list the final image

.. GENERATED FROM PYTHON SOURCE LINES 136-141

.. code-block:: Python

    images = [initial]
    for i in range(5):
        images.append(initial.copy())
    images.append(final)








.. GENERATED FROM PYTHON SOURCE LINES 142-144

Because the .copy() method does not copy the calculator, you need to set a
new one for the created images

.. GENERATED FROM PYTHON SOURCE LINES 144-149

.. code-block:: Python

    for image in images:
        image.calc = EMT()
        image.set_constraint(constraint)









.. GENERATED FROM PYTHON SOURCE LINES 150-151

Create a Nudged Elastic Band (:class:`~ase.mep import NEB`) object

.. GENERATED FROM PYTHON SOURCE LINES 151-155

.. code-block:: Python

    from ase.mep import NEB

    neb = NEB(images)





.. 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 156-158

Make a starting guess for the minimum energy path by performing a linear
interpolation from the initial to the final image

.. GENERATED FROM PYTHON SOURCE LINES 158-160

.. code-block:: Python

    neb.interpolate()








.. GENERATED FROM PYTHON SOURCE LINES 161-162

Perform the NEB calculation minimizing the force below 0.05 eV/A

.. GENERATED FROM PYTHON SOURCE LINES 162-167

.. code-block:: Python

    from ase.optimize import MDMin

    minimizer = MDMin(neb)
    minimizer.run(fmax=0.05)





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

 .. code-block:: none

           Step     Time          Energy          fmax
    MDMin:    0 09:41:46       18.237088        0.993544
    MDMin:    1 09:41:46       18.190945        0.836011
    MDMin:    2 09:41:46       18.105033        0.460778
    MDMin:    3 09:41:46       18.056201        0.122138
    MDMin:    4 09:41:46       18.035931        0.102285
    MDMin:    5 09:41:46       18.022022        0.118145
    MDMin:    6 09:41:46       18.010798        0.088852
    MDMin:    7 09:41:46       18.005186        0.127228
    MDMin:    8 09:41:46       18.003865        0.114049
    MDMin:    9 09:41:46       18.001954        0.088821
    MDMin:   10 09:41:46       18.000262        0.063304
    MDMin:   11 09:41:46       17.998668        0.053980
    MDMin:   12 09:41:46       17.997130        0.040716

    np.True_



.. GENERATED FROM PYTHON SOURCE LINES 168-170

Visualize the minimum energy path (MEP) in side view to see the motion.
Here, we look at the surface slab in yz direction.

.. GENERATED FROM PYTHON SOURCE LINES 170-175

.. code-block:: Python

    for image in images:
        fig, ax = plt.subplots()
        plot_atoms(image, ax, rotation=('-90x, 90y,0z'))
        ax.set_axis_off()




.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_004.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_004.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_005.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_005.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_006.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_006.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_007.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_007.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_008.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_008.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_009.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_009.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_010.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_010.png
         :class: sphx-glr-multi-img





.. GENERATED FROM PYTHON SOURCE LINES 176-177

Plot the variation of potential energy

.. GENERATED FROM PYTHON SOURCE LINES 177-191

.. code-block:: Python

    potential_energies = [image.get_potential_energy() for image in images]

    fig, ax = plt.subplots()
    plt.plot(
        range(len(potential_energies)),
        potential_energies - potential_energies[0],
        marker='+',
    )
    plt.xlabel('Image number')
    plt.ylabel('Potential energy (eV)')

    diff = np.max(potential_energies) - potential_energies[0]
    print(f'The energy barrier is {diff:.4f} eV.')




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


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

 .. code-block:: none

    The energy barrier is 0.1154 eV.




.. GENERATED FROM PYTHON SOURCE LINES 192-194

You can visualize the NEB path using ASE GUI after saving the
tajectory in a file

.. GENERATED FROM PYTHON SOURCE LINES 194-199

.. code-block:: Python


    from ase.io import write

    write('neb_path.traj', images, format='traj')








.. GENERATED FROM PYTHON SOURCE LINES 200-203

Otherwise, you can use ``ase gui neb_path.traj`` command in your terminal and
visualize the energy curve by plotting ``i, E[i] - E[1]``.
You now can answer those questions :

.. GENERATED FROM PYTHON SOURCE LINES 205-210

* How is the shape of the potential (symmetric/asymmetric) and does this make
  sense for this process (when looking at the moving adatom in the
  simulation)?
* What is the energy barrier?


.. GENERATED FROM PYTHON SOURCE LINES 213-220

Beyond your first NEB calculation
----------------------------------

You now can redo the same process to find the energy barrier
to cross one row.
The following code will produce the result (by making use of the previously
initialized code), though we encourage you to try by yourself.

.. GENERATED FROM PYTHON SOURCE LINES 222-227

.. code-block:: Python



    final = initial.copy()
    final.positions[-1, 0] += a








.. GENERATED FROM PYTHON SOURCE LINES 228-229

Plot the images

.. GENERATED FROM PYTHON SOURCE LINES 229-234

.. code-block:: Python

    for image in [initial, final]:
        fig, ax = plt.subplots()
        plot_atoms(image, ax, rotation=('-90x, 0y, 0z'))
        ax.set_axis_off()




.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_012.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_012.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_013.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_013.png
         :class: sphx-glr-multi-img





.. GENERATED FROM PYTHON SOURCE LINES 235-236

Construct a list of images:

.. GENERATED FROM PYTHON SOURCE LINES 237-242

.. code-block:: Python

    images = [initial]
    for i in range(5):
        images.append(initial.copy())
    images.append(final)








.. GENERATED FROM PYTHON SOURCE LINES 243-245

Make a mask of zeros and ones that select fixed atoms (the
two bottom layers):

.. GENERATED FROM PYTHON SOURCE LINES 245-254

.. code-block:: Python

    mask = initial.positions[:, 2] - min(initial.positions[:, 2]) < 1.5 * h
    constraint = FixAtoms(mask=mask)
    print(mask)

    for image in images:
        # Let all images use an EMT calculator:
        image.calc = EMT()
        image.set_constraint(constraint)





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

 .. code-block:: none

    [ True  True False False  True  True False False  True  True False False
      True  True False False  True  True False False  True  True False False
      True  True False False  True  True False False  True  True False False
      True  True False False  True  True False False  True  True False False
      True  True False False  True  True False False  True  True False False
      True  True False False False]




.. GENERATED FROM PYTHON SOURCE LINES 255-256

Relax the initial and final states:

.. GENERATED FROM PYTHON SOURCE LINES 256-259

.. code-block:: Python

    QuasiNewton(initial).run(fmax=0.05)
    QuasiNewton(final).run(fmax=0.05)





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

 .. code-block:: none

                    Step[ FC]     Time          Energy          fmax
    BFGSLineSearch:    0[  0] 09:41:49       17.881710       0.0269
                    Step[ FC]     Time          Energy          fmax
    BFGSLineSearch:    0[  0] 09:41:49       17.898108       0.2165
    BFGSLineSearch:    1[  2] 09:41:50       17.883509       0.0755
    BFGSLineSearch:    2[  4] 09:41:50       17.880775       0.0216

    np.True_



.. GENERATED FROM PYTHON SOURCE LINES 260-261

Create a Nudged Elastic Band:

.. GENERATED FROM PYTHON SOURCE LINES 261-263

.. code-block:: Python

    neb = NEB(images)








.. GENERATED FROM PYTHON SOURCE LINES 264-266

Make a starting guess for the minimum energy path (a straight line
from the initial to the final state):

.. GENERATED FROM PYTHON SOURCE LINES 266-268

.. code-block:: Python

    neb.interpolate()








.. GENERATED FROM PYTHON SOURCE LINES 269-270

Relax the NEB path:

.. GENERATED FROM PYTHON SOURCE LINES 270-274

.. code-block:: Python

    minimizer = MDMin(neb)
    minimizer.run(fmax=0.05)






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

 .. code-block:: none

           Step     Time          Energy          fmax
    MDMin:    0 09:41:50       21.182521        7.493573
    MDMin:    1 09:41:50       19.417286        3.689216
    MDMin:    2 09:41:50       18.645079        0.922870
    MDMin:    3 09:41:50       18.533896        0.503897
    MDMin:    4 09:41:50       18.525368        0.439764
    MDMin:    5 09:41:50       18.507383        0.263744
    MDMin:    6 09:41:50       18.491239        0.227551
    MDMin:    7 09:41:50       18.479659        0.160835
    MDMin:    8 09:41:50       18.468241        0.152303
    MDMin:    9 09:41:50       18.458050        0.091391
    MDMin:   10 09:41:50       18.450556        0.112507
    MDMin:   11 09:41:50       18.447753        0.095248
    MDMin:   12 09:41:50       18.447405        0.085207
    MDMin:   13 09:41:50       18.446629        0.058262
    MDMin:   14 09:41:50       18.445840        0.046679

    np.True_



.. GENERATED FROM PYTHON SOURCE LINES 275-277

Visualize the MEP in side view to see the motion.
We are now viewing the xz direction of the cell.

.. GENERATED FROM PYTHON SOURCE LINES 277-284

.. code-block:: Python


    for image in images:
        fig, ax = plt.subplots()
        plot_atoms(image, ax, rotation=('-90x, 0y, 0z'))
        ax.set_axis_off()





.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_014.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_014.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_015.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_015.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_016.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_016.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_017.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_017.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_018.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_018.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_019.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_019.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_020.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_020.png
         :class: sphx-glr-multi-img





.. GENERATED FROM PYTHON SOURCE LINES 285-286

Plot the variation of potential energy

.. GENERATED FROM PYTHON SOURCE LINES 286-301

.. code-block:: Python

    potential_energies = [image.get_potential_energy() for image in images]

    fig, ax = plt.subplots()
    plt.plot(
        range(len(potential_energies)),
        potential_energies - potential_energies[0],
        marker='+',
    )
    plt.xlabel('Image number')
    plt.ylabel('Potential energy (eV)')

    diff = np.max(potential_energies) - potential_energies[0]
    print(f'The energy barrier is {diff:.4f} eV.')





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


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

 .. code-block:: none

    The energy barrier is 0.5641 eV.




.. GENERATED FROM PYTHON SOURCE LINES 302-304

Finding the third mechanism
----------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 308-314

A third diffusion process can be found: Diffusion by an exchange process.
You can read more about it in the paper listed :mod:`here <ase.mep.dimer>`.
Find the barrier for this process, and compare the energy barrier
with the two other ones.
The following code will produce the result (by making use of the previously
initialized code), though we encourage you to try by yourself.

.. GENERATED FROM PYTHON SOURCE LINES 314-335

.. code-block:: Python


    a = 4.0614
    b = a / sqrt(2)
    h = b / 2
    initial = Atoms(
        'Al2',
        positions=[(0, 0, 0), (a / 2, b / 2, -h)],
        cell=(a, b, 2 * h),
        pbc=(1, 1, 0),
    )
    initial *= (2, 2, 2)
    initial.append(Atom('Al', (a / 2, b / 2, 3 * h)))
    initial.center(vacuum=4.0, axis=2)


    final = initial.copy()
    # move adatom to row atom 14
    final.positions[-1, :] = initial.positions[14]
    # Move row atom 14 to the next row
    final.positions[14, :] = initial.positions[-1] + [a, b, 0]








.. GENERATED FROM PYTHON SOURCE LINES 336-337

Visualize the initial and final images

.. GENERATED FROM PYTHON SOURCE LINES 337-342

.. code-block:: Python

    for image in [initial, final]:
        fig, ax = plt.subplots()
        plot_atoms(image, ax, rotation=('-60x, 10y,0z'))
        ax.set_axis_off()




.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_022.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_022.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_023.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_023.png
         :class: sphx-glr-multi-img





.. GENERATED FROM PYTHON SOURCE LINES 343-344

Construct a list of images:

.. GENERATED FROM PYTHON SOURCE LINES 344-349

.. code-block:: Python

    images = [initial]
    for i in range(5):
        images.append(initial.copy())
    images.append(final)








.. GENERATED FROM PYTHON SOURCE LINES 350-352

Make a mask of zeros and ones that select fixed atoms (the
two bottom layers):

.. GENERATED FROM PYTHON SOURCE LINES 352-356

.. code-block:: Python

    mask = initial.positions[:, 2] - min(initial.positions[:, 2]) < 1.5 * h
    constraint = FixAtoms(mask=mask)
    print(mask)





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

 .. code-block:: none

    [ True  True False False  True  True False False  True  True False False
      True  True False False False]




.. GENERATED FROM PYTHON SOURCE LINES 357-358

Let all images use an EMT calculator:

.. GENERATED FROM PYTHON SOURCE LINES 358-362

.. code-block:: Python

    for image in images:
        image.calc = EMT()
        image.set_constraint(constraint)








.. GENERATED FROM PYTHON SOURCE LINES 363-364

Relax the initial and final states:

.. GENERATED FROM PYTHON SOURCE LINES 364-367

.. code-block:: Python

    QuasiNewton(initial).run(fmax=0.05)
    QuasiNewton(final).run(fmax=0.05)





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

 .. code-block:: none

                    Step[ FC]     Time          Energy          fmax
    BFGSLineSearch:    0[  0] 09:41:53        4.639215       0.2538
    BFGSLineSearch:    1[  2] 09:41:53        4.622063       0.1580
    BFGSLineSearch:    2[  4] 09:41:53        4.613340       0.0622
    BFGSLineSearch:    3[  5] 09:41:53        4.611640       0.0336
                    Step[ FC]     Time          Energy          fmax
    BFGSLineSearch:    0[  0] 09:41:53        4.639215       0.2538
    BFGSLineSearch:    1[  2] 09:41:54        4.622063       0.1580
    BFGSLineSearch:    2[  4] 09:41:54        4.613340       0.0622
    BFGSLineSearch:    3[  5] 09:41:54        4.611640       0.0336

    np.True_



.. GENERATED FROM PYTHON SOURCE LINES 368-369

Create a Nudged Elastic Band:

.. GENERATED FROM PYTHON SOURCE LINES 369-371

.. code-block:: Python

    neb = NEB(images)








.. GENERATED FROM PYTHON SOURCE LINES 372-374

Make a starting guess for the minimum energy path (a straight line
from the initial to the final state):

.. GENERATED FROM PYTHON SOURCE LINES 374-376

.. code-block:: Python

    neb.interpolate()








.. GENERATED FROM PYTHON SOURCE LINES 377-378

Relax the NEB path:

.. GENERATED FROM PYTHON SOURCE LINES 378-382

.. code-block:: Python

    minimizer = MDMin(neb)
    minimizer.run(fmax=0.05)






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

 .. code-block:: none

           Step     Time          Energy          fmax
    MDMin:    0 09:41:54        5.433815        1.152831
    MDMin:    1 09:41:54        5.352153        1.037842
    MDMin:    2 09:41:54        5.179224        0.774873
    MDMin:    3 09:41:54        5.033743        0.467209
    MDMin:    4 09:41:54        4.952171        0.268936
    MDMin:    5 09:41:54        4.905037        0.159289
    MDMin:    6 09:41:54        4.877853        0.127328
    MDMin:    7 09:41:54        4.868893        0.137692
    MDMin:    8 09:41:54        4.868167        0.128909
    MDMin:    9 09:41:54        4.866425        0.108529
    MDMin:   10 09:41:54        4.864372        0.092540
    MDMin:   11 09:41:54        4.862194        0.087441
    MDMin:   12 09:41:54        4.859889        0.083497
    MDMin:   13 09:41:54        4.857550        0.078463
    MDMin:   14 09:41:54        4.855276        0.071014
    MDMin:   15 09:41:54        4.853106        0.058560
    MDMin:   16 09:41:54        4.851041        0.045468

    np.True_



.. GENERATED FROM PYTHON SOURCE LINES 383-384

Visualize the MEP in side view to see the motion

.. GENERATED FROM PYTHON SOURCE LINES 384-390

.. code-block:: Python


    for image in images:
        fig, ax = plt.subplots()
        plot_atoms(image, ax, rotation=('-60x, 10y,0z'))
        ax.set_axis_off()




.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_024.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_024.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_025.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_025.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_026.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_026.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_027.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_027.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_028.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_028.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_029.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_029.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_030.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_030.png
         :class: sphx-glr-multi-img





.. GENERATED FROM PYTHON SOURCE LINES 391-392

Plot the variation of potential energy

.. GENERATED FROM PYTHON SOURCE LINES 392-407

.. code-block:: Python

    potential_energies = [image.get_potential_energy() for image in images]

    fig, ax = plt.subplots()
    plt.plot(
        range(len(potential_energies)),
        potential_energies - potential_energies[0],
        marker='+',
    )
    plt.xlabel('Image number')
    plt.ylabel('Potential energy (eV)')

    diff = np.max(potential_energies) - potential_energies[0]
    print(f'The energy barrier is {diff:.4f} eV.')





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


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

 .. code-block:: none

    The energy barrier is 0.2394 eV.




.. GENERATED FROM PYTHON SOURCE LINES 408-415

.. hint::
  When opening a trajectory with :program:`ase gui` with calculated energies,
  the default plot window shows the energy versus frame number.
  To get a better feel of the energy barrier in an NEB calculation;
  choose  :menuselection:`Tools --> NEB`.
  This will give a smooth curve   of the energy as a function of the NEB path
  length, with the slope at each point estimated from the force.

.. GENERATED FROM PYTHON SOURCE LINES 418-433

Performing Dimer-method calculation
-----------------------------------

In the NEB calculations above we knew the final states, so all we had
to do was to calculate the path between the initial state
and the final state.
But in some cases we do not know the final state.
Then the :mod:`Dimer method <ase.mep.dimer>` can be used to
find the transition state.
The result of a Dimer calculation will hence not be the complete particle
trajectory as in the NEB output, but rather the configuration of
the transition-state image.

The following code will find the transition-state image of the
jump along the row.

.. GENERATED FROM PYTHON SOURCE LINES 434-474

.. code-block:: Python


    from ase.io import Trajectory
    from ase.mep import DimerControl, MinModeAtoms, MinModeTranslate

    a = 4.0614
    b = a / sqrt(2)
    h = b / 2
    initial = Atoms(
        'Al2',
        positions=[(0, 0, 0), (a / 2, b / 2, -h)],
        cell=(a, b, 2 * h),
        pbc=(1, 1, 0),
    )
    initial *= (2, 2, 2)
    initial.append(Atom('Al', (a / 2, b / 2, 3 * h)))
    initial.center(vacuum=4.0, axis=2)


    initial_copy = initial.copy()

    N = len(initial)  # number of atoms


    # Make a mask of zeros and ones that select fixed atoms - the two
    # bottom layers:
    mask = initial.positions[:, 2] - min(initial.positions[:, 2]) < 1.5 * h
    constraint = FixAtoms(mask=mask)
    initial.set_constraint(constraint)

    # Calculate using EMT:
    initial.calc = EMT()

    # Relax the initial state:
    QuasiNewton(initial).run(fmax=0.05)
    e0 = initial.get_potential_energy()

    # To save the trajectory file
    traj = Trajectory('dimer_along.traj', 'w', initial)
    traj.write()





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

 .. code-block:: none

                    Step[ FC]     Time          Energy          fmax
    BFGSLineSearch:    0[  0] 09:41:55        4.639215       0.2538
    BFGSLineSearch:    1[  2] 09:41:55        4.622063       0.1580
    BFGSLineSearch:    2[  4] 09:41:56        4.613340       0.0622
    BFGSLineSearch:    3[  5] 09:41:56        4.611640       0.0336




.. GENERATED FROM PYTHON SOURCE LINES 475-476

Making dimer mask list:

.. GENERATED FROM PYTHON SOURCE LINES 476-500

.. code-block:: Python

    d_mask = [False] * (N - 1) + [True]

    # Set up the dimer:
    d_control = DimerControl(
        initial_eigenmode_method='displacement',
        displacement_method='vector',
        logfile=None,
        mask=d_mask,
    )
    d_atoms = MinModeAtoms(initial, d_control)

    # Displacement settings:
    displacement_vector = np.zeros((N, 3))
    # Strength of displacement along y axis = along row:
    displacement_vector[-1, 1] = 0.001
    # The direction of the displacement is set by the a in
    # displacement_vector[-1, a], where a can be 0 for x, 1 for y and 2 for z.
    d_atoms.displace(displacement_vector=displacement_vector)

    # Converge to a saddle point:
    dim_rlx = MinModeTranslate(d_atoms, trajectory=traj, logfile=None)
    dim_rlx.run(fmax=0.001)






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

 .. code-block:: none


    np.True_



.. GENERATED FROM PYTHON SOURCE LINES 501-503

Visualize the Initial state and the saddle point in side view to see
the change

.. GENERATED FROM PYTHON SOURCE LINES 503-513

.. code-block:: Python


    for image in [initial, initial_copy]:
        fig, ax = plt.subplots()
        plot_atoms(image, ax, rotation=('-90x, 90y,0z'))
        ax.set_axis_off()


    diff = initial.get_potential_energy() - e0
    print(f'The energy barrier is {diff:.4f} eV.')




.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_032.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_032.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_033.png
         :alt: neb selfdiffusion
         :srcset: /examples_generated/tutorials/images/sphx_glr_neb_selfdiffusion_033.png
         :class: sphx-glr-multi-img


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

 .. code-block:: none

    The energy barrier is 0.1090 eV.




.. GENERATED FROM PYTHON SOURCE LINES 514-521

* Compare the transition-state images of the NEB and Dimer as viewed
  in the GUI. Are they identical?
* What is the energy barrier? How does it compare to the one found in the
  NEB calculation?
* Do the same as above for the jump across the row and the exchange
  process by copying and modifying the Dimer script, while remembering
  that you have to  give the relevant atoms a kick in a meaningful direction.


.. _sphx_glr_download_examples_generated_tutorials_neb_selfdiffusion.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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