8.3. Recording a Potential Energy Profile: The Intrinsic Reaction Coordinate (IRC)

A transition state search merely guarantees that one reaches a first-order saddle point on the PES; there is, however, no guarantee that it is really linked to the reactant and product. (For instance, there may be further transition states lying in between.) Although the visualisation of the imaginary vibrational mode may give some hint as to whether one has found a reasonable saddle point, this information is not sufficient to guarantee a connection to the reactant and product well. Following the negative Hessian eigenvalues from the transition state to the next local minimum, on the other hand, will immediately reveal whether such a connection exists. This Hessian-following in mass-weighted coordinates is referred to as a search along the Intrinsic Reaction Coordinate (IRC). This search will create a potential energy profile for the reaction, as well as a trajectory that visualises the lowest energy path from reactant to TS to product or reactant (depending on the direction that one chooses).

As psi4 automatically performs a frequency calculation we can verify whether the transition state that was computed is the proper transition state.

Warning

You do not need to run this notebook. The outputs forward.p and backward.p are provided to you in the next notebook. This is here for completeness. The calculation takes about 30 - 45 minutes depending on the direction you choose.

import psi4
import numpy as np

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

# Read the transition state coordinates from disk
with open('ts.xyz', 'r') as f:
    ts_string = f.read()  # this is the xyz file as a string
    print(ts_string)

ts = psi4.geometry(f"""
-1 1
symmetry c1
 C    0.169891741896   -0.303396542243   -0.246565081910
 C    0.191012418901   -0.532643890512    1.275916186757
CL    0.143804012854    1.566013819697   -0.686628413220
 H    1.120024392980   -0.636063203527   -0.672877892311
 C   -1.016755132582   -0.934703159138   -0.942141431016
 O    0.339535270596   -1.849728813088    1.472891387260
 H    1.010561961934    0.121411357267    1.686229845174
 H   -0.762913236261   -0.079611519921    1.679061412831
 H   -0.957983639113   -0.847359581313   -2.033892804890
 H   -1.031376458892   -1.990374831218   -0.649080196171
 H   -1.947496420408   -0.464044585762   -0.599503342745""")

8.4. Choose a direction

Start from the TS that we optimized and run the IRC calculation in one direction.

You can choose forward or backward for the irc_direction

irc_direction =  'forward' # 'backward'

psi4.core.set_output_file(f'irc.log', False)
psi4.set_options({
     "geom_maxiter":500,
     "full_hess_every":0,
     "opt_type":"irc",
     "irc_step_size":0.25,
     "ensure_bt_convergence":True,
     "irc_direction":irc_direction,
     })

E, history = psi4.optimize('b3pw91/6-31+G*',molecule=ts, return_history=True)

 55 displacements needed.
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
Warning: thermodynamics relations excluded imaginary frequencies: ['312.3714i']
Warning: used thermodynamics relations inappropriate for low-frequency modes: ['123.0283' '199.5642' '200.4680' '278.6674' '310.2934' '399.1543']
Optimizer: Optimization complete!

coordinates = history['coordinates']
energies = history['energy']


# Convert to kcal mol-1

# Convert to kcal mol-1
e = [energy * psi4.constants.hartree2kcalmol for energy in energies]
# set highest value to 0
e = np.array(e)-np.max(e)

import pickle
import json

pickle.dump((backward_e, 
             mol2traj(backward, coordinates_backward), 
             coord2traj_array(coordinates_backward)), open( f"{irc_direction}.p", "wb" ) )