From: Hwankyu Lee (leeh3_at_nhlbi.nih.gov)
Date: Wed Feb 25 2009 - 14:55:39 CST

Thank you soooo much for your comments. Now, everything works
fine! I'm using COM, because I want to build up coarse-grained
models. Again, Thanks,

best,
Hwankyu.

On Feb 25, 2009, at 3:43 PM, J T wrote:

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