Scripts: mapsum.py

# Report average density around CA atoms.
def report_atom_densities(atom_names = ('CA',), radii = (3,), grid_spacing = 0.3):

    from chimera import openModels as om, Molecule
    import VolumeViewer
    mols = om.list(modelTypes = [Molecule])
    maps = VolumeViewer.volume_list()
    grids = [grid_points(r, grid_spacing) for r in radii]
    for mol in mols:
        print 'Molecule %s, ' % mol.name,
        for map in maps:
            atoms = [a for a in mol.atoms if a.name in atom_names]
            rtext = ' '.join('%.4g'%r for r in radii)
            ngtext = ' '. join('%d'%len(g) for g in grids)
            print 'average density near %d atoms in map %s, ' % (len(atoms), map.name)
            print 'radii %s, sampling grid spacing %.4g, grid points %s' % (rtext, grid_spacing, ngtext)
            aves = [atom_densities(atoms, g, map) for g in grids]
            for i,a in enumerate(atoms):
                rid = a.residue.id
                atext = '   '.join('%8.3g'%r_ave[i] for r_ave in aves)
                print '%4d %4s %4s %s' % (rid.position, rid.chainId, a.name, atext)

# Average density around atoms.
def atom_densities(atoms, sampling_grid, map):

    centers = atom_coordinates(atoms, map.openState.xform)
    aves = average_density(centers, sampling_grid, map)
    return aves

# Compute atom coordinates in map coordinate system
def atom_coordinates(atoms, map_xf):

    xf = map_xf.inverse()
    axyz = [xf.apply(a.xformCoord()).data() for a in atoms]
    return axyz

# Compute average density inside sphere using grid centered on sphere.
def average_density(centers, sampling_grid, map):

    n = len(sampling_grid)
    aves = []
    for c in centers:
        v = map.interpolated_values(sampling_grid + c)
        aves.append(float(v.sum())/n)
    return aves

# Produce array of grid points inside sphere centered at 0,0,0.
def grid_points(radius, grid_spacing):

    g = []
    n = int(float(radius + grid_spacing)/grid_spacing)
    r2 = radius*radius
    for i in range(-n,n+1):
        for j in range(-n,n+1):
            for k in range(-n,n+1):
                dx,dy,dz = i*grid_spacing,j*grid_spacing,k*grid_spacing
                if dx*dx + dy*dy + dz*dz < r2:
                    g.append((dx,dy,dz))
    import numpy
    ga = numpy.array(g)
    return ga

# Report densities for open molecule and map.
report_atom_densities(('CA','CB'), (3,5,7))
from Accelerators import standard_accelerators as sa
sa.show_reply_log()