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