Making a Molecular editor/visualizer: Object oriented programming, data structures, and molecules

Question

I am new to programming and I am trying to solve my first big problem and write my first big program. I have looked for open source examples of code to learn from, but so far have only found code in languages I don't fully understand or that does related but still-too-distant things for me to really learn from. I am having trouble taking some conceptual steps here.

I want make a simple piece of software to build, modify, and later represent small organic molecules. This is mainly a learning exercise. A user would give a SMILES string or choose from a basic set of starter molecules, and then could build on that molecule either graphically or via a text input syntax. However, I am not even at that point of complexity yet. I can't even fully comprehend how to make classes/object to store the molecules. So, my question succinctly: How do I use classes/objects to build molecules while retaining all levels of information, and what data structures should I use as attributes to which objects? And, can objects be attributes to other objects?

Here is my train of thought so far: I was thinking to have a "Molecule" class, then a "Atom" class/subclass and a "Bond" subclass, and maybe a "FunctionalGroup" subclass also. That seems like a good place to start, but maybe I'm misunderstanding OOP and this is bad. But then my problem really gets confusing (for me). Even though I have all these concepts/ideas/classes, I don't fully grasp what datastructures are needed to represent the molecule. A list of atoms would be a good thing to have. Could that list be a list of Atom objects? I would also need to have a way to store connectivity. A 2D matrix seems like a good idea, with bond order as integers in the matrix positions.

At this point I'm starting to get overwhelmed with the task. Is everything I'm doing so far make sense? Attaching a display/drawing aspect on top of this may mean that I need to rewrite/rework a lot of these things, but I'm just trying to get to a point where I can at least store molecules with the relevant data and then access that data to check/modify. I was thinking of doing this in Python, so the code/classes would perhaps look something like this: http://pastebin.com/uUi1BMzr

Perhaps this is really a programming questions for StackOverflow, but I thought it was specific enough to go here. Any assistance would be greatly appreciated, even if you just point out where I have made conceptual blunders. Thanks in advance.

Answer 1

writing software is an iterative process - write code -> see how far you can get then plan next steps -> write code -> repeat. At this stage in learning the art I'd suggest you jump in soon just to test the waters. No need to plan the entire system up front. Yes python would be a good first language. For visualization try out MatPlotLib, also NumPy and SciPy are handy. Industrial scale software always relies on pulling in pre-build libraries, instead of writing everything yourself, yet simple solutions can and should be written yourself, especially when learning to program. Your OO layout seems OK for now. When/if your object relationships need to change later, the act of re-factoring your code itself is a worthy experience to master. Welcome aboard !

Answer 2

Your goal has a lot of challenges. I'll break them down into parts.

SMILES is not a trivial language to parse, and the rules for aromaticity perception are not well defined. The detailed grammar definition from the OpenSMILES project should help some.

SMILES defines the topology, but does not give 2D or 3D information. Making either is hard. (That is, if you want it to look nice.)

Really, you should look to the RDKit cheminformatics toolkit (or OpenBabel, but I prefer RDKit). It has a SMILES parser built in, as well as 2D layout, and I believe 3D conformation generation. OpenBabel does as well.

Then for display you'll have to figure out the GUI system. Actually, here the CDK cheminformatics toolkit, in Java, is the most advanced.

But you're in the basics of how to represent a molecule. There are differences between small-molecule and large-molecule (protein, DNA) data models, but since you're interested in SMILES that means you're small-molecule oriented.

You can look at the API docs for RDKit, OpenBabel, CDK, OEChem, and Indigo, among many others. That will give you an idea of the diverse ways that people develop their class API. Of these, I prefer OEChem the most, followed by RDKit. Even though OEChem is open source, the API is online and freely readable, along with examples of use.

In short, have a Molecule class, with a list of Atom and Bond instances. "mol.AddAtom(element number)" creates a new atom, with no bonds, "mol.AddBond(atom1, atom2, bond_type)" make a bond connection. Each bond needs to know the atoms its connected to, and each atom needs a list of bonds. This leads to a lot of cycles in the data structure, but it's needed so that various algorithms, like connectivity search, can be done in linear time.

Do not use a 2D matrix. While feasible for a small molecule, it doesn't scale that well, and there's no need for it. Few algorithms need the connectivity matrix, and it's easily generated if/when needed.

Don't have a "FunctionalGroup". It's too specialized. Use something like a "Subset" or "Fragment", which contains the list of atoms and bonds you're interested in. That way you can also handle things like "selected atoms" and "ring substructure" and "scaffold" by referring to the specific subset.

I looked at your pastebin. The parser shouldn't work that way. You should separate parsing from your actual molecule structure. Try something like this:

class Molecule(object):
    def __init__(self):
        self.atoms = []
        self.bonds = []
        self._atom_id = 0
        self._bond_id = 0
    def _next_atom_id(self):
        atom_id = self._atom_id
        self.atom_id += 1
        return atom_id
    def AddAtom(self, eleno):
        self.atoms.append(Atom(self, self._next_atom_id(), eleno))
    def AddBond(self, atom1, atom2, bondtype):
        assert atom1.molecule is atom2.molecule
        self.bonds.append(Bond(self, self._next_bond_id(),
                               atom1, atom2, bondtype))

class Atom(object):
    def __init__(self, molecule, id, eleno):
        self.molecule = molecule
        self.id = id
        self.eleno = eleno
        self.charge = 0
        self.isotope = 0
   ..

and then the parser for a simple linear chain like "C C O" is:

def parse_linear_chain(text):
   mol = Molecule()
   prev_atom = None
   for atom_symbol in text.split():
     eleno = lookup_symbol[atom_symbol]
     atom = mol.NewAtom(eleno)
     if pre_atom is not None:
       mol.AddBond(prev_atom, atom, 1)
     prev_atom = atom
   return mol

Of course, a full SMILES parser is much more complicated than this, and the full data model has to handle things like hydrogen counts, which are often implicit.

The OpenBabel, RDKit, and CDK mailing lists are also good places to go, should you decide to use one of those toolkits. There's also the "Blue Obelisk" Q&A site, hosted by Shapado.

Answer 3

Another approach to get started might be to have a look at some code that already does something related to your problem. In that case your exercise might even end up in some other program, and that would be nice, wouldn't it?

Programs that might be of interest for you are

• The MD simulation and molecular modeling package MMTK (as already suggested by Nate above)

• The visualization package PyMol

Answer 4

Learning the details of a molecular object system is a wonderful way for chemists to learn object oriented programming. You'll find that implementation of such a system will refine your molecular intuitions. You should do some hard thinking about atoms, molecules, and perhaps collections of molecules in terms of attributes and methods.

Here are some nice python (a bit older) slides that may help: http://www.wag.caltech.edu/home/rpm/python_course/Lecture_4.pdf

For checking your work: in addition to openbabel (which has python bindings!) and MMTK, there is ELBOW in phenix.

For the polyglot in you, there is also PerlMol (Perlmol.org). PerlMol is written in object oriented perl and can be downloaded from CPAN.