## VMD-L Mailing List

**From:** Hwankyu Lee (*leeh3_at_nhlbi.nih.gov*)

**Date:** Wed Feb 25 2009 - 14:55:39 CST

**Next message:**Peter Freddolino: "Re: namd-l: charm27 forcefield_topology"**Previous message:**Yinglong Miao: "to get angle from matrix of rotation around an axis"**In reply to:**J T: "Re: question about "measure dihed" with COM points"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

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

*>>>
*

*>>>
*

*>>
*

*>>
*

*>
*

**Next message:**Peter Freddolino: "Re: namd-l: charm27 forcefield_topology"**Previous message:**Yinglong Miao: "to get angle from matrix of rotation around an axis"**In reply to:**J T: "Re: question about "measure dihed" with COM points"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]