Scripts: xyzsdmap.py

# -----------------------------------------------------------------------------
# Create a volume data set from a text file of xyz point coordinates
# and standard deviations along each axis.
#
# Example Chimera command
#
#       run xyzsdmap.py points.xyz 1.5
#
# reads text file points.xyz that contains
#
#       #  x      y       z      sx   sy   sz
#       15.749  16.032 -11.160  1.3  2.1  3.1
#       16.540  14.853 -10.693  2.11 3.2  1.4
#       17.946  15.275 -10.267  1.4  1.3  3.8
#       ...
#
# creating a volume data set with grid spacing 1.5 Angstrom.
#
from math import sqrt, pi
def xyz_map(xyz_path,
            grid_spacing,
            pad = None,                    # default is 10 times grid spacing
            cutoff_range = 5,              # in standard deviations
            ):

    xyz_list = []
    sdev_list = []
    f = open(xyz_path, 'r')
    while True:
        line = f.readline()
        if line == '':
            break
        if line[0] != '#' and line.strip():
            fields = [float(x) for x in line.split()]
            if len(fields) >= 6:
                xyz_list.append(fields[0:3])
                sdev_list.append(fields[3:6])
    f.close()
    from numpy import array, float32, zeros
    xyz = array(xyz_list, float32)
    sdevs = array(sdev_list, float32)

    from _multiscale import bounding_box
    xyz_min, xyz_max = bounding_box(xyz)
    if pad is None:
        pad = 10*grid_spacing
    origin = [x-pad for x in xyz_min]
    step = grid_spacing
    ijk = (xyz - origin) / step
    from math import pow, pi, ceil
    normalization = pow(2*pi,-1.5)*pow(step,-3)
    weights = normalization/(sdevs[:,0]*sdevs[:,1]*sdevs[:,2])
    shape = [int(ceil((xyz_max[a] - xyz_min[a] + 2*pad) / step))
             for a in (2,1,0)]
    matrix = zeros(shape, float32)

    from _gaussian import sum_of_gaussians
    sum_of_gaussians(ijk, weights, sdevs, cutoff_range, matrix)
    matrix *= normalization

    from os.path import basename
    name = basename(xyz_path)
    from VolumeData import Array_Grid_Data
    grid = Array_Grid_Data(matrix, origin, (step,step,step), name = name)

    from VolumeViewer import volume_from_grid_data
    v = volume_from_grid_data(grid)

    return v

# -----------------------------------------------------------------------------
#
if 'arguments' in globals():
    args = globals()['arguments']
    args = list(args[:1]) + [float (a) for a in args[1:]]
    xyz_map(*args)