Consistent View/ Map Inspect

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