# How to obtain the system force on a fixed atom

From: chengh_at_ringo.chem.pitt.edu
Date: Wed Dec 16 2009 - 08:40:37 CST

Dear NAMD developers and users,

Does anyone know how to get the system force on a fixed atom. I used
tclForces "loadtotalforces". It returned a value of zero. To solve this
problem, I used a very large homonic restraint (K=100000.0) to mimic
system force. Any better solutions?

By the way, I noticed that "loadtotalforces" or "loadforces" did not work
using NAMD 2.7b1. I have to use an older version NAMD 2.6 to make it work.
I attached the script I used. Can anybody fix this problem in NAMD 2.7b1?

Many thanks,

Mary H. Cheng
Department of Chemistry
University of Pittsburgh

##############################################################################
tclForces on
tclForcesScript {
set aid1 168200
set Tclfreq 10
set k 100000.0
set x0 3.706
set y0 1.398
set z0 -26.642

proc calcforces {} {
global aid1 Tclfreq
global k z0 x0 y0

set t [getstep]
set rx [lindex \$r(\$aid1) 0]
set ry [lindex \$r(\$aid1) 1]
set rz [lindex \$r(\$aid1) 2]
print "COORDs: \$t \$r(\$aid1)"

\$k*(\$rx-\$x0)*(\$rx-\$x0)/2.0+\$k*(\$ry-\$y0)*(\$ry-\$y0)/2.0+\$k*(\$rz-\$z0)*(\$rz-\$z0)/2.0]
set force "[expr -\$k*(\$rx-\$x0)] [expr -\$k*(\$ry-\$y0)] [expr
-\$k*(\$rz-\$z0)]"

if { [array exists forc] } {
set f1 \$forc(\$aid1)
set f1x [lindex \$f1 0]
set f1y [lindex \$f1 1]
set f1z [lindex \$f1 2]
print "ADDforce: \$t \$f1x \$f1y \$f1z"

}
if
{
[array
exists
fortot]
}
{
set ftot \$fortot(\$aid1)
set ftotx [lindex \$ftot 0]
set ftoty [lindex \$ftot 1]
set ftotz [lindex \$ftot 2]
print "TOTforce: \$t \$ftotx \$ftoty \$ftotz"
print "SYSforce:\$t [expr \$ftotx-\$f1x] [expr \$ftoty-\$f1y] [expr
\$ftotz-\$f1z]" }

}

}

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:53:36 CST