Tuesday, 27 February 2018

Calling all students: Google Summer of Code and CSA Trust Grant both happening now!

Deadlines are fast approaching for the following:

Open Chemistry Google Summer of Code: If you're a student and interested in hacking on open source chemistry software (and get paid a bit for the privilege), then get on over to the Google Summer of Code (GSoC) page of the Open Chemistry project. A whole bunch of Open Source chemistry tools have gathered together (including Open Babel, natch) and come up with project ideas that hopefully will spark interest. If you've always wanted to get involved with Open Source but didn't know how, this is a good chance to do so.

CSA Trust Grant: I was a recipient of a CSA Trust Grant myself, and funny story, I'm now a CSA Trustee and on the Grant Committee. Applications are now invited - it's pretty straightforward. If you look at the details of previous recipients you can see the sorts of things they applied for. The success rate is pretty high, so I really do encourage you to apply. And if you don't get it this year, well, try again next year (it worked for me).

Saturday, 3 February 2018

Using those other processors

This year, I've decided it's time to make more use of those other processors that my PC has. At some point, it gets embarrassing to see my single CPU job trundling along at its max speed, while 80 or 90% of the processing power is just sitting idle.

My typical input is a SMILES file such as ChEMBL or PubChem, or a file I've created by processing these, or a set of reaction SMILES extracted from patents. In any case, the problem is the same - a large number of inputs that can be processed in parallel. Sure, you can split the file (see the 'split' command, though be careful to split on lines) but I tend to avoid that as it's a bit clunky. If on Linux, there's the magic "GNU parallel". But what I really should be doing is making use of multiple CPUs directly in my Python scripts.

And actually, that's not too hard. The multiprocessing module makes it pretty easy once you've figured it out. Here's my goto template for these calculations. I try to keep it as slimmed down as possible, because there's some magic going on behind the scenes and I don't want to have to debug any complex problems. (Note that there's also a multithreading module, but CPython cannot parallelize CPU-bound calculations using threads due to its Global Interpreter Lock.)

import multiprocessing as mp
import pybel

def calculate(data):
    return pybel.readstring("smi", data).molwt

if __name__ == "__main__":
    POOLSIZE = 4 # the number of CPUs
    CHUNKSIZE = 1000
    pool = mp.Pool(POOLSIZE)
    with open("output.txt", "w") as out:
        with open(r"C:\Tools\LargeData\chembl_23.smi", "r") as inp:
            for result in pool.imap(calculate, inp, CHUNKSIZE):
            # for result in pool.imap_unordered(calculate, inp, CHUNKSIZE):
            # for result in map(calculate, inp): # no multiprocessing
                out.write("%f\n" % result)

1. This blog is now Python 3.
2. By editing the commented-out code, you can choose one of three variations. During development, you can avoid multiprocessing entirely by just using a regular map (equivalent to itertools.imap in Python 2). If using multiprocessing, you can choose to have the output in the same order as the input, or just have it as it comes; the latter is presumably faster.
3. If you need to kill a multiprocessing job, CTRL+C just won't do it, as it spawns an additional process to replace it (if anyone knows of a way to make this work please let me know). You need to use the Process Manager, choose the master Python process and kill that one.

Tuesday, 23 January 2018

Which came first, the atom or the bond?

Let's suppose you need to iterate over the neighbours of an atom, and you want to know both the corresponding bond order and the atomic number of the neighbour. Given that toolkits typically provide iterators over the attached bonds or the attached atoms, but not both simultaneously, how exactly should you do this?

Should you:
  • (a) iterate over the neighbouring atoms, and then request the bond joining the two atoms?
  •     for nbr in ob.OBAtomAtomIter(atom):
            bond = atom.GetBond(nbr)
  • or (b) iterate over the attached bonds, and ask for the atom at the other end?
  •     for bond in ob.OBAtomBondIter(atom):
            nbr = bond.GetNbrAtom(atom)

Obviously, either way you get your atom and bond. But which is more efficient? Clearly, the answer to this depends on the internals of the toolkit. But if you assume that each atom knows its attached atoms and bonds, then it's only the second step that determines the relative efficiency. That is:
  • (a) given two atoms find the bond that joins them, versus
  • (b) given an atom and a bond find the atom at the other end

Since the implementation of (a) will probably involve the same test that (b) is doing plus additional work, it follows that (b) must be more efficient. I never really thought about this before until I was writing the kekulization code for Open Babel. It's the sort of thing that's useful to work out once and then apply in future without thinking. Sure, the speed difference may be minimal but given that you have to choose, you might as well write it the more efficient way.

Thursday, 18 January 2018

Implementing the Sayle tautomer hash with Open Babel

One of the consequences of last year's overhaul of handling of kekulization, aromaticity and implicit hydrogens is the ability to easily calculate structure hashes such as Roger Sayle's tautomer hash, which I wrote up on the NextMove blog a while ago.

I routinely use this hash to handle tautomers, particularly when dealing with R groups. It doesn't solve all tautomer issues (e.g. ones that involve carbon) but it can quickly bring you from having no support for handling tautomers to getting you 95% of the way. In fact, I've been thinking about adding this (and some of the other hashes that Roger has come up with) as cansmi output options. Anyhoo, here's an implementation in Python:
import pybel
ob = pybel.ob

def tautomerhash(smi):
    """Take a SMILES and return the Sayle tautomer hash:
    mol = pybel.readstring("smi", smi).OBMol
    mol.DeleteHydrogens() # just in case
    formalcharges = 0
    hcount = 0
    for atom in ob.OBMolAtomIter(mol):
        formalcharges += atom.GetFormalCharge()
        if atom.GetAtomicNum() != 6: # non-carbon
            hcount += atom.GetImplicitHCount()
    for bond in ob.OBMolBondIter(mol):
    mol.SetAromaticPerceived() # no point triggering perception
    return "%s_%d" % (pybel.Molecule(mol).write("can").rstrip(), hcount-formalcharges)

if __name__ == "__main__":
    smis = ["*c1c(c(C(=N)O)cc2nc([nH]c12)C(=O)[O-])N(=O)=O",
    for smi in smis:
The two SMILES in the example code above are those from the original blog post. Here they give the following two identical hashes (note to self: 'fix' the bracketed asterisk):

Sunday, 17 December 2017

Faster toolkit, faster! Part II

A while ago I described some work I was doing to improve the overall speed of the Open Babel toolkit. Here I want to focus on a particular use-case and where I've got to.

This use case is SMILES to SMILES conversion. Now, while this particular transformation might not sound very interesting, it does encompass both SMILES reading and SMILES writing in one handy package, and both of these operations are often important when dealing with databases or datasets of chemical structures. It also exercises several areas of the toolkit such as kekulization, handling of aromaticity, and stereo perception (or not, as we'll see). Canonicalization is also relevant here, but I didn't do any work on that (and it needs some).

To begin with, some timings. To convert 100K ChEMBL molecules from smi to smi took 10m7s with OB 2.4.1. With the current development version it takes 31s. One change to the defaults is that the dev version does not reperceive the stereo. If you turn on stereo perception (-aS), it takes 1m13s. You can speed things up if you also avoid reperceiving the aromaticity (-aa) and read it as provided in the input; then the conversion only takes 19s.

[21 days later] So I was originally going to describe the results from the Visual Studio profiler for this conversion. But then I said, hey, I might as well fix that one, and that one there, and, well, you know how it goes - this part is actually quite fun, when you make a small change and see the speed go up. Anyway, the conversion that used to take 19s, now takes 11.0s. If you're interested, the speedups included things like replacing std::endl by "\n", caching option values, avoiding string copies, avoiding use of stringstream, avoiding SSSR calculation, and using reserve() on vectors. It was often surprising what things appeared high on the list in the profiler. I can see a few more things that could be improved, but I'm going to leave it there for the moment.

So, in summary, this particular conversion has gone from slow to fast, with a speedup of 55x. There's always more that could be done, but it's respectable.

Open Babel in a snap II

This is a quick update on a previous post, where I mentioned that I had made Open Babel available as a snap. That snap is the released version, 2.4.1. I have now put in place a semi-automated procedure that builds a snap of the development version. You can only have one snap version installed at a time, but you can switch between them.

To install the stable version use:
sudo snap install openbabel

To install the development version use:
sudo snap install openbabel --channel=edge

You can switch between them with:
sudo snap refresh openbabel --channel=stable # or edge

To see which you have installed, use "snap list", or run "openbabel.obabel" and look at the version number and date.

Notes: I'm using Launchpad to do this. Rather than base it directly off the openbabel master (which would require me to check-in snapcraft specific files), I've set it up so that it runs off a branch (named "snaps") in my own repo. Every so often, I merge master into this and a new snap will be created. To fully automate it, I will need to have a cronjob to do that merge automatically.

Thursday, 2 November 2017

Announcing more closures

It is a good idea to avoid extending an existing format. If additional information needs to be stored, many formats provide a means to store it separately so that software that does not support the extension will still be able to read and process the core information. For example, with SMILES, ChemAxon have famously used the title field to store additional info.

However, this is not always possible if one runs up against a fundamental limitation of a format.

I recently ran into a situation where I needed to deal with the fact that the SMILES syntax only supports at most two digit ring closure numbers. This is the %nn notation. The problem is that %111 does not mean ring closure 111 but rather ring closures 11 and 1. If you need a three digit number you are out of luck - it simply can't be expressed as a SMILES string.

So, is it really necessary to have this ability in the first place? Daylight didn't think so. And they were mainly right - when writing SMILES, if you reuse the numbers and pick a traversal order that closes rings as soon as possible then it may be possible to avoid requiring three digits (though a diamond or graphene substructure may cause problems - a refutable hypothesis).

But the thing is, any particular implementation may not be reusing ring closure numbers, or may be favoring a different traversal order (e.g. nicer looking SMILES are made by following ring substituents first). Indeed, remembering that this syntax can be used to create a bond between any two atoms, you may wish to use the ring closure syntax as a means to construct a molecule, and there you don't want to be limited to just 100 bonds; e.g. by stitching together a large molecule by combining shorter SMILES strings (see SmiLib for example, and the note under Restrictions), or via this SMILES hack of Andrew's.

When Bob Hanson of Jmol had to deal with this issue he settled on the syntax "%(number)" (discussed here in 2012, and more recently in J. Cheminf. 2016, 8, 50). One alternative might be "%%", but again this would be limited, this time to three digits. Might as well come up with a general solution as Bob has done. It uses an extra character for the three-digit case but this is going to be rare in the first place, and these SMILES are going to be so big that having a few more characters is not going to break the bank.

Support for this syntax has recently been added to RDKit and Open Babel. For more info see the associated pull requests (here and here). Hopefully we will start to see support for this more generally in other toolkits.

1. According to the Jmol docs, "%(n)" supports any non-negative number (integer, presumably?). My implementation supports up to 5 digits.

Image credit:
image by InExtremiss on Flickr licensed CC-BY