On Jan 9, 2008, at 11:56 AM, Buz Barstow wrote:
Does anyone have a python routine that can be used to calculate the
moment
of inertia matrix and principal axes of a selection of atoms?
I compulsively write python code sometimes--its that fun of a
language...
If you can get the x, y, z, and mass of your atoms, then you can use
the little module below (moi.py). It approximates atomic nuclei as
dimensionless ;-). An example usage is at the bottom. "JTools" is my
own PDB library, but it gives an idea of use. Once you have the moi
tensor, I think Eigenvalue decomposition gets the principle axes, so
this becomes a 1-op with numpy. Please let me know if you find errors.
James
========================================================================
#! /usr/bin/env python
"""
moi.py
(c) 2008, James Stroud
This module is released under the GNU Public License 3.0.
See <http://www.gnu.org/licenses/gpl.html>.
A module to find the moment of inertia of some atoms.
See <http://en.wikipedia.org/wiki/Moment_of_inertia> for
a mathematical discussion.
Anywhere "ref" is optional defaults to the center of mass
except where indicated.
The "sel" argument is always a list of Atoms.
"""
import math
class Atom(object):
def __init__(self, x, y, z, mass):
self.x = float(x)
self.y = float(y)
self.z = float(z)
self.mass = float(mass)
def __str__(self):
return "%s @ (%s, %s, %s)" % (self.x, self.y, self.z, self.mass)
class Point(object):
def __init__(self, x, y, z):
self.x = x
self.y = y
self.z = z
def __str__(self):
return "(%s, %s, %s)" % (self.x, self.y, self.z)
def dist(atom, ref=None):
"""
Calculates distance between an atom and a ref.
If ref is not given, then the origin is assumed.
"""
if ref is None:
d = math.sqrt(atom.x**2 + atom.y**2 + atom.z**2)
else:
x = atom.x - ref.x
y = atom.y - ref.y
z = atom.z - ref.z
d = math.sqrt(x**2 + y**2 + z**2)
return d
def com(sel):
"""
Calculates the center of mass of a list of atoms.
"""
m = sum([atom.mass for atom in sel])
c = []
for i in 'xyz':
mr = []
for atom in sel:
mr.append(atom.mass * getattr(atom, i))
c.append(sum(mr)/m)
return Point(*c)
def kd(i, j):
"""
The Kronaker delta.
"""
return int(i==j)
def moi_element(i, j, sel, ref):
"""
Generalized moment of inertia element--both diagonal and off
diagonal. Coordinate axes i and j should be 'x', 'y', or 'z'.
The sel argument is a list of Atoms.
The ref arguemnt is a Point.
"""
delta = kd(i, j)
refi = getattr(ref, i)
refj = getattr(ref, j)
el = []
for atom in sel:
r2d = (dist(atom, ref)**2) * delta
ri = getattr(atom, i) - refi
rj = getattr(atom, j) - refj
el.append(atom.mass*(r2d - ri*rj))
return sum(el)
def tensor_moi(sel, ref=None):
"""
Calculates the moment of inertia tensor.
"""
if ref is None:
ref = com(sel)
indices = [(i,j) for i in 'xyz' for j in 'xyz']
tensor = [moi_element(i, j, sel, ref) for (i,j) in indices]
return [tensor[0:3], tensor[3:6], tensor[6:9]]
def scalar_moi(sel, ref=None):
"""
Calculates the scalar moment of inertia.
"""
if ref is None:
ref = com(sel)
i = []
for atom in sel:
i.append(atom.mass * dist(atom, ref)**2)
return sum(i)
def test():
"""
You will want to construct your own test.
"""
import JTools
import numpy
s = JTools.build_pdb('1ass.pdb')
sel = []
for p in s[0]:
for a in p:
sel.append(Atom(a.x, a.y, a.z, a.get_weight()))
print "Center of mass of chain 0:"
print com(sel)
print
print "The moi as tensor:"
tensor = numpy.array(tensor_moi(sel))
print tensor
print
print "The eigenvalue decomposition:"
decomp = numpy.linalg.eig(tensor)
print "The eigenvalues:"
print decomp[0]
print "The eigenvectors:"
print decomp[1]
print
print "The moi as scalar:"
print scalar_moi(sel)
print
print "The moi as scalar around ori:"
print scalar_moi(sel, Point(0, 0, 0))
print
if __name__ == "__main__":
test()
========================================================================
Thanks! and all the best,
--Buz
--
James Stroud
UCLA-DOE Institute for Genomics and Proteomics
611 Charles E. Young Dr. S.
Los Angeles, CA 90095
http://www.jamesstroud.com
