VMD-L Mailing List

From: Arneh Babakhani (ababakha_at_mccammon.ucsd.edu)
Date: Thu Jul 05 2007 - 14:57:03 CDT

Hi,

Is it possible to get a "Divide by Zero" error when running a tcl script
in VMD when you're not really dividing by zero?!

I'm running a script that loops through the frames of my traj. It
aborts with the "Divide by Zero" error when it gets to a frame (but not
always at the same frame, seems to be random where it aborts). In my
script, there's only two places that I'm actually dividing anything. So
I used a "puts" line to output the variable in the denominator, to check
it out. In neither of the two places, the variable in question is never
zero when it aborts, not even close.

Any ideas???? Script is pasted below if you'd like to look at it,

Thanks,

Arneh

-------------

proc CarbAngle {outputfile} {

  set outfile [open $outputfile w]
  set n [expr {1 * [molinfo top get numframes] } ]

  set Membrane [atomselect top "resname DMPC"]
  set Indole [atomselect top "resname TRP and sidechain"]

  for {set frame 0} {$frame < $n} {incr frame} {

    $Membrane frame $frame
    set MembraneCenter [measure center $Membrane]
    set MembraneCenterZ [lindex $MembraneCenter 2]

    $Indole frame $frame
    set IndoleCenter [measure center $Indole]
    set IndoleCenterZ [lindex $IndoleCenter 2]

    set Xmin [expr [lindex $IndoleCenter 0] - 5.00]
    set Xmax [expr [lindex $IndoleCenter 0] + 5.00]
    set Ymin [expr [lindex $IndoleCenter 1] - 5.00]
    set Ymax [expr [lindex $IndoleCenter 1] + 5.00]

    if {$IndoleCenterZ > $MembraneCenterZ} {
      set CarbonylOs [atomselect top "(name OH or name OF) and (x<$Xmax
and x>$Xmin) and (y<$Ymax and y>$Ymin) and (z>$MembraneCenterZ)"]
    } else {
      set CarbonylOs [atomselect top "(name OH or name OF) and (x<$Xmax
and x>$Xmin) and (y<$Ymax and y>$Ymin) and (z<$MembraneCenterZ)"]
    }

    set Oxygens [$CarbonylOs get index]
    set NumberOfOxygens [llength $Oxygens]

    set Angle 0

    for {set counter 0} {$counter < $NumberOfOxygens} {incr counter} {

      set CurrentOxygenIndex [lindex $Oxygens $counter]
      set CurrentOxygen [atomselect top "index $CurrentOxygenIndex"]
      set CurrentOxygenCoord [ lindex [$CurrentOxygen get {x y z}] 0]

      set CurrentCarbon [atomselect top "(name C1A or name C2A) and
within 2 of index $CurrentOxygenIndex"] ; $CurrentCarbon frame $frame
      set CurrentCarbonCoord [ lindex [$CurrentCarbon get {x y z}] 0]

      set Carbonyl [vecsub $CurrentOxygenCoord $CurrentCarbonCoord]
      set CarbonylLength [veclength $Carbonyl]
      set CarbonylDot [vecdot $Carbonyl {0 0 1}] ; set CarbonylDotLength
[veclength $CarbonylDot]
      puts $outfile "$frame $CarbonylLength"
      set CurrentAngleRadians [expr
{acos($CarbonylDotLength/$CarbonylLength)}]
      set Angle [expr $Angle + $CurrentAngleRadians]

    }
    puts $outfile "$frame $NumberOfOxygens"
    set Angle [expr $Angle / $NumberOfOxygens ]

    set Angle [expr {$Angle * 180/3.14159}]

  }

  close $outfile
}