Quick navigation:        Home   |    Site Map   ||    References   |    Biography   ||    Copyright   |    Other copyright   |    Contact us   |   
Protein structure
 

Re: [ccp4bb] Principal axes and moment of inertia matrix
 

Basic tutorials: 		 
 

CCP4bb navigation

CCP4bb <-- 2008 <-- January 2008 <-- 09 January 2008
Previous message:
Subject: Re: unbiased electron density map
From: William Scott wgscott {- at -} CHEMISTRY {- dot -} UCSC {- dot -} EDU
Date: 2008-01-09 		Next message:
Subject: Re: unbiased electron density map
From: Dale Tronrud det102 {- at -} UOXRAY {- dot -} UOREGON {- dot -} EDU
Date: 2008-01-09


Subject: Re: Principal axes and moment of inertia matrix
From: James Stroud jstroud {- at -} MBI {- dot -} UCLA {- dot -} EDU
Date: 2008-01-09
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 .

A module to find the moment of inertia of some atoms.

See 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
ProteinCrystallography.org: Copyright 2006-2008 by Quid United Ltd