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