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

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

sorry, didn't read to the end of your post.

yes, i used the same approach (harmonic constraints) to get an average
force from on the virtual spring extension.

On Wed, Dec 16, 2009 at 10:34 AM, Roman Petrenko <rpetrenko_at_gmail.com> wrote:
> 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
>> "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]"        }
>>
>>
>>  }
>>
>>
>> }
>>
>>
>>
>>
>
>
>
> --
> Roman Petrenko
> Physics Department
> University of Cincinnati
>

-- 
Roman Petrenko
Physics Department
University of Cincinnati

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