# AMBER Archive (2006)Subject: Re: AMBER: A BUG report, "ptraj corrplane"

From: Myunggi Yi (myunggi_at_gmail.com)
Date: Sun Dec 10 2006 - 19:32:22 CST

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_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu