Scripts: xyzmap.py

# -----------------------------------------------------------------------------
# Create a volume data set from a text file of xyz point coordinates.
#
# Example Chimera command
#
#       run xyzmap.py points.xyz 5 1
#
# reads text file points.xyz that contains
#
#       # coords from 1fav
#       15.749  16.032 -11.160
#       16.540  14.853 -10.693
#       17.946  15.275 -10.267
#       ...
#
# creating a volume data set with resolution (Gaussian width) 5 Angstroms
# and grid spacing 1 Angstrom.
#
from math import sqrt, pi
def xyz_map(xyz_path,
            resolution,
            grid_spacing = None,           # default is 1/3 resolution
            pad = None,                    # default is 3 times resolution
            cutoff_range = 5,              # in standard deviations
            sigma_factor = 1/(pi*sqrt(2))  # resolution / standard deviation
            ):

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

    from _multiscale import bounding_box
    xyz_min, xyz_max = bounding_box(xyz)
    if pad is None:
        pad = 3*resolution
    origin = [x-pad for x in xyz_min]
    if grid_spacing is None:
        step = resolution / 3.0
    else:
        step = grid_spacing
    ijk = (xyz - origin) / step
    weights = [1.0] * len(xyz)
    sdev = sigma_factor * resolution / step
    from numpy import zeros, float32
    sdevs = zeros((len(xyz),3), float32)
    sdevs[:] = sdev
    from math import pow, pi, ceil
    normalization = pow(2*pi,-1.5)*pow(sdev*step,-3)
    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)