.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_generated/tutorials/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_tutorials_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-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_selfdiffusion_001.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_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_selfdiffusion_002.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_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_selfdiffusion_003.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_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] 18:46:10 17.939914 0.2540 BFGSLineSearch: 1[ 2] 18:46:10 17.897422 0.1631 BFGSLineSearch: 2[ 4] 18:46:10 17.884347 0.0522 BFGSLineSearch: 3[ 5] 18:46:10 17.881710 0.0269 Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 18:46:10 17.939914 0.2540 BFGSLineSearch: 1[ 2] 18:46:10 17.897422 0.1631 BFGSLineSearch: 2[ 4] 18:46:10 17.884347 0.0522 BFGSLineSearch: 3[ 5] 18:46:10 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 /home/ase/.local/lib/python3.13/site-packages/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 18:46:10 18.237088 0.993544 MDMin: 1 18:46:10 18.190945 0.836011 MDMin: 2 18:46:10 18.105033 0.460778 MDMin: 3 18:46:10 18.056201 0.122138 MDMin: 4 18:46:10 18.035931 0.102285 MDMin: 5 18:46:10 18.022022 0.118145 MDMin: 6 18:46:10 18.010798 0.088852 MDMin: 7 18:46:10 18.005186 0.127228 MDMin: 8 18:46:10 18.003865 0.114049 MDMin: 9 18:46:10 18.001954 0.088821 MDMin: 10 18:46:10 18.000262 0.063304 MDMin: 11 18:46:10 17.998668 0.053980 MDMin: 12 18:46:10 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_selfdiffusion_004.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_004.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_005.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_005.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_006.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_006.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_007.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_007.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_008.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_008.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_009.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_009.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_010.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_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_selfdiffusion_011.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_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_selfdiffusion_012.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_012.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_013.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_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] 18:46:13 17.881710 0.0269 Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 18:46:13 17.898108 0.2165 BFGSLineSearch: 1[ 2] 18:46:13 17.883509 0.0755 BFGSLineSearch: 2[ 4] 18:46:13 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 18:46:13 21.182521 7.493573 MDMin: 1 18:46:13 19.417286 3.689216 MDMin: 2 18:46:13 18.645079 0.922870 MDMin: 3 18:46:13 18.533896 0.503897 MDMin: 4 18:46:13 18.525368 0.439764 MDMin: 5 18:46:13 18.507383 0.263744 MDMin: 6 18:46:13 18.491239 0.227551 MDMin: 7 18:46:13 18.479659 0.160835 MDMin: 8 18:46:13 18.468241 0.152303 MDMin: 9 18:46:13 18.458050 0.091391 MDMin: 10 18:46:13 18.450556 0.112507 MDMin: 11 18:46:14 18.447753 0.095248 MDMin: 12 18:46:14 18.447405 0.085207 MDMin: 13 18:46:14 18.446629 0.058262 MDMin: 14 18:46:14 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_selfdiffusion_014.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_014.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_015.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_015.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_016.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_016.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_017.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_017.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_018.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_018.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_019.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_019.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_020.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_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_selfdiffusion_021.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_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 `. 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_selfdiffusion_022.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_022.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_023.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_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] 18:46:16 4.639215 0.2538 BFGSLineSearch: 1[ 2] 18:46:16 4.622063 0.1580 BFGSLineSearch: 2[ 4] 18:46:16 4.613340 0.0622 BFGSLineSearch: 3[ 5] 18:46:16 4.611640 0.0336 Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 18:46:16 4.639215 0.2538 BFGSLineSearch: 1[ 2] 18:46:16 4.622063 0.1580 BFGSLineSearch: 2[ 4] 18:46:16 4.613340 0.0622 BFGSLineSearch: 3[ 5] 18:46:16 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 18:46:16 5.433815 1.152831 MDMin: 1 18:46:16 5.352153 1.037842 MDMin: 2 18:46:16 5.179224 0.774873 MDMin: 3 18:46:16 5.033743 0.467209 MDMin: 4 18:46:16 4.952171 0.268936 MDMin: 5 18:46:16 4.905037 0.159289 MDMin: 6 18:46:16 4.877853 0.127328 MDMin: 7 18:46:16 4.868893 0.137692 MDMin: 8 18:46:16 4.868167 0.128909 MDMin: 9 18:46:17 4.866425 0.108529 MDMin: 10 18:46:17 4.864372 0.092540 MDMin: 11 18:46:17 4.862194 0.087441 MDMin: 12 18:46:17 4.859889 0.083497 MDMin: 13 18:46:17 4.857550 0.078463 MDMin: 14 18:46:17 4.855276 0.071014 MDMin: 15 18:46:17 4.853106 0.058560 MDMin: 16 18:46:17 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_selfdiffusion_024.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_024.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_025.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_025.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_026.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_026.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_027.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_027.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_028.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_028.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_029.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_029.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_030.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_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_selfdiffusion_031.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_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 ` 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] 18:46:18 4.639215 0.2538 BFGSLineSearch: 1[ 2] 18:46:18 4.622063 0.1580 BFGSLineSearch: 2[ 4] 18:46:18 4.613340 0.0622 BFGSLineSearch: 3[ 5] 18:46:18 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_selfdiffusion_032.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_032.png :class: sphx-glr-multi-img * .. image-sg:: /examples_generated/tutorials/images/sphx_glr_selfdiffusion_033.png :alt: selfdiffusion :srcset: /examples_generated/tutorials/images/sphx_glr_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_selfdiffusion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: selfdiffusion.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: selfdiffusion.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: selfdiffusion.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_