Molecular Dynamics Introduction

Molecular dynamics (MD) is a computational method that simulates position as a function of time producing a "movie" of the system's behavior. Using computers to apply Newton's equations of motion many, many times, CHARMM MD software is able to model the movements and fluctuations of individual atoms. Consider a one-body system from classical mechanics

F = ma
^-δU/_δx = m(^δ2x/_δt²)

where F is force, m is mass, a is acceleration, U is potential energy, and t is time. These equations are deterministic: given the velocity and position for a time, the position and velocity can be known at any previous or later time.

For every atom, CHARMM applies this principle to calculate acceleration and then moves the atoms for a timestep. In order to model atom movement accurately, the difference in potential energy before and after the timestep needs to be small. A large shift in potential energy means that the motion of atoms during the applied timestep is unrealistic. The chemical properties of atoms determine the potential energy curves and therefore the appropriate step size. Hydrogen bonds vibrate rapidly and require a timestep on the order of a few femptoseconds (10^-15 s).

To derive the Verlet integration algorithm, add the following Taylor series expansions together:

x(t+δt) = x(t) + v(t)δt + ¹/₂ a(t)δt² + O(δ³)
x(t-δt) = x(t) - v(t)δt + ¹/₂ a(t)δt² + O(δ³)

x(t+δt) = 2x(t) - x(t-δt) + a(t)δt² + O(δ³)

where x(t) is position at time t, v(t) velocity, and a(t) acceleration. Note that by only expanding to 3 terms we have exchanged an exact solution for an approximation which makes an error on the order of δ³.

Once a timestep is chosen, hardware limitations and parallel efficiency determine the overall simulation timescale. A simulation of a lipid bilayer may contain from 30k to 80k atoms. On a commodity cluster (like our Rocks clusters), this simulation can reach 100 ns in a reasonable amount of time. Specialized supercomputers (such as Anton) are capable of producing 10 μs in a fraction of the processor time.

MD has several advantages. It gives us a computational microscope for investigation. The atomic scale of the simulation means that measurements have a very high resolution. You can directly measure the structure and properties of the system without perturbing it. You also have complete control over the system and simulation. There is also no contaminatin or sample degradation.

MD also has several concerns. How accurately can you represent the potential energy of real life molecules? Another concern is sampling. The MD approach makes the statistical argument that average properties of the microscopic states agree with the macroscopic properties measured in a conventional lab. Is the system well-equilibrated and is the simulation long enough to obtain a good representation of its properties? Another concern is timescale. Some physical phenomena occur on a timescale that traditional MD cannot reach. Thus MD is unable to provide a microscopic movie of these phenomena. The desired behavior of the system (Ex: protein-folding, energy equilibration, diffusion, mixing, ...) needs have an attainable timescale.

An important part of MD is an accurate potential energy model. Chemistry describes these potential energy functions. Energy potential can be separated into intra-molecular forces (the arrangement of bonds) and inter-molecular forces (effects from nearest neighbor, non-bonded). Bonded forces are modeled with terms describing bond stretching, angle bending, and rotation about a bond known as a dihedral or torsion angle.

U_bonded = U_bond + U_angle + U_dihe
U_bond = k_l(l - l₀)
U_angle = k_θ(θ - θ₀)
U_dihe = Σ [k_φ(cos(nφ + δ)+1))]

Bond and angle potential are modeled with a harmonic function that use a "spring" constant. They use force constants k_l and k_θ, and equilibrium values l₀ and θ₀. The dihedral term describes a periodic potential and is summed for the 1-4 atoms that define the dihedral. It uses force constant k_φ, multiplicity n, and offset δ. Non-bonded forces are modeled with terms describing electrostatic potential (strong force) and Van der Vaals potential (weak force).

U_nbonded = U_elec + U_vdw
U_elec = Σ [^q_iq_j/_{εD rij}]
U_vdw = Σ 4ε_ij[(^σ_ij/_rij)¹² - (^σ_ij/_rij)⁶]

The electric potential uses charges q_i and q_j for atoms (i,j), ε_D is the dielectric constant, and r_ij is distance. Van der Vaals interactions are described with a Lennard-Jones 6-12 potential: ε_ij is the potential energy minimum between atoms (i,j), σ_ij is the position of that minimum, and r_ij is the distance.

The summations for non-bonded forces include terms for all atom pairs (i,j). That is, every atom feels these two forces from every other atom. By their nature, bonded forces will only have a few terms to calculate but a system with n atoms will need to compute n² non-bonded terms. That is a very expensive calculation. In a traditional evaluation of potential energy terms, CHARMM will spend about 80% of processor time calculating non-bonded forces.

As you can see, these potantial energy terms require the use of parameters. Every atom and molecule needs to be parameterized to model their real-life counterpart. Quantum mechanical calculations provide a very accurate tool to parameterize the CHARMM force field. Additionally, simulations can be compared to real life samples to determine the validity of a force field scheme.

Here is an excellent CHARMM tutorial at EMbnet. It will go over MD and force fields as well as basic CHARMM principles. Throughout the tutorial you may find that you wish to look up the meaning of particular CHARMM commands and their usage. The CHARMM webiste provides documentation by version. The dynamc.doc (dynamics), corman.doc (coordinate manipulation), select.doc (atom selection), minimiz.doc (energy minimization), miscom.doc (miscellaneous commands), and nbonds.doc (treatment of non-bonded atom interaction) detail commonly used commands. Try looking up the commands from the tutorial in the CHARMM documentation.

See the quick guide on the Wabash Guide for the location of charmm executables. For the embnet tutorial you may need to use your own topology and parameter files, collectively known as toppar. These files will be in the toppar folder in the CHARMM installation. Try looking in /home/charmm/c34b2/toppar/.

The Charmm-GUI website is also helpful for beginners. You can build a lipid bilayer simulation from scratch. Navigate to "Input Generator" and then "Membrane Builder". The website will walk you through the creation of a lipid bilayer and generate input files. Some steps may take a while to complete and the whole process may take an hour or two. To keep things simple try a membrane only system. Note the synopsis of each of the five steps at the bottom of the page as you complete each step.

Back to home page.