VMD-L Mailing List

From: Andres Morales N (andresmoralesn2_at_hotmail.com)
Date: Tue May 19 2009 - 11:14:21 CDT

Hi!
 
I wrote a script to align and calculate average strcuture form a .dcd file.
I wait you can look it and give me some suggestions.
 
when I run this script appears the following mistake: domain error: argument not in valid range
 So this mistake does not appears everytime I run the script.
 
 
I wait somebody can help me with this.

##### Align
##move atoms of each strcuture from a dcd file to a position where alpha carbon of first residue will be located at point {0 0 0}
set O [atomselect top "protein and resid 1 and alpha"]
set nf [molinfo top get numframes]
set protein [atomselect top all]
for {set i 0 } {$i < $nf } { incr i } {
$O frame $i
$protein frame $i
set ocoord [lindex [$O get {x y z}] 0]
$protein moveby [vecscale -1 $ocoord]
}
## Align vectors (between CA of resid 1 and 15) of each frame
set P0 [atomselect top "protein and resid 15 and alpha" frame 0]
set P0coord [lindex [$P0 get {x y z}] 0]
set P [atomselect top "protein and resid 15 and alpha"]
for {set i 0 } {$i < $nf } { incr i } {
$P frame $i
$protein frame $i
set Pcoord [lindex [$P get {x y z}] 0]
global M_PI
set cosine [expr [vecdot $Pcoord $P0coord] / ( [veclength $Pcoord] * [veclength $P0coord])]
set angulo1 [expr acos($cosine)*(180.0/$M_PI)]
set A [veccross $Pcoord $P0coord]
set M1 [transabout $A $angulo1]
$protein move $M1
}
 
## Align plane (formed between CA of resid 1 15 and 50) of each frame
set Q0 [atomselect top "protein and resid 50 and alpha" frame 0]
set Q0coord [lindex [$Q0 get {x y z}] 0]
set Q [atomselect top "protein and resid 50 and alpha"]
set N0 [veccross $P0coord $Q0coord]
set O0 [atomselect top "protein and resid 1 and alpha" frame 0]
set O0coord [lindex [$O0 get {x y z}] 0]
for {set i 0 } {$i < $nf } { incr i } {
$P frame $i
$Q frame $i
$protein frame $i
set Pcoord [lindex [$P get {x y z}] 0]
set Qcoord [lindex [$Q get {x y z}] 0]
set N [veccross $Pcoord $Qcoord]
global M_PI
set cosine2 [expr [vecdot $N $N0] / ( [veclength $N] * [veclength $N0])]
set angulo2 [expr acos($cosine2)*(180.0/$M_PI)]
set B [veccross $N $N0]
set M2 [transabout $B $angulo2]
$protein move $M2
}
 
 ### Average strcuture
 
set sel [atomselect top all]
set nf [molinfo top get numframes]
set n [expr $nf - 1]
set newpos [measure avpos $sel first 0 last $n]
$sel set {x y z} $newpos
$sel writepdb "promedio.pdb"
 
 
Thanks for your suggestions
 

Andres

Invite your mail contacts to join your friends list with Windows Live Spaces. It's easy! Try it!
_________________________________________________________________
Invite your mail contacts to join your friends list with Windows Live Spaces. It's easy!
http://spaces.live.com/spacesapi.aspx?wx_action=create&wx_url=/friends.aspx&mkt=en-us