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

From: chengh_at_ringo.chem.pitt.edu
Date: Mon Dec 21 2009 - 15:18:47 CST

Dear Jerome and other NAMD developers,

I have tested the NAMD 2.7b1 and NAMD 2.7b2 istalled on psc supercomputing
center. "loadforces" "loadtotalforces", "loadforces" and "loadcoords" do
not work at all. The error message is something like "FATAL ERROR: can't
read "r(5941)": no such element in array".

I strongly suggest to implement a way to get the system force on a fixed
atom (not the "zero" value, but the force calculated due to molecular
interactions). This is extremely usefull for calculating the ion diffusion
constant inside an ion channel. Using harmonic restraint may introduce
artificial oscillation behariors for the random force autocorrelaton
function calculations.

Many thanks,

Mary

> 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:38 CST