***, HCN <---> NHC Isomerization Reaction Path
 basis=3-21G

 rcn=1.18282 ang      ! Starting geometry is transition state
 rnh=1.40745 ang
 alpha=55.05 degree

 symmetry,x           ! Cs Symmetry
 geometry={
           C
           N,1,rcn
           H,2,rnh,1,alpha}

 int
 rhf,gap_min=0     !set default shift to zero to avoid convergence to wrong state
 optg,root=2,saveact=hcn_ts,rewind            ! Find and store the TS
 {optg,method=qsdpath,dir=1, numhess=5,hesscentral,saveact=hcn_path}
                                              ! find IRC in positive direction

 readvar,hcn_ts.act                           ! Reset geometry to TS
 {optg,method=qsdpath,dir=-1,numhess=5,hesscentral,saveact=hcn_path,append}
                                              !find IRC in negative direction

 readvar,hcn_path.act

 alpha=alpha*pi/180    !convert angle to radian

 table,irc,rcn,rnh,alpha,e_opt   !tabulate results

 {table,irc,e_opt                !plot energy profile as function of irc
  plot,file='hcn_eopt.plot'}

 {table,irc,rcn,rnh,alpha        !plot distances and angle as function of irc
  plot,file='hcn_dist.plot'}