[Chimera-users] Re: Gaussians

[Thomas Goddard]

Hi Lars,

  I'm not sure what you mean by visualizing some Gaussians.  If you
mean you want to make a volume which contains some specified Gaussians
then I can give you a separate Python program that will do that.
I already wrote it for another user.  It is included at the end of
this email.  It outputs a netcdf format volume file.

     Tom

% ./spots.py
Make a volume array that is the sum of specified Gaussians.
The first line of the input file has the grid cell size.
Each subsequent line has 5 fields giving x, y, z, width, height.
A volume grid is made that includes each Gaussian out to 3 widths.
The resulting volume is written in netcdf format.
Syntax: spots.py <infile> <outfile>

% cat spotfile
.5
15.49100  7.05000 17.16500   0.38   1.00
14.94600  7.55000 15.88800   0.35   1.00
16.01100  7.60400 14.81400   0.34   1.00
...

% ./spots.py spotfile gspottest.nc

The spots.py code follows.  Note that it contains a path to your
Chimera installation because it needs Chimera modules to write the
netcdf volume.


#!/usr/bin/env python
# -----------------------------------------------------------------------------
# Make a volume array that is the sum of specified Gaussians.
# The resulting volume is written in netcdf format.
#
# Syntax: spots.py <infile> <outfile>
#
# The input file has a grid cell size on the first line followed by one
# line per gaussian with 5 fields specifying x, y, z, width, height.
#
import os.path
import sys

# Need Chimera path to use VolumeData and Numeric modules
CHIMERA_PATH = '/usr/local/chimera'

sys.path.insert(0, os.path.join(CHIMERA_PATH, 'lib'))
sys.path.insert(0, os.path.join(CHIMERA_PATH, 'share'))

import Numeric

# -----------------------------------------------------------------------------
#
def make_spot_grid(path):

    step_size, spot_parameters = read_spots(path)
    
    range = 3         # in spot widths
    xyz_min, xyz_max = bounding_box(spot_parameters, range)
    xyz_origin = xyz_min
    xyz_step = (step_size, step_size, step_size)
    xyz_size = map(lambda a,b: a - b, xyz_max, xyz_min)
    import math
    grid_size = map(lambda s, step, math=math: 1 + int(math.ceil(s/step)),
                    xyz_size, xyz_step)

    a = spot_array(spot_parameters, range, xyz_origin, step_size, grid_size)
    from VolumeData import Array_Grid_Data
    g = Array_Grid_Data(a, xyz_origin, xyz_step)

    return g

# -----------------------------------------------------------------------------
#
def read_spots(path):

    f = open(path, 'r')
    lines = f.readlines()
    f.close()

    step_size = float(lines[0])

    spot_params = []
    for line in lines[1:]:
        params = map(float, line.split())
        spot_params.append(params)

    return step_size, spot_params

# -----------------------------------------------------------------------------
#
def bounding_box(spot_parameters, range):

    xmin = ymin = zmin = None
    xmax = ymax = zmax = None
    for x, y, z, w, h in spot_parameters:
        pad = range * w
        if xmin == None or x - pad < xmin: xmin = x - pad
        if ymin == None or y - pad < ymin: ymin = y - pad
        if zmin == None or z - pad < zmin: zmin = z - pad
        if xmax == None or x + pad > xmax: xmax = x + pad
        if ymax == None or y + pad > ymax: ymax = y + pad
        if zmax == None or z + pad > zmax: zmax = z + pad

    return (xmin, ymin, zmin), (xmax, ymax, zmax)

# -----------------------------------------------------------------------------
#
def spot_array(spot_parameters, range, xyz_origin, step_size, grid_size):

    shape = list(grid_size)
    shape.reverse()             # array indices are z, y, x.
    array = Numeric.zeros(shape, Numeric.Float)
    for x, y, z, width, height in spot_parameters:
        add_spot((x,y,z), width, height, range, xyz_origin, step_size, array)
    return array

# -----------------------------------------------------------------------------
#
def add_spot(xyz, width, height, srange, xyz_origin, step_size, array):

    ijk = map(lambda a,b,s=step_size: (a - b)/s, xyz, xyz_origin)
    w = width / step_size
    w2 = w * w
    r = srange * width / step_size
    import math
    imin, jmin, kmin = map(lambda l, r=r, math=math: int(math.floor(l-r)), ijk)
    imax, jmax, kmax = map(lambda l, r=r, math=math: int(math.ceil(l+r)), ijk)
    for k in range(kmin, kmax+1):
        for j in range(jmin, jmax+1):
            for i in range(imin, imax+1):
                di, dj, dk = i-ijk[0], j-ijk[1], k-ijk[2]
                s2 = di*di + dj*dj + dk*dk
                e = height * math.exp(-s2/(2*w2))
                array[k,j,i] += e

# -----------------------------------------------------------------------------
#
def bfactor_backbone_spots(path, molecule = None):

    if molecule == None:
        import chimera
        molecule = chimera.openModels.list()[0]
    out = open(path, 'w')
    out.write('.5\n')     # step size
    import math
    for a in molecule.atoms:
        if a.name in ('C', 'N', 'CA') and hasattr(a, 'bfactor'):
            x, y, z = a.coord().xyz.data()
            w = math.sqrt(a.bfactor / (8 * math.pi**2))
            h = 1
            out.write('%8.5f %8.5f %8.5f %6.2f %6.2f\n' % (x,y,z,w,h))
    out.close()

# -----------------------------------------------------------------------------
#
def syntax():
  
  msg = ('Make a volume array that is the sum of specified Gaussians.\n' +
         'The first line of the input file has the grid cell size.\n' +
         'Each subsequent line has 5 fields giving x, y, z, width, height.\n' +
         'A volume grid is made that includes each Gaussian out to 3 widths.\n' +
         'The resulting volume is written in netcdf format.\n' +
         'Syntax: spots.py <infile> <outfile>\n')
  sys.stderr.write(msg)
  sys.exit(1)

# -----------------------------------------------------------------------------
#
def run_command():

    if len(sys.argv) != 3:
        syntax()

    inpath = sys.argv[1]
    outpath = sys.argv[2]

    grid = make_spot_grid(inpath)

    from VolumeData.netcdf.netcdf_grid import write_grid_as_netcdf
    write_grid_as_netcdf(grid, outpath)

# -----------------------------------------------------------------------------
#
if __name__ == '__main__':
    run_command()