VMD-L Mailing List

From: Lee, Hwan Kyu (NIH/NHLBI) [E] (leeh3_at_nhlbi.nih.gov)
Date: Wed Feb 25 2009 - 13:35:07 CST

I just posted my message, saying that I'm still have a trouble. I just realized that my script repeatedly calculates only first frame, not all frames. I loaded XXX.psf file, and then loaded CHARMM trj files (XXX.trj). Then, I executed my script, but it only calculates the first frame. My script is shown below, again. Could you tell me how to fix this problem? Thank you very much for your help.

best,
Hwankyu.
--------------------------
proc calc {atom1 atom2 atom3 atom4} {
   set 12 [vecsub $atom2 $atom1]
   set 21 [vecsub $atom1 $atom2]
   set 23 [vecsub $atom3 $atom2]
   set 32 [vecsub $atom2 $atom3]
   set 34 [vecsub $atom4 $atom3]
   set n1 [veccross $12 $23]
   set n2 [veccross $23 $34]
   set cosangle [expr [vecdot $n1 $n2] / ([veclength $n1] * [veclength $n2])]
   set angle [expr {180 * acos($cosangle) / 3.14159265358}]
   set sign [vecdot [veccross $n1 $n2] $23]
   if { $sign < 0 } {
     set angle [expr {0-$angle}]
   } elseif { $sign == 0 } {
     set angle 180.0
   }
   return $angle
}

set numFrames [molinfo top get numframes]

set m 1

for {set pp 1000} {$pp < $numFrames} {incr pp} {

        set n 1

        while {$n < 35} {

        set n1 [expr {$n+1}]
        set n2 [expr {$n+2}]
        set n3 [expr {$n+3}]

        set sel1 [atomselect top "name C1 O1 C2 and resid $n"]
        set sel2 [atomselect top "name C1 O1 C2 and resid $n1"]
        set sel3 [atomselect top "name C1 O1 C2 and resid $n2"]
        set sel4 [atomselect top "name C1 O1 C2 and resid $n3"]

        set sel11 [measure center $sel1 weight mass]
        set sel21 [measure center $sel2 weight mass]
        set sel31 [measure center $sel3 weight mass]
        set sel41 [measure center $sel4 weight mass]

        set angle [calc $sel11 $sel21 $sel31 $sel41]

        $sel1 delete
        $sel2 delete
        $sel3 delete
        $sel4 delete

        incr n
        set result [open "dihedral.xvg" a]
        puts $result "$m\t$angle"
        close $result
        incr m
        }
}
--------------------------------------

-----Original Message-----
From: J T [mailto:jtibbitt_at_odu.edu]
Sent: Tue 2/24/2009 10:37 PM
To: Lee, Hwan Kyu (NIH/NHLBI) [E]
Cc: vmd-l_at_ks.uiuc.edu
Subject: Re: vmd-l: question about "measure dihed" with COM points
 
VMDs 'measure dihed' command takes four atom indices as input. You
are trying to give 4 xyz coordinate sets as input instead of atom
indices. I've already written a small procedure script for
calculating the dihedral angle given 4 sets of coordinates. In your
script, you would use it like:

set angle [calc $sel11 $sel21 $sel31 $sel41]

# CALCULATE DIHEDRAL ANGLES
# n1 and n2 = normal vectors formed by atoms 1 2 3 and 2 3 4,
respectively
# n1 = (x1 y1 z1)
# n2 = (x2 y2 z2)
# angle = acos (n1 . n2)/(|n1|*|n2|)
#
# Given 3 points: (x1 y1 z1) (x2 y2 z2) (x3 y3 z3)
# Two vectors are: (x2-x1 y2-y1 z2-z1) (x3-x2 y3-y2 z3-z2)
# Normal Vector = | i j k |
# | x1-x2 y1-y2 z1-z2 |
# | x2-x3 y2-y3 z2-z3 |
#
proc calc {atom1 atom2 atom3 atom4} {
   set 12 [vecsub $atom2 $atom1]
   set 21 [vecsub $atom1 $atom2]
   set 23 [vecsub $atom3 $atom2]
   set 32 [vecsub $atom2 $atom3]
   set 34 [vecsub $atom4 $atom3]
   set n1 [veccross $12 $23]
   set n2 [veccross $23 $34]
   set cosangle [expr [vecdot $n1 $n2] / ([veclength $n1] *
[veclength $n2])]
   set angle [expr {180 * acos($cosangle) / 3.14159265358}]
   set sign [vecdot [veccross $n1 $n2] $23]
   if { $sign < 0 } {
     set angle [expr {0-$angle}]
   } elseif { $sign == 0 } {
     set angle 180.0
   }
   return $angle
}

I'm pretty sure it works well. Hope you can use it.
Jeff Tibbitt

On Feb 24, 2009, at 8:09 PM, Lee, Hwan Kyu (NIH/NHLBI) [E] wrote:

> Dear VMDer,
>
> I calculated center of mass for each residue, and then calculated
> dihedral angles for those four center of mass. I wrote tcl script
> like below, but got an error, "expected integer but got
> "11.965297699" measure dihed: bad atom index". It looks like
> $sel11 cannot be recognized as a coordinate of COM. How can I fix
> this problem?
> Thanks for your help in advance.
>
> -----------------
> set numFrames [molinfo top get numframes]
>
> for {set pp 300} {$pp < $numFrames} {incr pp} {
>
> set n 1
>
> while {$n < 35} {
>
> set n1 [expr {$n+1}]
> set n2 [expr {$n+2}]
> set n3 [expr {$n+3}]
>
> set sel1 [atomselect top "resid $n"]
> set sel2 [atomselect top "resid $n1"]
> set sel3 [atomselect top "resid $n2"]
> set sel4 [atomselect top "resid $n3"]
>
> set sel11 [measure center $sel1 weight mass]
> set sel21 [measure center $sel2 weight mass]
> set sel31 [measure center $sel3 weight mass]
> set sel41 [measure center $sel4 weight mass]
>
> set dih [measure dihed [list $sel11 $sel21 $sel31 $sel41]]
>
> $sel1 delete
> $sel2 delete
> $sel3 delete
> $sel4 delete
> $sel11 delete
> $sel21 delete
> $sel31 delete
> $sel41 delete
>
> incr n
> set result [open "dihedral.xvg" a]
> puts $result "$dih"
> close $result
> }
> }
> -----------------
> best,
> Hwankyu.
>
>