Thursday 4 November 2010

Automorphisms, isomorphisms, symmetry classes and canonical labelling

I've just been spending some time ensuring that Tim's new symmetry and isomorphism code in Open Babel can be accessed from Python (you can already access it from C++). It didn't make it into the 2.3.0 official release but I'll push it out into the Windows Python in the next while.

First the symmetry classes and canonical labelling. These can be calculated and accessed as follows (the test molecule is shown in the image above):
import pybel
ob = pybel.ob

mol = pybel.readstring("smi", "c1c(O)cccc1(O)")

gs = ob.OBGraphSym(mol.OBMol)
symclasses = ob.vectorUnsignedInt()
gs.GetSymmetry(symclasses)

print symclasses
# <openbabel.vectorUnsignedInt; proxy of <Swig Object of type 'std::vector< unsigned int > *' at 0x004DFE90> >
print list(symclasses)

# [4, 5, 1, 3, 2, 3, 5, 1]

canonlabels = ob.vectorUnsignedInt()
ob.CanonicalLabels(mol.OBMol, symclasses, canonlabels)
print list(canonlabels)

# [4, 2, 1, 3, 5, 7, 6, 8]

Does the molecule have any automorphisms?
automorphs = ob.vvpairUIntUInt()
ob.FindAutomorphisms(mol.OBMol, automorphs, symclasses)
for x in automorphs:
    print x

# ((0, 0), (1, 6), (6, 1), (2, 7), (7, 2), (3, 5), (5, 3), (4, 4))
# ((0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7))

Next we'll use the example molecule as a query to see how it matches another identical molecule which has the atoms in a different order:
query = ob.CompileMoleculeQuery(mol.OBMol)
mapper = ob.OBIsomorphismMapper.GetInstance(query)

mol_b = pybel.readstring("smi", "c1(O)cccc(O)c1")

isomorph = ob.vpairUIntUInt()
mapper.MapFirst(mol_b.OBMol, isomorph)
print list(isomorph)

# [(0, 7), (1, 5), (3, 4), (4, 3), (5, 2), (6, 0), (2, 6), (7, 1)]

isomorphs = ob.vvpairUIntUInt()
mapper.MapAll(mol_b.OBMol, isomorphs)
for x in isomorphs:
    print x

# ((0, 7), (1, 5), (3, 4), (4, 3), (5, 2), (6, 0), (2, 6), (7, 1))
# ((0, 7), (1, 0), (2, 1), (3, 2), (4, 3), (5, 4), (6, 5), (7, 6))

7 comments:

  1. Nice work! I should probably re-read your other posts on this, but a quick question; what are you using to find automorphisms and to canonically label? Is is nauty?

    ReplyDelete
  2. The underlying code is all by Tim Vandermeersch, I should have pointed out. Here is some text he has written:
    "Open Babel implements a sophisticated canonicalization algorithm that can handle molecules or molecular fragments. The atom symmetry classes are the initial graph invariants and encode topological and chemical properties. A cooperative labeling procedure is used to investigate the automorphic permutations to find the canonical code. Although the algorithm is similar to the original Morgan canonical code{ref Morgan65}, various improvements are implemented to improve performance. Most notably, the algorithm implements heuristics from the popular nauty package {ref McKay 81}. Stereochemistry is also handled by canonical code. Stereochemistry is not a property of an isolated atom since different labelings can lead to different parities. This is further complicated by the possibility of symmetry-equivalent stereocenters and stereocenters whose configuration is interdependent."...coming soon to a journal near you :-)

    ReplyDelete
  3. Thanks for demo!
    Sorry if my question is out of the topic, but has OpenBabel anything to deal with periodicity, e.g. polymers or 2D films like graphene?

    ReplyDelete
  4. @jam31: I think your question is better asked at openbabel-discuss@lists.sourceforge.net. The short answer is, I think, no, but others may have some ideas if you describe what you are interested in.

    ReplyDelete
  5. Hi Noel, I tried to run your code, but I had this error:
    Traceback (most recent call last):
    File "test.py", line 3, in
    mappings = pybel.ob.vvpairUIntUInt()
    AttributeError: 'module' object has no attribute 'vvpairUIntUInt'

    Could you let me know which ob version you are using? Thanks!
    Boyu

    ReplyDelete
  6. It's in the development version, so you should either use this, or wait for 2.3.1.

    ReplyDelete
  7. I used the version in trunk, it works with this particular one, thank you!

    ReplyDelete