Ganesh Swami

Named after Anton van Leeuwenhoek, the father of microbiology, this is the latest from D. E. Shaw Research. This paper describes the massively-parallel hardware part of their super software-hardware technology. I’ve previously described their software here.

Scheduled for completion at the end of 2008, Anton is built with the aim of doing millisecond simulations. In contrast with approaches like Folding@Home, they are interested in a single trajectory. Folding@Home patches together many small trajectories that are more feasible on a regular desktop PC. Their initial configuration will have 512 nodes.

Many of the ideas implemented in their software simulation engine Desmond were developed with Anton in mind. For long range electrostatic interactions, they use a method called k-space Gaussian split Ewald, reducing the computational burden from O(n^2) to O(n log(n)). Another advantage is that the same kernel is used for force interpolation and spreading saving real estate on the chip.

As any molecular dynamics aficionado knows, about 90% of the computational time is spent in computing non-bonded forces. With this is mind, they’ve done everything possible to speed up the process devoting much of the ASIC space. One of their ideas that struck me was the explicit computation between all pairs and then to go back subtract correction terms. This leads to a clear separation in the accounting between what particles are required for what forces – a major win for scalability.

They get major speedups by building a specialized hardware datapath and control logic tuned for common molecular dynamics communication patterns. They even have dedicated support to accumulate forces for reducing latency. I was previously excited by the addition of the dot product instruction haddps to recent generations of processors. Anton does dot products in hardware. Dot products can be used everytime you’d want to “squish” a vector.

Anton does all computation in fixed point arithmetic, which is faster and more accurate for the same silicon area (their claim.) Molecular dynamics is expensive in computation and communication, but not so much in memory. All of their data for a 25,000 particle system fits in L1 cache.

One of the courses in Engineering at SFU (ENSC 250) deals with designing a RISC processor from the ground up. I vaguely remember building branch prediction logic in VHDL, not really knowing if I’d ever use that knowledge again. If only I’d known…

Anyways, I’m missing out on some of the finer details of their technology. Go read the paper for the real deal.

I walked into a gold mine: a list of thirteen classic papers in applied math. The webpage has other key papers, but this list should keep me occupied for some time. I’m reproducing the list for my own reference:

• J.W. Cooley and J.W. Tukey, “An algorithm for the machine calculation of complex Fourier series,” Math. Comp., 19 (1965), pp. 297–301.
• R. Courant, K. Friedrichs, and H. Lewy, “On the partial difference equations of mathematical physics,” IBM J. Res. Develop., 11 (1967), pp. 215–234.
• A.S. Householder, “Unitary triangularization of a nonsymmetric matrix,” J. Assoc. Comput. Mach., 5 (1958), pp. 339–342.
• C.F. Curtiss and J.O. Hirschfelder, “Integration of stiff equations,” Proc. Nat. Acad. Sci. U.S.A., 38 (1952), pp. 235–243.
• C. de Boor, “On calculating with B-splines,” J. Approximation Theory, 6 (1972), pp. 50–62.
• R. Courant, “Variational methods for the solution of problems of equilibrium and vibrations,” Bull. Amer. Math. Soc., 49 (1943), pp. 1–23.
• G. Golub and W. Kahan, “Calculating the singular values and pseudo-inverse of a matrix,” J. Soc. Indust. Appl. Math. Ser. B Numer. Anal., 2 (1965), pp. 205–224.
• A. Brandt, “Multi-level adaptive solutions to boundary-value problems,” Math. Comp., 31 (1977), no. 138, pp. 333–390.
• M.R. Hestenes and E. Stiefel, “Methods of conjugate gradients for solving linear systems,” J. Research Nat. Bur. Standards, 49 (1952), pp. 409–436.
• R. Fletcher and M.J.D. Powell, “A rapidly convergent descent method for minimization,” Comput. J., 6 (1963/1964), pp. 163–168.
• G. Wanner, E. Hairer, and S.P. Nørsett, “Order stars and stability theorems,” BIT, 18 (1978), no. 4, pp. 475–489.
• N. Karmarkar, “A new polynomial-time algorithm for linear programming,” Combinatorica, 4 (1984), no. 4, pp. 373–395.
• L. Greengard and V. Rokhlin, “A fast algorithm for particle simulations,” J. Comput. Phys., 73 (1987), no. 2, pp. 325–348.

I came across a paper on the pre-print archive: “Introduction to protein folding for physicists.” Proteins are polymers of amino acids and the sequence of these amino acids completely determine the native three dimensional structure. The human genome project gives us the sequences of all proteins in the human body. In theory it should be possible to study diseases like Alzheimer’s that are caused by the misfolding. Right?

Unfortunately, the combinatorial nature of the problem makes it impossible to do a brute force search. This paper is a decent introduction to the vocabulary of molecular biologists. The abstract:

[...] we provide an exhaustive account of the reasons underlying the protein folding problem enormous relevance and summarize the present-day sta- tus of the methods aimed to solving it. We also provide an introduction to the particular structure of these biological heteropolymers, and we phys- ically define the problem stating the assumptions behind this (commonly implicit) definition. [...]

The author recognizes the inadequacies of implicit solvent models. He also asks a key question: is kinetic information important for determining which minima is chosen. This was a similar question that I had while doing my docking project. Most docking algorithms right now only use static information.

References:

• Introduction to protein folding for physicists, by Pablo Echenique. arXiv:0705.1845v1 [physics.bio-ph]

• Picture courtesy of Wikipedia

Another interesting paper came out last week on the preprint archive, “A matrix generalization of Euler identity.” The Euler identity reads,

and is traditionally derived from the Taylor expansion of sine and cosine. In this paper, the author generalizes the Euler identity to matrices. He does this by introducing a complex matrix, known as the imaginary unit matrix

Just like , .

He then proves the Euler identity for these matrices as:

for all real x.

A friend brought this paper(pdf) to my attention “A common, avoidable source of error in molecular dynamics integrators” recently published in the Journal of Chemical Physics by D. E. Shaw Research. I had previously written about their award winning paper at SC06.

In constrained molecular dynamics simulations using some of the most popular molecular dynamics codes, calculation of the velocities of constrained particles is based solely on the differences in particle positions during two successive time steps. This creates a numerical instability that the authors’ show to be significant in a typical single-precision floating-point simulation. They describe a simple modification that eliminates this source of instability and demonstrate that this change substantially reduces the energy drift of a sample single-precision NVE simulation.

The paper is really short (2 pages) and describes recipes for modifying constraint algorithms SHAKE, RATTLE and SETTLE. Below is a graph of the energy drift in double-precision, single-precision and modified versions of GROMACS and Desmond (that’s D.E. Shaw’s MD engine.) They simulated a 30A cube of 901 water molecules with periodic boundary conditions for one nanosecond with a one femtosecond timestep. See how bad single-precision GROMACS is for NVE simulations?