VMD-L Mailing List

From: John Stone (johns_at_ks.uiuc.edu)
Date: Fri Jan 25 2008 - 00:15:14 CST

Hi,
  If you're just trying to measure angles between three atoms, I'd
suggest using the built-in "measure angle" command in VMD 1.8.6 rather
than doing it with a script. The built-in command will be much faster:

measure angle {{<atomid1> ?<molid1>?} {<atomid2> ?<molid2>?} {<atomid3> ?<molid3>?}} ?molid <default molid>? [?frame <frame|all>? | ?first <first>? ?last <last>?] -- angle

Cheers,
  John Stone
  vmd_at_ks.uiuc.edu

On Thu, Jan 24, 2008 at 12:39:40PM -1000, Irene Newhouse wrote:
>
> I found a script for computing the angle between 2 vectors defined by atom centers using VMD on the AMBER mail reflector. I've been trying to modify it to compute the angle between 3 atoms. I keep getting this error:
>
> Info) Segments: 1Info) Fragments: 34 Protein: 30 Nucleic: 0dcdplugin) detected standard 32-bit DCD file of native endiannessdcdplugin) CHARMM format DCD file (also NAMD 2.1 and later)Info) Using plugin dcd for coordinates from file h5ltnw100.dcdInfo) Finished with coordinate file h5ltnw100.dcd.0380atomselect0atomselect1atomselect2can't read "1": no such variableInfo) VMD for LINUXAMD64, version 1.8.6 (April 6, 2007)Info) Exiting normally.
> I commented out all the lines that contain a "1" , but I still get that error. Here's the script as
> I last ran it:
>
> #modified from script found in AMBER archive
> # Monomer 1
> set outfile [open atheta.txt w]
> # load pdb file defining atoms & dcd file
> mol load pdb tanw.pdb dcd h5ltnw100.dcd
> set nf [molinfo top get numframes]
> set pt1 [atomselect top "resid 1449 and name C2"]
> set pt2 [atomselect top "resid 1448 and name C1"]
> set pt3 [atomselect top "resid 1447 and name C1"]
> #set conv [expr 180.0/acos(-1.0)]
> #calculate first the 2 vectors v1 & v2 & their lengths v1l & v2l
> #then the angle
> #get the coords for 2 points
> for {set i 0} {$i <$nf} {incr i} {
> $pt1 frame $i
> $pt2 frame $1
> set coor1 [measure center $pt1]
> set coor2 [measure center $pt2]
> #do the first vector
> set v1 [vecsub $coor1 $coor2]
> set v1l [veclength $v1]
>
> #get coords for next 2 points
> $pt2 frame $i
> $pt3 frame $i
> set coor3 [measure center $pt2]
> set coor4 [measure center $pt3]
> #do 2nd vector
> set v2 [vecsub $coor3 $coor4]
> set v2l [veclength $v2]
>
> #do dotproduct
> set dotp [vecdot $v1 $v2]
> #do angle
> set theta [expr $conv * acos($dotp/($v1l*$v2l))]
> #if {$theta > 90.0} {
> # set theta [expr 180.0 -$theta]
> #}
> #puts $outfile "[expr($i+1)*10] $theta"
> puts $outfile "$i $theta"
> }
> close $outfile
>
> Thanks for any help!
>
> Irene Newhouse
> _________________________________________________________________
> Connect and share in new ways with Windows Live.
> http://www.windowslive.com/share.html?ocid=TXT_TAGHM_Wave2_sharelife_012008 

-- 
NIH Resource for Macromolecular Modeling and Bioinformatics
Beckman Institute for Advanced Science and Technology
University of Illinois, 405 N. Mathews Ave, Urbana, IL 61801
Email: johns_at_ks.uiuc.edu                 Phone: 217-244-3349
  WWW: http://www.ks.uiuc.edu/~johns/      Fax: 217-244-6078