# measure tilt angle through the rms fitting with projection (+/- sign)
# using least-square plane
# decompose all three angles
source tilt-trj-setup.tk
# rms fitting whole trajectory to the first snap shot
proc rms_first { sel ref } {
set all [atomselect [$sel molid] "all"]
set nf [molinfo [$sel molid] get numframes]
for { set i 1 } { $i < $nf } { incr i } {
$sel frame $i
$all frame $i
$all move [measure fit $sel $ref]
}
}
# project a vector to a certain plane
# Usage: projection $v1 x ==> project vector v1 to y-z plane
proc projection { vec axis } {
if {$axis == "x"} {
set vec "0 [lindex $vec 1] [lindex $vec 2]"
} elseif {$axis == "y"} {
set vec "[lindex $vec 0] 0 [lindex $vec 2]"
} else {
set vec "[lindex $vec 0] [lindex $vec 1] 0"
}
return $vec
}
# return rotation matrix to align $v1 to $v2 around $sel
proc orient { sel vector1 vector2 } {
set com [measure center $sel]
set vec1 [vecnorm $vector1]
set vec2 [vecnorm $vector2]
# compute the angle and axis of rotation
set rotvec [veccross $vec1 $vec2]
set sine [veclength $rotvec]
set cosine [vecdot $vec1 $vec2]
set angle [expr atan2($sine,$cosine)]
# return the rotation matrix
return [trans center $com axis $rotvec $angle rad]
}# start main program
proc tilt_trj { sel file } {
set fout [open $file w]
set moli [$sel molid]
set nf [molinfo $moli get numframes]
set r2d [expr 180.0/acos(-1.0)]
# align plane normal to Z axis
set all [atomselect $moli "all"]
set ref [atomselect $moli "sta and name CA"]
$all frame 0
$ref frame 0
set zcur [lsqplane $ref]
# check the sign and rotate
if {[lindex $zcur 2] < 0.0} { set zcur [vecscale $zcur -1]}
set A [orient $ref $zcur {0 0 1}]
$all move $A
# move to origin
$all move [transoffset [vecinvert [measure center $ref]]]
# align the selected chain to X
set sel2 [atomselect $moli "[$sel text] and sta"]
$sel2 frame 0
set v1 [measure center $sel2]
set v1 [vecnorm [projection $v1 z]]
set a1 [expr acos([vecdot $v1 {1 0 0}])]
# check sign and rotate
set sig [lindex [veccross $v1 {1 0 0}] 2]
if { $sig < 0.0 } { set a1 [expr -1.0 * $a1]}
$all move [transaxis z $a1 rad]
# rms fitting
set ref [atomselect $moli "[$sel text]" frame 0]
rms_first $sel $ref
set ref [atomselect $moli "sta and name CA"]
for { set i 1 } { $i < $nf } { incr i } {
$all frame $i
$ref frame $i
$sel2 frame $i
set zcur [lsqplane $ref]
if {[lindex $zcur 2] < 0.0} { set zcur [vecscale $zcur -1]}
set zcur [vecnorm $zcur]
set xcur [vecsub [measure center $sel2] [measure center $all]]
# projection to zcur
set tmp [vecscale $zcur [vecdot $xcur $zcur]]
# plane projected xcur
set xcur [vecsub $xcur $tmp]
set xcur [vecnorm $xcur]
# calc radial tilt and rotation angle
set xpro [vecnorm [projection $xcur z]]
set radt [expr $r2d*acos([vecdot $xcur $xpro])]
if {[lindex $xcur 2] < 0.0} { set radt [expr -1.0*$radt]}
set rot [expr $r2d*acos([vecdot $xpro {1 0 0}])]
if {[lindex $xcur 1] < 0.0} { set rot [expr -1.0*$rot]}
# rotate around the subunit to remove the radial tilt
set A [orient $sel2 $xcur $xpro]
$all move $A
# selection update
$all frame $i
$sel2 frame $i
# rotate around Z to remove rotation
set A [orient $sel2 $xpro {1 0 0}]
$all move $A
# Now the pure lateral tilt (and translation) exist
# selection update
$sel2 frame $i
$ref frame $i
# recalc normal vector and tilt angle
set zcur [lsqplane $ref]
if {[lindex $zcur 2] < 0.0} { set zcur [vecscale $zcur -1]}
set zcur [vecnorm $zcur]
set latt [expr $r2d*acos([vecdot $zcur {0 0 1}])]
if {[lindex $zcur 1] > 0.0} { set latt [expr -1.0*$latt]}
puts $fout "[expr $i*0.01] $radt $rot $latt"
}
close $fout
}
set sel [atomselect 1 "name CA and chain A"]
tilt_trj $sel sS1.tilt
set sel [atomselect 1 "name CA and chain B"]
tilt_trj $sel sS2.tilt
set sel [atomselect 1 "name CA and chain C"]
tilt_trj $sel sS3.tilt
set sel [atomselect 1 "name CA and chain D"]
tilt_trj $sel sS4.tilt
set sel [atomselect 1 "name CA and chain E"]
tilt_trj $sel sS5.tilt