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

From: Roman Petrenko (rpetrenko_at_gmail.com)
Date: Wed Dec 16 2009 - 09:34:36 CST

Dear Mary,
why should it be non-zero? what does newton's first law say?
what you mean is probably constraining an atom around a certain
position and then calculating the average force on the atom. for that
use harmonic constraints (see NAMD tutorials and user guide).

On Wed, Dec 16, 2009 at 9:40 AM, <chengh_at_ringo.chem.pitt.edu> wrote:
> 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]"        }
>
>
>  }
>
>
> }
>
>
>
>

```--
Roman Petrenko
Physics Department
University of Cincinnati
```

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