# Procedure to draw principle axes of the ellipsoids
# based on ANISOU records from a PDB file.
# 
# VERSION 1.0
# 
# By Aaron Oakley 
# Research School of Chemistry
# Australian National University
# Email: oakley@rsc.anu.edu.au
# Complaints: /dev/null
# 
# Please cite:  Burden, C.J. and Oakley, A.J. (2007) 
#               Anisotropic Atomic Motions in High Resolution 
#               Protein Crystallography Molecular Dynamics 
#               Simulations. Physical Biology 4:79-90
# 
# you need to source the Hume linear algebra library (la.tcl)
# before use! See http://www.hume.com/la/la.html and
# licence conditions therein.
# 
# you need to source draw_anisou.tcl
# 
# This file is distributed WITHOUT ANY WARRANTY;
# without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE. Should the
# software prove defective YOU ASSUME THE COST OF
# ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
# 
# 
# Usage eg:
# draw_anisou file.pdb "segname A and noh"


proc draw_anisou {filename sel} {
   set pi 3.14159265

# the following swithches can be used to customize behavior:
   set npd_flag 0 ; # draw Non-Positive Definite atoms
                    # using b-factor yes=1 no=0
   set anisoub  0 ; # draw atoms without ANISOU record 
                    # using the isotropic B-factor. yes=1 no=0 

   set radius 0.07   ; # Radius for drawing axes
   set resolution 10 ; # Resolution of axis cylinders

   graphics 0 delete all

   set selindex [[atomselect top "$sel"] get serial] 

   set file [open $filename r]

   set eof 1

   # Start the read cycle:
   while {$eof >= 0} {
      # Get the next line in the input PDB file:
      set eof [gets $file line]


      set record [lindex $line 0]

      set atom_flag [string compare $record "ATOM"]
      set hetatm_flag [string compare $record "HETATM"]
      set anisou_flag [string compare $record "ANISOU"]

      if {$atom_flag == 0 || $hetatm_flag == 0} {
         
         set serial [string trim [string range $line 6 10]]
         lappend seriallist $serial
         set atom($serial) [string range $line 12 15]
         set altloc($serial) [string range $line 16 16]
         set resname($serial) [string range $line 17 19]
         set chain($serial) [string range $line 21 21]
         set resid($serial) [string range $line 22 26]

         set x($serial) [string range $line 30 37]
         set y($serial) [string range $line 38 45]
         set z($serial) [string range $line 46 53]
         set b($serial) [string range $line 60 65]
      }

      if {$anisou_flag == 0} {
         set serial [string trim [string range $line 6 10]]
         set u11($serial) [expr [string range $line 28 34] / 10000.0]
         set u22($serial) [expr [string range $line 35 41] / 10000.0]
         set u33($serial) [expr [string range $line 42 48] / 10000.0]
         set u12($serial) [expr [string range $line 49 55] / 10000.0]
         set u13($serial) [expr [string range $line 56 62] / 10000.0]
         set u23($serial) [expr [string range $line 63 69] / 10000.0]
      }
   } ; # end of read cycle

   # finished with PDB file
   close $file

   # Start the draw cycle:
   foreach r $selindex {
      set iso 0
      set draw 0

      if {[info exists atom($r)] == 1} {
       puts "Atom $r $atom($r) $altloc($r) $resname($r) $chain($r) $resid($r)"
       if {[info exists u11($r)] == 1} {

          # Use "mevsvd" procedure in la1.0 to find eigenvectors/eigenvalues
          # of ANISOU matrix by singular value decomposition
          set matrix "2 3 3 $u11($r) $u12($r) $u13($r) $u12($r) $u22($r) $u23($r) $u13($r) $u23($r) $u33($r)"
          La::mevsvd_br matrix evals
 
          #get the eiginvectors and eigenvalues
          set x1 [lindex $matrix 3] ; set y1 [lindex $matrix 6] ; set z1 [lindex $matrix 9]
          set x2 [lindex $matrix 4] ; set y2 [lindex $matrix 7] ; set z2 [lindex $matrix 10]
          set x3 [lindex $matrix 5] ; set y3 [lindex $matrix 8] ; set z3 [lindex $matrix 11]
          set s1 [lindex $evals 3]
          set s2 [lindex $evals 4]
          set s3 [lindex $evals 5]
 
          puts "Eigenvector, eigenvalue: $x1 $y1 $z1 , $s1"
          puts "Eigenvector, eigenvalue: $x2 $y2 $z2 , $s2"
          puts "Eigenvector, eigenvalue: $x3 $y3 $z3 , $s3"
 
          set draw 1
 
          if {$s1 < 0 || $s2 < 0 || $s3 < 0} {
             puts "THIS ANISOU MATRIX IS NON-POSITIVE DEFINITE"
             if {$npd_flag == 1} {
                set iso 1 ; set draw 1
             } else {set iso 0 ; set draw 0} 
          }
          puts " "
       } else {
          if {$anisoub == 1} {
             set iso 1
          }
       } ; # end if u11 exists
 
       # in the case that the ANISOU is NPD or there is no ANISOU info, use the B factor:
       if {$iso == 1} {
          puts "The isotropic B-factor will be used to plot the axes: $b($r)"
          set u [expr ($b($r) / (8.0 * $pi * $pi))]
          puts "u = $u"
          puts " "
          set x1 1 ; set y1 0 ; set z1 0
          set x2 0 ; set y2 1 ; set z2 0
          set x3 0 ; set y3 0 ; set z3 1
          set s1 $u ; set s2 $u ; set s3 $u
          set draw 1
       }


       if {$draw == 1} {
          set s1 [expr sqrt($s1)] ; set s2 [expr sqrt($s2)] ; set s3 [expr sqrt($s3)]
 
          set element [string range $atom($r) 1 1]
          set col green
          foreach e {H C N O S P} c {white cyan blue red yellow purple} {
             if {[string compare $element $e] == 0} {set col $c}
          }
 
          foreach {xx yy zz ss} {$x1 $y1 $z1 $s1 $x2 $y2 $z2 $s2 $x3 $y3 $z3 $s3} {
             #the value of ss is the unit variance
             #To draw an ellipsoid such that the probability
             #of finding the atom inside the ellipsoid is
             #50%,  the critical value is 1.5382*ss. See ORTEP-III manual
 
             set ss [expr $ss * 1.5382]
             set sx [expr $x($r) - $xx * $ss]
             set sy [expr $y($r) - $yy * $ss]
             set sz [expr $z($r) - $zz * $ss]
             set ex [expr $x($r) + $xx * $ss]
             set ey [expr $y($r) + $yy * $ss]
             set ez [expr $z($r) + $zz * $ss]
             set start "$sx $sy $sz" ; set end "$ex $ey $ez"
 
             graphics 0 color $col
             graphics 0 cylinder $start $end radius $radius resolution $resolution
          } ; # end foreach
       } ; # end if

      } ; # end if atom exists

     } ; # END OF DRAW CYCLE

} ; # END OF PROC