III. An Empirical Energy Function: Free Energy vs. Potential Energy

Finally, given a defined system and its initial atomic coordinates, we need a function describing the energy of the system for any configuration, $\vec{R}$, of the atomic coordinates. A functional form must be chosen for the energy as well as the associated numerical constants. For macromolecular simulation potentials, these parameters number in the thousands and include spring stiffnesses and equilibrium distances, torsional barriers and periodicities, partial charges, and Lennard-Jones coefficients. The energy function and its associated constants are contained in the Parameter File (.prm). Development of parameter sets is a laborious process. Both the functional form and numerical parameters require extensive optimization. A brute-force iteration of simulation and parameter modification is performed to improve agreement between simulations of model systems and information derived from ab initio calculations, small-molecule spectroscopy, and educated guessing.

Free Energy vs. Potential Energy

For a system held at constant $N\cal{V}\it {T}$, the Helmholtz free energy, $A \equiv U - TS = -kT \rm {ln} (Z)$, is a minimum at equilibrium, where $U = \langle E \rangle$ is the average total energy of the system (kinetic energy plus potential energy), T is the absolute temperature, S is the entropy, and Z is the partition function (eq 3). Suppose the pressure P is held constant instead of the volume (constant $NPT$, the `isothermal isobaric' ensemble). In this case, the Gibbs free energy, $G \equiv U + P\cal{V} \it {- TS}$, is minimized at equilibrium. Note that the enthalpy, $H = U + P\cal{V}$, is the quantity at constant pressure that corresponds to U at constant volume. Differences in G drive chemical reactions.

Some empirical energy functions are designed to approximate the Gibbs free energy G. For example, in Monte Carlo studies of protein structure prediction, the energy function may be based simply on the likelihood of residues of type i and j being within a certain distance of each other. The probabilities p are determined by counting the number of times that residues i and j are found close to each other in the protein structures deposited in the Protein Data Bank. They are then converted into $\Delta G$-like energies by: $p_{ij} \propto e^{-\Delta G_{ij}/k_BT}$. Because the p's are derived from structures at constant T and P determined experimentally, these energy functions account for entropic contributions to the Gibbs free energy in an approximate way.

In most molecular dynamics software packages, however, the empirical energy function, $V(\vec{R})$ (not to be confused with the volume $\cal{V}$), is developed to approximate the potential energy of the system. In general, it does not include entropic effects in any effective way. Many simulations have been performed at constant energy, E. That is, E is fixed and T fluctuates about an average value as energy is exchanged between the kinetic energy and the potential energy. In principle, simulations performed at constant T and P mimic experimental conditions better than simulations at constant E. Recently, an improved constant-$PT$ algorithm has been developed [11]. Constant E simulations have the advantage that they allow energy conservation to be checked. Any significant drifts in E indicate a problem that should be tracked down before continuing the simulation. Although it fluctuates, the temperature is still well defined at constant E, and differences between dynamics at constant T and constant E are generally not too significant on the time scales currently accessible to MD simulation (100's of ps to a few ns). However, the constant-$PT$ simulation may well become the standard as large solvated systems are simulated over longer time scales.