.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_generated/03-tutorials/neb_selfdiffusion.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_generated_03-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 ` 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-59 .. 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')) .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_001.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 60-61 Then, lets add the Al addatom which will be the one moving on the surface. .. GENERATED FROM PYTHON SOURCE LINES 61-63 .. code-block:: Python initial.append(Atom('Al', (a / 2, b / 2, 3 * h))) .. GENERATED FROM PYTHON SOURCE LINES 64-65 Center the cell in vacuum along the z axis .. GENERATED FROM PYTHON SOURCE LINES 65-67 .. code-block:: Python initial.center(vacuum=4.0, axis=2) .. GENERATED FROM PYTHON SOURCE LINES 68-69 Visualize the new atom in the cell .. GENERATED FROM PYTHON SOURCE LINES 69-73 .. code-block:: Python fig, ax = plt.subplots() plot_atoms(initial, ax, rotation=('-60x, 10y,0z')) .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_002.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 74-84 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 84-86 .. code-block:: Python final = initial.copy() .. GENERATED FROM PYTHON SOURCE LINES 87-89 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 89-91 .. code-block:: Python final.positions[-1, 1] += b .. GENERATED FROM PYTHON SOURCE LINES 92-93 Visualize the new atom in the cell .. GENERATED FROM PYTHON SOURCE LINES 93-97 .. code-block:: Python fig, ax = plt.subplots() plot_atoms(final, ax, rotation=('-60x, 10y,0z')) .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_003.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 98-102 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 102-110 .. 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 111-114 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 114-120 .. code-block:: Python initial.calc = EMT() initial.set_constraint(constraint) final.calc = EMT() final.set_constraint(constraint) .. GENERATED FROM PYTHON SOURCE LINES 121-124 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 124-130 .. 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] 15:27:24 17.939914 0.2540 BFGSLineSearch: 1[ 2] 15:27:24 17.897422 0.1631 BFGSLineSearch: 2[ 4] 15:27:24 17.884347 0.0522 BFGSLineSearch: 3[ 5] 15:27:24 17.881710 0.0269 Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 15:27:24 17.939914 0.2540 BFGSLineSearch: 1[ 2] 15:27:24 17.897422 0.1631 BFGSLineSearch: 2[ 4] 15:27:24 17.884347 0.0522 BFGSLineSearch: 3[ 5] 15:27:24 17.881710 0.0269 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 131-133 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 133-138 .. code-block:: Python images = [initial] for i in range(5): images.append(initial.copy()) images.append(final) .. GENERATED FROM PYTHON SOURCE LINES 139-141 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 141-146 .. code-block:: Python for image in images: image.calc = EMT() image.set_constraint(constraint) .. GENERATED FROM PYTHON SOURCE LINES 147-148 Create a Nudged Elastic Band (:class:`~ase.mep import NEB`) object .. GENERATED FROM PYTHON SOURCE LINES 148-152 .. code-block:: Python from ase.mep import NEB neb = NEB(images) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/ase/.local/lib/python3.14/site-packages/ase/mep/neb.py:410: 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 153-155 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 155-157 .. code-block:: Python neb.interpolate() .. GENERATED FROM PYTHON SOURCE LINES 158-159 Perform the NEB calculation minimizing the force below 0.05 eV/A .. GENERATED FROM PYTHON SOURCE LINES 159-164 .. 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 15:27:24 18.237088 0.993544 MDMin: 1 15:27:24 18.190945 0.836011 MDMin: 2 15:27:24 18.105033 0.460778 MDMin: 3 15:27:24 18.056201 0.122138 MDMin: 4 15:27:24 18.035931 0.102285 MDMin: 5 15:27:24 18.022022 0.118145 MDMin: 6 15:27:24 18.010798 0.088852 MDMin: 7 15:27:24 18.005186 0.127228 MDMin: 8 15:27:24 18.003865 0.114049 MDMin: 9 15:27:24 18.001954 0.088821 MDMin: 10 15:27:24 18.000262 0.063304 MDMin: 11 15:27:24 17.998668 0.053980 MDMin: 12 15:27:24 17.997130 0.040716 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 165-167 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 167-171 .. code-block:: Python for image in images: fig, ax = plt.subplots() plot_atoms(image, ax, rotation=('-90x, 90y,0z')) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_004.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_004.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_005.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_005.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_006.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_006.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_007.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_007.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_008.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_008.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_009.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_009.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_010.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_010.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 172-173 Plot the variation of potential energy .. GENERATED FROM PYTHON SOURCE LINES 173-187 .. 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/03-tutorials/images/sphx_glr_neb_selfdiffusion_011.png :alt: neb selfdiffusion :srcset: /examples_generated/03-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 188-190 You can visualize the NEB path using ASE GUI after saving the tajectory in a file .. GENERATED FROM PYTHON SOURCE LINES 190-195 .. code-block:: Python from ase.io import write write('neb_path.traj', images, format='traj') .. GENERATED FROM PYTHON SOURCE LINES 196-199 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 201-206 * 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 209-216 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 218-223 .. code-block:: Python final = initial.copy() final.positions[-1, 0] += a .. GENERATED FROM PYTHON SOURCE LINES 224-225 Plot the images .. GENERATED FROM PYTHON SOURCE LINES 225-229 .. code-block:: Python for image in [initial, final]: fig, ax = plt.subplots() plot_atoms(image, ax, rotation=('-90x, 0y, 0z')) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_012.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_012.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_013.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_013.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 230-231 Construct a list of images: .. GENERATED FROM PYTHON SOURCE LINES 232-237 .. code-block:: Python images = [initial] for i in range(5): images.append(initial.copy()) images.append(final) .. GENERATED FROM PYTHON SOURCE LINES 238-240 Make a mask of zeros and ones that select fixed atoms (the two bottom layers): .. GENERATED FROM PYTHON SOURCE LINES 240-249 .. 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 250-251 Relax the initial and final states: .. GENERATED FROM PYTHON SOURCE LINES 251-254 .. 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] 15:27:27 17.881710 0.0269 Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 15:27:27 17.898108 0.2165 BFGSLineSearch: 1[ 2] 15:27:27 17.883509 0.0755 BFGSLineSearch: 2[ 4] 15:27:27 17.880775 0.0216 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 255-256 Create a Nudged Elastic Band: .. GENERATED FROM PYTHON SOURCE LINES 256-258 .. code-block:: Python neb = NEB(images) .. GENERATED FROM PYTHON SOURCE LINES 259-261 Make a starting guess for the minimum energy path (a straight line from the initial to the final state): .. GENERATED FROM PYTHON SOURCE LINES 261-263 .. code-block:: Python neb.interpolate() .. GENERATED FROM PYTHON SOURCE LINES 264-265 Relax the NEB path: .. GENERATED FROM PYTHON SOURCE LINES 265-269 .. 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 15:27:27 21.182521 7.493573 MDMin: 1 15:27:27 19.417286 3.689216 MDMin: 2 15:27:28 18.645079 0.922870 MDMin: 3 15:27:28 18.533896 0.503897 MDMin: 4 15:27:28 18.525368 0.439764 MDMin: 5 15:27:28 18.507383 0.263744 MDMin: 6 15:27:28 18.491239 0.227551 MDMin: 7 15:27:28 18.479659 0.160835 MDMin: 8 15:27:28 18.468241 0.152303 MDMin: 9 15:27:28 18.458050 0.091391 MDMin: 10 15:27:28 18.450556 0.112507 MDMin: 11 15:27:28 18.447753 0.095248 MDMin: 12 15:27:28 18.447405 0.085207 MDMin: 13 15:27:28 18.446629 0.058262 MDMin: 14 15:27:28 18.445840 0.046679 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 270-272 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 272-278 .. code-block:: Python for image in images: fig, ax = plt.subplots() plot_atoms(image, ax, rotation=('-90x, 0y, 0z')) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_014.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_014.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_015.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_015.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_016.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_016.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_017.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_017.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_018.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_018.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_019.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_019.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_020.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_020.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 279-280 Plot the variation of potential energy .. GENERATED FROM PYTHON SOURCE LINES 280-295 .. 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/03-tutorials/images/sphx_glr_neb_selfdiffusion_021.png :alt: neb selfdiffusion :srcset: /examples_generated/03-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 296-298 Finding the third mechanism ---------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 302-308 A third diffusion process can be found: Diffusion by an exchange process. You can read more about it in the paper listed :mod:`here `. 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 308-329 .. 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 330-331 Visualize the initial and final images .. GENERATED FROM PYTHON SOURCE LINES 331-335 .. code-block:: Python for image in [initial, final]: fig, ax = plt.subplots() plot_atoms(image, ax, rotation=('-60x, 10y,0z')) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_022.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_022.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_023.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_023.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 336-337 Construct a list of images: .. GENERATED FROM PYTHON SOURCE LINES 337-342 .. code-block:: Python images = [initial] for i in range(5): images.append(initial.copy()) images.append(final) .. GENERATED FROM PYTHON SOURCE LINES 343-345 Make a mask of zeros and ones that select fixed atoms (the two bottom layers): .. GENERATED FROM PYTHON SOURCE LINES 345-349 .. 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 350-351 Let all images use an EMT calculator: .. GENERATED FROM PYTHON SOURCE LINES 351-355 .. code-block:: Python for image in images: image.calc = EMT() image.set_constraint(constraint) .. GENERATED FROM PYTHON SOURCE LINES 356-357 Relax the initial and final states: .. GENERATED FROM PYTHON SOURCE LINES 357-360 .. 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] 15:27:31 4.639215 0.2538 BFGSLineSearch: 1[ 2] 15:27:31 4.622063 0.1580 BFGSLineSearch: 2[ 4] 15:27:31 4.613340 0.0622 BFGSLineSearch: 3[ 5] 15:27:31 4.611640 0.0336 Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 15:27:31 4.639215 0.2538 BFGSLineSearch: 1[ 2] 15:27:31 4.622063 0.1580 BFGSLineSearch: 2[ 4] 15:27:31 4.613340 0.0622 BFGSLineSearch: 3[ 5] 15:27:31 4.611640 0.0336 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 361-362 Create a Nudged Elastic Band: .. GENERATED FROM PYTHON SOURCE LINES 362-364 .. code-block:: Python neb = NEB(images) .. GENERATED FROM PYTHON SOURCE LINES 365-367 Make a starting guess for the minimum energy path (a straight line from the initial to the final state): .. GENERATED FROM PYTHON SOURCE LINES 367-369 .. code-block:: Python neb.interpolate() .. GENERATED FROM PYTHON SOURCE LINES 370-371 Relax the NEB path: .. GENERATED FROM PYTHON SOURCE LINES 371-375 .. 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 15:27:31 5.433815 1.152831 MDMin: 1 15:27:31 5.352153 1.037842 MDMin: 2 15:27:31 5.179224 0.774873 MDMin: 3 15:27:31 5.033743 0.467209 MDMin: 4 15:27:31 4.952171 0.268936 MDMin: 5 15:27:31 4.905037 0.159289 MDMin: 6 15:27:31 4.877853 0.127328 MDMin: 7 15:27:31 4.868893 0.137692 MDMin: 8 15:27:31 4.868167 0.128909 MDMin: 9 15:27:31 4.866425 0.108529 MDMin: 10 15:27:31 4.864372 0.092540 MDMin: 11 15:27:31 4.862194 0.087441 MDMin: 12 15:27:31 4.859889 0.083497 MDMin: 13 15:27:31 4.857550 0.078463 MDMin: 14 15:27:31 4.855276 0.071014 MDMin: 15 15:27:31 4.853106 0.058560 MDMin: 16 15:27:31 4.851041 0.045468 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 376-377 Visualize the MEP in side view to see the motion .. GENERATED FROM PYTHON SOURCE LINES 377-382 .. code-block:: Python for image in images: fig, ax = plt.subplots() plot_atoms(image, ax, rotation=('-60x, 10y,0z')) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_024.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_024.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_025.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_025.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_026.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_026.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_027.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_027.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_028.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_028.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_029.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_029.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_030.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_030.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 383-384 Plot the variation of potential energy .. GENERATED FROM PYTHON SOURCE LINES 384-399 .. 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/03-tutorials/images/sphx_glr_neb_selfdiffusion_031.png :alt: neb selfdiffusion :srcset: /examples_generated/03-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 400-407 .. 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 410-425 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 ` 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 426-466 .. 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] 15:27:33 4.639215 0.2538 BFGSLineSearch: 1[ 2] 15:27:33 4.622063 0.1580 BFGSLineSearch: 2[ 4] 15:27:33 4.613340 0.0622 BFGSLineSearch: 3[ 5] 15:27:33 4.611640 0.0336 .. GENERATED FROM PYTHON SOURCE LINES 467-468 Making dimer mask list: .. GENERATED FROM PYTHON SOURCE LINES 468-492 .. 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 493-495 Visualize the Initial state and the saddle point in side view to see the change .. GENERATED FROM PYTHON SOURCE LINES 495-504 .. code-block:: Python for image in [initial, initial_copy]: fig, ax = plt.subplots() plot_atoms(image, ax, rotation=('-90x, 90y,0z')) diff = initial.get_potential_energy() - e0 print(f'The energy barrier is {diff:.4f} eV.') .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_032.png :alt: neb selfdiffusion :srcset: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_032.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_neb_selfdiffusion_033.png :alt: neb selfdiffusion :srcset: /examples_generated/03-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 505-512 * 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_03-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 ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: neb_selfdiffusion.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: neb_selfdiffusion.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_