Molecular Mechanics (Part 1: Bonded Interactions)

Molecular dynamics refers quite generally to the practice of evaluating equations of motion given a potential energy function—the negative gradient of which is the force—on the particles that make up a model of a molecular system. While an accurate model can be derived by using a quantum mechanical potential to describe the system, such an approach is prohibitively expensive for macromolecular studies. In this post, I will describe molecular mechanical force fields as a model for evaluating potential energies and forces on each atom in a molecular system. While force fields vary in their exact functional forms and parametrization schemes, they are broadly composed of two classes of interactions between particles—those between atoms separated by a small number (typically 1 to 4) chemical bonds and those between all other pairs of atoms in the system.

Chemical bonds are modeled as a simple analytical function of the distance between bonded atoms which are often fitted to high-level quantum mechanical potential energy scans and/or vibrational spectroscopic measurements. The most common functional forms used to model chemical bonds are even-order truncated Taylor series expansions—like simple quadratic or quartic functions—and the Morse potential. These functions are shown alongside an exact QM energy scan for the hydrogen molecule below.

Potential energy scan for Hydrogen molecule alongside 3 simple, analytical functions fitting the potential.

Potential energy scan for Hydrogen molecule alongside 3 simple, analytical functions fitting the potential.

It is clear from this figure that for large vibrational energies—i.e., where bond vibrations exhibit strong anharmonicity—the more complex quartic and Morse potentials are better models. For most bonds, however, a simple quadratic description suffices at low temperatures around 300K. Such descriptions are not appropriate for all simulations, however, and we may wish to develop models utilizing atypical functional forms that are better suited for our applications. [1]

In the past, such a project necessitated the development of a new program or expansion of an existing program to run the simulations followed by extensive validation to make sure the code works. The code would then require extensive optimization and potentially even radical redesign in order to be used for production-level calculations to generate a statistically sufficient amount of data. With the incredibly flexible support for custom forces unique to OpenMM, however, a researcher need only assign a functional form for the valence bond terms and assign parameters for each of the terms she wants to model. OpenMM will generate an efficient kernel to compute both energies and analytical forces given an equation for the bond term that can run directly on GPUs—the fastest computational hardware widely available for molecular simulations today. This remarkable technology allows researchers to create pre-tested, production-ready software for brand new molecular models in minutes, drastically reducing the months-long process of writing, optimizing, and validating the code—even if he is not a skilled programmer!

But I digress... back to the force fields. Force fields also typically model atoms separated by two bonds using a function of the angle between those two bonds (again, typically a simple quadratic function whose frequencies can be computed or determined by vibrational analysis). Atoms separated by three bonds—forming torsions—often require a functional form like Fourier series that permit modeling the periodic, highly flexible rotations around the central bond. Force fields also frequently have potentials designed to keep aromatic functional groups in a planar configuration using a so-called improper torsion since it is a torsion angle defined between 4 atoms that are not connected the same way (see the figure below for an illustration). The final potential energy function I will mention here that is used in many force fields is coupled-torsion correction map (CMAP). This is a potential evaluated on a discrete, interpolating grid that describes how the flexibilities of two torsions may be coupled together in complex ways.

Schematics of various bonded potential terms. A) Bond term. B) Angle and Urey-Bradley distance terms. C) Torsion (or dihedral) term. D) Newman projection of the same torsion as C. E) Improper torsion used for maintaining planarity. F) Coupled-torsion correction map (CMAP) term.

Schematics of various bonded potential terms. A) Bond term. B) Angle and Urey-Bradley distance terms. C) Torsion (or dihedral) term. D) Newman projection of the same torsion as C. E) Improper torsion used for maintaining planarity. F) Coupled-torsion correction map (CMAP) term.

This ends the first part of my discussion about typical biomolecular force fields, and hopefully shows how OpenMM can be an invaluable tool in your tool set for carrying out molecular simulations. My next post will focus on arguably the most important [2] part of a force field—the nonbonded part of the potential.

[1] Zhang, L. and Hioaki, T., Wear, 1997, 211 (1) p. 44 (http://dx.doi.org/10.1016/S0043-1648(97)00073-2)
[2] Okay, okay. All parts are very important and must be applied collectively for any force field to work.