From: Axel Kohlmeyer (akohlmey_at_cmm.chem.upenn.edu)
Date: Thu Jul 05 2007 - 14:53:16 CDT

On Thu, 5 Jul 2007, Arneh Babakhani wrote:

AB> Hi,

hi arneh,

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

you may not divide directly in your code, but implicitly
in some other procedure, that you are using.

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

well, if you could upload some usable chunk of trajectory, that
reproduces the error to the VMD biocore biofs, we could help trying
to track down the error. at the moment, the only obvious problem,
that i am seeing from a superficial look at the code, is that you have
a few memory leaks in there by creating plenty of atom selections
without deleting them. not sure whether this would generate an error,
but if you run out of memory (what platform are you running on?,
which version of VMD?), it could cause some memory corruption...

cheers,
   axel.

AB>
AB> Thanks,
AB>
AB> Arneh
AB>
AB> -------------
AB>
AB> proc CarbAngle {outputfile} {
AB>
AB> set outfile [open $outputfile w]
AB> set n [expr {1 * [molinfo top get numframes] } ]
AB>
AB> set Membrane [atomselect top "resname DMPC"]
AB> set Indole [atomselect top "resname TRP and sidechain"]
AB>
AB> for {set frame 0} {$frame < $n} {incr frame} {
AB>
AB> $Membrane frame $frame
AB> set MembraneCenter [measure center $Membrane]
AB> set MembraneCenterZ [lindex $MembraneCenter 2]
AB>
AB> $Indole frame $frame
AB> set IndoleCenter [measure center $Indole]
AB> set IndoleCenterZ [lindex $IndoleCenter 2]
AB>
AB> set Xmin [expr [lindex $IndoleCenter 0] - 5.00]
AB> set Xmax [expr [lindex $IndoleCenter 0] + 5.00]
AB> set Ymin [expr [lindex $IndoleCenter 1] - 5.00]
AB> set Ymax [expr [lindex $IndoleCenter 1] + 5.00]
AB>
AB> if {$IndoleCenterZ > $MembraneCenterZ} {
AB> set CarbonylOs [atomselect top "(name OH or name OF) and (x<$Xmax
AB> and x>$Xmin) and (y<$Ymax and y>$Ymin) and (z>$MembraneCenterZ)"]
AB> } else {
AB> set CarbonylOs [atomselect top "(name OH or name OF) and (x<$Xmax
AB> and x>$Xmin) and (y<$Ymax and y>$Ymin) and (z<$MembraneCenterZ)"]
AB> }
AB>
AB> set Oxygens [$CarbonylOs get index]
AB> set NumberOfOxygens [llength $Oxygens]
AB>
AB> set Angle 0
AB>
AB> for {set counter 0} {$counter < $NumberOfOxygens} {incr counter} {
AB>
AB> set CurrentOxygenIndex [lindex $Oxygens $counter]
AB> set CurrentOxygen [atomselect top "index $CurrentOxygenIndex"]
AB> set CurrentOxygenCoord [ lindex [$CurrentOxygen get {x y z}] 0]
AB>
AB> set CurrentCarbon [atomselect top "(name C1A or name C2A) and
AB> within 2 of index $CurrentOxygenIndex"] ; $CurrentCarbon frame $frame
AB> set CurrentCarbonCoord [ lindex [$CurrentCarbon get {x y z}] 0]
AB>
AB> set Carbonyl [vecsub $CurrentOxygenCoord $CurrentCarbonCoord]
AB> set CarbonylLength [veclength $Carbonyl]
AB> set CarbonylDot [vecdot $Carbonyl {0 0 1}] ; set CarbonylDotLength
AB> [veclength $CarbonylDot]
AB> puts $outfile "$frame $CarbonylLength"
AB> set CurrentAngleRadians [expr
AB> {acos($CarbonylDotLength/$CarbonylLength)}]
AB> set Angle [expr $Angle + $CurrentAngleRadians]
AB>
AB> }
AB> puts $outfile "$frame $NumberOfOxygens"
AB> set Angle [expr $Angle / $NumberOfOxygens ]
AB>
AB> set Angle [expr {$Angle * 180/3.14159}]
AB>
AB>
AB>
AB> }
AB>
AB> close $outfile
AB> }
AB>
AB>
AB>
AB>

-- 
=======================================================================
Axel Kohlmeyer   akohlmey_at_cmm.chem.upenn.edu   http://www.cmm.upenn.edu
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.