Re: Simple TCL (Forces) question - defining bond vector

From: Victor Ovchinnikov (ovchinnv_at_MIT.EDU)
Date: Thu Dec 06 2007 - 12:53:37 CST

Jim,

try this

set bond2 [vecsub $p($aid1) $p($aid2)]

# to do this by hand,
set bond2 {}
foreach {a b c} $p($aid1) {d e f } $p($aid2) {
 lappend bond2 [expr {$a-$d}] [expr {$b-$e}] [expr {$c-$f}]
}

note that the forces on atoms 1 and 2 should be opposite in direction if
your potential is k(r1-r2)^2
Victor

Jim Pfaendtner wrote:
> Dear NAMD-L,
>
> I am trying easily to define a vector between two atoms inside of a
> TCL forces script. Taking a cue from the tutorial I have written this:
>
> .. intro stuff defining aid1 and aid2 ...
>
> addatom $aid1
> addatom $aid2
> proc calcforces {} {
>
> global aid1 aid2 k
> loadcoords p
> set bond [getbond $p($aid1) $p($aid2)]
>
> ###
> set bond2 [please help me here :) ]
> ###
> addenergy [expr $k*$bond*$bond]
> set force [vecscale [expr -$k*2] $bond2]
> addforce $aid1 $force
> addforce $aid2 $force
> }
> }
>
> I know that bond2 needs to be a vector that is simply equal to:
> $p($aid1) - $p($aid2) but I can't get the syntax right and I struck
> out searching.
>
> Could someone please suggest the correct syntax for defining a new
> vector as the difference of two other vectors?
>
> Thanks very much,
> Jim
>
>

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:45:40 CST