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):
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 dθ, and then the magnitude of the force is -E'(θ) dθ/δ.
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)2and taking partial differentials with respect to ux yields
2r ∂r = 2 (ux - vx) ∂uxThe 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.
∂r/∂ux = (ux - vx) / r
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 dθ, and then the magnitude of the force is -E'(θ) dθ/δ.
No comments:
Post a Comment