Minimization - Vibrational Calculations

Contents

Application of Minimization to Vibrational Theory

The vibrational frequencies and modes of a molecule are strictly dynamic properties. However, it is possible to calculate the harmonic vibrational frequencies of a molecule from just information at its equilibrium geometry by expanding the potential energy surface as a Taylor series, truncating after the second term, and considering infinitesimal displacements. This harmonic approximation usually gives a good description of the true frequencies and normal modes and can be quite valuable for tasks ranging from evaluating the quality of a forcefield to understanding vibrational shifts induced by conformational changes or other interactions. The harmonic vibrational frequencies also can be used for zero-point vibrational corrections and for deriving vibrational free energy contributions. These effects can be important in comparing conformational energies and rotational barrier heights.

Following Wilson et al. (1955), the kinetic energy of the nuclei is:

Eq. 4-11:
where the coordinates represent displacements from an equilibrium structure. If the 3N Cartesian coordinates are replaced with 3N mass-weighted coordinates as follows:

Eq. 4-12:
where m is the mass of the atom associated with the coordinate, and x runs over the y and z coordinates, as well as x, then the kinetic energy has the following simple form:

Eq. 4-13:
When the potential energy of the system is expanded as a Taylor series in the same coordinates, it yields:

Eq. 4-14:
V0 is simply a constant--the energy scale can be chosen so that V0 = 0. The definition of an equilibrium structure is that the force on each atom is zero. The second term in Eq. 4-14 also is zero, leaving the following second-order approximation of the potential energy:

Eq. 4-15:
Using this approximation of the potential energy in Newton's equations of motion yields the following simultaneous second-order differential equations:

Eq. 4-16:
The solution to these equations can be of the form:

Eq. 4-17:
where the Ai are related to the relative amplitudes of the vibrational motion, 1/2 is proportional to the vibrational frequency, and is a phase. Substituting Eq. 4-17 in Eq. 4-16 yields a set of algebraic equations:

Eq. 4-18:
where ij is a Kronecker delta, which equals one if i = j and zero otherwise. This is an eigenvalue problem that is readily solved numerically by standard techniques. The second derivatives of the potential energy, often called the force constants, can be analytically evaluated for most energy expressions used in molecular mechanics, in terms of the Cartesian coordinates of the atoms. A simple transformation to the mass-weighted coordinates then gives the values needed in Eq. 4-18.

Vibrational Frequencies

The Discover program calculates the vibrational frequencies by calculating the second derivative matrix, mass weighting it, and then diagonalizing it to obtain the eigenvalues. These eigenvalues are then converted to vibrational frequencies in wavenumbers as follows:

Eq. 4-19:
where the conversion factor F converts the units from kcal mol-1 to wavenumbers. Of the 3N coordinates used to calculate the energy and vibrational frequencies, six correspond to net translations and rotations of the molecule (five for linear systems). These six modes have no restoring force and therefore have vibrational frequencies of zero for a minimized structure. Due to numerical inaccuracies, the Discover program reports these frequencies with small values, which are typically less than 0.1 cm-1. If the structure is not perfectly minimized, the first-order terms in the Taylor expansion of the potential surface in Eq. 4-14 do not vanish. In turn, they introduce terms that couple the net rotations of the molecule with the internal motions and both perturb the internal vibrational frequencies and give apparent frequencies for the three rotations. Thus, the magnitude of the six "zero" frequencies is a good indication of the quality of the calculation.

Transition States

For calculating vibrational frequencies, the structure must be minimized and the gradients must be zero. This does not mean that the configuration of the molecule must be at a minimum, but rather that it must be at a stable point on the surface. If the structure is not at a minimum, but rather is at a saddle point or transition state in one or more directions, this is reflected in the eigenvalues and the reported vibrational frequencies. For a saddle point, at least one eigenvalue is negative, which means that the curvature of the surface along at least one normal mode is negative. By convention, the imaginary frequencies for such modes are reported as negative. Therefore, the reported vibrational frequencies describe the character of the stable point on the surface. If all frequencies are real, it is a minimum; if one frequency is imaginary, the structure is in a simple transition state; if two or more frequencies are imaginary, a double or more complicated transition state is indicated. The normal modes corresponding to the frequencies can be analyzed to understand the reaction path going through such a transition state.

Thermodynamics

The quantum mechanical solution for the vibrational energy of a set of uncoupled harmonic oscillators, which corresponds to the classical treatment outlined above, is:

Eq. 4-20:
where the summation is over all the vibrational frequencies, ni is the vibrational quantum number for each vibration, h is Planck's constant, and i is the vibrational frequency. This leads to the following correction to the classical forcefield energy:

Eq. 4-21:
where k is the Boltzmann constant and T is the temperature. The first term, hi / 2, is the zero-point correction; the second term corrects for the average thermal population of vibrational levels at the temperature T. This leads to a vibrational free energy correction of:

Eq. 4-22:
and a vibrational entropy of:

Eq. 4-23:

General Methodology for Vibrational Calculations

Computation of the harmonic vibrational frequencies requires the storage and diagonalization of the second-derivative matrix, which has dimensions of 3N x 3N, where N is the number of atoms. The work involved in the diagonalization scales as N 3 and quickly becomes prohibitively expensive for more than a few hundred atoms. The frequencies are valid only for well minimized structures having maximum derivatives no greater than approximately 0.001 kcal mol-1 Å-1.

The forcefield is a major consideration. Most forcefield development has emphasized structures and energetics rather than vibrational frequencies. As a result, the frequencies calculated with forcefields such as AMBER, MM2, CHARMm and, to a lesser extent, CVFF, may often be in error by several hundred wavenumbers. The inclusion of cross terms such as bond-bond and bond-angle terms is crucial for the accurate reproduction of experimental frequencies. The CVFF forcefield includes such cross terms and was, in part, parameterized to reproduce experimental frequencies, which explains its moderately good performance for vibrational calculations. The recent Class II forcefields MM3 and CFF91 were explicitly designed to evenly weight vibrational frequencies as well as structural and energetic properties. Therefore, they provide the most reasonable and consistent results, usually within 50-100 cm-1. This error of up to approximately 100 cm-1 appears to be the current limit of general-purpose, transferable forcefields.

Beyond these considerations, the vibrational frequencies can be used for two classes of problems. The first is for determining the shape of the potential energy surface; that is, the characterization of stable points as minima, transition states, or other points. For this purpose, the question of forcefield accuracy is less important. The qualitative, rather than quantitative, shape of the surface is all that is important. The second use is for comparison with experimental results. The vibrational frequencies and thermodynamic corrections depend strongly on the forcefield as well as on the fundamental harmonic approximation invoked. By nature, low-frequency modes are less harmonic. The torsion and nonbond interactions, which dominate low-frequency modes, are fundamentally anharmonic; hence the interpretation of the calculated low-frequency modes should take this into account. Unfortunately, these low-frequency modes make the largest contributions to the vibrational entropy.

Example Vibrational Calculation

Hagler et al. (1979c) investigated the stability of the -helix and C7 conformations of short polypeptides, explicitly calculating the vibrational corrections to the enthalpy and free energy. Their results show that, for the specific cases they examined, the conformational difference in the vibrational contribution to the enthalpy was essentially negligible. This result is not too surprising, because the enthalpy contribution is dominated by the high-frequency modes such as the C-H stretching frequencies. These high-energy modes would not be expected to change much between conformations differing primarily in torsional angles. In contrast, the entropic contributions to the free energy differences were substantial and approximately the same magnitude as the total enthalpy differences. Essentially, the entropic contribution is significant for the same reason that the vibrational contribution to the enthalpy is small: the low-frequency modes, particularly the torsions, contribute heavily to the free energy. For instance, for Boc-Met6-OMe the forcefield energy of the minimized C7 conformation was 5.09 kcal mol-1 lower in energy than for the -helical form. Correcting the difference in enthalpies for the vibrations decreased the difference by only 0.05 kcal mol-1. The contribution to the change in entropy, however, led to a further decrease in the free energy of 3.99 kcal mol-1. Although the energy difference of the minimized structures appeared to be 5.09 kcal mol-1, the estimated difference of the free energies was only 1.05 kcal mol-1. Several other polypeptides of varying lengths and compositions were studied. In all cases, the results were similar to those for the example given here, which clearly indicates that the vibrational contributions to relative free energies can be important.

Main access page Theory/Methodology access.

Minimization access Example Minimization Calculation Molecular Dynamics

Copyright Biosym/MSI