.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_generated/tutorials/constraint_diffusion.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_constraint_diffusion.py: .. _constraints diffusion tutorial: ============================================================ Constrained Calculations - Surface diffusion energy barriers ============================================================ In this tutorial, we will calculate the energy barrier that was found using the :mod:`NEB ` method in the :ref:`diffusion tutorial` tutorial. Here, we use a simple :class:`~ase.constraints.FixedPlane` constraint that forces the Au atom to relax in the *yz*-plane only: .. GENERATED FROM PYTHON SOURCE LINES 12-24 .. code-block:: Python from ase.build import add_adsorbate, fcc100 from ase.calculators.emt import EMT from ase.constraints import FixAtoms, FixedPlane from ase.optimize import QuasiNewton # 2x2-Al(001) surface with 3 layers and an # Au atom adsorbed in a hollow site: 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 25-26 We can visualize the structure with ase visualize: .. GENERATED FROM PYTHON SOURCE LINES 26-38 .. code-block:: Python import matplotlib.pyplot as plt from ase.visualize.plot import plot_atoms fig, (ax1, ax2) = plt.subplots(1, 2) plot_atoms(slab, ax1) plot_atoms(slab, ax2, rotation='-90x') ax1.set_title('top view') ax2.set_title('side view') ax1.set_axis_off() ax2.set_axis_off() .. image-sg:: /examples_generated/tutorials/images/sphx_glr_constraint_diffusion_001.png :alt: top view, side view :srcset: /examples_generated/tutorials/images/sphx_glr_constraint_diffusion_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 39-40 Alternatively, you can use also use view directly: .. GENERATED FROM PYTHON SOURCE LINES 40-44 .. code-block:: Python # $ from ase.visualize import view # $ view(slab) .. GENERATED FROM PYTHON SOURCE LINES 45-46 We can now continue fixing the atoms in the slab: .. GENERATED FROM PYTHON SOURCE LINES 46-57 .. code-block:: Python # Fix second and third layers: mask = [atom.tag > 1 for atom in slab] # print(mask) fixlayers = FixAtoms(mask=mask) # Constrain the last atom (Au atom) to move only in the yz-plane: plane = FixedPlane(-1, (1, 0, 0)) slab.set_constraint([fixlayers, plane]) .. GENERATED FROM PYTHON SOURCE LINES 58-61 Now we can perform the calculation optimizing the displacement of the gold atom along the x-axis. We do structure optimization here using the EMT potential: .. GENERATED FROM PYTHON SOURCE LINES 61-83 .. code-block:: Python # Use EMT potential: slab.calc = EMT() for i in range(5): qn = QuasiNewton(slab, trajectory=f'mep{i}.traj') qn.run(fmax=0.05) # Move gold atom along x-axis: slab[-1].x += slab.get_cell()[0, 0] / 8 # Let's visualize the saved trajectory. # Here is code to visualize # a side-view of the path (unit cell repeated twice): from ase.io import read configs = [read(f'mep{i}.traj', '-1') for i in range(5)] # for easier visualization, let's repeat the structures configs_repeated = [config.repeat((2, 1, 1)) for config in configs] .. rst-class:: sphx-glr-script-out .. code-block:: none Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 10:57:13 3.323870 0.2462 BFGSLineSearch: 1[ 1] 10:57:13 3.314754 0.0378 Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 10:57:13 3.730166 1.6460 BFGSLineSearch: 1[ 2] 10:57:13 3.546230 0.2605 BFGSLineSearch: 2[ 4] 10:57:13 3.513795 0.1080 BFGSLineSearch: 3[ 6] 10:57:13 3.502279 0.0500 Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 10:57:13 3.764262 1.0444 BFGSLineSearch: 1[ 2] 10:57:13 3.708446 0.1203 BFGSLineSearch: 2[ 4] 10:57:13 3.695835 0.0668 BFGSLineSearch: 3[ 6] 10:57:13 3.691660 0.0618 BFGSLineSearch: 4[ 7] 10:57:13 3.690089 0.0300 Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 10:57:13 3.580034 0.6521 BFGSLineSearch: 1[ 1] 10:57:13 3.544894 0.1627 BFGSLineSearch: 2[ 3] 10:57:13 3.516850 0.1946 BFGSLineSearch: 3[ 4] 10:57:13 3.508545 0.1325 BFGSLineSearch: 4[ 5] 10:57:13 3.503316 0.1022 BFGSLineSearch: 5[ 6] 10:57:13 3.500392 0.0550 BFGSLineSearch: 6[ 8] 10:57:13 3.499584 0.0392 Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 10:57:13 3.463634 0.5332 BFGSLineSearch: 1[ 2] 10:57:13 3.356702 0.4028 BFGSLineSearch: 2[ 3] 10:57:13 3.333298 0.1536 BFGSLineSearch: 3[ 4] 10:57:13 3.322604 0.1003 BFGSLineSearch: 4[ 5] 10:57:13 3.318477 0.0658 BFGSLineSearch: 5[ 6] 10:57:13 3.316534 0.0483 .. GENERATED FROM PYTHON SOURCE LINES 84-85 We can visualize the structures with `ase.visualize.plot.animate`: .. GENERATED FROM PYTHON SOURCE LINES 85-95 .. code-block:: Python from ase.visualize.plot import animate animate( configs_repeated, ax=None, interval=500, # in ms; same default value as in FuncAnimation rotation=('-90x,0y,0z'), ) .. video:: /examples_generated/tutorials/images/sphx_glr_constraint_diffusion_002.mp4 :class: sphx-glr-single-img :height: 480 :width: 640 :autoplay: :loop: .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 96-97 Let's plot the energy and look at the barrier. .. GENERATED FROM PYTHON SOURCE LINES 97-107 .. code-block:: Python # get the potential energies of the structures energies = [config.get_potential_energy() for config in configs] # set last energy value to 0 for easier comparison energies = [energy - energies[-1] for energy in energies] plt.ylabel(r'$E(i) - E_{\mathrm{final}}$ (meV)') plt.xlabel('Image number') plt.plot(range(1, len(energies) + 1), energies) .. image-sg:: /examples_generated/tutorials/images/sphx_glr_constraint_diffusion_003.png :alt: constraint diffusion :srcset: /examples_generated/tutorials/images/sphx_glr_constraint_diffusion_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 108-115 The barrier is found to be 0.35 eV - exactly as in the :ref:`NEB ` tutorial. The result can also be analysed with the command :command:`ase gui mep?.traj -n -1` (choose :menuselection:`Tools --> NEB`). .. GENERATED FROM PYTHON SOURCE LINES 118-124 .. seealso:: * :mod:`ase.mep.neb` * :mod:`ase.constraints` * :ref:`diffusion tutorial` * :func:`~ase.build.fcc100` .. _sphx_glr_download_examples_generated_tutorials_constraint_diffusion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: constraint_diffusion.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: constraint_diffusion.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: constraint_diffusion.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_