About Molecular Simulation

Molecular simulation is a general term for the use of computer models to describe chemical systems at an atomic level of detail. In a molecular simulation, a chemical system is described using the individual positions and orientations of every atom or molecule. Different methods can be used to obtain both thermodynamic and kinetic properties of the system.

A molecular simulation specifies a collection of atoms and molecules in a simulation cell, interacting through a potential, and evolving according to some simulation algorithm.

Potentials

The central concept in molecular simulation is that of the interaction potential, also known as the force field, which governs the interactions between different molecules, and between different atoms on the same molecule.

Most molecular simulations use "empirical" potentials for describing the interactions between molecules, that consist of convenient functional forms with parameters fit to either experimental data or the results of high-level quantum mechanical calculations. These techniques model the system with purely classical mechanics so that isotope effects through zero-point energies are neglected. More seriously, it is difficult to model chemically reactive systems with classical models, because the breaking and forming of bonds involved electrons and is a fundamentally quantum-mechanical process.

The earliest potentials used simple hard spheres to model rare-gas atoms. Modern potentials break the interactions between atoms and molecules into different types, and model each type with different functions:

• Atom-atom repulsions are represented by either an exponential repulsion term, or an inverse power-law (usually r^-12) function.
• Induced-dipole - induced-dipole attractions, also known as London forces, dispersion forces, or van der Waals interactions, are modeled with an r^-6 function.
• Charge-charge interactions are modeled with Coulomb's Law potentials, and Charge-dipole and dipole-dipole interactions either using "point dipoles" and the corresponding classical elecrostatics, or by representing the dipoles as a collection of "partial charges" and using charge-charge interactions.
• Bonds are either considered as completely rigid, or modeled as stiff harmonic oscillators.
• Bond angles and bond torsions are modeled with harmonic and sinusoidal potentials
• Charge - induced-dipole interactions can be modeled either with fluctuating charge models, in which the total charge distribution between several sites on a molecule is allowed to adjust dynamically, or with polarizability where the dipoles and charges are determined self-consistently by the electric field they experience.

The most formally realistic potentials used are based on solving, using electronic structure methods (usually Density Functional Theory) the Schroedinger Equation for all of the electrons in the system, thus obtaining the electronic Potential Energy Surface (PES) upon which the nuclei move. In most simulations of this type, the only significant approximations are that the PES is not "exact" due to approximations in the electronic structure calculations, and that the nuclei are still treated as classical particles. (This last approximation will likely be lifted in the near future with more advanced methods.) Because DFT calculations are not especially good at describing phenomena due to electron correlation, these potentials do not do a great job on the dispersion forces mentioned above.

Simulation cells

The ususal workaround for this problem is to impose what are known as periodic boundary conditions (PBCs) on the system, by requiring that the (macroscopic) system of interest be composed of a very large number of identical "cells", replicated in space in the same fashion as the unit cell of a crystal. Therefore, we only need to keep track of the positions, velocities, etc., for the particles in one of these cells, since they are all identical. In practice, cells containing between a few hundred and a few thousand molecules are sufficiently realistic, except in special, and fairly obvious, cases.

Implementation of PBCs is easy; when a molecule moves out of one side of the cell, it re-enters on the other side (e.g., the "shared" side of the next cell, which is identical to the cell it just left.) That is, particle motion works exactly the same way as the old video games "asteroids" and "pac-man", among others!

Most simulations are done using cubic simulation cells, primarily for ease of implementation. One can also use truncated octahedral or rhombic dodecahedral cells, as well as somewhat more exotic solutions to this problem.

Simulation Algorithms

There are three broad categories of algorithms used in this type of study, which I will relatively arbitrarily categorize as molecular mechanics, molecular dynamics, and Monte Carlo.

Molecular mechanics

Molecular dynamics

It can be shown that, in the limit of short time-steps, this procedure samples states accessible in the microcanonical ensemble. Since most experimental work is done under conditions of constant temperature, and either constant volume or constant pressure, the simple simulation algorithm just described has been augmented with additional features that allow one to specify the "configurational" temperature, or allow the simulation to access a range of energies and/or pressures that correspond to either the canonical or isothermal-isobaric ensembles.

The main strengths of molecular dynamics are that it efficiently samples the given ensemble, and that it provides dynamical quantities, such as velocity autocorrelation functions, dynamic scattering factors, and diffusion constants. The main weakness of molecular dynamics is an inability to access very long time scales, on the order of one microsecond or greater.

Monte Carlo

The Monte Carlo simulation method is also easily adapted to "force" the system to explore unfavorable regions of phase-space, which is useful in calculating free energies and free energy differences between phases, among other things. (This kind of adaptation can also be done in MD simulations, but is considerably trickier to implement.)

The drawback, of course, is that one does not really obtain dynamical properties from this type of simulation, except in an approximate and somewhat roundabout way.