# [Chimera-users] Re: Gaussians

Thomas Goddard goddard at cgl.ucsf.edu
Wed Dec 10 23:36:11 PST 2003

```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):

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

# -----------------------------------------------------------------------------
#

f = open(path, 'r')
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:
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()

```