Scripts: torspace.py

import chimera
from chimera import selection, replyobj
from BondRotMgr import bondRotMgr

pdbfile = chimera.openModels.open('./kk.pdb')
mol = pdbfile[0]
mybonds = ['#0:2.C@CP3,CPG','#0:2.C@CPG,OPA','#0:2.C@CPG,C01','#0:2.C@C01,C02']
myatoms = ['#0:2.C@1H0','#0:2.C@1HPK','#0:2.C@2H0','#0:2.C@2HPK','#0:2.C@3H0','#0:2.C@C01','#0:2.C@C02','#0:2.C@C03','#0:2.C@CPG','#0:2.C@CPK','#0:2.C@HPG','#0:2.C@O02','#0:2.C@O03','#0:2.C@OPA']
# Eric's alternative:

mybondlist = []
for addbond in mybonds:
    a1, a2 = selection.OSLSelection(addbond).atoms()
    bond = (set(a1.bonds)&set(a2.bonds)).pop()
    mybondlist.append(bond)

torsionlist = []
for mybond in mybondlist:
    torsionlist.append(bondRotMgr.rotationForBond(mybond))

myatomlist = []
for addatom in myatoms:
    thisatom = selection.OSLSelection(addatom).atoms()[0]
    myatomlist.append(thisatom)

mycoordlist = []
for torsion in torsionlist:
    for angle in range(0,360,1):
        torsion.increment(angle)
        for atom in myatomlist:
            mycoordlist.append(atom.coord())

gridData = {}
from math import floor
from numpy import array, float32
from _contour import affine_transform_vertices

spacing = 0.1
pts = array([a for a in mycoordlist], float32)

# add a half-voxel since volume positions are
# considered to be at the center of their voxel
from numpy import floor, zeros
pts = floor(pts/spacing + 0.5).astype(int)

for pt in pts:
    center = tuple(pt)
    gridData[center] = gridData.get(center, 0) + 1

# generate volume
axisData = zip(*tuple(gridData.keys()))
minXyz = [min(ad) for ad in axisData]
maxXyz = [max(ad) for ad in axisData]

# allow for zero-padding on both ends
dims = [maxXyz[axis] - minXyz[axis] + 10 for axis in range(3)]
print minXyz, maxXyz, dims

from numpy import zeros, transpose
volume = zeros(dims, int)
for index, val in gridData.items():
        adjIndex = tuple([index[i] - minXyz[i] + 5 for i in range(3)])
        volume[adjIndex] = val

from VolumeData import Array_Grid_Data
# the "cushion of zeros" means d-5...
gd = Array_Grid_Data(volume.transpose(), [(d-5) * spacing for d in minXyz], [spacing] * 3)
gd.name = 'volumen'

# show volume
import VolumeViewer
dataRegion = VolumeViewer.volume_from_grid_data(gd)
vd = VolumeViewer.volumedialog.volume_dialog(create=True)
vd.message("Volume can be saved from File menu")