8.5. IRC analysis

In this notebook, we combine the both the forward and backward irc search into one trajectory and one list of energy values so that we can plot/analyze the minimum free energy path for converting the chloro-alcohol into the epoxide.

import pickle
import sys
sys.path.append("..")
from helpers import *

import psi4

import matplotlib.pyplot as plt 

8.5.1. Download the other part of the trajectory from your colleague

This load the pickled objects from disk

backward_e, traj_backward, traj_backward_np = pickle.load( open( "backward.p", "rb" ) )
forward_e, traj_forward, traj_forward_np= pickle.load( open( "forward.p", "rb" ) )

Look at the energy list and find out which you need to reverse.

Remember that we substracted the highest energy from all the values. Think about which point is the highest energy point on the minimum energy path connecting reactant and products.

You can use variablename.reverse() to reverse the list.

backward_e = list(backward_e)
forward_e = list(forward_e)

#list.reverse() #this is an inplace operation meaning you do not need to reassign the variable

# combine the energies 
# Note that we need exclude the ts point from one of the trajectories
# You can remove one index from the lists using list1[:-1] 
# Lists can be combined using lists = l1 + l2
energies = None 

# combine the trajectories which are lists
traj = None

8.5.2. Visualize the trajectory

drawXYZGeomSlider(traj)

Exercise 8

Is the stereochemistry at the carbon at which the reaction takes place retained?

8.5.3. Plot the minium energy path

fig, ax =  plt.subplots(1)

ax.plot(energies)
ax.set_xlabel('IRC step')
ax.set_ylabel('Pot Energy (kcal mol$^{-1}$)')

plt.show()

Exercise 7

Take a screenshot of the graph of the potential energy profile you recorded. Why is the barrier for the epoxide formation so low? Will this be the overall barrier for the reaction as depicted in the previous section?

8.5.4. Analyse Bonds

For the visualization we used a list of xyz file formated strings. For the analysis it is more convenient to work with a list of numpy arrays.

Let’s have a look at different bond distances involving different atoms in the molecule.

We can fist have a look at the first and last conformations, which correspond to the starting and ending point of the IRC that we performed.

drawXYZSideBySide_labeled(traj[0], traj[-1])

traj_array = traj_backward_np + traj_forward_np

bond_length_XX = []

for coord in traj_array:
    # collect bond length between atom 0 and 1 over the trajectory
    bond_length_XX.append(calculate_bond(coord[0], coord[1])) 

Exercise 9

How do the C-Cl and the two relevant C-O bond lengths change during the trajectory? Does the C-C bond in the ring contract as the epoxide is formed? Show a graph depicting the evolution of these parameters as the reaction progresses.

fig, ax = plt.subplots(1)

# add a call to ax.plot for each distance that you want to plot 
ax.plot(bond_length_XX)

# add an entry to the list of labels for each distance that you want to plot
plt.legend(["label"]) 

ax.set_xlabel('IRC step')
ax.set_ylabel('bond [A]')
plt.show()

Exercise 10

What is happening to the methyl group as the reaction proceeds? Find a suitable parameter (angle, dihedral) to describe and characterise possible changes you observe, change the code below. Explain in your report what atoms you considered and take and include the evolutions of the chosen parameters during the IRC procedure.

#Modify the code below, chosing proper atoms to monitor the changes of the methyl group 
methyl_angle = []
methyl_dihedral = []

for coord in traj_array:
    #modify here
    methyl_angle.append(calculate_angle(coord[#], coord[#], coord[#])) 
    methyl_dihedral.append(calculate_dihedral(coord[#], coord[#], coord[#], coord[#])) 

fig, ax = plt.subplots(1)
ax.plot(methyl_angle)
ax.plot(methyl_dihedral)

plt.legend(["angle (#-#-#)", "dihedral (#-#-#-#)"])
        
ax.set_xlabel('IRC step')
ax.set_ylabel('angle [deg]')
plt.show()