``````

In :

#!/usr/bin/python

'''
Daniel J. Sindhikara
sindhikara@gmail.com
Handles 3D volumetric data
'''

from __future__ import division
import numpy as np
import os

``````
``````

In :

#import sys
#isnotebook ='py' not in sys.argv
#if isnotebook:
#    print "Running as IPython notebook"

``````
``````

In :

class Grid:
'''
Contains volumetric data
'''

def __init__(
self,
distribution,
origin,
gridcount,
deltas,
concentration=-1.0,
):
if type(distribution) is list:
self.distribution = np.array(distribution)
elif type(distribution) is np.ndarray:
self.distribution = distribution
self.origin = origin
self.gridcount = gridcount
self.deltas = deltas
self.concentration = concentration  # in molar

def getvalue(self, coord):
return linearinterpolatevalue(self.distribution, self.origin,
self.deltas, coord)

def writedx(self, filename):
printdxfrom3d(self.distribution, self.origin, self.deltas,
self.gridcount, filename)
def writegrd(self, filename):
printgrd(self.distribution, self.origin, self.deltas,
self.gridcount, filename)

def nearestgridindices(self, coord):
'''
Given a 3D cartesian coordinate, return nearest grid indices
'''

gridindices = [int(round((coord[i] - self.origin[i])
/ self.deltas[i])) for i in range(3)]
return gridindices

def coordfromindices(self, indices):
return [indices[i] * self.deltas[i] + self.origin[i] for i in
range(3)]

def coarseRDF(self, coord):
'''
Given a 3D cartesian coordinate, return 1D distribution as list, spaced by gridspacing
'''

dist = []
myindices = self.nearestgridindices(coord)
for shellindex in shellindices:
avg = 0.0
try:
for indexonshell in shellindex:
avg += self.distribution[myindices
+ indexonshell][myindices
+ indexonshell][myindices
+ indexonshell]
avg = avg / len(shellindex)
dist.append(avg)
except IndexError:
break
return dist

def interpRDF(
self,
coord,
delta,
limit,
concentration = False # If not False, print coordination number
):
'''
Given 3D cartesian coordinate (list-like object of floats),
delta (float), and limit (float),

return RDF

by averaging linear interpolated g(r_vec) on points on spherical shell.
'''

spherefilename = \
os.path.join(os.path.dirname(os.path.realpath(__file__)),
'data', 'points', '200.pts')
spherepoints = [[float(element) for element in line.split()]
rdf = []
volints = [0.]
for radius in np.arange(0, limit, delta):
mysum = 0
for point in spherepoints:
+ coord[dim] for dim in range(3)])
rdf.append(float(mysum / len(spherepoints)))
if not concentration:
return rdf
else:
coordnums = [volint * concentration for volint in volints]
return rdf, coordnums

``````
``````

In :

def rdf23dgrid(
rdf,
rdfdelta,
griddelta,
gridorigin,
shellindices,
):
'''
Convert 1d RDF to 3D grid.
'''

gridcount = [len(shellindices) * 2 + 1] * 3  # + 1 is to make it odd, let
mydist = np.array([[[1.0 for i in range(gridcount)] for j in
range(gridcount)] for k in
range(gridcount)])

indexmultiplier = int(griddelta / rdfdelta)  # from 3D grid index to 1D rdf index
for (i, shell) in enumerate(shellindices):  # loop through shells
myrdf = rdf[i * indexmultiplier]
for shellindex in shell:
distindex = tuple([shellindex[dim] + int(gridcount[dim]
/ 2) for dim in range(3)])

# above will put shell 0 at the int(gridcount/2) aka numshells/2 + 1

mydist[distindex] = myrdf
return Grid(mydist, gridorigin, gridcount, [griddelta] * 3)

``````
``````

In :

def dx2Grid(dxfilename):
'''
Reads a dx file into a Grid class object
'''

# for now Ghetto rigged to use old readdx function

(distribution, origin, deltas, gridcount) = readdx(dxfilename)
return Grid(distribution, origin, gridcount, deltas)

``````
``````

In :

'''
Use gzip.open if it's a .gz file.
'''
from gzip import open as gopen
if 'gz' in filename:
f = gopen(filename,'rb')
else:
f = open(filename,'r')
return(f)

``````
``````

In :

'''
Returns grid data and single 4D array containing dx data [index][xindex][yindex][zindex]
and the origin, deltas and gridcounts.
'''
dxlines = []
for i in range(10):  # only need the first few lines to get grid data
dxfile.close()
gridcounts = []
origin = []
deltas = [0, 0, 0]
startline = 0
for (i, line) in enumerate(dxlines):
splitline = line.split()
if len(splitline) > 2:
if splitline == '1':
gridcounts.append(int(splitline))
gridcounts.append(int(splitline))
gridcounts.append(int(splitline))

# print "# gridcounts ",gridcount

if splitline == 'origin':
origin.append(float(splitline))
origin.append(float(splitline))
origin.append(float(splitline))
if splitline == 'delta':
if float(splitline) > 0:
deltas = float(splitline)
if float(splitline) > 0:
deltas = float(splitline)
if float(splitline) > 0:
deltas = float(splitline)
if splitline == '3':
numpoints = int(splitline)

# print "# Total number of gridpoints is ",numpoints

startline = i + 1
if startline > 1:
break

dxfile.close()
splittext = dxtext.split()
del splittext[0:splittext.index('follows') + 1]  # get rid of header text, last word is "follows"
floats = []
for element in splittext:
if len(element) > 0:
try:
floats.append(float(element))
except ValueError:
pass

# Assign to 3D numpy array
assert len(floats) == gridcounts*gridcounts*gridcounts
import numpy as np
distribution = np.array(floats)
distribution = np.reshape(distribution, gridcounts)

return (distribution, origin, deltas, gridcounts)

``````
``````

In :

'''
Reads 'UxDATA' (UVDATA or UUDATA) files into memory.
Returns dictionary of 3D numpy arrays, and the origin, deltas and gridcounts.
Dictionary keys refer to <speciesname>.<distributiontype>

Since these files contain many distribution types, the type
of distributions to keep should be specified. Default 'g' -> g(r)

First line should look something like:
## <version> <gridcounts> <gridspacing> <number of species>
## 20130619        64     0.5     2

'''
import numpy as np
linesplit = filelines.split()
if len(linesplit) != 5:
exit("Error! Cannot parse first line of UxDATA file")
gridcounts = [int(linesplit) for dim in range(3)]
deltas = [float(linesplit) for dim in range(3)]
numspecies = int(linesplit)

distcolumns = {}
origin = False
for line in filelines[:10]: # Should contain at least one label line
if 'g(r)' in line:
labels = line.split()
for disttype in disttypes:
for lindex, label in enumerate(labels):
if disttype in label and 'k' not in label:
distcolumns[disttype] = lindex - 1
if not origin and '#' not in line:
origin = [float(element) for element in line.split()[:3]]

dists = {}
for linenum, line in enumerate(filelines):
splitline = line.split()
if '##<' in line: #New species
speciesname = splitline
for key in distcolumns.keys():
#print key
distname = speciesname+'.'+key
dists[distname] = []
#print distname
if "#" not in line and len(splitline) > 2:
# Probably a data line
for key, colnum in distcolumns.iteritems():
distname = speciesname+'.'+key
dists[distname].append(float(splitline[colnum]))

for key, value in dists.iteritems():
assert len(value) == gridcounts*gridcounts*gridcounts
myarray = np.array(value)
dists[key] = np.reshape(myarray, gridcounts, order='F')
return dists, origin, deltas, gridcounts

``````
``````

In :

def data2Grids(uxdatafilename, disttypes=['g']):
'''
Reads a UxDATA into Grid class objects
'''
(distributions, origin, deltas, gridcounts) = readUxDATA(uxdatafilename, disttypes=disttypes)
grids = {}
for name, dist in distributions.iteritems():
grids[name] = Grid(dist, origin, gridcounts, deltas)
return grids
#dists = data2Grids('/Users/sindhikara/Dropbox/scripts/modules/Grid/grid/tests/data/UxDATAfiles/UVDATA.sample', disttypes=['g','c'])

``````
``````

In :

'''
Returns lists of 3D numpy and the origin, deltas and gridcounts.

First two lines should look something like:
16.00000000000000         16.00000000000000         16.00000000000000
32           32           32

'''
print "Reading TINKER style guv file"
print "Warning, origin set to default 0.0, 0.0, 0.0"
print "Distributions are unlabeled."
print "Assuming distributions start on line 8"

sidelengths = [float(element) for element in filelines.split()]
gridcounts = [int(element) for element in filelines.split()]
deltas = [sidelengths[dim] / gridcounts[dim] for dim in range(3)]
origin = [0.0]*3
numspecies = int(filelines)
print "Found %d species." % (numspecies)
datalines = [line.split() for line in filelines[7:]]
assert len(datalines) == gridcounts*gridcounts*gridcounts
dists = []

for specnum in range(numspecies):
dists.append(np.reshape(np.array([float(dataline[specnum]) for dataline in datalines]),
gridcounts))
return dists, origin, deltas, gridcounts

``````
``````

In :

def TKRguv2Grids(filename):
'''
Reads a TINKER guv into Grid class objects
'''

(distributions, origin, deltas, gridcounts) = readTKRguv(filename)
grids = []

for dist in distributions:
grids.append(Grid(dist, origin, gridcounts, deltas))
return grids
#grids = TKRguv2Grids('/Users/sindhikara/Dropbox/scripts/modules/Grid/grid/tests/data/TKRguv/h2o.guv.sample')
#grids.writedx('delme.dx')

``````
``````

In :

def h5ToGrids(h5filename):
'''
Reads distributions from an MDF hdf5 (h5) file into a dictionary
of Grid class objects.
Returns list of dictionaries:

gridlist[speciesindex]["distributiontype"] = GridObject

'''
import tables
h5file = tables.openFile(h5filename)
print "Warning. Box will be centered on 0,0,0 (not read)"

#Get universal parameters
gridcounts = h5file.root.parameters_uv.num_grid_points[:]
origin = [-deltas[dim]*gridcounts[dim]*0.5 for dim in range(3)]

speciesnames = h5file.root.parameters_vv.solvent.element_0.names[:]
grids = dict((speciesname, {}) for speciesname in speciesnames)
for key, value in h5file.root._v_children.iteritems():
#
#print "key (in childen) is ",key
if "uv" in key[1:] and len(key) < 5:
# This leave is an array of distributions
# key[1:] is to avoid uvv
# len(key) < 5 is to avoid /parameters_uv/
value = value[:].T
for specnum, dist in enumerate(value):
assert len(dist) == gridcounts * gridcounts * gridcounts
threeddist = np.reshape(dist,gridcounts).T
grids[speciesnames[specnum]][key] = Grid(threeddist, origin, gridcounts, deltas)
return grids

#gridlist = h5ToGrids('/Users/sindhikara/Dropbox/scripts/modules/Grid/grid/tests/data/MDF_H5/tip3p.uv.h5')

``````
``````

In :

``````
``````

In :

def printdxfrom1dzfast(
values,
origin,
delta,
gridcounts,
filename,
):
''' Print a dx file'''
print "Deprecated Function"
f = open(filename, 'w')
f.write("#DX file from Dan's program\n")
f.write('object 1 class gridpositions counts {0} {1} {2}\n'.format(gridcounts,
gridcounts, gridcounts))
f.write('origin {0} {1} {2}\n'.format(origin, origin,
origin))
f.write('delta {0} 0 0\n'.format(delta))
f.write('delta 0 {0} 0\n'.format(delta))
f.write('delta 0 0 {0}\n'.format(delta))
f.write('object 2 class gridconnections counts {0} {1} {2}\n'.format(gridcounts,
gridcounts, gridcounts))
f.write('object 3 class array type double rank 0 items {0} data follows\n'.format(gridcounts
* gridcounts * gridcounts))
for value in values:
f.write('{0}\n'.format(value))
f.write('object {0} class field\n'.format(filename))
f.close()

``````
``````

In :

def printdxfrom3d(
distribution,
origin,
delta,
gridcounts,
filename,
):
''' Print a dx file given a 3d list
This function is used by Grid objects.
'''
f = open(filename, 'w')
f.write("#DX file from Dan Sindhikara's 'grid' program.\n")
f.write('object 1 class gridpositions counts {0} {1} {2}\n'.format(gridcounts,
gridcounts, gridcounts))
f.write('origin {0} {1} {2}\n'.format(origin, origin,
origin))
f.write('delta {0} 0 0\n'.format(delta))
f.write('delta 0 {0} 0\n'.format(delta))
f.write('delta 0 0 {0}\n'.format(delta))
f.write('object 2 class gridconnections counts {0} {1} {2}\n'.format(gridcounts,
gridcounts, gridcounts))
f.write('object 3 class array type double rank 0 items {0} data follows\n'.format(gridcounts
* gridcounts * gridcounts))
for i in range(gridcounts):
for j in range(gridcounts):
for k in range(gridcounts):
f.write('{0}\n'.format(distribution[i][j][k]))
f.write('object {0} class field\n'.format(filename))
f.close()

``````
``````

In :

def printgrd(
distribution,
origin,
deltas,
gridcounts,
filename,
):
''' Print a .grd file readable by Accelrys DS given a 3d list.
Note: This is not UHBD .grd format!!
This function is used by Grid objects.
'''
boxdims = [(gridcounts[dim]-1)*deltas[dim] for dim in range(3)]
f = open(filename, 'w')
f.write("Volumetric data written by grid.py.\n")
f.write('(1p,e12.5)\n')
f.write('%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n' % (boxdims,
boxdims, boxdims,
90.0, 90.0, 90.0))
f.write('%5d %5d %5d\n' % (gridcounts-1, gridcounts-1, gridcounts-1))

#pleft and pright are number of grid points left or right of 0,0,0
pleft = [ origin[dim] / deltas[dim] for dim in range(3)]
pright = [gridcounts[dim] + pleft[dim] - 2 for dim in range(3)]

# First number (1 or 3) represents fast variable (x or z)
f.write('    1 %6d %6d %6d %6d %6d %6d\n' % (pleft, pright, pleft, pright, pleft, pright))
for z in range(gridcounts):
for y in range(gridcounts):
for x in range(gridcounts):
f.write('%12.5E\n' % (distribution[x][y][z]))
f.close()
# Output is unreadable by VMD and Chimera. Maybe a mistake in the instructions?
#griddata = dx2Grid('/Users/sindhikara/Desktop/Data/RismMap/1L2Y/1L2Y.O.1.dx')
#griddata.writegrd('DanTest.grd')

``````
``````

In :

def getcoordfromindices(indices, origin, deltas):
'''
Returns coordinates as a length 3 list of floats.
# Optimization on June 10, 2013
'''
return [float(indices[i]) * deltas[i] + origin[i] for i in range(3)]

``````
``````

In :

def getindicesfromcoord(coord, origin, deltas):
indices = []
for i in range(3):
indices.append(int((coord[i] - origin[i]) / deltas[i] + 0.5))
return indices

``````
``````

In :

def precomputeshellindices(maxindex):
'''return a list of 3d lists containing the indices in successive search shells
i.e. 0 = [0,0,0]
1 = [[-1,-1,-1],[-1,0,-1],..
essentially how much to shift i,j,k indices from center to find grid point on shell
This will make evacuation phase faster
'''

from math import sqrt
shellindices = [[[0, 0, 0]]]
for index in range(1, maxindex):

# range

indicesinthisshell = []
for i in range(-index, index + 1):
for j in range(-index, index + 1):
for k in range(-index, index + 1):

# print "math.sqrt(i*i + j*j + k*k))=",math.sqrt(i*i + j*j + k*k),"index = ",index

if int(sqrt(i * i + j * j + k * k)) == index:  # I think this will miss some
indicesinthisshell.append((i, j, k))
shellindices.append(tuple(indicesinthisshell))
return tuple(shellindices)

``````
``````

In :

def createprecomputedindicesjson(numshells=40):
'''
stores a local file called shells.json
'''

shellindices = precomputeshellindices(numshells)
from json import dump
import os
outfile = os.path.join(os.path.dirname(__file__), "data", 'shells.json')
f = open(outfile, 'w')
dump(shellindices, f)
f.close()

``````
``````

In :

import os
infile = os.path.join(os.path.dirname(__file__),"data", 'shells.json')
f = open(infile, 'rb')
return shellindices

``````
``````

In :

def getlinearweightsandindices(origin, deltas, coord):
'''
Given one 3D coordinate, return 8 corner indices and weights
that would allow direct linear interpolation
This subroutine is separated from linearinterpolatevalue to allow
precomputation of weights and indices.

This is using the formula for Trilinear interpolation
Significantly optimized (given typical Python constraints).
'''

#output = np.empty(indices.shape)
x0i = int((coord - origin) / deltas)
y0i = int((coord - origin) / deltas)
z0i = int((coord - origin) / deltas)

x0 = x0i * deltas + origin
y0 = y0i * deltas + origin
z0 = z0i * deltas + origin

x1 = x0 + 1
y1 = y0 + 1
z1 = z0 + 1

#Check if xyz1 is beyond array boundary:
#x1[np.where(x1==input_array.shape)] = x0.max()
#y1[np.where(y1==input_array.shape)] = y0.max()
#z1[np.where(z1==input_array.shape)] = z0.max()

dx = coord - x0
dy = coord - y0
dz = coord - z0

normweights = ((1-dx)*(1-dy)*(1-dz),\
dx*(1-dy)*(1-dz),\
(1-dx)*dy*(1-dz),\
(1-dx)*(1-dy)*dz,\
dx*(1-dy)*dz,\
(1-dx)*dy*dz,\
dx*dy*(1-dz),\
dx*dy*dz)

cindices = ((x0i,y0i,z0i),
(x0i+1,y0i,z0i),
(x0i,y0i+1,z0i),
(x0i,y0i,z0i+1),
(x0i+1,y0i,z0i+1),
(x0i,y0i+1,z0i+1),
(x0i+1,y0i+1,z0i),
(x0i+1,y0i+1,z0i+1))
return (normweights, cindices)

``````
``````

In :

def linearinterpolatevalue(
distribution,
origin,
deltas,
coord,
):
'''given a 3d coordinate, using a linear interpolation from the 8 nearest gridpoints,
estimate the value at that coordinate
'''

(weights, cindices) = getlinearweightsandindices(origin, deltas,
coord)
value = 0
for (i, mycindices) in enumerate(cindices):
try:
value += distribution[mycindices] * weights[i]
except:
print 'Failed to find gridpoint at', mycindices
print 'coordinate=', coord
return False
return value

``````
``````

In :

def calcrdf(
distribution,
origin,
deltas,
coord,
deltar=0.1,
numsumgrids=20,
):
'''
Calculates the radial distribution function about a point using the 3d distribution
'''

sum = 0.0
count = 0.0
for i in range(int(2.0 * radius / delta)):
for j in range(int(2.0 * radius / delta)):
for k in range(int(2.0 * radius / delta)):
x = float(i) * delta - radius
y = float(j) * delta - radius
z = float(k) * delta - radius
thisrad = (x ** 2 + y ** 2 + z ** 2) ** 0.5
< radius + delta / 2.0:
mycoord = [center + x, center + y,
center + z]
sum += linearinterpolatevalue(distribution,
origin, deltas, mycoord)
count += 1.0
sum = sum / count
return sum

/ deltar))]
gr = []
subdelta = deltar * (float(i) + 0.5) / float(numsumgrids)

``````

Notebook Testing below

``````

In :

#if isnotebook: