VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_cmm.chem.upenn.edu)
Date: Fri Sep 07 2007 - 17:27:26 CDT

On Fri, 7 Sep 2007, jo hanna wrote:

JH> Hi

hi jo,

before going into detail, i strongly recommend that you first
test and debug each fragment and VMD feature before writing
a whole script and being confused.

JH>
JH> I am trying to calculate the angle between specific atoms in my files (pdb
JH> and xtc) at each frame. The script I'm using is as follows:
JH>
JH> **
JH>
JH> set outfile [open angle.dat w]
JH> set frame0 [atomselect top "all within 3 of index 2" frame 0]

frame0 is not use below, so why define it here?

JH> set nf [molinfo top get numframes]

JH> for {set i 0} {$i<$nf} {incr i} {
JH> set A [atomselect top "index 2" frame $i]
JH> set B [atomselect top "index 0" frame $i]
JH> set C [atomselect top "index 1" frame $i]

the preceding code is one of the most common errors in
VMD scripts and a certain way to create a memory leak.

each atomselect command creates a new tcl procedure and
returns its (unique) name. since those names depend on
how many calls to atomselect were made before, you have
to store its name in a variable to have a generic script.

you can more easily handle this by putting the following
code fragment _outside_ the loop

set A [atomselect top "index 2"]
set B [atomselect top "index 0"]
set C [atomselect top "index 1"]

and then using within the loop

$A frame $i
$B frame $i
$C frame $i

to advance to the frame of interest.

JH> global M_PI
JH>
JH> # initialize arrays
JH> set a() 0; set b() 0; set c() 0; set ab() 0; set ac() 0;
JH>

to get the xyz triple for a give coordinate you have to use

[$A get {x y z}] which will return the coordinates. however,
since selections can span multiple atoms, you will receive a
list of lists as a result. so to get the x value you'd need to do.

set a(x) [lindex [$A get {x}] 0]

but it would actually be more easier to use a different strategy
(note this is completely off the top off my head and untested,
but i hope you can see where it is leading):

set atoms [atomselect top {index 0 1 2}]

and then

$atoms frame $i

and

set coords [$atoms get {x y z}]

and

set av [lindex $coords 2]
set bv [lindex $coords 0]
set cv [lindex $coords 1]

set ab [vecsub $av $bv]
set ac [vecsub $av $cv]

put $outfile "$i [expr acos([vecdot $ab $ac]/([veclength $ab] * [veclength $ac]) *(180.0/$M_PI)]"

this would make use of a lot of vmd internal functions
the whole programming less error prone (every line you
do not have to write cannot contain a bug).

cheers,
   axel.

JH> # split coordinates
JH> set a(x) [lindex $A 0]; set a(y) [lindex $A 1]; set a(z) [lindex $A 2];
JH> set b(x) [lindex $B 0]; set b(y) [lindex $B 1]; set b(z) [lindex $B 2];
JH> set c(x) [lindex $C 0]; set c(y) [lindex $C 1]; set c(z) [lindex $C 2];
JH>
JH> # setup vectors
JH> set ab(x) [expr $b(x) - $a(x)];
JH> set ab(y) [expr $b(y) - $a(y)];
JH> set ab(z) [expr $b(z) - $a(z)];
JH>
JH> set ac(x) [expr $c(x) - $a(x)];
JH> set ac(y) [expr $c(y) - $a(y)];
JH> set ac(z) [expr $c(z) - $a(z)];

JH>
JH> # compute cosine
JH> set cosine [expr \
JH> ($ab(x)*$ac(x) + $ab(y)*$ac(y) + $ab(z)*$ac(z)) / \
JH> (sqrt(pow($ab(x),2) + pow($ab(y),2) + pow($ab(z),2)) * \
JH> sqrt(pow($ac(x),2) + pow($ac(y),2) + pow($ac(z),2)))];
JH>
JH> puts $outfile "set angle [expr acos($cosine)*(180.0/$M_PI)]"
JH> }
JH> close $outfile
JH>
JH> **
JH>
JH> Using this gives me an error that I don't understand and therefore can't
JH> fix:
JH>
JH> syntax error in expression "atomselect10 - atomselect9": variable references
JH> resuire preceeding $
JH>
JH> Can someone please advise?
JH>
JH> Thanks
JH> Jo
JH> 

-- 
=======================================================================
Axel Kohlmeyer   akohlmey_at_cmm.chem.upenn.edu   http://www.cmm.upenn.edu
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.