.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_generated/03-tutorials/diffusion_neb.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_diffusion_neb.py: .. _diffusion tutorial: ============================================================================= Surface diffusion energy barriers using the Nudged Elastic Band (NEB) method ============================================================================= .. GENERATED FROM PYTHON SOURCE LINES 10-21 .. code-block:: Python import matplotlib.pyplot as plt from ase.build import add_adsorbate, fcc100 from ase.calculators.emt import EMT from ase.constraints import FixAtoms from ase.io import read from ase.mep import NEB from ase.optimize import BFGS, QuasiNewton from ase.parallel import world from ase.visualize.plot import plot_atoms .. GENERATED FROM PYTHON SOURCE LINES 22-25 First, set up the initial and final states: 2x2-Al(001) surface with 3 layers and an Au atom adsorbed in a hollow site: .. GENERATED FROM PYTHON SOURCE LINES 27-31 .. code-block:: Python slab = fcc100('Al', size=(2, 2, 3)) add_adsorbate(slab, 'Au', 1.7, 'hollow') slab.center(axis=2, vacuum=4.0) .. GENERATED FROM PYTHON SOURCE LINES 32-33 Make sure the structure is correct: .. GENERATED FROM PYTHON SOURCE LINES 35-38 .. code-block:: Python fig, ax = plt.subplots() plot_atoms(slab, ax) .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_diffusion_neb_001.png :alt: diffusion neb :srcset: /examples_generated/03-tutorials/images/sphx_glr_diffusion_neb_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 39-40 Fix second and third layers: .. GENERATED FROM PYTHON SOURCE LINES 42-46 .. code-block:: Python mask = [atom.tag > 1 for atom in slab] print(mask) slab.set_constraint(FixAtoms(mask=mask)) .. rst-class:: sphx-glr-script-out .. code-block:: none [np.True_, np.True_, np.True_, np.True_, np.True_, np.True_, np.True_, np.True_, np.False_, np.False_, np.False_, np.False_, np.False_] .. GENERATED FROM PYTHON SOURCE LINES 47-48 Use EMT potential: .. GENERATED FROM PYTHON SOURCE LINES 50-52 .. code-block:: Python slab.calc = EMT() .. GENERATED FROM PYTHON SOURCE LINES 53-54 Initial state: .. GENERATED FROM PYTHON SOURCE LINES 56-59 .. code-block:: Python qn = QuasiNewton(slab, trajectory='initial.traj') qn.run(fmax=0.05) .. rst-class:: sphx-glr-script-out .. code-block:: none Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 15:25:49 3.323870 0.2462 BFGSLineSearch: 1[ 1] 15:25:49 3.314754 0.0378 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 60-61 Final state: .. GENERATED FROM PYTHON SOURCE LINES 63-67 .. code-block:: Python slab[-1].x += slab.get_cell()[0, 0] / 2 qn = QuasiNewton(slab, trajectory='final.traj') qn.run(fmax=0.05) .. rst-class:: sphx-glr-script-out .. code-block:: none Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 15:25:49 3.320051 0.1208 BFGSLineSearch: 1[ 1] 15:25:49 3.316117 0.0474 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 68-71 .. note:: Notice how the tags are used to select the constrained atoms Now, do the NEB calculation: .. GENERATED FROM PYTHON SOURCE LINES 73-92 .. code-block:: Python initial = read('initial.traj') final = read('final.traj') constraint = FixAtoms(mask=[atom.tag > 1 for atom in initial]) images = [initial] for i in range(3): image = initial.copy() image.calc = EMT() image.set_constraint(constraint) images.append(image) images.append(final) neb = NEB(images, method='improvedtangent') neb.interpolate() qn = BFGS(neb, trajectory='neb.traj') qn.run(fmax=0.05) .. rst-class:: sphx-glr-script-out .. code-block:: none Step Time Energy fmax BFGS: 0 15:25:49 4.219952 3.520752 BFGS: 1 15:25:49 3.937039 2.176460 BFGS: 2 15:25:49 3.719814 0.435144 BFGS: 3 15:25:49 3.709652 0.230134 BFGS: 4 15:25:49 3.708878 0.244118 BFGS: 5 15:25:49 3.706087 0.257725 BFGS: 6 15:25:49 3.698531 0.213294 BFGS: 7 15:25:49 3.692120 0.246238 BFGS: 8 15:25:49 3.692278 0.187553 BFGS: 9 15:25:49 3.693494 0.172952 BFGS: 10 15:25:49 3.692670 0.151720 BFGS: 11 15:25:49 3.690813 0.073884 BFGS: 12 15:25:49 3.690200 0.071211 BFGS: 13 15:25:49 3.690380 0.078122 BFGS: 14 15:25:49 3.690427 0.103668 BFGS: 15 15:25:49 3.689897 0.100381 BFGS: 16 15:25:49 3.689034 0.055057 BFGS: 17 15:25:49 3.688737 0.028898 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 93-94 Visualize the results with:: .. GENERATED FROM PYTHON SOURCE LINES 94-97 .. code-block:: Python fig, ax = plt.subplots() plot_atoms(final, ax) .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_diffusion_neb_002.png :alt: diffusion neb :srcset: /examples_generated/03-tutorials/images/sphx_glr_diffusion_neb_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 98-103 or from the command line: $ ase gui neb.traj@-5: and select Tools->NEB. .. GENERATED FROM PYTHON SOURCE LINES 106-114 You can also create a series of plots like above, that show the progression of the NEB relaxation, directly at the command line:: $ ase nebplot --share-x --share-y neb.traj For more customizable analysis of the output of many NEB jobs, you can use the :class:`ase.mep.NEBTools` class. Some examples of its use are below; the final example was used to make the figure you see above. .. GENERATED FROM PYTHON SOURCE LINES 116-125 .. code-block:: Python import matplotlib.pyplot as plt from ase.io import read from ase.mep import NEBTools images = read('neb.traj@-5:') nebtools = NEBTools(images) .. GENERATED FROM PYTHON SOURCE LINES 126-127 Get the calculated barrier and the energy change of the reaction. .. GENERATED FROM PYTHON SOURCE LINES 129-131 .. code-block:: Python Ef, dE = nebtools.get_barrier() .. GENERATED FROM PYTHON SOURCE LINES 132-133 Get the barrier without any interpolation between highest images. .. GENERATED FROM PYTHON SOURCE LINES 135-137 .. code-block:: Python Ef, dE = nebtools.get_barrier(fit=False) .. GENERATED FROM PYTHON SOURCE LINES 138-139 Get the actual maximum force at this point in the simulation. .. GENERATED FROM PYTHON SOURCE LINES 141-143 .. code-block:: Python max_force = nebtools.get_fmax() .. 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 144-145 Create a figure like that coming from ASE-GUI. .. GENERATED FROM PYTHON SOURCE LINES 147-149 .. code-block:: Python fig = nebtools.plot_band() .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_diffusion_neb_003.png :alt: $E_\mathrm{f} \approx$ 0.374 eV; $E_\mathrm{r} \approx$ 0.373 eV; $\Delta E$ = 0.001 eV :srcset: /examples_generated/03-tutorials/images/sphx_glr_diffusion_neb_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 150-151 Create a figure with custom parameters. .. GENERATED FROM PYTHON SOURCE LINES 153-157 .. code-block:: Python fig = plt.figure(figsize=(5.5, 4.0)) ax = fig.add_axes((0.15, 0.15, 0.8, 0.75)) nebtools.plot_band(ax) .. image-sg:: /examples_generated/03-tutorials/images/sphx_glr_diffusion_neb_004.png :alt: $E_\mathrm{f} \approx$ 0.374 eV; $E_\mathrm{r} \approx$ 0.373 eV; $\Delta E$ = 0.001 eV :srcset: /examples_generated/03-tutorials/images/sphx_glr_diffusion_neb_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none
.. GENERATED FROM PYTHON SOURCE LINES 158-172 .. note:: For this reaction, the reaction coordinate is very simple: The *x*-coordinate of the Au atom. In such cases, the NEB method is overkill, and a simple constraint method should be used like in this tutorial: :ref:`constraints diffusion tutorial`. .. seealso:: * :mod:`ase.mep.neb` * :mod:`ase.constraints` * :ref:`constraints diffusion tutorial` * :func:`~ase.build.fcc100` .. GENERATED FROM PYTHON SOURCE LINES 174-178 Restarting NEB ============== Restart NEB from the trajectory file: .. GENERATED FROM PYTHON SOURCE LINES 180-181 read the last structures (of 5 images used in NEB) .. GENERATED FROM PYTHON SOURCE LINES 183-192 .. code-block:: Python images = read('neb.traj@-5:') for i in range(1, len(images) - 1): images[i].calc = EMT() neb = NEB(images, method='improvedtangent') qn = BFGS(neb, trajectory='neb_restart.traj') qn.run(fmax=0.005) .. rst-class:: sphx-glr-script-out .. code-block:: none Step Time Energy fmax BFGS: 0 15:25:49 3.688737 0.028898 BFGS: 1 15:25:49 3.688734 0.026547 BFGS: 2 15:25:49 3.688731 0.021515 BFGS: 3 15:25:49 3.688730 0.021379 BFGS: 4 15:25:49 3.688718 0.009577 BFGS: 5 15:25:49 3.688718 0.008389 BFGS: 6 15:25:49 3.688718 0.009122 BFGS: 7 15:25:49 3.688718 0.009290 BFGS: 8 15:25:49 3.688716 0.007568 BFGS: 9 15:25:49 3.688714 0.005563 BFGS: 10 15:25:49 3.688715 0.002547 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 193-205 Parallelizing over images with MPI ================================== Instead of having one process do the calculations for all three internal images in turn, it will be faster to have three processes do one image each. In order to be able to run python with MPI you need a special parallel python interpreter, for example gpaw python (see `GPAW parallel runs `_) and set ``parallel=True`` in the NEB calculation. The example below can then be run with ``gpaw -P3 python diffusion_parallel.py``: .. GENERATED FROM PYTHON SOURCE LINES 207-227 .. code-block:: Python initial = read('initial.traj') final = read('final.traj') constraint = FixAtoms(mask=[atom.tag > 1 for atom in initial]) n_images = 1 # Set to number of processes you use with mpiexec images = [initial] j = world.rank * n_images // world.size # my image number for i in range(n_images): image = initial.copy() if i == j: image.calc = EMT() image.set_constraint(constraint) images.append(image) images.append(final) neb = NEB(images, parallel=True, method='improvedtangent') neb.interpolate() qn = BFGS(neb, trajectory='neb.traj') qn.run(fmax=0.05) .. rst-class:: sphx-glr-script-out .. code-block:: none Step Time Energy fmax BFGS: 0 15:25:49 4.219952 3.520752 BFGS: 1 15:25:49 3.937039 2.176459 BFGS: 2 15:25:49 3.730687 0.621235 BFGS: 3 15:25:49 3.710767 0.211379 BFGS: 4 15:25:49 3.707518 0.163942 BFGS: 5 15:25:50 3.704890 0.187028 BFGS: 6 15:25:50 3.699676 0.197511 BFGS: 7 15:25:50 3.694802 0.156084 BFGS: 8 15:25:50 3.691892 0.081850 BFGS: 9 15:25:50 3.691165 0.050493 BFGS: 10 15:25:50 3.690969 0.048196 np.True_ .. _sphx_glr_download_examples_generated_03-tutorials_diffusion_neb.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: diffusion_neb.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: diffusion_neb.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: diffusion_neb.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_