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

[Ralf W. Grosse-Kunstleve]

Hi Buz,

> 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?


If you have Phenix, CCI Apps or just cctbx installed, try the script below.
Example:

  mmtbx.python script.py 1NSF.pdb 'resname GLU and resid 735'

See also: http://www.phenix-online.org/download/

Ralf


from mmtbx.monomer_library import pdb_interpretation
from mmtbx.monomer_library import server
from scitbx.math import principal_axes_of_inertia
from scitbx.array_family import flex
from libtbx.utils import Usage, plural_s
import sys, os

def run(args):
  if (len(args) < 2 or not os.path.isfile(args[0])):
    raise Usage('pdb_file "atom_selection" [...]')
  mon_lib_srv = server.server()
  ener_lib = server.ener_lib()
  processed_pdb_file = pdb_interpretation.process(
    mon_lib_srv=mon_lib_srv,
    ener_lib=ener_lib,
    file_name=args[0],
    log=sys.stdout)
  print
  acp = processed_pdb_file.all_chain_proxies
  selection_cache = acp.stage_1.selection_cache()
  aal = acp.stage_1.atom_attributes_list
  for selection_string in args[1:]:
    print selection_string
    isel = acp.iselection(string=selection_string, cache=selection_cache)
    print "  %d atom%s selected" % plural_s(isel.size())
    points = flex.vec3_double()
    for i_seq in isel:
      print "    %s" % aal[i_seq].pdb_format()
      points.append(aal[i_seq].coordinates)
    pai = principal_axes_of_inertia(points=points)
    print "  center of mass:", pai.center_of_mass()
    print "  inertia tensor:", pai.inertia_tensor()
    es = pai.eigensystem()
    print "  eigenvalues:", list(es.values())
    print "  eigenvectors:", list(es.vectors())
    print

if (__name__ == "__main__"):
  run(sys.argv[1:])