Dear Amber users,
I had modified the source code of the function "lsqplane" in
"action.c" of ptraj.
I think I converted the codes of subroutines "solve_cubic_eq" and
"lsqplane" to tcl code correctly.
After calculation, I've got strange results.
The normal vectors that I got was about 90 degrees different from that
I expected.
What are the "a b c" in the code?
Is this vector perpendicular to a plane?
The following is my tcl code.
===============================
proc solve_cubic_eq { a b c d } {
#
# Solves a cubic equation
# ax^3 + bx^2 + cx + d = 0
# using "Cardan's formula"
# (see: Bronstein, S.131f)
#
set PI 3.141592654
set one3 [expr 1.0 / 3.0]
set one27 [expr 1.0 / 27.0]
# Coeff. for normal form x^3 + rx^2 + sx + t = 0
set r [expr $b / $a]
set s [expr $c / $a]
set t [expr $d / $a]
# Coeff. for red. eq. y^3 + py + q = 0 with y = x + r/3 bzw. (x = y - r/3)
set p [expr $s - $r * $r * $one3]
set q [expr 2.0 * $r * $r * $r * $one27 - $r * $s * $one3 + $t]
# Dummy variables
set rho [expr sqrt(-1.0 * $p * $p * $p * $one27)]
set phi [expr acos(-1.0 * $q / (2.0 * $rho))]
# Discriminante(?)
set D [expr pow(($p * $one3),3.0) + $q * $q * 0.25]
# x real -> one real solution
# three real solutions (d < 0) | one real solution + one real double solution
# or one real triple solution (d = 0)
if { $D > 0 } {
set u [expr pow((-1.0 * $q * 0.5 + sqrt($D)), $one3)]
set v [expr -1.0 * $p / $u * $one3]
return [expr ($u + $v) - $r * $one3]
} else {
lappend dtmp [expr 2.0 * pow($rho, $one3) * cos($phi * $one3) - $r * $one3]
lappend dtmp [expr 2.0 * pow($rho, $one3) * cos(($phi + 2.0 * $PI)
* $one3) - $r * $one3]
lappend dtmp [expr 2.0 * pow($rho, $one3) * cos(($phi + 4.0 * $PI)
* $one3) - $r * $one3]
lsort -real $dtmp
return [lindex $dtmp 0]
}
unset dtmp
}
# main calculation
proc lsqplane { sel file } {
set fout [open $file w]
set nf [molinfo 0 get numframes]
set n [$sel num]
set r2d [expr 180.0/acos(-1.0)]
for { set i 0 } { $i < $nf } { incr i } {
$sel frame $i
set crd [$sel get { x y z }]
set cg [measure center $sel]
for { set j 0 } { $j < $n } { incr j } {
lappend cx [expr [lindex [lindex $crd $j] 0] - [lindex $cg 0]]
lappend cy [expr [lindex [lindex $crd $j] 1] - [lindex $cg 1]]
lappend cz [expr [lindex [lindex $crd $j] 2] - [lindex $cg 2]]
}
#
# if { $n == 3 } {
# set x1 [expr [lindex $cx 1] - [lindex $cx 0]]
# set y1 [expr [lindex $cy 1] - [lindex $cy 0]]
# set z1 [expr [lindex $cz 1] - [lindex $cz 0]]
# set x2 [expr [lindex $cx 2] - [lindex $cx 1]]
# set y2 [expr [lindex $cy 2] - [lindex $cy 1]]
# set z2 [expr [lindex $cz 2] - [lindex $cz 1]]
# set a [expr $y1 * $z2 - $z1 * $y2]
# set b [expr $z1 * $x2 - $x1 * $z2]
# set c [expr $x1 * $y2 - $y1 * $x2]
# } else {
#
set dSumXX 0.0
set dSumYY 0.0
set dSumZZ 0.0
set dSumXY 0.0
set dSumXZ 0.0
set dSumYZ 0.0
for { set k 0 } { $k < $n } { incr k } {
set dSumXX [expr $dSumXX + [lindex $cx $k] * [lindex $cx $k]]
set dSumYY [expr $dSumYY + [lindex $cy $k] * [lindex $cy $k]]
set dSumZZ [expr $dSumZZ + [lindex $cz $k] * [lindex $cz $k]]
set dSumXY [expr $dSumXY + [lindex $cx $k] * [lindex $cy $k]]
set dSumXZ [expr $dSumXZ + [lindex $cx $k] * [lindex $cz $k]]
set dSumYZ [expr $dSumYZ + [lindex $cy $k] * [lindex $cz $k]]
}
# Calc coeff. for -l^3 + o * l^2 + p * l + q = 0
set o [expr $dSumXX + $dSumYY + $dSumZZ]
set p [expr $dSumXY * $dSumXY + $dSumXZ * $dSumXZ + $dSumYZ * \
$dSumYZ - ($dSumXX * $dSumYY + $dSumXX * $dSumZZ + $dSumYY * \
$dSumZZ)]
set q [expr $dSumXX * $dSumYY * $dSumZZ + 2.0 * $dSumXY * $dSumXZ \
* $dSumYZ - ($dSumXX * $dSumYZ * $dSumYZ + $dSumYY * $dSumXZ \
* $dSumXZ + $dSumZZ * $dSumXY * $dSumXY)]
# Solve cubic eq.
set root [solve_cubic_eq -1.0 $o $p $q]
# Calc determinantes
set a [expr ($dSumYY - $root) * $dSumXZ - $dSumXY * $dSumYZ]
set b [expr ($dSumXX - $root) * $dSumYZ - $dSumXY * $dSumXZ]
set c [expr $dSumXY * $dSumXY - ($dSumYY - $root) * ($dSumXX - $root)]
# }
# Normalize
set dnorm [expr 1.0 / sqrt( $a * $a + $b * $b + $c * $c)]
set a [expr $a * $dnorm]
set b [expr $b * $dnorm]
set c [expr $c * $dnorm]
lappend vec $a $b $c
set a0 [expr $r2d * acos($a/[veclength $vec])]
set a1 [expr $r2d * acos($b/[veclength $vec])]
set a2 [expr $r2d * acos($c/[veclength $vec])]
puts $fout "[expr ($i+1.0)/1000.0] $vec $a0 $a1 $a2"
unset cx cy cz
unset vec
}
close $fout
}
=============================
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Wed Dec 13 2006 - 05:21:32 PST
