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
"fixed". Then using "loadforces" and "loadtotalforces" together to get the
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
        addatom $aid1
        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]
        loadcoords r
        set rx [lindex $r($aid1) 0]
        set ry [lindex $r($aid1) 1]
        set rz [lindex $r($aid1) 2]
        print "COORDs: $t $r($aid1)"

        addenergy [expr
$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)]"
        addforce $aid1 $force

        loadtotalforces fortot
        loadforces forc

        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 - 05:22:36 CST