Hi JD, Take a look at the new Axes tool (under Tools... Structure Analysis in recent builds of Chimera): http://www.cgl.ucsf.edu/chimera/docs/ContributedSoftware/ structuremeas/structuremeas.html#axes It can show protein helix axes, but you can also calculate an axis for any set of selected atoms. Then you can perform other measurements such as distance to axis, crossing angle between two axes, etc. However, I am not sure whether this computed axis is what you need - it is not exactly an axis of inertia (doesn't depend on the atomic weights of the atoms, for example). Best, Elaine ----- Elaine C. Meng, Ph.D. meng@cgl.ucsf.edu UCSF Computer Graphics Lab and Babbitt Lab Department of Pharmaceutical Chemistry University of California, San Francisco http://www.cgl.ucsf.edu/home/meng/index.html On Mar 19, 2008, at 8:40 AM, Jean-Didier Maréchal wrote:

Hi everyone, I was wondering if there is a method available to calculate the principal axis of inertia of a protein (or a part of a protein) in chimera?

If not what would be the easiest way to move forward? All the best, JD

Just wondering why you weighted by atomic number rather than mass (a.element.mass)? --Eric On Mar 19, 2008, at 11:42 AM, Tom Goddard wrote:

Hi JD,

Attached is a script that computes the inertia axes of molecules. Open a molecule, then open the script. It will print the axes in the reply log and show an inertia ellipsoid. It is weighting the atoms by their atomic number, not by the average isotopic mass. Tested in Chimera 1.2498. Should work in older Chimera versions.

Tom

Jean-Didier Maréchal wrote:

...

Hi everyone,

I was wondering if there is a method available to calculate the principal axis of inertia of a protein (or a part of a protein) in chimera?

If not what would be the easiest way to move forward?

All the best,

JD

# ---------------------------------------------------------------------- ------- # Compute inertia tensor principle axes for molecule. # def inertia_ellipsoid(m):

atoms = m.atoms n = len(atoms) from _multiscale import get_atom_coordinates xyz = get_atom_coordinates(atoms) anum = [a.element.number for a in atoms] from numpy import array, dot, outer, argsort, linalg wxyz = array(anum).reshape((n,1)) * xyz mass = sum(anum)

c = wxyz.sum(axis = 0) / mass # Center of mass v33 = dot(xyz.transpose(), wxyz) / mass - outer(c,c) eval, evect = linalg.eigh(v33)

# Sort by eigenvalue size. order = argsort(eval) seval = eval[order] sevect = evect[:,order]

return sevect, seval, c

# ---------------------------------------------------------------------- ------- # def ellipsoid_surface(center, axes, lengths, color = (.7,.7,.7,1)):

from Icosahedron import icosahedron_triangulation varray, tarray = icosahedron_triangulation(subdivision_levels = 3, sphere_factor = 1.0) from numpy import dot, multiply es = dot(varray, axes) ee = multiply(es, lengths) ev = dot(ee, axes.transpose()) ev += center

import _surface sm = _surface.SurfaceModel() sm.addPiece(ev, tarray, color) return sm

# ---------------------------------------------------------------------- ------- # def show_ellipsoid(axes, d2, center, model):

from math import sqrt d = [sqrt(e) for e in d2] sm = ellipsoid_surface(center, axes, d) sm.name = 'inertia ellipsoid for %s' % m.name from chimera import openModels as om om.add([sm], sameAs = model)

# ---------------------------------------------------------------------- ------- # def print_axes(axes, d2, m):

from math import sqrt paxes = ['\tv%d = %6.3f %6.3f %6.3f d%d = %6.3f' % (a+1, axes[a][0], axes[a][1], axes[a][2], a+1, sqrt(d2[a])) for a in range(3)] from chimera.replyobj import info info('Inertia axes for %s\n%s\n' % (m.name, '\n'.join(paxes))) from Accelerators.standard_accelerators import show_reply_log show_reply_log()

# ---------------------------------------------------------------------- ------- # from chimera import openModels as om, Molecule mlist = om.list(modelTypes = [Molecule]) for m in mlist: axes, d2, center = inertia_ellipsoid(m) print_axes(axes, d2, m) show_ellipsoid(axes, d2, center, m) _______________________________________________ Chimera-users mailing list Chimera-users@cgl.ucsf.edu http://www.cgl.ucsf.edu/mailman/listinfo/chimera-users

Oops. That script I posted for calculating the inertia of a molecule had some problems. The printed matrix of axes was transposed though the correct axes were used for the displayed ellipsoid. Also I was using inertia matrix I_ij = sum_over_k(m_k*r_ki*r_kj) omitting the diagonal term. Should have been I_ij = sum_over_k(m_k*(|r_k| - r_ki*r_kj)). And finally I incorrectly computed the ellipsoid size along each axis. Everything was a little bit wrong. Here's the fixed version. Tested in Chimera version 1.2500 but should work with earlier versions. Tom # ----------------------------------------------------------------------------- # Compute inertia tensor principle axes for molecule. # def inertia_ellipsoid(m): atoms = m.atoms n = len(atoms) from _multiscale import get_atom_coordinates xyz = get_atom_coordinates(atoms) amass = [a.element.mass for a in atoms] from numpy import array, identity, dot, outer, argsort, linalg wxyz = array(amass).reshape((n,1)) * xyz mass = sum(amass) c = wxyz.sum(axis = 0) / mass # Center of mass i = (xyz*wxyz).sum()*identity(3) - dot(xyz.transpose(), wxyz) # Adjust inertia tensor to be about center of mass i -= mass*dot(c,c)*identity(3) - mass*outer(c,c) i /= mass eval, evect = linalg.eigh(i) # Sort by eigenvalue size. order = argsort(eval) seval = eval[order] sevect = evect[:,order] # Make rows of 3 by 3 matrix the principle axes. return sevect.transpose(), seval, c # ----------------------------------------------------------------------------- # def ellipsoid_surface(center, axes, lengths, color = (.7,.7,.7,1)): from Icosahedron import icosahedron_triangulation varray, tarray = icosahedron_triangulation(subdivision_levels = 3, sphere_factor = 1.0) from numpy import dot, multiply es = dot(varray, axes.transpose()) ee = multiply(es, lengths) ev = dot(ee, axes) ev += center import _surface sm = _surface.SurfaceModel() sm.addPiece(ev, tarray, color) return sm # ----------------------------------------------------------------------------- # def show_ellipsoid(axes, d2, center, model): from math import sqrt d = [sqrt(e) for e in d2] elen = [sqrt(d[1]*d[2]), sqrt(d[0]*d[2]), sqrt(d[0]*d[1])] sm = ellipsoid_surface(center, axes, elen) sm.name = 'inertia ellipsoid for %s' % m.name from chimera import openModels as om om.add([sm], sameAs = model) # ----------------------------------------------------------------------------- # def print_axes(axes, d2, m): from math import sqrt paxes = ['\tv%d = %6.3f %6.3f %6.3f d%d = %6.3f' % (a+1, axes[a][0], axes[a][1], axes[a][2], a+1, sqrt(d2[a])) for a in range(3)] from chimera.replyobj import info info('Inertia axes for %s\n%s\n' % (m.name, '\n'.join(paxes))) from Accelerators.standard_accelerators import show_reply_log show_reply_log() # ----------------------------------------------------------------------------- # from chimera import openModels as om, Molecule mlist = om.list(modelTypes = [Molecule]) for m in mlist: axes, d2, center = inertia_ellipsoid(m) print_axes(axes, d2, m) show_ellipsoid(axes, d2, center, m)