#
# readcube, a script to read gaussian "cube" files in vmd
# $Id: readcube.tcl,v 1.2 2005/01/12 00:55:23 johns Exp $
# First released version
#
# Gaussian "cube" format described here:
#   http://www.gaussian.com/00000430.htm 
#
# Pease note that normally the cutoff to be used
# with cubefiles is much smaller than the standard
# vmd cutoff (e.g. either you multiply your cube values
# by 1E3 or you set the cutoff manually)
#
# Example usage:
#  source readcube.tcl
#  mol pdbload 1bnf
#  readcube name.cube
#    [then make new isosurface and volume_slice representations]
#
# Authors
# -------
# Written by Francesco Luigi Gervasio
# CSCS, Manno, Switzerland
# 
# Small fixes and improvements by John Stone
# Theoretical Biophysics Group
# University of Illinois
#


proc readcube {fileName} {
  if [ catch {open $fileName r} in] {
    puts stderr "Cannot open the file $fileName"

    return;  # exit early, no need to close files etc.
  } else {
    set topLabel [gets $in ]  
    puts "topLabel = >$topLabel<"

    set topLabel [gets $in ]  
    puts "Second line = >$topLabel<"

    set natoms [read $in 6]
    if {$natoms > 0} {
      # unformatted file format, what do we do specially?
    } else {
      # formatted file format, what do we do specially?
      # fix sign if negative
      set natoms [expr -$natoms]
    }
    puts "Number of atoms: $natoms"

    ##
    ## Read in gaussian's origin (non scaled yet..)
    ##
    set xOrigin [read $in 12]
    set yOrigin [read $in 12]
    set zOrigin [read $in 12]

    ##
    ## abs(N1) is number of points, sign indicates the type of
    ## units being used, whether Bohr or Angstroms. 
    ##
    set xgrid [read $in 6]
    if {$xgrid > 0} {
      # This is the conversion factor between a.u. and angstroems
      set convfact 0.529177;     # Grid units are in angstroms already 
    } else {
      # This is the conversion factor between a.u. and angstroems
      set convfact 0.529177;     # what should this really be???
      set xgrid [expr -$xgrid];  # fix sign if negative
    }

    ##
    ## Calculate the X/Y/Z box side vectors for VMD
    ## 
    set xVec [list [expr ([read $in 12] * $convfact * $xgrid)] \
                   [expr ([read $in 12] * $convfact * $xgrid)] \
                   [expr ([read $in 12] * $convfact * $xgrid)]]

    set ygrid [read $in 6]
    set yVec [list [expr ([read $in 12] * $convfact * $ygrid)] \
                   [expr ([read $in 12] * $convfact * $ygrid)] \
                   [expr ([read $in 12] * $convfact * $ygrid)]]

    set zgrid [read $in 6]
    set zVec [list [expr ([read $in 12] * $convfact * $zgrid)] \
                   [expr ([read $in 12] * $convfact * $zgrid)] \
                   [expr ([read $in 12] * $convfact * $zgrid)]]

    ##
    ## Convert from Gaussian to VMD volume origin
    ##
    set xOrigin [expr $xOrigin * $convfact]
    set yOrigin [expr $yOrigin * $convfact]
    set zOrigin [expr $zOrigin * $convfact]
    puts "Origin: $xOrigin $yOrigin $zOrigin"
 
    puts " xgrid and vectorX: $xgrid $xVec"
    puts " ygrid and vectorY: $ygrid $yVec"
    puts " zgrid and vectorZ: $zgrid $zVec"

    ## 
    ## Read in the atom data (atomic number, charge, x/y/z coords)
    ##
    for {set j 0} {$j < $natoms} {incr j} {
      ##
      ## test code to draw spheres to verify that the 
      ## coordinate systems match up correctly.
      ##
      # set anum [read $in  5]
      # set qq   [read $in 12]
      # set xx   [expr [read $in 12] * $convfact]
      # set yy   [expr [read $in 12] * $convfact]
      # set zz   [expr [read $in 12] * $convfact]
      # draw sphere [list $xx $yy $zz] radius 0.5

      set atoms [gets $in];  # eat any remaining characters
    }

    #############################
    # Here it reads the grid data
    #############################
    set valList "" 
    set valCount [expr $xgrid*$ygrid*$zgrid]
    while { [gets $in line] >=0} {
      set linelength [llength $line]
      #puts "lunghezza: $linelength"
      set readList $line
      append valList " "
      append valList $readList 
    }   

    puts "Length of list is [llength $valList], wanted total= $valCount"
  }

  close $in
  mol volume [molinfo top] "Gaussian cube" [list $xOrigin $yOrigin $zOrigin] $xVec $yVec $zVec $xgrid $ygrid $zgrid $valList
}