Skip to main content

How to Perform IRC Calculations

  • Updated

IRC stands for Intrinsic Reaction Coordinate. It is a method used to generate a minimum energy reaction path by starting from a transition state (TS), found through methods like NEB (Nudged Elastic Band), and optimizing the structure such that the energy decreases toward the reactant and product directions, respectively.

In this article, we will explain how to perform IRC calculations using Sella, a library developed by the Zádor group.

You can install Sella via pip:

pip install Sella

As of the time of updating this article (Jan 2026), the operation has been verified using following versions:

  • Python: 3.11

  • ASE: 3.27.0

  • Sella: 2.3.5

     

1. setup

import the required libraries.

import numpy as np
from IPython.display import Image
import matplotlib.pyplot as plt

import ase
from ase.constraints import FixAtoms
from ase import Atoms

from sella import IRC
from sella import Sella, Constraints

from pfcc_extras import view_ngl

from pfp_api_client import ASECalculator, Estimator
estimator = Estimator()
calculator = ASECalculator(estimator)

 

2. Loading the TS Structure

Please perform an NEB calculation or similar method to generate a TS structure. From those calculation results, the structure with the highest energy and decreasing force is your TS structure.

(Rerefence) Atomistic Simulation Tutorial 5-2, Reaction analysis on heterogeneous catalysts, 4. Check NEB calculation results, 

images = ase.io.read("NEB_images.xyz", index="::")
for i in images:
    i.calc = calculator
    energies = [image.get_total_energy() for image in images]

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(energies)
ax.set_xlabel("replica")
ax.set_ylabel("energy")
ax.grid(True)
ax.xaxis.set_ticks(np.array(range(0,20,2)))

fig.show()

 

def calc_max_force(atoms):
    return ((atoms.get_forces() ** 2).sum(axis=1).max()) ** 0.5

mforces = [calc_max_force(image) for image in images]
plt.plot(range(len(mforces)), mforces)
plt.xlabel("replica")
plt.ylabel("max force [eV]")
plt.xticks(np.arange(0, len(mforces), 2))
plt.grid(True)
plt.show()

 

In the case of the above figures, replica No. 8 would be the TS structure.

 

3. TS Optimization

The TS structure obtained from NEB does not undergo operations to strictly converge it to a saddle point. If you perform vibrational analysis or an IRC calculation directly from this TS structure, you may end up with an incorrect number of imaginary frequencies (it should be exactly one) or generate an abnormal pathway. Here, we will use the sella library to properly converge the TS structure.

TS = images[8].copy()

c = FixAtoms(indices=[atom.index for atom in TS if atom.position[2] <= 2])
TS.set_constraint(c)

TS.calc = calculator
print((((TS.get_forces()**2).sum(1))**0.5).max())

TSopt = Sella(TS)
TSopt.run(fmax=0.03)
print((((TS.get_forces()**2).sum(1))**0.5).max())

 

4. IRC Calculation

Here we execute the IRC calculation using the structure obtained after the TS optimization. This must be run in both directions of the reaction, so we execute direction='forward' and direction='reverse' separately.

TS_f = TS.copy()
TS_f.calc = calculator
TS_r = TS.copy()
TS_r.calc = calculator

irc_f = IRC(TS_f, trajectory='irc_f.traj', dx=0.1, eta=1e-4, gamma=0.4)
irc_f.run(fmax=0.01, steps=9999, direction='forward')
irc_r = IRC(TS_r, trajectory='irc_r.traj', dx=0.1, eta=1e-4, gamma=0.4)
irc_r.run(fmax=0.01, steps=9999, direction='reverse')

Then, concatenate the generated trajectories.

IRC_f_images = ase.io.read("irc_f.traj", index="::")
IRC_f_images.reverse()
IRC_r_images = ase.io.read("irc_r.traj", index="::")
IRC_images = IRC_f_images + IRC_r_images

You can also display the combined trajectory in a viewer.

v = view_ngl(IRC_images[87], representations=["ball+stick"], replace_structure=True)
v.view.control.spin([1,0,0], np.deg2rad(-90))
v

 

 

It is also possible to output it as a GIF or APNG video.

from pfcc_extras.visualize.povray import traj_to_apng
traj_to_apng(IRC_images[::4], rotation="0x,0y,0z", clean=True)

anim.png

 

Finally, let's review the history of the energy and forces.

energies = [image.get_total_energy() for image in IRC_images]

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(energies)
ax.set_xlabel("replica")
ax.set_ylabel("energy")
ax.grid(True)
fig.show()
def calc_max_force(atoms):
    return ((atoms.get_forces() ** 2).sum(axis=1).max()) ** 0.5

mforces = [calc_max_force(image) for image in IRC_images]
plt.plot(range(len(mforces)), mforces)
plt.xlabel("replica")
plt.ylabel("max force [eV]")
plt.xticks(np.arange(0, len(mforces), 20))
plt.grid(True)
plt.show()

 

 

5. Tips

If the TS optimization or the IRC calculation fails, please re-evaluate whether the initial TS structure generated by the NEB calculation is reasonable.

 

 

 

 

 

Related to

Related articles