Orbital Interaction Mechanism Along a Reaction Pathway#

Output video showing the orbital mechanism, overlap, and orbital energy gaps for the oxidative addition reaction of palladium into the C-H bond of methane. The dashed verticle line shows the C-H distance of the transition state.

A popular analysis method in quantum chemical modelling involves performing fragment analyses along a reaction coordinate to track how certain properties change during a reaction. This often involves some kind of partitioning of the total electronic energy, such as the activation strain model (ASM) and/or the energy decomposition analysis (EDA). For a comprehensive guide, see Vermeeren, P., et al. Nature Protoc. 2020.

../_images/OxAdd_rxn.png — The oxidative addition reaction of palladium into the methane C-H bond.#

ASM/EDA has been successfully applied to many different chemical systems. Subsequent KS-MO analysis has also been successfully performed, however, this is usually done at a single consistent geometry. With PyOrbb we can now perform KS-MO analysis along the full reaction coordinate. This example shows how this could be a usefull tool to aid the understanding of chemical reactions at the orbital interaction scale.

We will apply this idea to the oxidative addition of a palladium atom into the methane C-H bond. Palladium catalysts are often used in cross-coupling reactions and the oxidative addition is usually the rate-determining step as it involves breaking a relatively strong bond (C-H in this case). Because it is the rate-determining step this reaction is the target of numerous computational and experimental studies.

We have calculated the intrinsic reaction coordinate path (IRC) from the transition state (TS) of the oxidative addition reaction. We then used PyFrag to easily perform fragment analyses at each converged structure along the IRC path. We split the molecular system into the CH₄ and Pd fragments.

../_images/OxAddOrbint.png — The dominant orbital interaction pattern for the oxidative addition reaction.#

We know from past research (Diefenbach, A., et al. J. Phys. Chem. A 2004, Diefenbach, A., et al. J. Chem. Theory and Comp. 2005) that at the start of the oxidative addition reaction the stabilizing interaction between the empty Pd(5S) orbital and the filled CH₄(𝜎_C-H) dominates. As the reaction progresses towards the TS the orbital interaction between the filled Pd(4d_x2-y2) orbital and the empty (𝜎^*_C-H) orbital becomes more important. The video generated by the script shows the evolution of the orbital interaction mechanism of the oxidative addition reaction of Pd into CH₄, as well as the overlaps and orbital energy gaps of the Pd(5s)/CH₄(𝜎_C-H) and Pd(4d_x2-y2)/CH₄(𝜎^*_C-H) interactions. It shows that as the reaction proceeds the overlaps for both interactions increase significantly, with the interaction between Pd(5s)/CH₄(𝜎_C-H) having the highest overlap across all C-H bond distances, making it the dominant interaction at the start of the reaction. Note however, that the energy gap of the Pd(4d_x2-y2)/CH₄(𝜎^*_C-H) interaction quickly decreases as the reaction proceeds and approaches and crosses 0 eV. This means that the stabilizing orbital interaction for this interaction becomes more and more important as the reaction progresses. The generated video confirms the findings from the literature and is a quick and easy way to understand the orbital interactions given any Pyfrag analysis.

Download OxAdd_rkfs.zip

Download tracking.py

'''
This example shows how to generate a video of the evolution of the 
orbital interaction mechanism of the oxidative addition of a 
palladium atom into the methane C-H bond using PyOrbb.

This example uses ADF fragment calculations along the IRC path of 
the oxidative addition reaction. All data was computed at ZORA-BP86/TZ2P 
with Good numerical quality.
'''

import pyorbb
import os
import matplotlib.pyplot as plt
import numpy as np
import moviepy.editor as mvp


# obtain all adf.rkf files
files = [os.path.join('OxAdd_rkfs', file) for file in os.listdir('OxAdd_rkfs') if file.endswith('.adf.rkf')]

# then load the Orbitals objects
orbs = [pyorbb.Orbitals(file) for file in files]

# we want to also have the key distance, the C-H distance
# we will sort on this distance later on
molecules = [orb.molecule for orb in orbs]
# C and H are atoms 1 and 3 in the molecule
distances = [mol[1].distance_to(mol[3]) for mol in molecules]
# we then get the order in which to generate the frames
distance_order = np.argsort(distances)

# sort the orbs and distances by the order
orbs = [orbs[i] for i in distance_order]
distances = [distances[i] for i in distance_order]

# make the directory to store the frames in
os.makedirs('frames', exist_ok=True)

# for each Orbitals object generate the diagram and save it to a file
# we want to also show the energy gaps and overlaps
frames = []
overlaps1 = []
overlaps2 = []
egaps1 = []
egaps2 = []
for i, orb in enumerate(orbs):
    fig, axes = plt.subplots(1, 3, figsize=(6.4*2, 4.8))
    # the first interaction
    Pd_5s = orb.sfos['Pd(5S)']
    sub_lumo = orb.sfos['Substrate(HOMO)']
    mos = [orb.mos['21AA'], orb.mos['17AA']]
    mix1 = pyorbb.analysis.mixing.Mixing(
        orb,
        sfos=[Pd_5s, sub_lumo], 
        mos=mos,
        connection_colors='b')

    overlaps1.append(abs(Pd_5s @ sub_lumo))
    egaps1.append(abs(Pd_5s.energy - sub_lumo.energy))


    # the second interaction
    Pd_4d = orb.sfos['Pd(2D:x2-y2)']
    sub_homo = orb.sfos['Substrate(LUMO)']
    mos = [orb.mos['22AA'], orb.mos['18AA']]
    mix2 = pyorbb.analysis.mixing.Mixing(
        orb,
        sfos=[Pd_4d, sub_homo], 
        mos=mos,
        connection_colors='r')

    overlaps2.append(abs(Pd_4d @ sub_homo))
    egaps2.append(abs(Pd_4d.energy - sub_homo.energy))


    # add the two mixings together
    main_mix = mix1 + mix2

    # draw the diagrams
    main_mix.draw_diagram(ax=axes[0], ylim=(-13, 2))
    plt.suptitle(f'C-H distance: {distances[i]: .2f} Å')

    # and then plot the overlaps
    axes[1].plot(distances[:i+1], overlaps1, label=r'$\langle 5S | 4AA \rangle$', c='b')
    axes[1].plot(distances[:i+1], overlaps2, label=r'$\langle 2D:x2-y2 | 5AA \rangle$', c='r')

    # and energy gaps
    axes[2].plot(distances[:i+1], egaps1, label=r'$| \varepsilon_{5S} - \varepsilon_{4AA} |$', c='b')
    axes[2].plot(distances[:i+1], egaps2, label=r'$| \varepsilon_{2D:x2-y2} - \varepsilon_{5AA} |$', c='r')

    # set some limits and labels
    axes[1].set_xlim(min(distances), max(distances))
    axes[1].set_ylim(0, 0.65)
    axes[1].vlines(1.614, 0, 0.6, colors='k', alpha=0.5, linestyle='dashed')
    axes[1].set_ylabel('Overlap')
    axes[1].set_xlabel('C-H distance / Å')

    axes[2].set_xlim(min(distances), max(distances))
    axes[2].set_ylim(0, 5)
    axes[2].vlines(1.614, 0, 5, colors='k', alpha=0.5, linestyle='dashed')
    axes[2].set_ylabel('Orbital Energy Gap / eV')
    axes[2].set_xlabel('C-H distance / Å')

    axes[1].legend(frameon=False)
    axes[2].legend(frameon=False)

    # save the frame and store the file location for later use
    frame_file = f'frames/{i}.png'
    frames.append(frame_file)

    plt.savefig(frame_file)
    plt.close()


# take the file locations of the frames we saved and make a movie
clips = [mvp.ImageClip(frame).set_duration(1/14) for frame in frames]
concat_clip = mvp.concatenate_videoclips(clips, method="compose")
concat_clip.write_videofile('orbint.mp4', fps=14)
concat_clip.write_gif('orbint.gif', fps=14)