Sunday, July 31, 2011

Molecule construction and visualization website

For a couple of months now, I've been at work on a website for constructing and visualizing molecules. In-browser molecular dynamics are done with a molecular mechanics modeler based on Norman Allinger's MM2 as described in Eric Drexler's Nanosystems. This is the same set of mathematics I used in an earlier effort in the same vein called NanoCAD in 1997. Unfortunately my knowledge of chemistry and molecular modeling hasn't grown very quickly in that time. I know a bit more from my time with Nanorex, where molecular modeling was mixed with gadgets to supply external forces ("jigs" in the parlance of our program), an idea that I believe is crucial to nanotechnology design software and also to scaling molecular simulations to much larger scales. I hope to use the code from this website as a starting point for working in that area.

I have expenses like everybody else, and I'm trying to think of ways to use this website to make a little money without tarnishing its educational potential or scientific credibility. I want the website to be readily available for use in schools and universities. If I end up putting ads on the website, I hope to make them tastefully small and out-of-the-way. I've noticed that the HTML5 canvas I'm relying upon for graphics doesn't work on iPhones, iPads, or Android phones, so there's an opportunity to sell mobile apps for those platforms, using their native graphics canvases. I'm very open to ideas to bring in a little revenue without being tacky. I've put a lot of work into this, and plan a lot more.

Longer term plans include adding jigs as discussed above, maybe an interface for a force-feedback joystick so that you can find out what Brownian motion feels like, using the website to get access to much better and faster simulators, and storing your own private library of molecules. This stuff started out as Java code with the website's JavaScript being co-developed, but at some point the JavaScript development took off and I didn't make time to keep the Java up to date. So I need to do that, and then the JAR file can become a useful computational chemistry tool.

Tuesday, July 26, 2011

Molecular dynamics and force fields

In the web's early days, I wrote a Java applet to do a little bit of molecular modeling in your web browser. I had picked up a copy of Eric Drexler's book Nanosystems, read the section on MM2, and understood it well enough to implement it in code. That was a lot of fun to work on, and as I've watched web technology progress, I've occasionally thought about taking another stab at it. I am now in the process of doing that, as some parts have gotten easier, others have gotten faster, and some have gotten just plain interesting. JavaScript, despite a few warts, is charming, and widespread deployment of HTML5 is a big help too.

Molecular modeling approaches like MM2 generally work by computing potential energy as a function of the relative positions of atoms within a molecule. This is done by decomposing the potential energy into a series of terms, relating to the lengths of chemical bonds, the angles between bonds that share a single atom, or the dihedral angle between two bonds linked by a third bond. These terms are parameterized based on the elements and hybridizations of the atoms involved. Force contributions of these terms can be computed independently and summed together to find the forces acting on each of the atoms. This is necessarily a simplification, and more accurate approaches exist involving solutions to the Schrodinger wave equation describing probabilistic locations of electrons and nuclei. But the simple mechanical approach is adequate for getting a sense of the molecule's general shape and how it perturbs over familiar temperature ranges.

For a potential energy function E(p) relating to some geometric parameter p, we can use the chain rule to get the force on an atom at position (x,y,z):
(fx, fy, fz) = -E'(p) (∂p/∂x, ∂p/∂y, ∂p/∂z)
where ∂p/∂x is the notation for the partial derivative of p in the x direction. The forms of the potential energy functions are pretty straightforward (1, 2), and taking their derivatives is not difficult. The remaining trick is to determine the partial derivatives of parameters with respect to x, y, and z. Let's consider a simple case where the parameter is the distance between two atoms at positions (ux, uy, uz) and (vx, vy, vz). Then the distance is r where
r2 = (ux - vx)2 + (uy - vy)2 + (uz - vz)2
and taking partial differentials with respect to ux yields
2r ∂r = 2 (ux - vx) ∂ux

∂r/∂ux = (ux - vx) / r
The same operations with uy and uz tell us that the force vector acting on the first atom is f = -E'(r) r / r where the direction of r is from the v atom to the u atom, and the inverse for the second atom. Here, boldface denotes a vector quantity.

It's usually easy to apply some geometric intuition and determine a unit vector in the direction of greatest change for a parameter p (the gradient of p with respect to a particular atom's position). For an angle θ involving three atoms at positions u, v and w, with v being the vertex, the gradient for u lies in the plane and is perpendicular to (u - v), pointing away from the third atom. Then the vector (∂θ/∂x, ∂θ/∂y, ∂θ/∂z) is in the same direction as the gradient unit vector, and it's necessary only to determine a scaling factor. That can be obtained by doing a little trigonometry to determine that a teeny move of distance δ in that direction will produce a parameter change , and then the magnitude of the force is -E'(θ) dθ/δ.