.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_generated/tutorials/neb_dissociation.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_neb_dissociation.py: .. _dissociation: NEB: Dissociation of a molecule using the NEB method ==================================================== In this tutorial we provide an illustrative example of a nudged-elastic band (NEB) calculation. For more information on the NEB technique, see :mod:`ase.mep.neb`. We consider the dissociation of a nitrogen molecule on the Cu (111) surface. The first step is to find the relaxed structures of the initial and final states. .. GENERATED FROM PYTHON SOURCE LINES 15-30 .. code-block:: Python import matplotlib.animation as animation import matplotlib.pyplot as plt import numpy as np from ase import Atoms from ase.build import add_adsorbate, fcc111 from ase.calculators.emt import EMT from ase.constraints import FixAtoms from ase.io import read, write from ase.mep import NEB from ase.optimize import QuasiNewton from ase.optimize.fire import FIRE from ase.visualize.plot import plot_atoms .. GENERATED FROM PYTHON SOURCE LINES 31-33 First, we create the initial state: An N\ :sub:`2`\ molecule on a Cu(111) slab .. GENERATED FROM PYTHON SOURCE LINES 33-44 .. code-block:: Python # Set up a (4 x 4) two layer Cu (111) slab slab = fcc111('Cu', size=(4, 4, 2)) slab.set_pbc((1, 1, 0)) # Add the N2 Molecule, oriented at 60 degrees: d = 1.10 # N2 bond length n2_mol = Atoms( 'N2', positions=[[0.0, 0.0, 0.0], [0.5 * 3**0.5 * d, 0.5 * d, 0.0]] ) .. GENERATED FROM PYTHON SOURCE LINES 45-47 Then we put the adsorbate onto the surface .. GENERATED FROM PYTHON SOURCE LINES 47-50 .. code-block:: Python add_adsorbate(slab, n2_mol, height=1.0, position='fcc') .. GENERATED FROM PYTHON SOURCE LINES 51-52 And add a calculator for the forces and energies: .. GENERATED FROM PYTHON SOURCE LINES 52-55 .. code-block:: Python slab.calc = EMT() .. GENERATED FROM PYTHON SOURCE LINES 56-58 We don't want to worry about the Cu degrees of freedom, so fix these atoms: .. GENERATED FROM PYTHON SOURCE LINES 58-62 .. code-block:: Python mask = [atom.symbol == 'Cu' for atom in slab] slab.set_constraint(FixAtoms(mask=mask)) .. GENERATED FROM PYTHON SOURCE LINES 63-64 Then we relax the structure and write the trajectory into a file .. GENERATED FROM PYTHON SOURCE LINES 64-70 .. code-block:: Python relax = QuasiNewton(slab) relax.run(fmax=0.05) print('initial state:', slab.get_potential_energy()) write('N2.traj', slab) .. rst-class:: sphx-glr-script-out .. code-block:: none Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 12:40:36 14.459066 9.8897 BFGSLineSearch: 1[ 1] 12:40:36 12.244808 4.3275 BFGSLineSearch: 2[ 2] 12:40:36 11.343194 2.1646 BFGSLineSearch: 3[ 3] 12:40:36 11.077854 0.6379 BFGSLineSearch: 4[ 4] 12:40:36 11.056067 0.2425 BFGSLineSearch: 5[ 6] 12:40:36 11.049129 0.2142 BFGSLineSearch: 6[ 8] 12:40:36 11.034808 0.1607 BFGSLineSearch: 7[ 9] 12:40:36 11.032950 0.0397 initial state: 11.032950397754657 .. GENERATED FROM PYTHON SOURCE LINES 71-73 Now the final state. Move the second N atom to a neighboring hollow site and relax. .. GENERATED FROM PYTHON SOURCE LINES 73-81 .. code-block:: Python slab[-1].position[0] = slab[-2].position[0] + 0.25 * slab.cell[0, 0] slab[-1].position[1] = slab[-2].position[1] relax.run() print('final state: ', slab.get_potential_energy()) write('2N.traj', slab) .. rst-class:: sphx-glr-script-out .. code-block:: none BFGSLineSearch: 8[ 11] 12:40:36 11.023212 0.0971 BFGSLineSearch: 9[ 13] 12:40:36 11.021738 0.1053 BFGSLineSearch: 10[ 15] 12:40:36 11.017702 0.0733 BFGSLineSearch: 11[ 16] 12:40:36 11.017285 0.0323 final state: 11.017284554450073 .. GENERATED FROM PYTHON SOURCE LINES 82-89 Having obtained these structures we set up an NEB calculation with 9 images. Using :func:`~neb.interpolate()` provides a guess for the path between the initial and final states. We perform the relaxation of the images and obtain the intermediate steps. First, we read the previous configurations .. GENERATED FROM PYTHON SOURCE LINES 89-93 .. code-block:: Python initial = read('N2.traj') final = read('2N.traj') .. GENERATED FROM PYTHON SOURCE LINES 94-96 Then, we make 9 images (note the use of copy) and, as before, we fix the Cu atoms .. GENERATED FROM PYTHON SOURCE LINES 96-104 .. code-block:: Python configs = [initial.copy() for i in range(8)] + [final] constraint = FixAtoms(mask=[atom.symbol != 'N' for atom in initial]) for config in configs: config.calc = EMT() config.set_constraint(constraint) .. GENERATED FROM PYTHON SOURCE LINES 105-107 Next, we make the NEB object, and call its interpolate() method to guess the intermediate steps .. GENERATED FROM PYTHON SOURCE LINES 107-111 .. code-block:: Python band = NEB(configs) band.interpolate() .. 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 112-116 Then we relax the configurations in the NEB object. We can use same FIRE optimization class we might apply to optimize an Atoms object; NEB presents appropriate degrees of freedom and gradients for this to optimise multiple configurations simultaneously. .. GENERATED FROM PYTHON SOURCE LINES 116-120 .. code-block:: Python relax = FIRE(band) relax.run() .. rst-class:: sphx-glr-script-out .. code-block:: none Step Time Energy fmax FIRE: 0 12:40:36 11.066750 0.523744 FIRE: 1 12:40:36 11.063989 0.462236 FIRE: 2 12:40:36 11.059699 0.348578 FIRE: 3 12:40:36 11.055743 0.199917 FIRE: 4 12:40:36 11.053732 0.053868 FIRE: 5 12:40:36 11.054186 0.126908 FIRE: 6 12:40:36 11.054140 0.123746 FIRE: 7 12:40:36 11.054054 0.117517 FIRE: 8 12:40:36 11.053934 0.108416 FIRE: 9 12:40:36 11.053791 0.096753 FIRE: 10 12:40:36 11.053638 0.082986 FIRE: 11 12:40:36 11.053488 0.067787 FIRE: 12 12:40:36 11.053353 0.052238 FIRE: 13 12:40:36 11.053230 0.039264 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 121-122 Finally, we compare intermediate steps to the initial energy .. GENERATED FROM PYTHON SOURCE LINES 122-136 .. code-block:: Python energy_initial = initial.get_potential_energy() n2_distances = [] energy_differences = [] for i, config in enumerate(configs): n2_distance = np.linalg.norm(config[-2].position - config[-1].position) delta_energy = config.get_potential_energy() - energy_initial n2_distances.append(n2_distance) energy_differences.append(delta_energy) print( f'{i:>3}\td = {n2_distance:>4.3f}AA ' f'\tdelta energy = {delta_energy:>5.3f}eV' ) .. rst-class:: sphx-glr-script-out .. code-block:: none 0 d = 1.661AA delta energy = 0.000eV 1 d = 1.739AA delta energy = 0.003eV 2 d = 1.829AA delta energy = 0.012eV 3 d = 1.936AA delta energy = 0.019eV 4 d = 2.055AA delta energy = 0.020eV 5 d = 2.183AA delta energy = 0.015eV 6 d = 2.319AA delta energy = 0.004eV 7 d = 2.460AA delta energy = -0.010eV 8 d = 2.612AA delta energy = -0.016eV .. GENERATED FROM PYTHON SOURCE LINES 137-141 After the calculation is complete, the energy difference with respect to the initial state is given for each image, as well as the distance between the N atoms. .. GENERATED FROM PYTHON SOURCE LINES 141-176 .. code-block:: Python energy_differences = [e * 1.0e03 for e in energy_differences] # in meV fig, axs = plt.subplots( nrows=2, gridspec_kw={'height_ratios': [0.5, 1.0]}, ) axs[0].plot(n2_distances, energy_differences) scat = axs[0].scatter(n2_distances[0], energy_differences[0], s=10.0**2.0) axs[0].set_ylabel(r'$E(i) - E_{\mathrm{initial}}$ (meV)') axs[0].set_xlabel(r'$d_{\mathrm{N}-\mathrm{N}} (\AA)$') # Plot the atomic structure, focussing on the two Nitrogen atoms plot_atoms(config, axs[1], rotation='0x', show_unit_cell=0) axs[1].set_axis_off() axs[1].set_ylim(0.0, 5.642) axs[1].set_xlim(0.0, 14.833) def animate(i): scat.set_offsets((n2_distances[i], energy_differences[i])) # Remove the previous atomic plot [p.remove() for p in axs[1].patches] plot_atoms(configs[i], axs[1], rotation='0x', show_unit_cell=0) axs[1].set_xlim(0.0, 14.833) axs[1].set_ylim(0.0, 5.642) axs[1].set_axis_off() return (scat,) ani = animation.FuncAnimation( fig, animate, repeat=True, frames=len(configs) - 1, interval=200 ) .. video:: /examples_generated/tutorials/images/sphx_glr_neb_dissociation_001.mp4 :class: sphx-glr-single-img :height: 480 :width: 640 :autoplay: :loop: .. _sphx_glr_download_examples_generated_tutorials_neb_dissociation.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: neb_dissociation.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: neb_dissociation.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: neb_dissociation.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_