This is a read-only mirror of pymolwiki.org
Difference between revisions of "Consistent View/ Map Inspect"
Jump to navigation
Jump to search
m |
|||
| Line 1,116: | Line 1,116: | ||
R_prime = (kr.inverse() * std_view.rotmat) | R_prime = (kr.inverse() * std_view.rotmat) | ||
| − | print "delta",std_view.origin_of_rotation - residue.CA | + | if verbose: print "delta",std_view.origin_of_rotation - residue.CA |
test = residue.CA + (kr.inverse() * std_view.difference_to_this_CA) | test = residue.CA + (kr.inverse() * std_view.difference_to_this_CA) | ||
Revision as of 00:10, 9 April 2010
DESCRIPTION
This is a toolkit for rapidly inspecting multiple maps and models. The right & left arrow keys navigate sequentially through the amino acid chain, but alternating between two structures (could be NCS-related chains or models solved in different space groups). Each view is rendered in a consistent orientation (the default is centered on alpha carbon, with nitrogen right, carbonyl left & beta carbon up). The view can be customized. It is necessary to edit the script to define the behavior for the problem at hand.
INSTALLATION
Three components are needed:
- A Python module, support.py, giving matrix algebra & Kabsch-style structure alignment (source code given below). This module must be in the current working directory, and PyMol started from the command line, e.g., on MacOSX:
/Applications/MacPyMol.app/Contents/MacOS/MacPyMOL
- Atomic coordinates must be loaded first, and if appropriate, maps and isomesh levels. An example .pml file is shown here, consisting of two structural models solved in different space groups:
load sand2b_001_2mFo-DFc_map.xplor, map2
isomesh msh2,map2,2.0
color blue, msh2
load sand3_rotfix_001_2mFo-DFc_map.xplor, map3
isomesh msh3,map3,2.0
color marine, msh3
load sand2b_001.pdb, mod02
load sand3_rotfix_001.pdb, mod03
set antialias, 2
set line_width, 3
- A final pml is then loaded (@consistent.pml) to define the viewing sequence for the problem at hand. An example is shown here. (The python / python end construct will not work in PyMOL versions below 1.0.)
python
from support import matrix,residue,std_residue,kabsch_rotation
from support import view_matrix, new_set_view
# This section is the generic preset view--DO NOT EDIT #########################
# Centered on CA, CB up, N left, C=O right
preset_view = view_matrix((
0.295077115, 0.440619737, -0.847830832,
-0.408159286, 0.860427082, 0.305112839,
0.863915384, 0.256014407, 0.433732271,
0., 0., -30.,
0., 0., 0.,
26., 33., 0.
))
preset_view.difference_to_this_CA = matrix.col((0.,0.,0.,))
preset_std_residue = std_residue()
preset_std_residue.N = matrix.col((16.1469993591, 33.0660018921, -20.0760002136))
preset_std_residue.CA = matrix.col((15.4309997559, 34.2599983215, -20.4640007019))
preset_std_residue.C = matrix.col((15.2889995575, 34.1189994812, -21.9699993134))
preset_std_residue.CB = matrix.col((16.2469997406, 35.516998291, -20.1380004883))
# END GENERIC VIEW #############################################################
def view_preset():
cmd.current_std_view = preset_view
cmd.current_std_residue = preset_std_residue
view_preset()
# EDIT THIS SECTION TO REFLECT THE PROBLEM AT HAND #############################
# In this example, there are two objects: structure mod02 & mod03
# We will be looking at chain B in each structure, each having 262 residues,
# for a total of 524 amino acids to view.
# For each of 524 views, we select the appropriate map and isomesh
cmd.no_of_structures = 2
cmd.global_ptr = 1 * cmd.no_of_structures
cmd.last_residue = preset_std_residue
def seqwalk():
chain_id = "B"
if cmd.global_ptr<524:
resid = cmd.global_ptr/2
this_object = cmd.global_ptr%2
if this_object==0:
atoms = cmd.get_model("object mod03 and chain %s"%chain_id)
obj_id = "mod03"
map = "map3"
mesh = "msh3"
else:
atoms = cmd.get_model("object mod02 and chain %s"%chain_id)
obj_id = "mod02"
map = "map2"
mesh = "msh2"
this = residue(chain_id,resid,atoms)
cmd.last_residue = this #remember this residue in case reset command is issued
print "Centering on /%s//%s/%s %d/CA; waiting for right mouse key"%(
obj_id,chain_id,this.residue_type,resid,)
cmd.center("/%s//%s/%d/CA"%(obj_id,chain_id,resid))
cmd.isomesh(
name=mesh,map=map,level= 2.0, selection="object %s and chain %s and resid %d"%(obj_id,chain_id,resid),buffer=8)
# No more edits below this line
kr = kabsch_rotation(cmd.current_std_residue.sites(),this.sites())
view_matrix = new_set_view(cmd.current_std_view,kr,this,verbose=False).view_matrix()
cmd.set_view(view_matrix)
cmd.global_ptr=cmd.global_ptr+1
# END OF PROBLEM-SPECIFIC SECTION TO EDIT #####################################
def seqwalk2():
cmd.global_ptr-=2
seqwalk()
cmd.set_key("right",seqwalk)
cmd.set_key("left",seqwalk2)
def view_reset():
cmd.current_std_view = view_matrix(cmd.get_view()).set_diff_to_CA(cmd.last_residue)
cmd.current_std_residue = cmd.last_residue
def view_goto(number):
cmd.global_ptr = cmd.no_of_structures*number
cmd.extend("view_reset",view_reset) # After customizing the view with GUI mouse & clicks, make the
# view persistent by typing "view reset"
cmd.extend("view_preset",view_preset)# Forget the user-customized view, go back to the generic view
# defined at the beginning of the program
cmd.extend("view_goto",view_goto) # Example: view_goto(36) --resets at residue 36
python end
USAGE
Once the code is set up as described above, you must mouse-click in PyMOL's 3D display window to enable the right & left arrow keys. These keys will navigate through the amino acid sequence, displaying each residue in a predefined orientation. Maps, if defined, will be redrawn for each view.
Three commands can be issued at the prompt:
- view_goto(N) -- will go to residue N. The view will not actually change until you mouse click in the 3D display window and hit the right arrow key.
- view_reset -- will accept the current view (in relation to the present alpha carbon) as the default. This is extremely useful for preparing views for presentation, which compare sidechains or ligands in two homologous structures. For example, if you are interested in a Mg++ ion adjacent to residue 56, you will first use the arrow keys to center on the residue 56 alpha carbon. Then recenter on the Mg++ ion, and rotate to the exact view desired. Typing view_reset will then produce the same view of the Mg++ ion in both structures.
- view_preset -- revert back to the generic view centered on the alpha carbon.
METHODOLOGY
Four atoms of the current residue (N, CA, CB, and C) are used for a Kabsch-style alignment against a reference orientation. For glycines, a hypothetical beta-carbon is modelled (but not shown) based on tetrahedral geometry. A known limitation of the approach is that the alignment is very local, i.e., neighboring residues may not precisely align between structures.
Adaptation of PyMOL's set_view command to display aligned views of separate structures was suggested by Herb Klei.
SUPPORTING MODULE
The following Python module, support.py, must be placed in the current working directory. Code is based on CCTBX (cctbx.sf.net) and is released under the cctbx license.
import math
#####################################################################
# This code has been adapted from cctbx.sf.net, file scitbx/matrix.py
flex = None
numpy = None
class rec(object):
container_type = tuple
def __init__(self, elems, n):
assert len(n) == 2
if (not isinstance(elems, self.container_type)):
elems = self.container_type(elems)
assert len(elems) == n[0] * n[1]
self.elems = elems
self.n = tuple(n)
def n_rows(self):
return self.n[0]
def n_columns(self):
return self.n[1]
def __neg__(self):
return rec([-e for e in self.elems], self.n)
def __add__(self, other):
assert self.n == other.n
a = self.elems
b = other.elems
return rec([a[i] + b[i] for i in xrange(len(a))], self.n)
def __sub__(self, other):
assert self.n == other.n
a = self.elems
b = other.elems
return rec([a[i] - b[i] for i in xrange(len(a))], self.n)
def __mul__(self, other):
if (not hasattr(other, "elems")):
if (not isinstance(other, (list, tuple))):
return rec([x * other for x in self.elems], self.n)
other = col(other)
a = self.elems
ar = self.n_rows()
ac = self.n_columns()
b = other.elems
if (other.n_rows() != ac):
raise RuntimeError(
"Incompatible matrices:\n"
" self.n: %s\n"
" other.n: %s" % (str(self.n), str(other.n)))
bc = other.n_columns()
if (ac == 0):
# Roy Featherstone, Springer, New York, 2007, p. 53 footnote
return rec((0,)*(ar*bc), (ar,bc))
result = []
for i in xrange(ar):
for k in xrange(bc):
s = 0
for j in xrange(ac):
s += a[i * ac + j] * b[j * bc + k]
result.append(s)
if (ar == bc):
return sqr(result)
return rec(result, (ar, bc))
def __rmul__(self, other):
"scalar * matrix"
if (isinstance(other, rec)): # work around odd Python 2.2 feature
return other.__mul__(self)
return self * other
def transpose_multiply(self, other=None):
a = self.elems
ar = self.n_rows()
ac = self.n_columns()
if (other is None):
result = [0] * (ac * ac)
jac = 0
for j in xrange(ar):
ik = 0
for i in xrange(ac):
for k in xrange(ac):
result[ik] += a[jac + i] * a[jac + k]
ik += 1
jac += ac
return sqr(result)
b = other.elems
assert other.n_rows() == ar, "Incompatible matrices."
bc = other.n_columns()
result = [0] * (ac * bc)
jac = 0
jbc = 0
for j in xrange(ar):
ik = 0
for i in xrange(ac):
for k in xrange(bc):
result[ik] += a[jac + i] * b[jbc + k]
ik += 1
jac += ac
jbc += bc
if (ac == bc):
return sqr(result)
return rec(result, (ac, bc))
def __div__(self, other):
return rec([e/other for e in self.elems], self.n)
def __truediv__(self, other):
return rec([e/other for e in self.elems], self.n)
def __floordiv__(self, other):
return rec([e//other for e in self.elems], self.n)
def __mod__(self, other):
return rec([ e % other for e in self.elems], self.n)
def __call__(self, ir, ic):
return self.elems[ir * self.n_columns() + ic]
def __len__(self):
return len(self.elems)
def __getitem__(self, i):
return self.elems[i]
def as_float(self):
return rec([float(e) for e in self.elems], self.n)
def as_int(self, rounding=True):
if rounding:
return rec([int(round(e)) for e in self.elems], self.n)
else:
return rec([int(e) for e in self.elems], self.n)
def each_abs(self):
return rec([abs(e) for e in self.elems], self.n)
def min(self):
result = None
for e in self.elems:
if (result is None or result > e):
result = e
return result
def max(self):
result = None
for e in self.elems:
if (result is None or result < e):
result = e
return result
def min_index(self):
result = None
for i in xrange(len(self.elems)):
if (result is None or self.elems[result] > self.elems[i]):
result = i
return result
def max_index(self):
result = None
for i in xrange(len(self.elems)):
if (result is None or self.elems[result] < self.elems[i]):
result = i
return result
def sum(self):
result = 0
for e in self.elems:
result += e
return result
def product(self):
result = 1
for e in self.elems:
result *= e
return result
def trace(self):
assert self.n_rows() == self.n_columns()
n = self.n_rows()
result = 0
for i in xrange(n):
result += self.elems[i*n+i]
return result
def norm_sq(self):
result = 0
for e in self.elems:
result += e*e
return result
def round(self, digits):
return rec([ round(x, digits) for x in self.elems ], self.n)
def __abs__(self):
assert self.n_rows() == 1 or self.n_columns() == 1
return math.sqrt(self.norm_sq())
length_sq = norm_sq # for compatibility with scitbx/vec3.h
length = __abs__
def normalize(self):
return self / abs(self)
def dot(self, other=None):
result = 0
a = self.elems
if (other is None):
for i in xrange(len(a)):
v = a[i]
result += v * v
else:
assert len(self.elems) == len(other.elems)
b = other.elems
for i in xrange(len(a)):
result += a[i] * b[i]
return result
def cross(self, other):
assert self.n in ((3,1), (1,3))
assert self.n == other.n
a = self.elems
b = other.elems
return rec((
a[1] * b[2] - b[1] * a[2],
a[2] * b[0] - b[2] * a[0],
a[0] * b[1] - b[0] * a[1]), self.n)
def is_r3_rotation_matrix_rms(self):
if (self.n != (3,3)): raise RuntimeError("Not a 3x3 matrix.")
rtr = self.transpose_multiply()
return (rtr - identity(n=3)).norm_sq()**0.5
def is_r3_rotation_matrix(self, rms_tolerance=1e-8):
return self.is_r3_rotation_matrix_rms() < rms_tolerance
def unit_quaternion_as_r3_rotation_matrix(self):
assert self.n in [(1,4), (4,1)]
q0,q1,q2,q3 = self.elems
return sqr((
2*(q0*q0+q1*q1)-1, 2*(q1*q2-q0*q3), 2*(q1*q3+q0*q2),
2*(q1*q2+q0*q3), 2*(q0*q0+q2*q2)-1, 2*(q2*q3-q0*q1),
2*(q1*q3-q0*q2), 2*(q2*q3+q0*q1), 2*(q0*q0+q3*q3)-1))
def r3_rotation_matrix_as_unit_quaternion(self):
# Based on work by:
# Shepperd (1978), J. Guidance and Control, 1, 223-224.
# Sam Buss, http://math.ucsd.edu/~sbuss/MathCG
# Robert Hanson, jmol/Jmol/src/org/jmol/util/Quaternion.java
if (self.n != (3,3)): raise RuntimeError("Not a 3x3 matrix.")
m00,m01,m02,m10,m11,m12,m20,m21,m22 = self.elems
trace = m00 + m11 + m22
if (trace >= 0.5):
w = (1 + trace)**0.5
d = w + w
w *= 0.5
x = (m21 - m12) / d
y = (m02 - m20) / d
z = (m10 - m01) / d
else:
if (m00 > m11):
if (m00 > m22): mx = 0
else: mx = 2
elif (m11 > m22): mx = 1
else: mx = 2
invalid_cutoff = 0.8 # not critical; true value is closer to 0.83
invalid_message = "Not a r3_rotation matrix."
if (mx == 0):
x_sq = 1 + m00 - m11 - m22
if (x_sq < invalid_cutoff): raise RuntimeError(invalid_message)
x = x_sq**0.5
d = x + x
x *= 0.5
w = (m21 - m12) / d
y = (m10 + m01) / d
z = (m20 + m02) / d
elif (mx == 1):
y_sq = 1 + m11 - m00 - m22
if (y_sq < invalid_cutoff): raise RuntimeError(invalid_message)
y = y_sq**0.5
d = y + y
y *= 0.5
w = (m02 - m20) / d
x = (m10 + m01) / d
z = (m21 + m12) / d
else:
z_sq = 1 + m22 - m00 - m11
if (z_sq < invalid_cutoff): raise RuntimeError(invalid_message)
z = z_sq**0.5
d = z + z
z *= 0.5
w = (m10 - m01) / d
x = (m20 + m02) / d
y = (m21 + m12) / d
return col((w, x, y, z))
def unit_quaternion_product(self, other):
assert self.n in [(1,4), (4,1)]
assert other.n in [(1,4), (4,1)]
q0,q1,q2,q3 = self.elems
o0,o1,o2,o3 = other.elems
return col((
q0*o0 - q1*o1 - q2*o2 - q3*o3,
q0*o1 + q1*o0 + q2*o3 - q3*o2,
q0*o2 - q1*o3 + q2*o0 + q3*o1,
q0*o3 + q1*o2 - q2*o1 + q3*o0))
def axis_and_angle_as_unit_quaternion(self, angle, deg=False):
assert self.n in ((3,1), (1,3))
if (deg): angle *= math.pi/180
h = angle * 0.5
c, s = math.cos(h), math.sin(h)
u,v,w = self.normalize().elems
return col((c, u*s, v*s, w*s))
def axis_and_angle_as_r3_rotation_matrix(self, angle, deg=False):
uq = self.axis_and_angle_as_unit_quaternion(angle=angle, deg=deg)
return uq.unit_quaternion_as_r3_rotation_matrix()
def rt_for_rotation_around_axis_through(self, point, angle, deg=False):
assert self.n in ((3,1), (1,3))
assert point.n in ((3,1), (1,3))
r = (point - self).axis_and_angle_as_r3_rotation_matrix(
angle=angle, deg=deg)
return rt((r, self-r*self))
def ortho(self):
assert self.n in ((3,1), (1,3))
x, y, z = self.elems
a, b, c = abs(x), abs(y), abs(z)
if c <= a and c <= b:
return col((-y, x, 0))
if b <= a and b <= c:
return col((-z, 0, x))
return col((0, -z, y))
def rotate_around_origin(self, axis, angle, deg=False):
assert self.n in ((3,1), (1,3))
assert axis.n == self.n
if deg: angle *= math.pi/180
n = axis.normalize()
x = self
c, s = math.cos(angle), math.sin(angle)
return x*c + n*n.dot(x)*(1-c) + n.cross(x)*s
def rotate(self, axis, angle, deg=False):
import warnings
warnings.warn(
message=
"The .rotate() method has been renamed to .rotate_around_origin()"
" for clarity. Please update the code calling this method.",
category=DeprecationWarning,
stacklevel=2)
return self.rotate_around_origin(axis=axis, angle=angle, deg=deg)
def vector_to_001_rotation(self,
sin_angle_is_zero_threshold=1.e-10,
is_normal_vector_threshold=1.e-10):
assert self.n in ((3,1), (1,3))
x,y,c = self.elems
xxyy = x*x + y*y
if (abs(xxyy + c*c - 1) > is_normal_vector_threshold):
raise RuntimeError("self is not a normal vector.")
s = (xxyy)**0.5
if (s < sin_angle_is_zero_threshold):
if (c > 0):
return sqr((1,0,0,0,1,0,0,0,1))
return sqr((1,0,0,0,-1,0,0,0,-1))
us = y
vs = -x
u = us / s
v = vs / s
oc = 1-c
return sqr((c + u*u*oc, u*v*oc, vs, u*v*oc, c + v*v*oc, -us, -vs, us, c))
def outer_product(self, other=None):
if (other is None): other = self
assert self.n[0] == 1 or self.n[1] == 1
assert other.n[0] == 1 or other.n[1] == 1
result = []
for a in self.elems:
for b in other.elems:
result.append(a*b)
return rec(result, (len(self.elems), len(other.elems)))
def cos_angle(self, other, value_if_undefined=None):
self_norm_sq = self.norm_sq()
if (self_norm_sq == 0): return value_if_undefined
other_norm_sq = other.norm_sq()
if (other_norm_sq == 0): return value_if_undefined
d = self_norm_sq * other_norm_sq
if (d == 0): return value_if_undefined
return self.dot(other) / math.sqrt(d)
def angle(self, other, value_if_undefined=None, deg=False):
cos_angle = self.cos_angle(other=other)
if (cos_angle is None): return value_if_undefined
result = math.acos(max(-1,min(1,cos_angle)))
if (deg): result *= 180/math.pi
return result
def accute_angle(self, other, value_if_undefined=None, deg=False):
cos_angle = self.cos_angle(other=other)
if (cos_angle is None): return value_if_undefined
if (cos_angle < 0): cos_angle *= -1
result = math.acos(min(1,cos_angle))
if (deg): result *= 180/math.pi
return result
def is_square(self):
return self.n[0] == self.n[1]
def determinant(self):
assert self.is_square()
m = self.elems
n = self.n[0]
if (n == 1):
return m[0]
if (n == 2):
return m[0]*m[3] - m[1]*m[2]
if (n == 3):
return m[0] * (m[4] * m[8] - m[5] * m[7]) \
- m[1] * (m[3] * m[8] - m[5] * m[6]) \
+ m[2] * (m[3] * m[7] - m[4] * m[6])
if (flex is not None):
m = flex.double(m)
m.resize(flex.grid(self.n))
return m.matrix_determinant_via_lu()
return determinant_via_lu(m=self)
def co_factor_matrix_transposed(self):
n = self.n
if (n == (0,0)):
return rec(elems=(), n=n)
if (n == (1,1)):
return rec(elems=(1,), n=n)
m = self.elems
if (n == (2,2)):
return rec(elems=(m[3], -m[1], -m[2], m[0]), n=n)
if (n == (3,3)):
return rec(elems=(
m[4] * m[8] - m[5] * m[7],
-m[1] * m[8] + m[2] * m[7],
m[1] * m[5] - m[2] * m[4],
-m[3] * m[8] + m[5] * m[6],
m[0] * m[8] - m[2] * m[6],
-m[0] * m[5] + m[2] * m[3],
m[3] * m[7] - m[4] * m[6],
-m[0] * m[7] + m[1] * m[6],
m[0] * m[4] - m[1] * m[3]), n=n)
assert self.is_square()
raise RuntimeError("Not implemented.")
def inverse(self):
assert self.is_square()
n = self.n
if (n[0] < 4):
determinant = self.determinant()
assert determinant != 0
return self.co_factor_matrix_transposed() / determinant
if (flex is not None):
m = flex.double(self.elems)
m.resize(flex.grid(n))
m.matrix_inversion_in_place()
return rec(elems=m, n=n)
if (numpy is not None):
m = numpy.asarray(self.elems)
m.shape = n
m = numpy.ravel(numpy.linalg.inv(m))
return rec(elems=m, n=n)
return inverse_via_lu(m=self)
def transpose(self):
elems = []
for j in xrange(self.n_columns()):
for i in xrange(self.n_rows()):
elems.append(self(i,j))
return rec(elems, (self.n_columns(), self.n_rows()))
def _mathematica_or_matlab_form(self,
outer_open, outer_close,
inner_open, inner_close, inner_close_follow,
label,
one_row_per_line,
format,
prefix):
nr = self.n_rows()
nc = self.n_columns()
s = prefix
indent = prefix
if (label):
s += label + "="
indent += " " * (len(label) + 1)
s += outer_open
if (nc != 0):
for ir in xrange(nr):
s += inner_open
for ic in xrange(nc):
if (format is None):
s += str(self(ir, ic))
else:
s += format % self(ir, ic)
if (ic+1 != nc): s += ", "
elif (ir+1 != nr or len(inner_open) != 0): s += inner_close
if (ir+1 != nr):
s += inner_close_follow
if (one_row_per_line):
s += "\n"
s += indent
s += " "
return s + outer_close
def mathematica_form(self,
label="",
one_row_per_line=False,
format=None,
prefix="",
matrix_form=False):
result = self._mathematica_or_matlab_form(
outer_open="{", outer_close="}",
inner_open="{", inner_close="}", inner_close_follow=",",
label=label,
one_row_per_line=one_row_per_line,
format=format,
prefix=prefix)
if matrix_form: result += "//MatrixForm"
result = result.replace('e', '*^')
return result
def matlab_form(self,
label="",
one_row_per_line=False,
format=None,
prefix=""):
return self._mathematica_or_matlab_form(
outer_open="[", outer_close="]",
inner_open="", inner_close=";", inner_close_follow="",
label=label,
one_row_per_line=one_row_per_line,
format=format,
prefix=prefix)
def __repr__(self):
n0, n1 = self.n
e = self.elems
if (len(e) <= 3):
e = str(e)
else:
e = "(%s, ..., %s)" % (str(e[0]), str(e[-1]))
return "matrix.rec(elems=%s, n=(%d,%d))" % (e, n0, n1)
def __str__(self):
return self.mathematica_form(one_row_per_line=True)
def as_list_of_lists(self):
result = []
nr,nc = self.n
for ir in xrange(nr):
result.append(list(self.elems[ir*nc:(ir+1)*nc]))
return result
def as_sym_mat3(self):
assert self.n == (3,3)
m = self.elems
return (m[0],m[4],m[8],
(m[1]+m[3])/2.,
(m[2]+m[6])/2.,
(m[5]+m[7])/2.)
def as_mat3(self):
assert self.n == (3,3)
return self.elems
def as_flex_double_matrix(self):
assert flex is not None
result = flex.double(self.elems)
result.reshape(flex.grid(self.n))
return result
def as_flex_int_matrix(self):
assert flex is not None
result = flex.int(self.elems)
result.reshape(flex.grid(self.n))
return result
def extract_block(self, stop, start=(0,0), step=(1,1)):
assert 0 <= stop[0] <= self.n[0]
assert 0 <= stop[1] <= self.n[1]
i_rows = range(start[0], stop[0], step[0])
i_colums = range(start[1], stop[1], step[1])
result = []
for ir in i_rows:
for ic in i_colums:
result.append(self(ir,ic))
return rec(result, (len(i_rows),len(i_colums)))
def __eq__(self, other):
if self is other: return True
if other is None: return False
if issubclass(type(other), rec):
return self.elems == other.elems
for ir in xrange(self.n_rows()):
for ic in xrange(self.n_columns()):
if self(ir,ic) != other[ir,ic]: return False
return True
def resolve_partitions(self):
nr,nc = self.n
result_nr = 0
for ir in xrange(nr):
part_nr = 0
for ic in xrange(nc):
part = self(ir,ic)
assert isinstance(part, rec)
if (ic == 0): part_nr = part.n[0]
else: assert part.n[0] == part_nr
result_nr += part_nr
result_nc = 0
for ic in xrange(nc):
part_nc = 0
for ir in xrange(nr):
part = self(ir,ic)
if (ir == 0): part_nc = part.n[1]
else: assert part.n[1] == part_nc
result_nc += part_nc
result_elems = [0] * (result_nr * result_nc)
result_ir = 0
for ir in xrange(nr):
result_ic = 0
for ic in xrange(nc):
part = self(ir,ic)
part_nr,part_nc = part.n
i_part = 0
for part_ir in xrange(part_nr):
i_result = (result_ir + part_ir) * result_nc + result_ic
for part_ic in xrange(part_nc):
result_elems[i_result + part_ic] = part[i_part]
i_part += 1
result_ic += part_nc
assert result_ic == result_nc
result_ir += part_nr
assert result_ir == result_nr
return rec(elems=result_elems, n=(result_nr, result_nc))
class mutable_rec(rec):
container_type = list
def __setitem__(self, i, x):
self.elems[i] = x
class row_mixin(object):
def __init__(self, elems):
super(row_mixin, self).__init__(elems, (1, len(elems)))
class row(row_mixin, rec): pass
class mutable_row(row_mixin, mutable_rec): pass
class col_mixin(object):
def __init__(self, elems):
super(col_mixin, self).__init__(elems, (len(elems), 1))
def random(cls, n, a, b):
uniform = random.uniform
return cls([ uniform(a,b) for i in xrange(n) ])
random = classmethod(random)
class col(col_mixin, rec): pass
class mutable_col(col_mixin, mutable_rec): pass
class sqr(rec):
def __init__(self, elems):
l = len(elems)
n = int(l**(.5) + 0.5)
assert l == n * n
rec.__init__(self, elems, (n,n))
################################################################################
# This code has been adapted from cctbx.sf.net, file scitbx/matrix/eigensystem.h
class real_symmetric:
def __init__(self,symmat3):
m = symmat3
self.relative_epsilon = 1.e-10
self.absolute_epsilon = 0
self.a = [symmat3[0],symmat3[3],symmat3[1],symmat3[4],symmat3[5],symmat3[2]]
self.vectors_ = [0.,0.,0.,0.,0.,0.,0.,0.,0.]
self.values_ = [0.,0.,0.];
self.min_abs_pivot_ = self.real_symmetric_given_lower_triangle(
self.a,
3,
self.vectors_,
self.values_,
self.relative_epsilon,
self.absolute_epsilon);
def real_symmetric_given_lower_triangle(self,
a, # size of memory pointed to by a must be n*(n+1)/2
n,
eigenvectors,
eigenvalues,
relative_epsilon,
absolute_epsilon):
assert (relative_epsilon >= 0);
assert (absolute_epsilon >= 0);
if (n == 0): return 0;
# The matrix that will hold the results is initially = I.
for x in xrange(n*n):
eigenvectors[x] = 0.0;
for x in xrange(0,(n*n),(n+1)):
eigenvectors[x] = 1.0;
# Setup variables
#std::size_t il, ilq, ilr, im, imq, imr, ind, iq;
#std::size_t j, k, km, l, ll, lm, lq, m, mm, mq;
#FloatType am, anorm, anrmx, cosx, cosx2, sincs, sinx, sinx2, thr, x, y;
# Initial and final norms (anorm & anrmx).
anorm=0.0;
iq=0;
for i in xrange(n):
for j in xrange(i+1):
if (j!=i): anorm+=a[iq]*a[iq];
iq+=1
anorm=math.sqrt(2.0*anorm);
anrmx=relative_epsilon*anorm/n;
if (anrmx < absolute_epsilon): anrmx = absolute_epsilon;
if (anorm>0.0):
# Compute threshold and initialise flag.
thr=anorm;
while (thr>anrmx): # Compare threshold with final norm
thr/=n;
ind=1;
while (ind):
ind=0;
l=0;
while (l != n-1): # Test for l beyond penultimate column
lq=l*(l+1)/2;
ll=l+lq;
m=l+1;
ilq=n*l;
while (m != n): # Test for m beyond last column
# Compute sin & cos.
mq=m*(m+1)/2;
lm=l+mq;
if (a[lm]*a[lm]>thr*thr):
ind=1;
mm=m+mq;
x=0.5*(a[ll]-a[mm]);
denominator=math.sqrt(a[lm]*a[lm]+x*x);
assert (denominator != 0);
y=-a[lm]/denominator;
if (x<0.0): y=-y;
sinx=y/math.sqrt(2.0*(1.0+(math.sqrt(1.0-y*y))));
sinx2=sinx*sinx;
cosx=math.sqrt(1.0-sinx2);
cosx2=cosx*cosx;
sincs=sinx*cosx;
# Rotate l & m columns.
imq=n*m;
for i in xrange(n):
iq=i*(i+1)/2;
if (i!=l and i!=m):
if (i<m): im=i+mq;
else: im=m+iq;
if (i<l): il=i+lq;
else: il=l+iq;
x=a[il]*cosx-a[im]*sinx;
a[im]=a[il]*sinx+a[im]*cosx;
a[il]=x;
ilr=ilq+i;
imr=imq+i;
x = (eigenvectors[ilr]*cosx) \
- (eigenvectors[imr]*sinx);
eigenvectors[imr] = (eigenvectors[ilr]*sinx) \
+ (eigenvectors[imr]*cosx);
eigenvectors[ilr] = x;
x=2.0*a[lm]*sincs;
y=a[ll]*cosx2+a[mm]*sinx2-x;
x=a[ll]*sinx2+a[mm]*cosx2+x;
a[lm]=(a[ll]-a[mm])*sincs+a[lm]*(cosx2-sinx2);
a[ll]=y;
a[mm]=x;
m+=1;
l+=1;
# Sort eigenvalues & eigenvectors in order of descending eigenvalue.
k=0;
for i in xrange(n-1):
im=i;
km=k;
am=a[k];
l=0;
for j in xrange(n):
if (j>i and a[l]>am):
im=j;
km=l;
am=a[l];
l+=j+2;
if (im!=i):
a[km]=a[k];
a[k]=am;
l=n*i;
m=n*im;
for jj in xrange(n):
am=eigenvectors[l];
eigenvectors[l] = eigenvectors[m];
l+=1;
eigenvectors[m] = am;
m+=1
k+=i+2;
# place sorted eigenvalues into the matrix_vector structure
k = 0
for j in xrange(n):
eigenvalues[j]=a[k];
k+=j+2;
return anrmx;
def vectors(self):
return matrix.sqr(self.vectors_)
def values(self):
return matrix.col(self.values_)
class module: pass
eigensystem = module()
eigensystem.real_symmetric = real_symmetric
matrix = module()
matrix.sqr = sqr
matrix.col = col
matrix.rec = rec
######################################################################
#This code is adapted from cctbx.sf.net, file scitbx/math/superpose.py
def kabsch_rotation(reference_sites, other_sites):
"""
Kabsch, W. (1976). Acta Cryst. A32, 922-923.
A solution for the best rotation to relate two sets of vectors
Based on a prototype by Erik McKee and Reetal K. Pai.
"""
assert len(reference_sites) == len(other_sites)
sts = matrix.sqr(other_sites.transpose_multiply(reference_sites))
sym_mat3_input = (sts * sts.transpose()).as_sym_mat3()
eigs = eigensystem.real_symmetric(sym_mat3_input)
vals = list(eigs.values())
vecs = list(eigs.vectors())
a3 = list(matrix.col(vecs[:3]).cross(matrix.col(vecs[3:6])))
a = matrix.sqr(list(vecs[:6])+a3)
b = list(a * sts)
for i in xrange(3):
d = math.sqrt(math.fabs(vals[i]))
if (d > 0):
for j in xrange(3):
b[i*3+j] /= d
b3 = list(matrix.col(b[:3]).cross(matrix.col(b[3:6])))
b = matrix.sqr(b[:6]+b3)
return b.transpose() * a
#######################################
#This code is new for this application:
class residue:
def __init__(self,chain,resid,atoms):
self.C = None
self.CA = None
self.CB = None
self.N = None
for atom in atoms.atom:
if chain != str(atom.chain): continue
if str(resid) != str(atom.resi): continue
if atom.name=="N":
self.N = matrix.col((atom.coord[0],atom.coord[1],atom.coord[2]))
if atom.name=="CA":
self.residue_type = atom.resn
self.CA = matrix.col((atom.coord[0],atom.coord[1],atom.coord[2]))
if atom.name=="C":
self.C = matrix.col((atom.coord[0],atom.coord[1],atom.coord[2]))
if atom.name=="CB":
self.CB = matrix.col((atom.coord[0],atom.coord[1],atom.coord[2]))
assert self.N != None
assert self.CA!= None
assert self.C != None
if self.CB==None: #generate a CB position for Glycine
self.CB = self.constructed_CB()
def constructed_CB(self):
#refer to the documentation for geometrical construction
unit_N = (self.N - self.CA).normalize()
assert abs(unit_N.length() - 1.0) < 0.0001
unit_C = (self.C - self.CA).normalize()
on_bisector = (unit_N + unit_C).normalize()
unit_rotation_axis = (unit_N - unit_C).normalize()
expected_angle_bisector_to_CB = math.acos (-1./math.sqrt(3.))
#Use Euler-Rodrigues formula for rotation
unit_on_CB = self.rotation_formula(on_bisector,unit_rotation_axis,
expected_angle_bisector_to_CB)
O_to_CB = 1.53 * unit_on_CB
return self.CA + O_to_CB
def rotation_formula(self,vector,axis,theta):
return math.cos(theta)*vector + \
math.sin(theta)*(axis.cross(vector)) + \
(1.-math.cos(theta))*((axis.dot(vector))*axis)
def sites(self):
diff_N = self.N - self.CA
diff_C = self.C - self.CA
diff_CB= self.CB- self.CA
all = []
for site in [diff_N,diff_C,diff_CB]:
for coord in [0,1,2]:
all.append(site[coord])
return matrix.rec(all,(3,3))
class std_residue(residue):
def __init__(self): pass
class view_matrix:
def __init__(self,getview_output):
#3x3 rotation matrix which transforms model to camera space
self.rotmat = matrix.sqr(getview_output[0:9])
#camera position in model space and relative to the origin of rotation
self.camera_position = matrix.col(getview_output[9:12])
#origin of rotation in model space
self.origin_of_rotation = matrix.col(getview_output[12:15])
#front plane distance from camera
self.front_plane = getview_output[15]
self.back_plane = getview_output[16]
self.orthoscopic_flag = getview_output[17]
def set_diff_to_CA(self,residue):
self.difference_to_this_CA = self.origin_of_rotation - residue.CA
return self
class new_set_view:
def __init__(self,std_view,kr,residue, verbose=True):
R_prime = (kr.inverse() * std_view.rotmat)
if verbose: print "delta",std_view.origin_of_rotation - residue.CA
test = residue.CA + (kr.inverse() * std_view.difference_to_this_CA)
if verbose:
print "set_view ( ",
for x in R_prime.elems: print "%10.4f,"%x,
for x in std_view.camera_position.elems: print "%10.4f,"%x,
for x in residue.CA.elems: print "%10.7f,"%x,
for x in [std_view.front_plane,std_view.back_plane,std_view.orthoscopic_flag]:
print "%10.4f,"%x,
print ")"
self.result = list(R_prime.elems)+\
list(std_view.camera_position.elems)+\
list(test.elems)+\
[std_view.front_plane,std_view.back_plane,std_view.orthoscopic_flag]
def view_matrix(self): return self.result
Author: Nick Sauter