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