Next: H-bonds and N-bonds, Previous: Geometric Constants, Up: PARAM1
Having established these geometric parameters, the next step was
the definition of force constants. Force constants describing simple
hydrocarbons and amides were obtained by least squares fitting to the
vibrational spectra of propane, butane, pentane, and N-methyl acetamide
(Snyder et al and Warshel, Levitt, and Lifson). This fitting was carried
out in stages, first choosing parameters and evaluating the resulting
normal modes until a set was obtained with correctly ordered vibration
frequencies. These parameters were then refined by least squares fitting
to vibrational spectra of all the isotopic variants available. This
fitting was performed with a user subroutine linked to CHARMM that
called the CHARMM energy and vibrational analysis routines for normal
mode analysis, and used an IMSL minimization routine (ZXSSQ).
This fitting algorithm uses finite differences to construct an
approximate first derivative vector and second derivative matrix,
requiring a rapidly increasing number of function evaluations as the
number of free variables rises. To reduce the complexity of this
fitting, parameters were grouped (bonds vs. angles vs. dihedrals and
impropers) and minimized separately until good agreement was obtained
over the entire range of observed frequencies. The experimental and
calculated normal modes for N-methyl acetamide are shown in the first
table. The most prominent differences between the current and the
previous parameterization are seen in the low frequency and out of plane
bending modes.
The dihedral angle potential for alkanes was found to be 1.6 kcal/mole by this procedure in good agreement with the value reported in the dipeptide and acetylcholine work of Gelin (1.8 kcal/mole). It should be noted that the previous parameter sets used a much lower value (0.2 kcal/mole). The old value resulted in essentially unhindered rotation between trans and gauche isomers and a very small gauche well depth (0.03 kcal/mole). The new parameters are in much better agreement with experimental observations as can be seen in the second table.
For the aromatic systems it was not possible to perform this least squares fitting because the simple potential energy form that we are using is not able to correctly order the various normal modes (Varsanyi and Dollish et al). This is apparently a result of the strong valence interactions between internal coordinates, and presumably could be overcome using a potential form that includes valence cross terms. The large number of modes (30) in benzene and the strong interaction of carbon and hydrogen motions also complicated the analysis of this system. For these reasons, benzene was analyzed by manually adjusting parameters for a model including all hydrogens. The heavy atom parameters obtained from this system agree with those used previously, and were kept fixed in modelling more elaborate aromatic systems.
Using the CHARMM potential energy function, aromatic out of plane bending behavior is determined by both the arrangement and parameterization of the improper dihedrals. For histidine, phenylalanine, and tyrosine the current parameterization and arrangement gave good agreement for the intrinsic ring modes (based on the benzene modelling), and they were retained. The ring substituent out of plane modes were not well described, and it was necessary to parameterize the system using the lowest out of plane bending modes for toluene (217 cm-1) and phenol (241 cm-1) to fix new constants for the ring CB and ring hydroxyl bends respectively. The present improper dihedral was retained, but in both cases the force constants required large increases (from 25 kcal/mol/rad2 to 120 and 150 respectively).
Tryptophan required both rearrangement and reparameterization of the improper dihedrals to obtain reasonable agreement with the skeletal out of plane bend (400 to 500 cm-1 in purines and napthalene) without creating very high frequency modes involving the ring intersection. The new arrangement replaces the improper dihedrals centered on the atoms at the intersection of the rings (CD2 and CE2) with terms that are trans two fold dihedrals beginning on one ring, running across the bridging atoms and ending on the opposite ring as shown below:
CZ2 NE1
* * * *
* * * *
CH2 CE2 *
| | *
| | CD1
| | *
| | *
CZ3 CD2 *
* * * *
* * * *
CE3 CG
|
|
CB
PARAM1 PARMFIXn
CZ2 CZ2----------NE1
* * *
* * *
CE2 * CE2
| *
| *
CD2 CD2
*
*
CG
NE1
*
*
CE2 CE2
| *
| *
CD2 * CD2
* * *
* * *
CE3 CE3-----------CG
In addition, two improper dihedrals were added as shown below.
CH2--------CE2 CE2 *
| | *
| * CD1 | * CD1
| * |
CD2 * CZ3-------CD2
This new arrangement distributes the out of plane bending forces evenly across the ring in a manner analogous to the delocalized electronic structure of the system, and does so without excessively elevating the frequency of CD2/CE2 out of plane modes. With the PARMFIX10 parameters these modes had frequencies of about 1000 cm-1, but the skeletal out of plane bend occurred at 120 cm-1. Simply raising the force constants without rearranging the dihedrals resulted in a bridge atoms out of plane mode at 3000 cm-1! Aside from being unrealistic, this high frequency would force the use of small step sizes to retain accuracy in the dynamics integration.