Re: How to obtain the system force on a fixed atom

From: Jérôme Hénin (jhenin_at_ifr88.cnrs-mrs.fr)
Date: Wed Dec 16 2009 - 10:38:55 CST

Dear Mary,
Does the fixedAtomsForces option make a difference in the behavior you see?
The stiff harmonic restraint is not a problem per se, as long as the
simulation remains stable, but the bad interaction with fixed atoms
could be a fixable bug. I don't remember whether the case was
considered when those TclForces functions were added.

As far as I know, you are the first to report the problem with
load(total)forces in NAMD 2.7b1. That's interesting, and should
certainly be addressed. Can you check that it still happens with
version 2.7b2, and describe the unexpected behavior precisely?

Best,
Jerome

2009/12/16 <chengh_at_ringo.chem.pitt.edu>:
> 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 - 15:53:37 CST