VMD-L Mailing List

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

Hi, thanks all. It was actually the NumberOfOxygens variable that was
going to zero, per the previous selection. I was able to circumvent it
with an if statement. So the "Divide by zero" error really does mean
divide by zero!

Arneh Babakhani wrote:
> 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
> }
>
>
>
>