Script for following the Intrinsic Reaction Co-ordinate

Initiate the running of GRACE and apply the CHARMM QM/MM system set up procedure.

source $env(GRACE_LIB)/grace.tcl

set GLOBAL(scratch) $env(GRACE_SCRATCH)

init i

Bomb -1

CHARMM QM/MM SYSTEM SET UP

Instruct GRACE to wait until the CHARMM QM/MM set-up has been achieved.

Prnl 0

Wrnl 0

wait

The frelax procedure is defined in much the same was as for the TS refinement script. It relaxes the environment after moving the core atoms in the IRC calculation.

proc frelax xvec\ gvec {

upvar $xvec xparam $gvec grad

global X Y Z WMAI DX DY DZ PARA QMSTART

qmmm silent unmake_pos=xparam

Scal xref stor 1

Scal yref stor 2

Scal zref stor 3

Cons fix sele core end

Scal zref reca 3

Scal yref reca 2

Scal xref reca 1

Prnl 5

Mini abnr nstep 1000 tolg 0.005 step 0.0001 nprint 1

Prnl 0

Scal xref stor 1

Scal yref stor 2

Scal zref stor 3

Cons fix sele none end

Scal zref reca 3

Scal yref reca 2

Scal xref reca 1

set X(1)

# there appears to be an unresolved bug in some cases

# it may have been sorted, but to be sure, hit the

# X cache in ensure upload

set X(1)

set f [fcharmm xparam grad]

xyz>p p=xyz x=X y=Y z=Z

array>list array=xyz list=lxyz

work< name=lats2s1_ircb.xyz var=lxyz

Quick 2 16

set d1 $PARA(DIST)

Quick 2 326

set d2 $PARA(DIST)

puts "+++ Distance 2-16 = $d1"

puts "+++ Distance 2-328 = $d2"

flush stdout

return $f

}

As with the TS refinement script, the first task is to declare the atoms described by other QM codes. All atoms used in CHARMM are defined as MM atoms (fmm), which in this case is the entire system.

qmmm \

fmm=all

Select the atoms in the system to be included in the core and print information on the selected atoms to standard output.

Prnl 5

Defi core sele (SEGID neu5) end

Prnl 0

acquire_selection selection=core var=lcore

puts "Selected atoms and types"

foreach i $lcore {

puts "$IUPC($i) $i $RESI($i) $RESN($i)"

}

flush stdout

Now read in the refined TS structure file.

work> name=lats2s1_final.xyz var=lxyz

list>array array=xyz list=lxyz

p>xyz p=xyz x=X y=Y z=Z

Fix the core of the refined TS structure and optimise the environment. The amount of optimisation required is likely to be very small if the refined TS structure is taken from the .xyz file.

Scal xref stor 1

Scal yref stor 2

Scal zref stor 3

Cons fix sele core end

Scal zref reca 3

Scal yref reca 2

Scal xref reca 1

Prnl 5

Mini abnr nstep 1000 tolg 0.005 step 0.0001 nprint 1

Prnl 0

Scal xref stor 1

Scal yref stor 2

Scal zref stor 3

Cons fix sele none end

Scal zref reca 3

Scal yref reca 2

Scal xref reca 1

set X(1)

The qmmm command declares the positions of the core atoms; it is used to change their positions (make_pos) in the frelax routine.

qmmm \

mm=core \

make_pos=pos \

start

The hessian matrix for the refined TS is obtained. Up until this point, no such calculation has been done: the TS refinement calculation generates the hessian for the approximate TS and then updates the hessian after each refinement.

gethess \

func=fcharmm \

position=pos \

restart=$GLOBAL(scratch)/lats2s1_ircf_start_hess.dump \

hessian=h

For the IRC calculation, the hessian is mass weighted and the resulting eigenvalues are printed to the output file.

acquire_selection selection=core var=sele_list

catch {unset masses}

set c 0

foreach i $sele_list {

incr c

set masses($c) $MASS($i)

}

la mass mass=masses matrix=h result=mh

la eign matrix=mh

pmatrix matrix=values

# make velocity the first root

list>array list=roots(1) array=velocity

The reaction path is now calculated. Two IRC calculations are done: one leading to reactants, the other leading to the product well. The backward reaction is achieved by scaling the velocities by –1. The first line of the next part of the script achieves this and is not required for the forward IRC calculation. It should be noted that the hessian for the refined TS structure is the same in both IRC calculations and doesn’t need to be re-evaluated. The hessian can be read in by placing RESTART in place of dump_name below the gethess command. The gradient tolerance of 0.1 kcal mol^-1 Å^-1 (0.42 kJ mol^-1 Å^-1) is rarely achieved. This is because optimising the environment structure can occasionally result in an increase in the total energy of the system and a premature exit from the IRC calculation. This is why a further optimisation of both the IRC reactant and product structures is carried out (in a separate calculation) using CHARMM.

la mbys matrix=velocity scale=-1.0 result=velocity

reaction_path \

func=frelax \

nfunc=100000 \

position=pos \

velocity=velocity \

masses=masses \

print=2 \

accuracy=15 \

tstep=1.0 \

toler=0.1 \

dump_name=$GLOBAL(scratch)/lats2s1_ircb.dump \

half=0.0

The reactant or product structure is then saved to a .xyz file or a .pdb file.

xyz>p x=X y=Y z=Z p=pos

array>list array=pos list=lpos

work< name=lats2s1_ircb.xyz var=lpos

# save end structure in pdb format

Open unit 67 write form name lats2s1_ircb.pdb

Write coor pdb unit 67

* pdb file of backward irc path

* lactone break (react)

Close unit 67

Finally, the running of the script is finished.

exit