tcl force

From: Seungho Choe (sec53_at_pitt.edu)
Date: Fri Oct 16 2009 - 16:34:58 CDT

Dear all,

I'm trying to restrain the center of mass (only z position) of a
molecule. Below is my tcl script. What I'm wondering is whether a new
energy term from the added force is applied correctly to the total
energy. Do I need to explicitly add the energy term in the script ? Is
there any way to check this ? Any comments are welcome. Thank you in
advance.

Best,
Seungho

************************************************************************
set numatoms 63230

set atoms {}
for {set i 63207} { $i <= $numatoms} { incr i} {
    lappend atoms $i
}

print "atoms:$atoms"

foreach atom $atoms {
    addatom $atom
}

set groupid [addgroup $atoms]

# Set spring constant
set k 7.0
set z0 -3

proc calcforces {} {

global atoms groupid numatoms k z0
loadcoords p

#print "$p($groupid)"

set z1 [lindex $p($groupid) 2]

print "z comp:$z1"

set r01 [expr $z1-$z0]
#print "r01:$r01"

if {$r01 > 0} {
    set n01 {0 0 -1}
} elseif {$r01 < 0} {
    set n01 {0 0 +1}
}

# Calculate the "force"
set force [expr abs($k*$r01)]

set forcez [vecscale $force $n01]

print "forcez:$forcez"

addforce $groupid $forcez

}

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:51:35 CST