Tracking Molecular Orbitals Along a Reaction Pathway#

It can often be insightful to investigate the evolution of orbitals along the potential energy surface of a chemical system. However, when investigating systems with large geometric deformations, there may be cases where the ordering of the molecular orbitals changes from one step to the next. In this example we investigate the MO energies along the stretching of the Cr-Cr bond of the chromium(I) hydride dimer.

../_images/CrH_dimer.png — The chromium(I) hydride dimer we investigate in this example. The Cr–Cr bond distance is varied between 1.2 and 3.0 angstrom over 30 steps.#

Without tracking#

To illustrate the problem we will obtain the MO energies of the 12Σ, 13Σ, 14Σ, and 15Σ MOs along the stretching coordinate. Selecting means that we simply take the energies of the MOs named “12Σ”, “13Σ”, “14Σ”, and “15Σ” from the adf.rkf files. In the figure below we see that the MOs 12Σ (blue line) and 13Σ (orange line) switch ordering around 1.6 Å when the energy of 13Σ dips below that of 12Σ. Similarly, around 2.8 Å we see that 15Σ should become more stable than 14Σ.

../_images/untracked.png — Energies of the 12Σ, 13Σ, 14Σ, and 15Σ orbitals along the bond stretching coordinate of the chromium(I) hydride dimer (HCr–CrH).#

Orbital tracking#

To obtain the correct energies of the MOs in question we rely on calculating and comparing the overlap populations of the MOs and FMOs between each step. We compute for each step (\(T\)) an order-3 tensor containing the overlap population of each MO \(\Psi_i^T\) and to each pair of FMOs \(\psi_j^T\) and \(\psi_k^T\) of step \(T\):

\[P_{ijk}^T = \langle \psi_j^T | \Psi_i^T \rangle \langle \psi_k^T | \Psi_i^T \rangle \langle \psi_j^T | \psi_k^T \rangle,\]

where \(\langle \psi_j^T | \Psi_i^T \rangle\) is the mixing coefficient of \(\psi_j^T\) into \(\Psi_i^T\), and \(\langle \psi_j^T | \psi_k^T \rangle\) is the overlap integral between \(\psi_j^T\) and \(\psi_k^T\). Starting from the matrix slice corresponding to the inital MO (e.g., \(P_{12Σ,jk}^1\) for MO 12Σ of the first step) we calculate the difference between the initial slice with all slices in the order-3 tensor of the next step. The MO corresponding to the slice with the smallest difference should be the most similar MO. We then continue this process for all given steps to “track” the MOs along the bond stretching coordinate.

\[\mathrm{argmin}_i \sum_j\sum_k |P_{\mathrm{initial}, jk}^T - P_{ijk}^{T+1}|,\]

where \(P_{\mathrm{initial}, jk}^T\) is the overlap population matrix for the initial MO at step \(T\).

PyOrbb includes an easy-to-use implementation for this algorithm (pyorbb.analysis.sequential.MOTracker). Using this implementation we obtain the correctly tracked MOs for our systems (see Figure below).

../_images/tracked.png — Energies of the 12Σ, 13Σ, 14Σ, and 15Σ orbitals tracked along the bond stretching coordinate of the chromium(I) hydride dimer (HCr–CrH).#

After tracking the MOs correctly we see that the MOs that previously switched order now cross each other. E.g., 12Σ now correctly keeps increasing in energy and becomes higher in energy than 13Σ around 1.6 Å. The same happens for 14Σ and 15Σ around 2.8 Å. We can also generate movies of the evolution of the orbitals along the bond stretching. We see that the MOs stay consistent along the complete coordinate. Conversely, the movies showing the untracked orbitals show clearly that the orbitals are not consistent across the whole coordinate, with large changes in their shapes.

Visualization#

12Σ	13Σ
14Σ	15Σ

Script and Resources#

The rkf files have been split into two separate zip files. Make sure to extract them into the same folder.

Download rkfs.1.zip

Download rkfs.2.zip

Download sequential_example.py

import matplotlib.pyplot as plt
import pyorbb
import os
import tcmu
from math import pi

# the systems described are fragment analyses of two HCr fragments
# each rkf file describes the system at different bond distances 
# between the Cr---Cr atoms

# the folder is structured like this:
# ./rkfs
# |-- step1.rkf
# |-- step2.rkf
# |-- ...
# |-- step30.rkf

# Obtain all the adf.rkf files we want to analyse
# sort them by the number in their name
rkfs_folder = './rkfs'
rkf_file_names = [f for f in os.listdir(rkfs_folder) if os.path.isfile(os.path.join(rkfs_folder, f)) and f != '.DS_Store']
rkf_file_names = sorted(rkf_file_names, key=lambda f: int(f.removeprefix('step').removesuffix('.rkf')))

# construct the full paths of the sorted rkf file names
adf_rkf_files = [os.path.join(rkfs_folder, f) for f in rkf_file_names]

# construct an MOTracker object with the obtained adf.rkf files
T = pyorbb.analysis.MOTracker(rkf_files=adf_rkf_files)

# energies of the 12SIGMA and 13SIGMA MOs tracked along the steps
tracked_E_12SIGMA = T.energy('12SIGMA')
tracked_E_13SIGMA = T.energy('13SIGMA')
tracked_E_14SIGMA = T.energy('14SIGMA')
tracked_E_15SIGMA = T.energy('15SIGMA')

# against the bond distance between the two Cr atoms (atoms 1 and 3)
coordinate = T.geometry(1, 3)

plt.figure(figsize=(5,3))
plt.title('Tracked MO energies')
plt.plot(coordinate, tracked_E_12SIGMA, label=r'12$\Sigma$')
plt.plot(coordinate, tracked_E_13SIGMA, label=r'13$\Sigma$')
plt.plot(coordinate, tracked_E_14SIGMA, label=r'14$\Sigma$')
plt.plot(coordinate, tracked_E_15SIGMA, label=r'15$\Sigma$')
plt.xlabel('Cr---Cr distance / Å')
plt.ylabel('MO energy / eV')
plt.legend(frameon=False)
plt.gca().spines[['right', 'top']].set_visible(False)
plt.tight_layout()
plt.savefig('./tracked.png', dpi=500)
plt.close()

# energies of the 12SIGMA and 13SIGMA MOs untracked
untracked_E_12SIGMA = [orbs.mos['12SIGMA'].energy for orbs in T.orbital_objects]
untracked_E_13SIGMA = [orbs.mos['13SIGMA'].energy for orbs in T.orbital_objects]
untracked_E_14SIGMA = [orbs.mos['14SIGMA'].energy for orbs in T.orbital_objects]
untracked_E_15SIGMA = [orbs.mos['15SIGMA'].energy for orbs in T.orbital_objects]
plt.figure(figsize=(5,3))
plt.title('Untracked MO energies')
plt.plot(coordinate, untracked_E_12SIGMA, label=r'12$\Sigma$')
plt.plot(coordinate, untracked_E_13SIGMA, label=r'13$\Sigma$')
plt.plot(coordinate, untracked_E_14SIGMA, label=r'14$\Sigma$')
plt.plot(coordinate, untracked_E_15SIGMA, label=r'15$\Sigma$')
plt.xlabel('Cr---Cr distance / Å')
plt.ylabel('MO energy / eV')
plt.legend(frameon=False)
plt.gca().spines[['right', 'top']].set_visible(False)
plt.tight_layout()
plt.savefig('./untracked.png', dpi=500)
plt.close()


# we generate videos of the tracked MOs to show that they indeed are correctly tracked
transforms = tcmu.geometry.Transform()
transforms.rotate(y=90/180*pi)
pyorbb.analysis.movie.make_orbital_movie('tracked_12SIGMA.mp4', T.track_mo('12SIGMA'), transform=transforms, fps=15)
pyorbb.analysis.movie.make_orbital_movie('tracked_13SIGMA.mp4', T.track_mo('13SIGMA'), transform=transforms, fps=15)
pyorbb.analysis.movie.make_orbital_movie('tracked_14SIGMA.mp4', T.track_mo('14SIGMA'), transform=transforms, fps=15)
pyorbb.analysis.movie.make_orbital_movie('tracked_15SIGMA.mp4', T.track_mo('15SIGMA'), transform=transforms, fps=15)

pyorbb.analysis.movie.make_orbital_movie('untracked_12SIGMA.mp4', [orbs.mos['12SIGMA'] for orbs in T.orbital_objects], transform=transforms, fps=15)
pyorbb.analysis.movie.make_orbital_movie('untracked_13SIGMA.mp4', [orbs.mos['13SIGMA'] for orbs in T.orbital_objects], transform=transforms, fps=15)
pyorbb.analysis.movie.make_orbital_movie('untracked_14SIGMA.mp4', [orbs.mos['14SIGMA'] for orbs in T.orbital_objects], transform=transforms, fps=15)
pyorbb.analysis.movie.make_orbital_movie('untracked_15SIGMA.mp4', [orbs.mos['15SIGMA'] for orbs in T.orbital_objects], transform=transforms, fps=15)