[molpro-user] optg/irc problem

[Jussi Eloranta]

Hi,

The following example is given in the molpro manual:
---

!examples/hcn_isomerization.com $Revision: 2009.1 $
  ***, HCN<--->  NHC Isomerization Reaction Path
  memory,1,m
  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
  optg,root=2,saveact=hcn_ts                   ! Find 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'}
---

This example works correctly only because it is given the correct transition state geometry as input. This
fist led me to think that I can combine both TS search and IRC calculation in one file. Of course,
in this case the first optg would take more than one iteration and then readvar,hcn_ts.act would
result in the initial (not the TS) geometry to be restored (note that variable refs without array indices
will pick the first element; the optimal TS is at the last array idx). So IMHO this example is a bit misleading
(even though the comment line in it states that it is the TS geometry)... On the other hand, I don't know what would be
a geometry/variable independent way of resting the TS in the input file?

Jussi Eloranta
Cal State Northridge