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()