Computational methods to determine the
structure and behavior of molecules

This discussion is restricted to efforts to determine the positions of the atoms (even more specifically the positions of the nuclei of these atoms) of a molecule. These atomic positions define what is commonly called the structure of a molecule, but excludes the shape and nature of the distribution of electrons around the nuclei. It is this electron distribution which largely determines the molecule's chemical reactivity, about which we will thus have little more to say.

There are two main methods that are commonly used, but they are both based on electrostatic-mechanical models for the molecule. Atoms bonded to each other are considered to be point masses connected by springs. Interactions between non-bonded atoms are represented by electrostatic forces due to non-integer charges on the atoms. There are many specific models that have been proposed by different groups for different classes of chemical compounds. The models are called Force Fields, even though they actually define energy and not force. The force fields differ from each other mainly in the complexity of the spring network and on the compressional and torsional spring constants. We now present two distinct, but closely related, methods in which these force fields are used.

Mechanics: The goal here is to determine a static structure, more specifically the structure that has the lowest energy as defined by the force field. It is implied that this is the structure of the molecule. Since the force field explicitly defines energy as a function of atomic positions, conventional algorithms for minimization of a function produce the local minimum.

Dynamics: The goal here is to determine the motion of the atoms of the molecule relative to each other, (the center of mass of the molecule is fixed). The terms of the force field are now interpreted as forces, i.e. the vector derivatives of the energy with respect to position are computed. These forces and the dynamics equation, F=ma (or here a=F/m) give the acceleration, and then the velocity, and finally the position of each atom. Often the temperature of the molecule is kept constant. This can be accomplished by periodically normalizing the average of the magnitudes of the velocities to a fixed value, which simulates energy flowing to and from a solvent bath of a fixed temperature in which the molecule is immersed.

The Mechanics approach always produces an answer, a structure. However, this structure is only a local minimum in multidimensional configuration space. This is presumedly the structure the molecule would end up in if it started at the starting structure and were immersed in a bath at absolute zero. If one is actually interested in the structure at a higher temperature, say room temperature, who knows what structures the minimization algorithms bypassed as they searched for the minimum energy. There are ways to search for a global minimum, but this becomes more and more difficult as the size of the molecule increases, since the dimension of the structure space in which you are trying to find the energy minimum is proportional to the number of atoms in the molecule.

The Dynamics approach produces a realistic (hopefully) time trajectory, i.e. history of structure versus time. However, it can't produce a "final" definite structure, if only for the fact that the atoms have to vibrate and rotate about equilibrium positions to represent a molecule at the simulated temperature, which is not usually absolute zero. However, can be more disturbing is a sudden change during the simulation into another structure, and then perhaps a change back to the original structure. Which one is the correct structure? The answer may be both, and thus the "correct structure" may need to be represented as a sequence of structures and the probability of the molecule being in each of them.

A more fundamental problem in molecular dynamics simulations is the fact that only a very short time can be simulated, in the range of a few nsec at most. Thus, one could start with an improbable structure that was stable for a usec, and not see the transition to the structure that you would find the molecule in the lab. This is the analog to not being able to search the structure space thoroughly enough in the mechanics approach to find the stable structure.