# read_dx:         An OpenDX-format grid reader for VMD
# Primary author:  Dave Sept - Washington University
# Hacked by:       Nathan Baker - Washington University
# Version:         $Id: my_functions,v 1.8 2005/12/07 16:54:20 todd_dolinsky Exp $
#
proc read_dx {fileName molNum} {
#
# This procedure reads in a dx grid and associates it with
# the molecule <molNum>
#
    puts "read_dx:  A program to read APBS OpenDX-format data into VMD"
    puts "Authors:  Dave Sept and Nathan Baker"
    puts "WARNING!  THIS PROGRAM ASSUMES YOU HAVE THE DATA ORDERED SUCH"
    puts "THAT Z INCREASES MOST QUICKLY, Y INCREASES NdnaEXT MOST QUICKLY,"
    puts "AND X INCREASES MOST SLOWLY."
    puts ""
    puts "Reading header information..."
    set valList ""
    set count 0
 
# the line below opens the filename for reading purposes only
    set in [open $fileName r] 
#---------------------------------------
# Gotcha... hehehehe
    set TitleLine [gets $in] 
    puts "Grid title:  $TitleLine"
    set InputLine [gets $in]
    set InputLine [gets $in]
    set InputLine [gets $in]
    set InputLine [gets $in]
#---------------------------------------- 
   scan $InputLine "object 1 class gridpositions counts %i %i %i" xgrid ygrid zgrid
    puts "Grid Points per Dimension = x:$xgrid y:$ygrid z:$zgrid"
    set InputLine [gets $in]
    scan $InputLine "origin %e %e %e" xOrigin yOrigin zOrigin
    set origin [list $xOrigin $yOrigin $zOrigin ]
    puts "Origin at: $xOrigin $yOrigin $zOrigin" 
    set InputLine [gets $in]
    scan $InputLine "delta %e %e %e" xdel dum2 dum1
    set InputLine [gets $in]
    scan $InputLine "delta %e %e %e" dum1 ydel dum2
    set InputLine [gets $in]
    scan $InputLine "delta %e %e %e" dum2 dum1 zdel
    set xVec [list [expr $xdel * $xgrid] 0 0 ]
    set yVec [list 0 [expr $ydel * $ygrid] 0 ]
    set zVec [list 0 0 [expr $zdel * $zgrid] ]
    puts "Grid Spacing = x:$xdel y:$ydel z:$zdel"
    # two more object lines - not needed
    set InputLine [gets $in]
    set InputLine [gets $in]
    set total [expr int($xgrid* $ygrid * $zgrid / 3)]
    puts "Reading data values..."

    while {$count < $total} {
        set InputLine [gets $in]
	  scan $InputLine "%e %e %e" v1 v2 v3
        lappend valList $v1 $v2 $v3
        incr count
    }

    if {[expr $xgrid * $ygrid * $zgrid - 3 * $total] == 2} {
	  scan $InputLine "%e %e" v1 v2
        lappend valList $v1 $v2 $v3
#the above line originally had only $v1 and $v2 and I added $v3
    }

    if {[expr $xgrid * $ygrid * $zgrid - 3 * $total] == 1} {
	  scan $InputLine "%e"  v1
        lappend valList $v1 $v2 $v3
# the above line originally had only $v1 and I added $v2 and $v3
    }

    # The value list we read has z increasing most quickly; however, VMD wants
    # x to increase most quickly

    puts "Reordering data values..."

    set dataList ""
    for {set k 0} {$k<$zgrid} {incr k} {
        for {set j 0} {$j<$ygrid} {incr j} {
            for {set i 0} {$i<$xgrid} {incr i} {
                set u [expr $i * $zgrid * $ygrid + $j * $zgrid + $k]
                lappend dataList [lindex $valList $u]
            }
        }
    }

      
    puts "Length of original list is [llength $valList]"
    puts "Length of reordered list is [llength $dataList]"
    puts "Passing volume data for molecule $molNum to VMD..."
    mol volume $molNum $TitleLine  $origin $xVec $yVec $zVec $xgrid $ygrid $zgrid $dataList

    close $in      
}

proc add_contour {molNum value color} {
#
# This procedure adds an isocontour to molecule <molNum>
# with the value <value> and colors it with colorId <color>
#
	set repnum [molinfo $molNum  get numreps]
	mol color colorid $color
	mol material Transparent
	mol addrep $molNum
	mol selection all
	mol modstyle $repnum $molNum Isosurface $value 0.0 0.0 0.0
}

proc load_pot {} {
    mol load pdb sphere.pqr
    mol modstyle 0 0 vdw
    rock x by 2
    display update
    read_dx pot.dx 0
    add_contour 0 -0.5 1
    add_contour 0 -1.0 3
    add_contour 0 -1.5 4
    add_contour 0 -2.0 7
    add_contour 0 -2.5 0
}

proc load_qdens {} {
    mol load pdb sphere.pqr
    mol modstyle 0 0 vdw
    rock x by 2
    display update

    read_dx qdens.dx 0
    add_contour 0 2.000 1
    add_contour 0 1.000 3
    add_contour 0 0.500 4
    add_contour 0 0.250 7
    add_contour 0 0.125 0

    display eyesep       0.065000
    display focallength  2.000000
    display height       6.000000 display distance     -2.000000
    display projection   Orthographic
    display nearclip set 0.500000
    display farclip  set 10.000000
    display depthcue   off
    display cuestart   0.500000
    display cueend     10.000000
    display cuedensity 0.400000
    display cuemode    Exp2
    display update

}