VMD-L Mailing List

From: J T (JTibbitt_at_odu.edu)
Date: Wed Feb 25 2009 - 14:43:01 CST

Hi. Right inside the for loop, you need to put this inside your loop:
molinfo top set frame $pp

Hopefully the calc script will work after yours is debugged. Just
for the heck of it, I attached the script I use for looping through
the frames and residues to calculate the backbone dihedrals for each
frame. Inside the double loop, it calls the calc procedure three
times for psi, omega and phi, once it has the coordinates it needs.
Just curious, why are you doing center of masses of residues?

Jeff

# Retrieve 4 Sets of Cartesian Coordinates Then
# Call Above Calculation Procedure to Calculate Dihedral Angles
(trajectory)
proc dihe {molnum} {
   set nf [molinfo $molnum get numframes]
   set nr [llength [lsort -integer -unique [[atomselect top all] get
residue]]]
   puts "frame resid psi omega phi"
   for { set i 0 } { $i < 10 } { incr i } {
     molinfo $molnum set frame $i
     for { set rnum 1 } { $rnum < $nr } { incr rnum } {
       set xyz1 [lindex [[atomselect $molnum "resid $rnum and name
N"] get {x y z}] 0 ]
       set xyz2 [lindex [[atomselect $molnum "resid $rnum and name
CA"] get {x y z}] 0 ]
       set xyz3 [lindex [[atomselect $molnum "resid $rnum and name
C"] get {x y z}] 0 ]
       set xyz4 [lindex [[atomselect $molnum "resid [expr {$rnum +
1}] and name N"] get {x y z}] 0 ]
       set xyz5 [lindex [[atomselect $molnum "resid [expr {$rnum +
1}] and name CA"] get {x y z}] 0 ]
       set xyz6 [lindex [[atomselect $molnum "resid [expr {$rnum +
1}] and name C"] get {x y z}] 0 ]
       set psi [calc $xyz1 $xyz2 $xyz3 $xyz4]
       set omega [calc $xyz2 $xyz3 $xyz4 $xyz5]
       set phi [calc $xyz3 $xyz4 $xyz5 $xyz6]
       puts [format "%4s %6s %10.1f %7.1f %7.1f" $i $rnum $psi $omega
$phi]
     }
   }
}

On Feb 25, 2009, at 2:35 PM, Lee, Hwan Kyu (NIH/NHLBI) [E] wrote:

> 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.
>>
>>
>
>