Transition-state searching

The script below is for the location of the transition-state for the bimolecular solvolysis reaction
of iso-propyl chloride. On the left handside are the GRACE commands with an explanation on the right handside.

A copy of the script can be found in the downloaded grace.tar file.

source $env(GRACE_LIB)/grace.tcl	Set location of GRACE files
set GLOBAL(scratch) ~/propyl	Set scratch disk location
init i	Initialise GRACE
PRNL 6 WRNL 6	Set CHARMM print and warning levels
!--------- BEGIN RTF ----------	Introduction of your CHARMM commands
Read rtf card * a simple rtf for propylcl + OH2 + TIP3P 24 1	Note that at the start of a line, uppercase refers to CHARMM commands and lowercase refer to GRACE commands
MASS 4 HT 1.000 ! H-bonding hydrogen MASS 12 CT 12.010 ! Carbon MASS 28 CL 35.000 ! Cl MASS 75 OT 15.999 ! TIP3P oxygen	Information on atom types used in simulation
! TIP3P RESI TIP3 0.000GROUP ATOM OH2 OT -0.83400 ATOM H1 HT 0.41700 ATOM H2 HT 0.41700 BOND OH2 H1 OH2 H2 H1 H2 ANGLE H1 OH2 H2 ACCEPTOR OH2 PATCHING FIRST NONE LAST NONE	TIP3P water information
RESI PRO 0.0 GROUP ATOM CR1 CT -0.1 ! HM1 HM2 HM3 ATOM XCL CL -0.2 ! \ | / ATOM CM1 CT -0.2 ! HO1 CM1 ATOM CM2 CT -0.2 ! \ | ATOM HD1 HT 0.1 ! OT1---- CR1 ----XCL ATOM HM1 HT 0.1 ! / | \ ATOM HM2 HT 0.1 ! HO2 HD1 CM2 ATOM HM3 HT 0.1 ! / | \ ATOM HM4 HT 0.1 ! HM4 HM5 HM6 ATOM HM5 HT 0.1 ! ATOM HM6 HT 0.1 ! ATOM OT1 OT -0.5 ! ATOM HO1 HT 0.25 ! ATOM HO2 HT 0.25 ! END Read para card * Simple parameters propylcl + H2O * BONDS HT OT 450.0 0.9572 HT HT 0.0 1.5139 ANGLES HT OT HT 55.0 104.52 DIHEDRALS	Propyl Chloride + nucleophilic water information such as: Name of residue and total charge on quantum system Atom labelling, atom type, specific atom charge and connectivity data for quick reference
NONBONDED nbxmod 5 atom cdiel shift vatom vdistance vswitch CUTNB 40.0 ctofnb 38.0 ctonnb 36.0 eps 1.0 e14fac 1.0 wmin 1.0 CT 0.000000 -0.080000 2.060000 ! All carbons HT 0.000000 -0.046000 0.224500 ! H bond H CL 0.000000 -0.126620 2.144580 ! Chlorine OT 0.000000 -0.152100 1.768200 ! TIP3P OXYGEN END	End of CHARMM structure and parameter data
Read sequence card * simple * 1 PRO Generate PRT setup	Read in sequence information The system contains only two segments PRT (QM) and SOLV (MM) made up from the PRO and TIP3P residues, respectively
Read sequence TIP3 494 Generate SOLV setup noangle nodihedral	Now read in the sequence for the 494 TIP3P water molecules
Open unit 10 read form name propyl_ts_inp.pdb Read coor pdb unit 10 Close unit 10	Now read in the co-ordinates of the MM & QM atoms from an external pdb file
Nbonds atoms switch cutnb 40.00 ctofnb 39.00 ctonnb 35.00 Updat inbfrq 1000 ihbfrq 0 cutnb 999.0 ctonnb 997.0 ctofnb 998.0 Cdie	Set non-bonding interactions and proximity warnings
Quan sele resi PRO end am1 charge 0 Prnl 0 Wrnl 0	Select residue that will be the QM system, set the charge (0) and computaional method (AM1)
wait	Wait for CHARMM to do it bit before we move on, or the output of grace and CHARMM get confused
proc frelax xvec\ gvec {	Set up procedure to perform relaxed environment functions. This procedure computes the f and g from fcharmm and then minimises the rest of the system around it, you must make the defined ‘core’ selection to represent the part of the system that is in the hessian, for this procedure to work
upvar $xvec xparam $gvec grad	Update the position and gradient vectors to xparam and 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 5000 tolg 0.005 step 0.0001 nprint 1000 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)	Always relax environment totally first set up core positions use the silent option so it does not print out lots of useless information. Start the initial environment optimisation using the ABNR optimiser.
set f [fcharmm xparam grad]	Get energy and gradient, the energy will be the total energy, but the gradient vector will only span the sub set defined by the QM/MM statement
Quick 1 12 set d1 $PARA(DIST) Quick 1 2 set d2 $PARA(DIST) puts "+++ Distance 1-12 = $d1" puts "+++ Distance 1-2 = $d2" flush stdout return $f }	Ask GRACE to instruct CHARMM to print out the bond lengths between bonds of interest
qmmm \ fmm=all	Define the parts of the QM/MM system. Grace does not ‘know’ if you have parts of the CHARMM system AM1 or not, so first make everything fixed MM
Prnl 5 Defi core sele resi PRO 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	Here we make the ‘core’ definition and set Prnl high so we can see if it has worked
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 5000 tolg 0.001 step 0.0001 nprint 20 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)	Now make the core bit movable MM and then generate a position vector for it, so GRACE 'internally' knows the parts that are movable. Carry out an energy minimisations using the ABNR optimiser
qmmm \ mm=core \ make_pos=pos \ start
gethess \ func=fcharmm \ position=pos \ dump_name=$GLOBAL(scratch)/hess.dump \ hessian=h
la symm matrix=h result=h la eign matrix=h	Diagonalise martix
pmatrix matrix=values	Print out roots values
optimise \ position=pos \ func=frelax \ nfunc=100000 \ omin=0.800 \ osmin=0.001 \ rmin=0.90 \ rmax=1.10 \ step=0.001 \ maxstep=0.01 \ mode=1 \ lock \ dump_name=$GLOBAL(scratch)/ts.dump \ method=ef \ hessian=h \ toler=0.3	Perform an EF search for a transition-state Follow mode 1, hopefully this is the right mode other wise we can stop the simulation, change the mode to the correct one and then restart it. ‘dump’ can be used to restart the simulation or used to view the mode that is being followed.
xyz>p x=X y=Y z=Z p=final_pos array>list array=final_pos list=lp array>list array=h list=lh la eign matrix=h array>list array=roots list=lr
Open unit 67 write form name propyl_ts.pdb Write coor pdb unit 67 * pdb of the transition-state for propylcl + H2O * Close unit 67 exit	Write out a pdb file of the final transition-state structure All done.