Script for G94/98 Hessian Generation

# tclsh gethess.grc

# load up grace

source $env(GRACE_LIB)/grace.tcl

# start grace without charmm

init i n

# open ts.txt as punch file

set file [open bu_sn1_rc.FChk]

# parse out bumph

while 1 {

gets $file in

puts $in

case $in in Nuclear* break

}

# acquire nuclear charges for atom types

while 1 {

gets $file in

puts $in

case $in in Current* break

append nuch "$in "

}

# get cartesian coordinates

while 1 {

gets $file in

puts $in

case $in in Cart* break

append lpos "$in "

}

# parse out gradients

while 1 {

gets $file in

puts $in

case $in in Cart* break

# get hessian

}

while 1 {

gets $file in

case $in in "" break

append hess "$in "

}

puts "NATOM = [llength $nuch] NHESS = [llength $hess]"

# Map for nuclear charge into atom type

set lmap " H he li be b c n o f ne na mg al p si s cl ar br"

# Conver input lists into tcl arrays

list>array list=lmap array=map

list>array list=nuch array=anuch

list>array list=hess array=h

list>array list=lpos array=pos

# scale carts into angstrom

la mbys matrix=pos scale=0.529177249 result=pos

# scale hessian into ?

!la mbys matrix=h scale=0.529177249 result=h

# No nuclear charge to type mapping

set i 0

foreach c $nuch {

incr i

# use lots of 'int' to get indexes integer

set atsy($i) $map([expr int($anuch($i))])

}

parray atsy

# make list of carst into x y and z

p>xyz p=pos x=x y=y z=z

# format hessian into camvib style hessian

set c 0

set c1 0

set list ""

foreach i [array names h] {

incr c

incr c1

append list [format %12.8f $h($c)]

if {$c1 == 6} {append list \n ; set c1 0}

}

# put camvib style heesian onto disk

work< var=list name=bu_sn1_rc.hess

# compute least distance connection patter

ic connect distance x=x y=y z=z var=connect

# compute valence coordiantes

ic valence x=x y=y z=z connect=$connect l_func non_redu torsion2 var=valence

# make camvib deck

ic camvib name=bu_sn1_rc.vib x=x y=y z=z atsy=atsy pre= "VIB

VALE SYMM IFCS VECT DISO

G94 punch deck " hessian=bu_sn1_rc.hess valence=$valence