OB, Avogadro and Molecular Modelling

CMake, CTest and CDash: Creating portable c++ software

2010-10-11T19:09:00.001-07:00

CMake is a powerful build system for C++ applications (and other languages) used by many software projects including KDE. The main difference between CMake and other build systems is that it generates native unix make files, Microsoft Visual C++ solutions or XCode projects. Before we adopted cmake, we had to maintain separate build files for windows and unix/mac platforms. A developer using a specific platform would not update the other platform's build files which often resulted in compile errors when another developer checked out the new code. While fixing these errors is usually a trivial task, it's still an inconvenience. Using CMake solves this problem without adding any additional dependencies. CMake is small and binaries for a large number of platforms are provided.

CTest is part of CMake and can be used to test software. It doesn't provide any code but functions as a way to run the tests. For example, in OpenBabel we have a header where we define some functions we can use to write unit tests: OB_ASSERT, OB_REQUIRE and OB_COMPARE. Implementations of these functions are simple but any unit testing framework (code library, headers) can be used. Our functions print "FAIL" to the standard output and we use this as regular expression to allow CTest to determine if a test has failed. The advantage of using CTest is that we can now use the tests to improve the quality of our software. This becomes possible because CMake/CTest can submit the results to CDash.

The OpenBabel dashboard is available here

CDash is a web application containing the status of building and testing your software. KitWare provides CDash hosting but you still need to set up build "sites". These sites are just computers that run a CTest script as scheduled task or cronjob. For example, several developers and users from OpenBabel run these scripts every night and these results are submitted to CDash for all developers to see. We have multiple configurations to cover the most frequently used platforms and compilers (win, linux, mac, gcc, cygwin, MSVC, 32/64 bit, ...) but this depends on the platforms you want to support for your project. CTest synchronizes the Nightly builds so all sites will use the same revision which makes it easy to spot differences between platforms.

We also use Continuous builds which means that every time a developer checks in changes, an (incremental) build will be done, the tests will run and the results are submitted to CDash. This is complementary to the Nightly builds. For example, if for a Nightly run, a unit test failed on all platforms and there were 20 commits. The continuous results will identify which commit caused the test to start failing.

CDash can also be used for other information beyond compile warnings/errors and test results. For example, it is easy to include code code coverage and memory leak detection. For code coverage, CTest has support for GNU GCov and a commercial tool Bullseye. For OpenBabel, we split the coverage in two to differentiate between the core library and plugins. This is possible by using subprojects but it currently takes two builds and I'm still experimenting with the script. I'm told support for this will be improved in future versions but if the number of subprojects is small, the multiuiple build approach works. It runs cmake and defines some variables that enable code coverage for the library or formats (i.e. -fprofile-ars -ftest-coverage for gcov). The building directory is emptied between builds to ensure the gcov files are removed before analysing the new coverage results. The script can be found here and the changes to the OpenBabel CMakeLists.txt files here.

Another supported feature is memory checking to find memory leaks. In linux, valgrind is supported. This takes some time to run if the number of tests is large but running this Nightly is enough. It is also possible to submit additional files.

Acknowledgements: I'd like to thank Marcus Hanwell and Kitware for providing CDash hosting and other contributions to OpenBabel and Avogadro.

The SMILES stereochemistry enigma

2010-09-16T17:37:00.000-07:00

Almost every cheminformatics system implementing smiles limits the support for stereochemistry to the tetrahedral and double bond cases. Some (e.g. OpenBabel, Jmol, ...) implement the square planar SP tags which only has 3 possible values. However, the daylight smiles page has very limited documentation when it comes to trigonal-bipyramidal and octahedral stereochemistry. All strictly needed information is there though and this blog post will give an overview of how to obtain the meaning of the TB3-TB20 and OH3-OH30 tags.

Since the trigonal-bipyramidal case is slightly easier it will be discussed first. As can be seen in the image below, not all ligands are equivalent. There are two axial ligands (S and F) and three equatorial ligands (Cl, Br, I).

The TB1 tag is documented on the daylight website meaning viewing from the first ligand atom to the last, the three remaining ligands are ordered counter clockwise (i.e. @). A valid smiles for the depicted molecule would be F[X@TB1](I)(Br)(Cl)S. TB2 means viewing along the same axis, the three remaining are clockwise (i.e. @@). With this information it is possible to write any molecule by hand since you can reorder the atoms as needed. However, when writing out a canonical smiles, the order of the atoms is determined by the canonicalization algorithm and there are 5! = 120 possible ways to order 5 ligands. Since the ligands are not equivalent, there must be a way to encode the axis ligand positions using the tags. This hypothesis is supported by the fact that there are 10 ways to position the axis atoms (a-e, a-d, a-c, a-b, b-e, b-d, b-c, c-e, c-d, d-e), 10 * 2 (@|@@) = 20 tags and each tag has 3! = 6 permutations associated with it resulting in 120 permutations total.

To find the order of these axis, the more subtle information from the daylight smiles specification can be used. The first clue is that TH1 and TH2 are two cases of a generic chiral specification method. The other clues are that tags are actually "chiral permutation designators" and "it's table driven". Now all that is left to do is enumerate all permutations lexicographically and assign the designators when a new axis is found. The first column is the permutation number, the second is the permutation itself, the 3th column is are the axis positions, the 4th is the winding and finally the designator. (note: TB1/TB2 are designated first, skipping permutation 2)

1: 01234 0-4 @ TB1
2: 01243 0-3 @ TB3
3: 01324 0-4 @@ TB2
4: 01342 0-3 @@ TB4
5: 01423 0-2 @ TB5
6: 01432 0-2 @@ TB6
7: 02134 0-4 @@ TB2
8: 02143 0-3 @@ TB4
9: 02314 0-4 @ TB1
10: 02341 0-3 @ TB3
11: 02413 0-2 @@ TB6
12: 02431 0-2 @ TB5
13: 03124 0-4 @ TB1
14: 03142 0-3 @ TB3
15: 03214 0-4 @@ TB2
16: 03241 0-3 @@ TB4
17: 03412 0-2 @ TB5
18: 03421 0-2 @@ TB6
19: 04123 0-1 @ TB7
20: 04132 0-1 @@ TB8
21: 04213 0-1 @@ TB8
22: 04231 0-1 @ TB7
23: 04312 0-1 @ TB7
24: 04321 0-1 @@ TB8
25: 10234 1-4 @ TB9
26: 10243 1-3 @ TB10
27: 10324 1-4 @@ TB11
28: 10342 1-3 @@ TB12
29: 10423 1-2 @ TB13
30: 10432 1-2 @@ TB14
31: 12034 2-4 @ TB15
32: 12043 2-3 @ TB16
33: 12304 3-4 @ TB17
34: 12340 3-4 @@ TB18 (40 @ -> @@)
35: 12403 2-3 @@ TB19 (40 @ -> @@
36: 12430 2-4 @@ TB20 (40 @ -> @@)

This enumeration can be done by hand or using a simple program. However, the information is easier to read when rearranged:

TB1/TB2: axis 0-4, @/@@
TB3/TB4: axis 0-3, @/@@
TB5/TB6: axis 0-2, @/@@
TB7/TB8: axis 0-1, @/@@
TB9/TB11: axis 1-4, @/@@
TB10/TB12: axis 1-3, @/@@
TB13/TB14: axis 1-2, @/@@
TB15/TB20: axis 2-4, @/@@
TB16/TB19: axis 2-3, @/@@
TB17/TB18: axis 3-4, @/@@

The octahedral designators can be obtained in the same way. Although there seem to be three axes to choose from, the requirement to list the remaining 4 ligands clockwise or counter clockwise restricts the number of axes to 1. The additional ligand results in 15 possible axis combinations giving rise to 30 designators. A c++ example program generating these designators can be found here. For OH1-OH30 the following designators are obtained:

OH1/OH2: axis 0-5, @/@@
OH3/OH4: axis 0-4, @/@@
OH5/OH6: axis 0-3, @/@@
OH7/OH8: axis 0-2, @/@@
OH9/OH10: axis 0-1, @/@@
OH11/OH13: axis 1-5, @/@@
OH12/OH14: axis 1-4, @/@@
OH15/OH16: axis 1-3, @/@@
OH17/OH18: axis 1-2, @/@@
OH19/OH21: axis 2-5, @/@@
OH20/OH22: axis 2-4, @/@@
OH23/OH24: axis 2-3, @/@@
OH25/OH30: axis 3-5, @/@@
OH26/OH29: axis 3-4, @/@@
OH27/OH28: axis 4-5, @/@@

Cheminformatics & Javascript: Part 2 (Ring perception)

2010-06-12T18:09:00.001-07:00

In the second post of the javascript series ring perception will be used as an example to illustrate some of the difficulties to overcome when developing a javascript cheminformatics library. Readers only interested in ring perception can skip the first part. The ring perception section contains some useful (algorithm) tips for speeding up ring perception not only in javascript.

The topics discussed in this post result from working on jchemhub. It's still a very young project soon to be renamed to kemia.

Some javascript benchmarks

Like I said in the previous post, firefox & chrome javascript performance is impressive Unfortunately, at the other end of the spectrum there is internet explorer. The way simple things like iterating over an array are done can make a huge difference. To illustrate, here are 4 code samples for iterating over the elements of an array. The preferred way is highlighted in green.


// example 1: check arr.length each iteration
for (var i = 0; i < arr.length; i++) {
  var tmp = arr[i] + arr[i];
}


// example 2: store arr.length in a variable
for (var i = 0, li = arr.length; i < li; i++) {
  var tmp = arr[i] + arr[i];
}

// example 3: for ... in ...
for (var i in arr) {
  var tmp = arr[i] + arr[i];
}

// example 4: forEach
function doStuff(element, index, arr) {
  var tmp = arr[index] + arr[index];
}
goog.array.forEach(doStuff);

It should be obvious that both for ... in ... and forEach should be avoided at all cost. To non-javascript developers, the for ... in .. syntax might be surprising. It doesn't give you the elements but the keys to the elements. As a result, you still need the var element = elements[i]; line and there is not much benefit over a classic c-style for loop. No optimization is done when running javascript and accessing an arr.length property for each iteration, will be slower than storing the property in a local variable.

The second benchmark shows the overhead of a function call for accessing an object property.


var obj = { x: [], getX: function(index) {return this.x[index]; } };
for (var i = 0; i < 10000; i++) {
  obj.x.push({});
}

// example 1: direct access
for (var i = 0; i < 10000; i++) {
  obj.x[i];
}

// example 2: wrapper function
for (var i = 0; i < 10000; i++) {
 obj.getX(i);
}

A final benchmark shows that common knowledge from using other languages may not always apply. When creating an array of known size, using push() to expand the array is the fastest method in chrome. This is not true for internet explorer and the best compromise are the first two examples.



// example 1: Use goog.array.repeat.
var arr1 = goog.array.repeat(2, 100000);

// example 2: Use access by index
var arr2 = [];
for (var i = 0; i < 100000; i++) {
  arr2[i] = 2;
}

// example 3: Allocate using access by index.
var arr3 = [];
arr3[99999] = 2;
for (var i = 0; i < 100000; i++) {
  arr3[i] = 2;
}

// example 4: Use push
var arr4 = [];
for (var i = 0; i < 100000; i++) {
  arr4.push(2);
}

More of these benchmarks will appear in the future on this page.

Ring perception

There haven't been any posts on this blog about ring perception before. The aim of ring perception is simple, it finds the rings in a molecule. Ring perception is also required in various other fields like electrical engineering. While it may be easy to get the concept, actually implementing it (efficiently) is not straightforward. There are various algorithms to compute the rings in a graph and the results are not the same. The algorithms can be classified by their output ring sets. In this post, exhaustive ring search and the Smallest Set of Smallest Rings (SSSR) will be discussed. An interesting approach combining both algorithms is also presented.

Exhaustive ring search generates a large ring set containing all possible rings. This includes many combinations of rings. The Hanser algorithm is the fastest exhaustive algorithm. However, finding all rings in a reasonable large molecule will be slow. Fortunately, a depth or maximum ring size can be specified. For many structures, rings up to size 6 are often enough and the algorithm is really fast in these cases. Increasing the maximum ring size quickly results in unacceptable performance though. Using other programming languages such as c++ or java, it might be possible to increase the maximal ring size to make it usable for depiction. The implementation in jchemhub is ported from the MX java cheminformatics library.

The Smallest Set of Smallest Rings (SSSR) is a popular kind of ring set. The number of sssr rings in a molecule (single connected fragment) is given by Cauchy's formulla:


m - n + 1       m: number of bonds
               n: number of atoms

If the molecule contains more disconnected fragments, 1 can be replaced with the fragment count. Knowing the number of rings in the sssr set in advance can be used for various optimizations. The simplest optimization is to conclude there are no rings if the formula is equal to 0. SSSR rings have their disadvantages though. For example, SSSR rings are not canonical, the found rings will depend on the atom order of the input molecule. This would be a good topic for an additional blog post but for now the SSSR will be good enough.

There are several algorithms to compute the SSSR such as Berger's and Figueras' algorithm. The algorithm in jchemhub implements a recent algorithm that makes use of path-included distance matrices. The algorithm is based on the Floyd-Warhsall algorithm for computing the distance matrix. The element i,j in a distance matrix is the length of the shortest path (i.e. number of bonds) connecting atoms with index i and j. Two additional matrices are created at the same time containing the shortest path and the shortest path+1 between atoms i and j. From these path-included matrices, the ring candidates can be computed. The method has an O(n^3) runtime where n is the number of atoms. This isn't bad for finding the SSSR but using it on large molecules can be time consuming.

The ideal algorithm would be fast when there are only small rings but it should also detect large rings. To solve this problem, both algorithms can be combined. First the Hanser algorithm detects all rings up to size 6. From this ring set, an attempt is made to construct the SSSR using a XOR function to test if a ring is already in the SSSR. If the full SSSR is found (i.e. the number of rings in the set after selection is equal to the calculated number of rings), the molecule doesn't contain any rings larger than 6. However, if the full SSSR set is not found the real SSSR algorithm is used to find all the rings of the SSSR including the larger rings.

To evaluate the combined ring search approach, the use case is important. When doing substructure search, the average time is the most important. If an individual molecule takes some more time to process, thi is not a problem as long as the whole process is fast. The graph below shows the average time for each algorithm. A rapid decrease in performance for the Hanser algorithm can be seen. Also, the combined algorithm is 3.4 times faster than the SSSR algorithm. These tests are done in nodejs.

A second graph shows the worst case for each algorithm. This is relevant for use cases like on-the-fly depiction. Here it doesn't matter if all molecules take 50ms more as long as there are no cases that seem to take forever and make the browser unresponsive. The SSSR algorithm can take up to 6 seconds limiting it's use. The same applies for the Hanser algorithm with larger maximal ring sizes.

In the previous graphs, a simple optimization is used. All terminal atoms are progressively removed until there are no more terminal atoms. This reduces the molecule's size and therefore improves the performance of the ring perception (regardless of the algorithm). This worked great on firefox and chrome but wasn't enough for internet explorer 7. Internet explorer is still widely used and it should be supported. Combining the algorithms and the simple atom removal brought IE7's time for perceiving rings in a 534 atom molecule down to around 10 seconds. To get it down to the current 500ms another trick was needed.

Ring membership for atoms can be computed using an O(n) algorithm. A breath-first iteration is done and visited atoms and bonds are tracked together with the atom depth. When a exploring the bonds from an atom to continue iterating, there will be cases where the bond isn't visited yet but the atom it connects to is. In these cases, the bond is a ring-closure. When a ring-closure is found a backtrack can be performed to set isInCycle properties for the atoms in the ring until the paths converge. Using this atom ring membership information, the molecule can be divided in merged ring systems (i.e. rings with at least one atom in common are part of the same ring system). Running ring perception on these smaller molecules makes a huge difference when using an O(n^3) algorithm. This together with some IE7 specific performance benchmarks allowed the algorithm to be efficient in IE7 too. The worst case for the 23000 molecules in nodejs is now only 49ms.

In conclusion, it is possible to make cheminformatics software using javascript, even for Internet Explorer. Doing so requires a clever combination of algorithms and programming though. The next bottleneck is rendering, VML in IE is slow...

Cheminformatics & Javascript: Part 1

2010-05-21T05:34:00.000-07:00

Richard Apodaca already made several posts on this topic on his blog. After experimenting with javascript for the past few weeks, I agree with him that javascript could have a role in future web applications. Exciting things happening in browser land. For example, javascript performance has increased significantly and HTML5 features like canvas are becoming more mainstream (i.e. Firefox, Chrome, ...). Even webGL is already supported by bleeding edge firefox/chrome.

This is the first post in a new series about cheminformatics using javascript. Since javascript is not java, I thought it would be a good idea to introduce some features of the javascript language. More interesting cheminformatics examples will follow.

Object literals & JSON
In javascript, you can create objects using object literals. This is also called JavaScript Object Notation (JSON). To create an empty object:


var a = {};

We can also create an object with properties and functions:


var atom = { element: "C", x: 4, y: 5, bonds: [3, 4], getSomething: function(){return 42;} };

atom.element // = "C"
atom.bonds // = [3, 4] <- this is an Array atom.getSomething(); // = 42

Javascript is a very dynamic language, so you can always add new properties and functions later.


atom.atomicNumber = 6;
atom.getAtomicNumber = function() {
 return this.atomicNumber;
}

The object literal or JSON format is also used to encode data when doing AJAX. While the X in AJAX stands for XML, the latter can be replaced by JSON. Many web API's return data in JSON (e.g. Google APIs).

Namespaces
There are no namespaces in javascript. However, to avoid name conflicts with other libraries, it is always a good idea to hide private parts of your code. This can be accomplished by using an anonymous function and directly calling it:


(function() {

 var hiddenVariable = 0;

 function hiddenFunction() {
   return 34;
 }

 window.visibleFunction = function() {
   return hiddenFunction();
 };

})();

Cross-domain AJAX
A normal AJAX request (i.e. same-domain) is done using the XMLHtmlRequest object. For security reasons, the browser doesn't allow the XMLHtmlRequest object to do cross-domain requests. To get around this, a technique called JSONP is used:


function myCallback(data) {
 ...do something with data...
}

function xAjax(url) {
 var head = document.getElementsByTagName("head")[0];
 var script = document.createElement("script");
 script.type = "text/javascript";
 script.src = url;
 head.appendChild(script);
}

xAjax("http://otherdomain.com/app1?smiles=CCC&callback=myCallback");

The xAjax function creates a new script element inside the head element. The src property is set to the url to load. The callback parameter in the url is used by the server to generate valid javascript:


myCallback({ ...JSON data... });

As a result, the browser will execute the myCallback function when the new script element is loaded and we can do something with the data.

Closure
Closure is a very cool javascript feature. In the example below, the private1 function has access to the outer function variables, even after returning from the outer function. This allows private1 to access the foo property. If self.foo is replaced by this.foo, the property is not found since this refers to the private1 function. This usage of closure (i.e self = this) is often found in javascript code.


function outer() {
 this.foo = 0;
 var self = this;

 function private1() {
   return self.foo;
 }

 this.public1 = function() {
   return private1();
 };
}

Javascript libraries
As for most languages, there are many libraries for javascript. Some of the better ones for general purposes are Prototype, jQuery and Dojo. Google also released their javascript closure tools. The tools include a compiler to minify javascript code. JChemHub, a new javascript cheminformatics library uses it extensively.

Detecting stereogenic units: An alternative method

2010-03-21T08:22:00.000-07:00

In case you already saw part of this post: I accidentally pressed publish before finishing. Sorry...

The detection of stereogenic units in a molecule isn't always a trivial task. This is illustrated by two previous blog posts (post 1 post 2) where the stereogenic units do not have 4 topological different ligands. The posts also list the set of rules that enable the detection of these stereogenic units as can be found in the referenced paper[1]. While there are only 3 rules, the amount of code needed to implement them is big. The rest of this blog post will describe an alternative method to detect these stereogenic units from a paper by the same authors[2].

A potential stereogenic center is any atom with at least 3 heavy atom neighbors. For some atoms (e.g. nitrogen), the barrier of inversion is low and these atoms are not stereogenic centers. The definition for a tetrahedral stereogenic atom is now:

A potential stereogenic centre really is a stereogenic centre if there exists no automorphic permutation causing an inversion of the configuration of only the potential stereogenic centre under consideration.

If there exists at least one automorphic permutation causing the inversion of the configuration of only the stereogenic centre under consideration, then the potential stereogenic centre can be a stereogenic centre if the number of topologically equivalent neighbours (ligands) of potential stereogenic centre is less than or equal to the number of configurations of these ligands.

At first sight, the new definition is much simpler compared to the 3 rules. However, when it comes to implementing this, there is a catch. The second part of the definition requires the number of configurations of a topologically equivalent ligand to be computed. This could eventually leads to code comparable to the previous rules. Nevertheless, the definition is brief and easy to remember. The need for the automorphic permutations limits it's manual use though.

To illustrate how the definition works, a few examples will be covered next. The first example is the most simple case where there are 4 different ligands.

For these cases, there exists no automorphims that only invert the center. In fact, for molecules like these there is only the identity permutation. The next example illustrates where an automorphism exists that inverts the potential stereogenic center only.

In both cases, there will exist an automorphism that will invert the center (i.e. exchange the two equivalent methyl groups). Therefore, both potential stereo centers are not stereo centers.

A more interesting example is 1,4-dihydroxy-cyclohexane. In this case, there are 4 automorphisms but none of them invert only center 1 or 4. 2 of the 4 permutations invert both center 1 and 4. Atoms 1 and 4 are therefore stereogenic centers.

The second part of the definition correctly identifies the central stereo centers in the following molecules. In all 3 cases, the number of configurations for the ligand is greater than the number of equivalent ligands.

The number of equivalent ligands are 3 for the two structures on top and 2 for the bottom structure. For the top two, the number of configurations is 4 (>= 3). Likewise, the number of configurations for the bottom structure is 2 (>=2). The central stereogenic center is successfully found in the 3 cases.

Finally, two examples where the number of configurations in the equivalent ligands is lower than the number of equivalent ligands. In both cases, there are 3 equivalent ligands and only 2 configurations for the ligand (2<=3). The central atoms are therefore not stereogenic since at least one configuration will be duplicated.

From the last two figures, it should be clear that counting the number of configurations in a ligand is not the same as the number of stereo centers times two. The paper doesn't really go into the details on how to do this. However, another concept introduced by the paper might be useful here: a set of interdependent stereo centers. These sets can be found be inspecting the automorphisms. A set can contain one or more stereo centers and these sets are the independent stereogenic units in a molecule. Some examples showing the sets of interdependent stereo centers in different colors should make this easier to understand.

For all examples given above, the number of configurations is the same as the number of interdependent sets of stereo centers in the ligand times 2. Unfortunately, this is not true in all cases. Below is an examples where this would not work.

Until I find a good way to count the number of configurations, the new definition is not preferable over the old rules. I'll keep looking though and will post my findings. Suggestions are always welcome of course.

[1] M. Razinger, K. Balasubramanian, M. Perdih, M. E. Munk,
Stereoisomere Generation in Computer-Enhanced Structure Elucidation,
J. Chem. Inf. Comput. Si. 1993, 33, 812-825
pdf
[2] M. Perdih, M. Razinger, Stereochemistry and Sequence Rules: A Proposal for Modification of Cahn-Ingold-Prelog System, Tetrahedron: Asymmetry, 1994, Vol. 5, No. 5, 835-861 http://dx.doi.org/10.1016/S0957-4166(00)86237-0

2010-03-09T17:14:00.001-08:00

It has been 6 months since my last post on stereochemistry. At the end of the post, the automorphisms for inositol were listed. These will now be used to enumerate all possible stereoisomers for a structure. Most of this post is only a summary, for details the Razinger paper can be consulted.

The first use of the automorphisms is to construct stereoindex vectors. These vectors contain 1 for all true stereocenters. Para-stereocenters are assigned -1 if the permutation interchanges equivalent atoms and 1 otherwise (see this post for examples of para-stereocenters). The automorphisms, or permutations in general can be represented as a matrix which can be multiplied by the stereoindex vectors. The resulting matrices are called signed permutation matrices and they are the key to revealing identical stereoisomers. For performance reasons only the stereogenic atoms are considered when constructing the signed permutation matrices. For example, inositol containing 6 tetrahedral stereocenters will have 6x6 signed permutation matrices.

Before the enumeration can begin, another matrix called the stereoparity matrix is needed. This matrix has n columns where n is equal to the number of stereocenters. There are 2^n rows in the matrix, one for each possible combination of the values 1 and -1. This is best illustrated with this example for 3 stereocenters.

The enumeration now begins by selecting row 1 from the stereoparitymatrix. This row is multiplied by each of the signed permutation matrices. The resulting row vectors are compared with the rows from the stereoparity matrix. Matching rows are redundancies of the same stereoisomer and are not considered in further steps. The identified stereoisomer can still become an enantiomer or diastereomer. This is determined by inverting row 1 and multiplying this with the signed permutation matrices. If any of the resulting vectors match a non redundant row, an enantiomers pair is found. If all resulting vectors match already redundant rows, a diastereomer is found. The same steps are repeated until all rows are assigned to a stereoisomer.

Determining the number and kind of stereoisomers alone isn't of much use though. A much more useful feature would be to generate each stereoisomer. For each stereoisomer, any of the rows assigned to it can be used. These rows can now be considered as parities for the stereogenic units in the molecule. The input molecule's data structures can be modified to match the parities in the row. Using inositol as example once more, the 9 possible stereoisomer smiles are:

The smiles are in fact canonical smiles. The information obtained by enumerating stereoisomers can also be used for this purpose. However, this time we are only interested in the stereoisomer specified by the input molecule. The redundant set of parities associated with it can be considered to be canonical candidates. Since the stereoparity matrix contains 2^n rows, every possible candidate for a given stereoisomer is guaranteed to be found. Canonical smiles for molecules like the ones below now becomes possible.

The current code can be found here.

PS: Sorry for the delay...

Avogadro's packmol plugin

2010-03-03T13:01:00.000-08:00

When running molecular dynamics (MD) simulations, one often needs to create an initial input structure. It is important that none of the atoms are too close to prevent the coordinates from "blowing up". Many MD programs provide their own utilities (e.g. gromacs' genbox) to create basic systems. For example, a protein with solvent (water) and counter ions. The counter ions are important to make the system charge neutral for algorithms like Ewald summation or the faster PME. While these utilities are often enough, creating complex input files is not always easy.

Luckily, there is a program dedicated to this task: Packmol. It uses simple text input files to define a system. The syntax is very easy to learn but writing the files by hand is a tedious task. Many of the required calculations can be automated and I started writing the Avogadro packmol plugin to do just that. The plugin should allow novice users to create complex systems and should also allow experienced users to make small adjustments (e.g. adding keywords, parameters).

The work flow of the plugin is as follows:
- Configure parameters in a wizard:

The plugin allows you to browse for all required files. This allows all files to be copied to a temporarily location to make sure all files are found. The screen shot shows the wizard for the most basic systems containing solvent and optionally solute.

- Clicking "Generate input file" shows the generated input file

If a new structure is added in the Text Mode tab, the plugin will ask for the file's location when the user clicks run. The syntax highlighting is there to prevent typing errors, it really helps :-)

- Clicking "Run Packmol" starts the program and shows the progress in the Output tab:

Once the program finishes, the resulting file is opened in Avogadro.

The result from the input file seen in the screen shot looks like:

One of the most important calculations required is the number of molecules required to fill a given volume. Given the density, molecular weight and volume, the number of molecules can easily be calculated using Avogadro's number. (This is the same calculation the volume guesser on packmol's site uses AFAIK)

The next step towards more complex systems would be to add solute and counter ions. I'm not completely sure how to subtract the volume taken by solute from the solvents volume. Computing the solute's volume might be an option. Ideas are welcome. :-)

Lipid bilayers are the next kind of system that can be created with packmol. This requires some more knowledge of the packmol syntax but it is still human readable. Before continuing, it might be good to review the bilayer input example from the packmol website. It has comments explaining how to the constraints are used to define the system. The wizard dialog currently looks like:

Currently only 1 "polar solvent" type and 1 "lipid" type are allowed but the idea is that this could allow for mixtures of lipids and so on.

To orient the lipids in the correct way, constraints are needed on specific atoms. This is currently a limitation since a user still has to specify the atom indexes for the polar head and lipophilic tail. In the screen shot below, the indexes have already been filled in.

The required indexes can easily be determined by opening the lipid structure file and enabling the label display type in Avogadro though.

The result from the input above (rendered with POV-Ray):

The next logical step would be to insert a membrane protein between the lipids. This requires the orientation and position of the protein in the membrane to be known. From the protein database files, there is no easy way to determine this. In many cases, the orientations of proteins in membrane (OPM) database provides exactly this. From their website:

The original PDB coordinates are transformed so that the calculated membrane normal coincides with the Z axis, and the origin of coordinates corresponds to the middle of membrane.

This is how the plugin generates the lipid bilayer which allows the protein to be inserted with only a few lines added:

Again, the volume of the protein isn't considered yet. Using slightly lower numbers should work for now. After waiting a while the result should look like:

The source code can be found here.

(The files from the OPM database contain dummy atoms marking where the membrane should be. Make sure to remove these if you try to replicate this)

The automorphism group & bliss

2009-09-29T10:15:00.000-07:00

In the previous post, an algorithm for computing the automorphism group was given. Although it drastically reduced the number of permutations to be checked, it was still much to slow for use in cheminformatics systems. A real world example is inositol where there are 2 orbits of 6 atoms. This example was also used for a stereochemistry-puzzler on Depth-First (answer). Since it is a good example (i.e. it contains "large" orbits) to test the performance and correctness (there are only 9 stereoisomers) of an implementation, I'll be using it for the following posts too.

While 6! * 6! = 518400 is already way faster better than 12! = 479001600, computing half a million matrix equations still takes a considerable amount of time. The actual time (seconds to minutes) depends on implementation details but the method is generally inefficient.

A faster set of algorithms to solve this problem are known as backtracking algorithms. Some well known examples include VFLib (for isomorphisms in general) and nauty (automorphisms). Another, relatively new algorithm called bliss (paper available from their site, MIT license) came to my attention. It is a very small library which computes a set of generators of the automorphism group. From these generators, the full automorphism group can easily be constructed. Using the bliss API was easy and straightforward (see permutation.cpp below). Using inositol as example again, bliss gives 2 generators:

Generator: 1 6 5 4 3 2 7 12 11 10 9 8
Generator: 2 3 4 5 6 1 12 7 8 9 10 11

Next the generators' inverses and all possible products can be added. This is repeated until all automorphisms are found. [note: the matrix equation doesn't have to be used, the inverse and products are always part of the automorphism group too -- don't forget the identity permutation :-) ]

The final 12 automorphisms constructed from the given generators are :

1 2 3 4 5 6 7 8 9 10 11 12
1 6 5 4 3 2 7 12 11 10 9 8
2 3 4 5 6 1 12 7 8 9 10 11
6 1 2 3 4 5 8 9 10 11 12 7
2 1 6 5 4 3 12 11 10 9 8 7
6 5 4 3 2 1 8 7 12 11 10 9
5 4 3 2 1 6 9 8 7 12 11 10
4 3 2 1 6 5 10 9 8 7 12 11
3 2 1 6 5 4 11 10 9 8 7 12
5 6 1 2 3 4 9 10 11 12 7 8
4 5 6 1 2 3 10 11 12 7 8 9
3 4 5 6 1 2 11 12 7 8 9 10

This whole computation takes .019s here :-D

The code for the example can be found here: permutation.h & permutation.cpp (requires OB trunk + libbliss)

In the next post, these automorphisms will be used to construct the signed permutation matrix to finally identify identical stereoisomers.

Molsketch: Reaction & Mechanism arrows

2009-09-16T11:14:00.001-07:00

Here is a short screencast showing some of the new features in Molsketch. First, there are the reaction and mechanism arrows. While there are only two major types of arrows, double clicking an arrow in move mode will bring up a dialog that allows you to change the arrow's appearance.

The second new feature is the .msk file format. Since it is actually cml, it will likely be renamed to cml later. Apart from the molecules, the format also saves the atom/bond colours and both types of arrows.

If I find some time in the near future, I will make a release and build a windows installer so windows users can try it out without going through the difficult build process.

More para-stereocenters, permutation groups and the automorphism group

2009-09-16T03:41:00.000-07:00

Rich Apodaca recently blogged about this stereochemistry related series of posts. He correctly pointed out that the automorphism group lies at the heart of the Razinger paper on stereoismer generation. However, before the automorphism group can be used, all stereocenters have to be found.

In a previous post I briefly described rule 1 and gave some examples. Below are more examples of these para-stereocenters identified by applying rule 1. As you will notice, most of these structures contain no stereocenter and are used to avoid false positives in the unit tests.
While rule 1 handles interdependant stereocenters for cycles, both rule 2 and rule 3 handle individual stereocenters and can be applied directly (rule 1 is applied recursively). The rules can schematically be summarized as:

2) tetracoordinate carbon is para-stereocenter if it attached:
2a) 1 or 2 pairs of identical symmetry classes (e.g. 1123, 1122), each pair contains at least :
2a') >= 1 true-stereocenter (have 4 different symmetry classes) -- or --
2a") >= 2 para-stereocenters (from rule 1)
2b) 3 or 4 identical symmetry classes (e.g. 1112, 1111), duplicated symmetry class ligand contains:
2b') >= 2 true-stereocenters -- or --
2b") >= 2 separate assemblies (merged rings) with each >= 2 para-stereocenters

3) double bond is para-stereocenter if terminal with identical symmetry classes contains at least:
3') >= 1 true-stereocenter
3") >= 2 para-stereocenters

Some examples of rule 2 & 3:

Now that all stereocenters are found, a (very) short introduction to permutations and permutation groups will help us get to the automorphism group (This can be skipped if you are familiar with it). A permutations is simply a reordering of a set of numbers. Although there are several ways to represent permutations, the matrix representation is most useful here. A permutation group is unsurprisingly a group or collection of Permutations. (see Razinger paper + wikipedia Permutation Group & Permutation Matrix)

The automorphism group is a subgroup of the Sn group containing all permutations. Since elements are not allowed to repeat, there will be n! permutations in Sn where n is the number of atoms. In general, the term "auto" means self and permutations belonging to the automorphism group change the atom labeling in such a way that the molecule isn't changed (i.e. no bonds broken or neighbors changed, the molecule remains itself). Mathematically this can be expressed using the adjacency matrix (A) and the permutation matrix (P) in the equation:

While a brute force check of all n! permutations would work for small molecules, it is not efficient enough for use with most drug-like molecules. To limit the number of automorphism candidates, symmetry classes can be used to provide an initial partitioning. Atoms with identical symmetry class are said to be in the same orbit. Only permutations within orbits are allowed (i.e. to be part automorphism group). These intra-orbit permutations can easily be computed in c++ using std::next_permutation. These intra-orbit permutations are all candidates for the automorphism group (G) and will be added to G if the equation holds (duplicates are ignored, especially the identity permutation will be duplicated many times).

Next, all inter-orbit permutations still need to be considered. This is achieved by "multiplying" the permutations (I used matrix multiplication but other methods are possible) and this gives rise to a new set of candidates which are also checked using the equation. All permutations in G now make up the automorphism group.

Below are 3 examples to illustrate how these "obrits" limit the number of permutations. The structure on the left, has 2 orbits, each containing 2 atoms. This structures' Sn group has 7! or 5040 permutations. Using the orbits we can reduce the number of candidates to 2! x 2! = 4. The structure in the middle also has 2 orbits but one orbit has 2 atoms and the other has 4 (6! = 720 -> 2! x 4! = 48). The structure on the right has 3 orbits, each containing 2 atoms (9! = 362880 -> 2! x 2! x 2! = 8).

While OpenBabel doesn't have code to compute the automorphism group yet, here is some proof of concept code which will be integrated later. It also contains a method to compute all permutations of the Sn group which can be used in combination with small molecules to validate the automorphism group generating algorithm.

Once the stereocenters are identified and the automorphism group has been found, this information can be used to eventually create a stereoparity matrix. The details and code will follow in a later post but knowing that duplicate rows in this matrix represent the same stereoisomer should be motivating :-) (for both stereoisomer generation & canonicalization)

2009-09-11T13:52:00.000-07:00

As promised, here are some molecules with stereogenic units that are considered "hard" to detect. As can be seen from the image, the previous definition of atoms having four different neighbours (by comparing symmetry classes) doesn't always work.

The term para-stereocenter is used and adopted for these cases.[1] In general para means ressemble, so para-stereocenters ressemble true-stereocenters but there are differences. The most important difference is that para-stereocenters depend on each other to be stereogenic. In fact this is rule 1 from the Razinger paper:

rule1: para-stereocenters are assigned to be stereogenic if there are at least two of them in a ring assemblage (single ring, sharing atom/bond, bridged, ...). In OpenBabel code this is called mergedRings but it basically means the same thing.

This simple rule allows all stereocenters to be found in the given examples. There is also a rule 2 and 3 but these will be covered in the following posts. Still need to commit the code for rule1 later tonight.

[1] M. Razinger, K. Balasubramanian, M. Perdih, M. E. Munk,
Stereoisomere Generation in Computer-Enhanced Structure Elucidation,
J. Chem. Inf. Comput. Si. 1993, 33, 812-825
http://www.mcs.csueastbay.edu/~kbalasub/reprints/282.pdf

Symmetry classes: How to get them and what to do with them

2009-09-11T06:38:00.000-07:00

For a long time, I've wanted to make a series of blog posts about my stereochemistry adventures. This is post will mainly be about symmetry but we'll get to it's relationship to stereochemistry at the end. In this post (and the following), symmetry will always mean topological symmetry (or graph symmetry). This doesn't imply any 3D orientation or special symmetry. For example, below is an image with symmetry classes assigned to the atoms (click to enlarge):

Perhaps, a better title would have been "... and what to do with them in the first place" but this seemed too long. However, to motivate my readers, I'll start by answering the second question. In OpenBabel we currently use the symmetry classes to identify stereogenic units (i.e. tetrahedral atoms, cis/trans double bonds, ...). For example, if an atom contains neighbors with 4 distinct symmetry classes, the atom is considered to be chiral or stereogenic. There are exceptions to this which will be discussed in later posts.

A second and perhaps more useful everyday example is to determine how many peaks there should be in an 13C-NMR spectrum. Since carbon atoms in the same symmetry class will be in a similar environment, their chemical shifts will be similar. In the example below, you can clearly see that there are 8 different carbon symmetry classes. Two classes have more than one atom and correctly show the symmetry (6 & 8). The spectrum taken from wikipedia doesn't hve a hight enough resolution but the peaks should be there :-)

For H-NMR this is slightly more complex due to stronger coupling.

Computing the symmetry classes is actually not that hard. The most common method use is the Morgan's Extended Connectivity (EC) algorithm:

Assign each atom it's degree as starting point

Iterator over the atoms and assign the summed values for each atom's neighbors to it.

Repeat untill the number of distinct classes remains the same for two iterations.

The degree is the Graph Invariant (GI) here. In reality OpenBabel uses another GI with more discriminating power to avoid isospectral points (i.e. points with the same symmetry class which are not symmetric) but the principle remains the same. The symmetry classes are also renumbered so they are all within 1 ... n (# symmetry classes).

Next post will be about the identification of stereocenters for symmetrical molecules.

Avogadro 0.9.8

2009-09-11T05:56:00.000-07:00

Avogadro 0.9.8 was released yesterday and the windows installer including the python bindings and dependencies is available from sourceforge. There are no big new features but several bugs are fixed including some problems with translations. Accents were not displayed correctly and it had something to do with the encoding (UTF-8) of the files but we found a solution and Avogadro looks fine in French.

Molsketch

2009-09-07T10:53:00.000-07:00

The last 2 months Nicola Zonta and I have been working on improving Molsketch. While development is still at an early stage and bugs may be present, great progress has been made. Nicola released a version called Zem to go together with the 3D Zodiac program.

If you want to try it out, a windows installer is available on zeden.org. Linux users can get the current development code for molsketch at my github page. NOTE: requires OpenBabel trunk: svn co https://openbabel.svn.sourceforge.net/svnroot/openbabel/openbabel/trunk.

The zem release doen't include the latest features yet but here is a screencast made with the current github version of molsketch. The big difference with zem 0.3 is the ItemPlugin concept which is illustrated in the screencast. The idea is that more of these plugins can be added. While I don't have an example of this yet I intend to make it possible to use REST API's in this way.

Feel free to make suggestions, propose ideas, file bugs, contribute, etc.

SourceForge's Community Choise Awards

2009-07-19T18:05:00.000-07:00

I'm pleased with Avogadro's nomination for SourceForge's Community Choise Awards. First of all, thanks to all people who have voted to get us nominated. Secondly, thanks to those who have already voted for the final. For those who haven't voted yet, voting ends July 21st. This is your last chance! :-)

Vote Now!

We've also made windows installer packages for the 0.9.7 release. These are already on sourceforge. The official realease mail etc. will follow soon.

Using the Avogadro library from python

2009-05-16T14:58:00.000-07:00

As we get closer to Avogadro 1.0, the python bindings are becoming more and more powerful. For example, take this small script (65 lines) which allows large files to be read in a separate thread. Once the file is read, clicking a title in the list will open a 3D subwindow with the molecule. Scripts like this could be used to preview various molecules, setup a very basic application for educational purposes etc.

Note: Requires Avogadro developer version (next release will probably be in two weeks)

SxsTrace, making Manifest-Hell more bearable

2009-04-28T08:52:00.000-07:00

For people comming from unix, Microsoft's side-by-side (SxS) assemblies can be annoying to say the least. The term Manifest-Hell, DLL-Hell's successor, is often used to describe various SxS problems.

It's not that we are complaining just because it is a microsoft product. The main problem is the cryptic error messages you get when you first encounter it.

The good news is that there is a tool called SxsTrace to debug these issues. It generates a detailed logfile containing all info you need. It even gives the linenumber in the manifest file where the error is found. Using SxsTrace, creating your own assemblies becomes easy :-) This makes me wonder why this tool isn't mentioned more in docs... Unfortunually there is also bad news, the SxsTrace tool is only found on windows vista (newer sxs.dll). Altough I'm happy I have vista now, I still wonder how XP users are supposed to fix these problems.

Avogadro gets new python wrappers

2008-12-06T12:55:00.000-08:00

The past week I've been very busy creating new python wrappers for Avogadro. About 70% of the classes have been wrapped and have associated unit tests to make sure the wrappings work.

We use boost::python to expose our classes to python and also use it to help with the embedding. However, there were several hurdles to overcome before the python bindings allowed us to do everything from python. In this post I'll summarize each of them and I hope this information may be useful for other people in a similar situation.

To start simple, I added to/from python array converters for Eigen, the template library we use for all our math. In python, the users can just use arrays from the numpy package to get/set atom positions, do rotations, ....

example:
>>> import Avogadro
>>> from numpy import *
>>> molecule = Avogadro.molecule
>>> atom = molecule.newAtom()
>>> atom.pos = array([1.0, 2.0, 3.0])
>>> print atom.pos
[1. 2. 3.]

Now, it gets more interesting, we also use Qt and although using SIP is an obvious option, I recently started to love boost:python and thought there had to be a way. Searching the web I found this, while the described situation may apply to some cases, the approach uses both boost::python and SIP to do the wrapping and makes SIP aware of boost::python. For Avogadro we need the reverse, all wrapping needs to be done by boost::python and our functions have to able to return and accept PyQt objects. Another search result was the documentation page of HippoDraw. They seem to have had the same problem but just used SIP. The page also mentions a discussion on the Python sig-c++ mailing list. Since these are the top results, I assume this an unexplored area. I would love to see how other projects have solved this, please let me know if you find any.

Anyway, after some hacking, Avogadro's python module is now able to return and accept real PyQt wrapped objects. For example, the Tool class is wrapped using::boost python, when you call

QWidget* Tool::settingsWidget();

you get a real PyQt object with all it's functions available etc...

>>> print tool
<Avogadro.Tool object at 0x94e06f4>
>>> print tool.settingsWidget
<PyQt4.QtGui.QWidget object at 0xb5b73d6c>

The following is also possible:

void GLWidget::setUndoStack( QUndoStack *undoStack )

>>> glwidget = Avogadro.GLWidget()
>>> print glwidget
<Avogadro.GLWidget object at 0x93f8d24>
>>> undostack = QUndoStack()
>>> print undostack
<PyQt4.QtGui.QUndoStack object at 0x956daec>
>>> glwidget.undoStack = undostack

So far so good, but there are still some things missing. For example how would you resize the GLWidget (our main view widget)? One option would be to just add boost::python definitions for the ones you think will be needed. However, seen how big Qt is, this does not scale very well. While I'd like to keep everything as simple as possible for the python users, I added a conversion function that allows us to convert a boost::python wrapped class to the relevant PyQt wrapped object. For our GLWidget class, this would be QGLWidget.

another example:
>>> print glwidget
<Avogadro.GLWidget object at 0x9da6d24>
>>> print Avogadro.toPyQt(glwidget)
<PyQt4.QtOpenGL.QGLWidget object at 0x9f167ac>
>>> Avogadro.toPyQt(glwidget).resize(640, 480)
>>> Avogadro.toPyQt(glwidget).show()

The toPyQt(...) function does not only allow us to use PyQt defined functions on our own classes, it also gives us signals/slots for free:

from PyQt4.Qt import *
import Avogadro
import sys

class MyClass(QObject):
def __init__(self):
QObject.__init__(self)

@pyqtSignature("")
def printHelloWorld(self):
print "Hello World"

app = QApplication(sys.argv) # PyQt wrapped object

button = QPushButton() # PyQt wrapped object
button.show()

molecule = Avogadro.Molecule() # boost:python wrapped object

myClass = MyClass() # real python class

# connect the signals
QObject.connect(button, SIGNAL("clicked()"), Avogadro.toPyQt(molecule), SIGNAL("updated()"))
QObject.connect(Avogadro.toPyQt(molecule), SIGNAL("updated()"), myClass, SLOT("printHelloWorld()"))

sys.exit(app.exec_())

Everytime you press the button, the signal is routed trhough a boost::python wrapped class and "Hello World" is printed :-)

That's all for now. More posts will follow demonstrating more cool python stuff and don't forget to check out the eye-candy Marcus has recently added to Avogadro.

Upcomming OpenBabel 2.2 release

2008-04-30T12:39:00.000-07:00

OpenBabel 2.2 will be released soon. In the past 9 months since the 2.1.1 release, the force field code has really improved.

Two new force fields are available in 2.2. The first is MMFF94, a force field that can be used for drug like molecules and ligand-protein interactions. Although this force field could be used for proteins themselves, protein specific ones like AMBER would be more suited, especially for performance reasons. The implementation is validated using the MMFF94 validation suite and the results are available here. Some special cases remain problematic such as bridged rings and emperical rules when parameters are missing. However, the errors are so small and rare that it is suited for most applications.

The second new force field, implemented by Geoff, is UFF. This force field has parameters for the entire periodic system and is therefore a valuable addition to the list. Bob Hanson has also ported UFF from openbabel to JMol. Towards OB 3.0 a biomolecule force field like AMBER will probably be added.

Another great improvement is the performance. The numbers on the updated wiki page speak for themselves. Minimizations run about 100 times faster!

This will be the first release with the OBBuilder class. This allows openbabel to generate 3D coordinates with a fragment based approach. The database with fragments was derived from Zinc database by Geoff. Stereochemistry will not be taken into account for now, but it's a great start.

2008-01-01T13:33:00.000-08:00

Hi,

First of all, I would like to welcome you to my blog. I'm a student farmaceutical science at the University of Antwerp. In my free time I work on the open source projects OpenBabel (http://openbabel.sourceforge.net/) and avogadro (http://avogadro.sourceforge.net/). This blog will mainly be about these two projects and their progress.