[Editor’s note: Piotr Janik (janikpiotrek@gmail.com) was a Google Summer of Code 2012 student at the Concord Consortium and is now a consultant working on our Next-Generation Molecular Workbench.]

Some time ago we described the core engine used in Molecular Workbench and our attempts to speed it up. At that time we focused mainly on the low-level optimization connected with reducing the number of necessary multiplications. This promising early work encouraged us to think even more about performance.

We next reviewed existing algorithms in the core of the molecular dynamics engine. To make a long story short, atoms interact with each other using two kinds of forces:

• Lennard-Jones forces (repulsion and short-range attraction)
• Coulomb forces (electrostatic and long-range attraction)

Atomic interactions are pairwise, meaning that we have to calculate forces between each pair of atoms while using the basic, naive algorithm. Having n atoms, we must perform about n^2/2 calculations. “The Big O” notation can be used and the computational complexity can be described as O(N^2), which means that the execution time of calculations grows very fast as the number of atoms used in the simulation increases. This is definitely an unwanted effect, but fortunately there are ways to reduce the complexity.

Solutions are different for short-range and long-range forces, so let’s start with short-range. “Short-range” means that atoms interact only while they are quite close to each other. Let’s use rCut as a symbol for the interaction maximum distance. So, one obvious optimization would be to limit calculations to pairs of atoms that are closer to each other than rCut. How? There are two popular approaches—cell lists and Verlet (neighbor) list algorithms.

The cell lists algorithm is based on the concept that we can divide the simulation area into smaller boxes or cells. Each cell dimension is equal to the maximum range of interaction between atoms—rCut. So, while calculating interactions for a given atom, it’s enough to take into account only atoms from the same box and its closest neighbors. Atoms in other boxes are too far to interact with this atom. This is both simple and effective, reducing computational complexity to O(N)! Note that it’s C * O(N) with a pretty significant C, unfortunately.

However, while calculating interactions between atoms in neighboring cells, still only 16% of atoms that we take into account are interacting! This is a waste of resources and where we find room for further optimizations. So, what about creating a list for each atom, which contains only atoms actually interacting with it? This Verlet or neighbor list algorithm as it’s called works well. The only problem is that we have to be smart about updating these lists, as atoms constantly change their position and, thus, their “neighborhood.” We can slightly extend these lists to also include some atoms outside the area of interaction. So each list should include atoms closer than rCut + d from the given atom, where d defines a buffer area size. Because of that, lists need to be updated only when the maximum displacement of some atom, measured since the moment of the previous lists update, is bigger than d. If it’s smaller, neighbor lists are still valid. Lists can be updated using the normal, naive algorithm (which still leaves the complexity O(N^2)), or even better, using the cell lists algorithm presented above! This ensures complexity O(N) and greatly reduces inefficiencies of the cell lists approach.

We’re also working on long-range forces optimization. Since we can no longer use the assumption that atoms interact only when they are close to each other, we can’t rely on the optimization strategies above. The algorithms are now more complicated. The problem of the electrostatic interaction is akin to a problem of gravitational interactions (called N-body problem), popular in astrophysics. One of the most common algorithms for speed-up of such calculations is the Barnes-Hut algorithm. We tried to implement it, but the overhead connected with creating additional data structures was bigger than potential performance gains. The reason is that the number of charged atoms we use in our models is too small to see the advantage of such an approach. As a result, we left our naive algorithms for long-range interactions, which perform better due to their simplicity.

However, we successfully implemented both short-range optimizations in Next-Generation Molecular Workbench and the results are spectacular. The speed-up varies from 20% for really small models (where the number of atoms is less than 50) to 700% for bigger ones (where the number of atoms is about 250). This is the really significant improvement and made complex models usable. As you can see, conceptual, algorithmic optimizations really matter!

We’re still thinking about further optimizations, both low level and algorithmic. Stay tuned as the Next-Generation MW is getting more and more computational power!