HCN → HNC: Transition State and Barrier¶
Hydrogen cyanide can isomerize to hydrogen isocyanide by moving the hydrogen from carbon to nitrogen. This example locates the transition state (TS) of that reaction, proves that it is a TS with a frequency calculation, and computes the barrier and the reaction energy from the two minima.
Step 1: find the transition state¶
Start from a bent guess with the hydrogen roughly between the two heavy atoms:
dft/pbe0/6-31g* ts(S0)
geom="""
C 0.000000 0.000000 0.000000
N 0.000000 0.000000 1.190000
H 1.100000 0.000000 0.500000
"""
openqp hcn-ts.oqp
ts(S0) searches for a first-order saddle point on the ground-state surface:
the optimizer climbs along one direction and minimizes along all the others.
The run takes about 4 seconds and converges in 9 steps:
Step Energy Shift RMSD Step Max Step RMSD Grad Max Grad
----------------------------------------------------------------------------------------------------------------------
1 -93.22942672 -93.22942672 0.000000 0.000000 0.006977 0.012290
...
8 -93.22809165 -0.00000149 0.000781 0.001613 0.000152 0.000291
9 -93.22809179 -0.00000014 0.000378 0.000829 0.000044 0.000085
| Quantity | Reference value |
|---|---|
| TS energy | −93.2280918 hartree |
| r(C–N) | 1.189 Å |
| r(C–H) | 1.182 Å |
| r(N–H) | 1.394 Å |
A TS search is more sensitive to the starting structure than a minimization. A guess far from the TS can take many more steps; it helps to place the moving atom near its expected position.
Step 2: confirm it is a transition state¶
A true TS has exactly one imaginary frequency. Run an analytic Hessian at the
converged geometry from opt.xyz:
dft/pbe0/6-31g* hess(S0,type=analytical)
geom="""
C 0.025858 0.000000 0.022436
N -0.050362 0.000000 1.208979
H 1.124504 0.000000 0.458585
"""
Mode Symmetry Frequency(cm-1) IR(km/mol) Raman(bohr^4/amu)
1 a' -1158.14 78.483570 59.221228
2 a' 2098.48 7.946939 242.971350
3 a' 2676.92 121.212840 491.302998
OpenQP prints imaginary frequencies as negative numbers. Mode 1, −1158 cm⁻¹,
is the only one, and its normal mode (shown above) moves the hydrogen from
carbon toward nitrogen. This is the reaction coordinate. (This table is from GitHub main
commit 0b75a06. OpenQP 1.3.1 prints the same frequencies with an IR(km/mol)
column 23 times too small; see the
frequency example.)
Step 3: the two minima¶
dft/pbe0/6-31g* opt
geom="""
C 0.000000 0.000000 0.000000
N 0.000000 0.000000 1.160000
H 0.000000 0.000000 -1.070000
"""
dft/pbe0/6-31g* opt
geom="""
C 0.000000 0.000000 0.000000
N 0.000000 0.000000 1.170000
H 0.000000 0.000000 2.170000
"""
Download hcn-min.oqp · Download hnc-min.oqp
| Structure | Energy (hartree) | Geometry |
|---|---|---|
| HCN | −93.3072282 | r(C–N) 1.155 Å, r(C–H) 1.071 Å |
| TS | −93.2280918 | see above |
| HNC | −93.2812843 | r(C–N) 1.175 Å, r(N–H) 0.999 Å |
With 1 hartree = 627.5095 kcal/mol:
- Barrier HCN → TS: 0.0791364 hartree = 49.7 kcal/mol
- Reaction energy HCN → HNC: 0.0259439 hartree = 16.3 kcal/mol (HNC is less stable)
These are electronic energies without zero-point corrections. To include them,
run hess(S0,type=analytical) at each of the three geometries and use the
zero-point energies from the thermochemistry block, as in the
frequency example.
Step 4: follow the reaction path (IRC)¶
The intrinsic reaction coordinate (IRC) is the steepest-descent path from the TS down to each minimum. It proves that this TS connects HCN with HNC. Run it in both directions from the TS geometry:
dft/pbe0/6-31g* irc(S0,hessian=analytical,maxit=120)
geom="""
C 0.025858 0.000000 0.022436
N -0.050362 0.000000 1.208979
H 1.124504 0.000000 0.458585
"""
Download hcn-irc-forward.oqp ·
hcn-irc-backward.oqp (first line
dft/pbe0/6-31g* irc(S0,direction=backward,hessian=analytical,maxit=120))
hessian=analytical starts the path along the imaginary mode of an analytic
Hessian, and maxit=120 allows up to 120 points (the default of 30 stops
before the minimum). Each direction takes about 3.5 minutes. Both runs end
with a line such as:
PyOQP: IRC Bracketed A Minimum Basin (49 points)
The forward path ends at −93.2811338 hartree after 49 points, next to the HNC
minimum (−93.2812843), and the backward path at −93.3072173 after 37 points,
next to the HCN minimum (−93.3072282). The log reports the step length as
step=0.100.
Try next¶
- Compare the barrier with a larger basis set, for example
def2-tzvp, at geometries reoptimized with that basis.
See also: Transition-State Searches, Optimization.