4020
J . Phys. Chem. 1985,89. 4020-4026
mechanism for two-photon ionization appears, however, to be different for tetracene in a boric acid glass at room temperature from that for tetracene in solid argon at 20 K. High-pressure mercury arc photolysis of 3 in boric acid produced a good yield of the 3 cation,” but a like irradiation of 3 in solid argon gave no products. The straightforward conclusion is that the ionic medium facilitates intersystem crossing to the TI state, which is sufficiently long lived to serve as the intermediate state for a biphotonic ionization process, but in solid argon intersystem crossing is less efficient. This is in accord with substantial changes in the relative quantum yields of fluorescence and phosphorescence for naphthalene between hydrocarbon and noble gas matrices.25 Robinson and Frosch also point out that nonradiative transition probabilities are affected by solvent-solute electrostatic interactions. Vibrational coupling of the SI guest state and the host matrix will clearly be greater for molecular matrices. The markedly different behavior for 1 and 2 in solid argon upon mercury arc irradiation was demonstrated in an anthracene experiment with about 10% phenanthrene impurity (based on the sharp SIorigin absorbance at 27 OOO and 29200 cm-’, respectively, for 1 and 2) and CH2C12electron trap. Irradiation produced strong 2 cation and no 1 cation absorptions.
Conclusions Argon resonance photoionization of five condensed ring aro(25) Robinson, G. W.; Frosch, R. P. J . Chem. Phys. 1963, 38, 1187.
matic precursors during condensation with excess argon at 20 K produced and trapped the radical cations as characterized by electronic absorption spectra. The argon matrix spectra are in agreement with glassy matrix spectra and band positions are in agreement with predictions from photoelectron spectra and HMO calculations. Vibrational fine structure in the absorption bands correlates with symmetric fundamentals from Raman spectra of the precursors. Irradiation of phenanthrene, 1,2-benzanthracene, and chrysene samples containing CHzClzwith the high-pressure mercury arc affected resonance two-photon ionization, but similar experiments with the more symmetrical isomers anthracene and tetracene failed to produce cations. Stronger precursor SI band intensities for the latter precursors suggest shorter SI lifetimes, which prevent the two-photon ionization process from involving SI as an intermediate state. The argon matrix provides a more inert environment for condensed ring aromatic cations than glassy matrices, which results in smaller red matrix shifts and better resolved vibrational structure in the more inert matrix. Acknowledgment. We gratefully acknowledge financial support from the National Science Foundation, a preprint from T. Bally on matrix spectra of anthracene and phenanthrene cations, and the assistance of R. T. Arlinghaus with phenanthrene experiments. Registry No. 1, 120-12-7; I+, 34512-28-2; 2, 85-01-8; 2+, 34504-68-2; 3, 92-24-0; 3+, 34512-31-7; 4, 56-55-3; 4+, 42299-25-2; 5, 218-01-9; 5+, 34474-38-9; Ar, 7440-37-1.
On Finding Stationary States on Large-Molecule Potential Energy Surfaces Dzung T. Nguyen and David A. Case* Department of Chemistry, University of California, Davis, California 9561 6 (Received: April IO, 1985)
We present a procedure for locating stationary points on molecular potential energy surfaces based on a constrained Newton-Raphson algorithm. This procedure is based on work by Simons et al. and by Cerjan and Miller but which has not previously been applied to systems with more than a few degrees of freedom. Attention is paid to practical aspects of efficiency and stability in finding both transition states and minima. For cyclooctane, using the AMBER empirical potential energy function, we can “walk” from the global (boat-chair) minimum to find three other low-energy conformers and nine low-energy transition states that connect them. We also report stationary points for N-methylglycylacetamide(”glycyl dipeptide”) and N-methylalanylacetamide (“alanyl dipeptide”). The results indicate that it is possible to find multiple conformational minima in a nearly automatic fashion for molecules of this size. Prospects for using these techniques on larger systems are discussed.
Introduction Many problems in conformational analysis require one to locate multiple (local) minima on molecular potential energy surfaces and to characterize the pthways of interconversion between them. Most common energy refinement techniques are not well suited to this task, since they rely upon a Taylor series expansion about an assumed starting point,’ and hence generally will find only that minimum energy conformer that is closest to the initial point. Global methods of searching for conformational minima exist, but these systematic searches typically become impractical for systems with more than a few degrees of freedoma2 Here we consider a class of methods that “walk” along potential energy surfaces, moving uphill from one local minimum to a nearby transition state (if one exists), and then downhill toward another local minimum. While such methods cannot guarantee that all (1) See, e&: Fletcher, R. ‘Practical Methods of Optimization”; Wiley: New York, 1980. (2) For an example, see: Millner, 0. E.,Jr.; Andcrsen, J. A. Biopolymers 1975, 14, 2159.
local minima will be found, they form an important addition to the usual energy minimization techniques of molecular mechanics. Algorithms of this sort have been used before for molecules with a small number of atoms; here we examine ways in which these can be ”scaled up” to larger systems. The idea of finding stationary points by walking uphill from one minimum, then downhill toward another one, is not new. In 1971, Crippen and Scheraga proposed a “method of gentlest a ~ c e n t ” in , ~ which a short step along the lowest normal mode direction is followed by energy minimization in the subspace orthogonal to the search direction. Iteration of this sequence moves away from a minimum along a low-energy pathway. Although this procedure does not locate transition states, it can tell when such a point has been passed, such that ordinary minimization techniques should lead to a new minimum. A similar algorithm has been described by HilderbrandtS4 ~
(3) Crippen, G. M.; Scheraga, H. A. Arch. Biochem. Biophys. 1971,144, 462. (4) Hilderbrandt, R. L. Comput. Chem. 1977, 1, 179.
0022-3654/8S/2089-4020$01 .SO10 0 1985 American Chemical Society
The Journal of Physical Chemistry, Vol. 89, No. 19, 1985 4021
Large-Molecule Potential Energy Surfaces In 1981, Cerjan and Miller5 proposed a modified NewtonRaphson technique in which step-length constraints are added by means of a Lagrange multiplier. For proper choices of the desired step length, this can lead to an automatic procedure for moving up a “stream-bed”, i.e., for increasing energy along one direction while simultaneously minimizing along other degrees of freedom. Simons et a1.6 suggested determining the step length by establishing a “trust radius” inside of which a local quadratic approximation to the potential is sufficiently accurate. The numerical examples considered in these papers were very encouraging, but contained at most a half dozen degrees of freedom. Here we discuss modifications we have found to be useful in dealing with systems having 30-70 degrees of freedom, using an empirical force field to represent the molecular potential energy surface.
Computational Method In order to establish notation, a brief outline of the constrained Newton-Raphson method will be given here. We begin with a Taylor series expansion of the potential energy about a point whose Cartesian coordinates are xo: V(x
+ xo) = V(x,) + g x + Y2xTHx
(1)
Here g and H a r e the first-derivative (gradient) and second-derivative (Hessian) arrays, respectively. Setting dV/dx = 0 in eq 1 yields the ordinary Newton-Raphson update x = -Rig
(2)
provided that H has an inverse. Aside from the difficulty of calculating and inverting the Hessian matrix, there are two problems with the Newton-Raphson algorithm for the present problem: (i) the step length (Le., the norm of x) may be large enough that the truncation of the Taylor series expansion in eq 1 is no longer sufficiently accurate; (ii) if the Hessian matrix is positive definite (as it will be near a minimum), the NewtonRaphson step will move toward a minimum energy point, rather than away from it as desired. Both of these problems can be ameliorated by constraining the step length to some predetermined value A. This is most easily carried out with the use of a Lagrange multiplier. Consider the Lagrangian function L:
that the eigenvalues are ordered such that bi < b,+,),then the step x will have the characteristic that it will be directed along the gradient, or uphill, in the v1 direction (since (X - b , ) is positive) and will oppose the gradient, or go downhill, along all other eigenvector directions (since (X - bi) is negative for i > 1). If, furthermore, X is chosen20 that 1/2b1< X < II2b2,then eq 6 shows that the contribution to V - V(x,,) will be negative for all directions except along the first eigenvector. This corresponds to a step of a fixed length that proceeds “uphill” along the softest eigenvector mode, while decreasing the energy along all other directions. As one nears a transition state, bl will become negative and X may be chosen to be zero; this yields the ordinary Newton-Raphson step, which is known to converge to transition states if the starting point is sufficiently good. Although this algorithm has given good results for test cases with up to six degrees of f r e e d ~ m ,we ~ , ~found that several practical matters had to be addressed before we were able to locate transition states reliably for larger systems. These matters are discussed in the following paragraphs. Removal of Translations and Rotations. For large molecules, especially those containing ring systems, it is inconvenient to work in internal coordinates, and Cartesian coordinates are commonly used. In Cartesian coordinates, the Hessian matrix at stationary points has six zero eigenvalues corresponding to overall translation and rotation. We have found it useful to employ a “level-shifting” procedure to modify H: H’ = H
+ cD.DT
(7)
< X < b2 (we assume
Here D is a 3 N X 6 matrix whose columns are unit vectors corresponding to translations and infinitesimal rotations,’ and c determines the magnitude of the level shift. We choose c to place the translation and rotation eigenvalues slightly above the highest internal motions. These shifted eigenvalues then do not interfere with the selection of X as outlined above. Other methods may be used to eliminate translations and rotations from Newton-Raphson-like methods. Thomas and Emersons proposed a Lagrange multiplier method which creates a 3N+6 dimensional matrix made up of H and D. Recently, Taylor and Simons9 have proposed a related method in which an auxiliary 6 X 6 linear problem is solved to obtain the Lagrange multipliers related to the constraints limiting translations and rotations. However, the ”extra” eigenvalues of the matrices used in these methods (Le., those beyond the 3N-6 internal modes) may be interspersed with the internal motions, complicating the choice of A. Since the level-shifting method has worked well for us, we have not pursued these other methods. Choice of an Optimal Step Length. The modified NewtonRaphson algorithm as outlined above required as input a value of A, the required step length. For too small a value there may be no useful solution of eq 4b, whereas too large a value may lead to steps outside the range of validity of the Taylor series expansion upon which our procedures are based. Cerjan and Miller proposed taking the smallest value of A for which eq 4b yields a solution having X between b , and b2, Le., they take the shortest step consistent with the algorithm outlined above. An alternative approach is to establish a “trust radius” inside of which the Taylor series expansion, truncated to second order, is accurate to some tolerance. This radius may change as one moves around on a surface, and a method proposed by Fletcher’ to update the trust radius has been used by Simons et aL6 There are several accuracy parameters that need to be chosen to use the Fletcher algorithm, and we have not made an exhaustive search of all implementations. Nevertheless, in a large number of trials, we were unable to achieve satisfactory performance, and we believe this may be related to the nature of the anharmonic portion of our potential. The nonbonded terms (which dominate many of the conformational changes we are interested in) rise very
(5) Cerjan, C. J.; Miller, W. H. J . Chem. Phys. 1981, 75, 2800. ( 6 ) Simons, J.; Jsrgenson, P.; Taylor, H.; Ozment, J. J . Phys. Chem. 1983, 87, 2745. See also: ONeal, D.; Taylor, H.; Simons, J. J. Phys. Chem. 1984, 88, 1510.
(7) Eckart, C. Phys. Rev. 1935, 47, 552. (8) Thomas, M. W.; Emerson, D. J . Mol. Struct. 1973, 16, 473. See also: van de Graaf, B.; Baas, J. M. A. J . Comput. Chem. 1984, 5 , 314. (9) Taylor, H.; Simons, J. J. Phys. Chem. 1985, 89, 684.
L(x,X) = V(X + xg)
+ ‘/ZX(A2 - X*X)
(3)
Setting d L / d x = dL/dX = 0 yield@ x = ( XI - H)-Ig
(4a)
A’ = gT(XI- H)-’g
(4b)
where I represents a unit matrix. Equation 4b may be solved for X (given A), and then (4a) solved for the step x. In general, there will be more than one solution for X in eq 4b, and the steps associated with these different solutions will have quite different properties. An easy way to evaluate the consequences of different choices of A is to perform a unitary transformation to a basis that diagonalizes H. Letting b, and vi be the eigenvalues and eigenvectors of H, eq 4 becomes x = C(X - b,)-yvi.g)g
(Sa)
I
The new potential energy, to second order, is
P:
P = V(X,) + C(X - bj)-*(vj-g)*(X- !/zbj) i If a solution to eq 5b exists such that b ,
(6)
4022
The Journal of Physical Chemistry, Vol. 89, No. 19, 1985
Nguyen and Case e)
00
bl
112b2
02
b2
.I
Figure 1. Plot of step length (A in A) vs. A for a case where b, < 1 / 2 b 2 . The values of X chose! by the Cerjan-Miller algorithm and this work are ~, n marked by XCM and A, respectively. A is defined as ( ~ . x / n ) ’ ’where is the number of atoms.
steeply, with an r-I2 dependence inside the van der Waals minimum. In such regions the difference between the true energy and P m a y be very large. Other regions of conformational space are smoother and more typical of the surfaces commonly used as test cases for iterative optimization algorithms. This variety makes it hard to choose a good trust radius. Even if a reliable trust radius were to be found, however, its use would still require the solution of eq 4b to determine A. We have chosen instead to follow a different path outlined in the next paragraph. Our approach begins by choosing A rather than A as our input parameter. Since the resulting step length may be too large, we have used the following ad hoc procedure for step length corrections: if the proposed step (obtained from eq 4a) is larger than b,,,, the step vector x is scaled to give it a length of b,,,; then if V(x, x) - V(x,)C e,,,, the step is taken; if not, x is repeatedly scaled by a factor a until the energy change becomes less than emax.In this way, both the step length and the energy change have fixed upper limits. For our surfaces this very simple method works quite well. (In practice, no corrections for steps that are too short have been necessary.) For the calculations below, we use a value of b,,, corresponding to a root-mean-square motion of 0.07 A per atom, and values of 0.1 kcal/mol and 0.7 for emaxand a. Choice of the Lagrange Multiplier A. As we discussed above, choosing X in the range b, C X C ‘/2b2(for bl positive) produces a step that increases Palong the first eigenvector direction, and decreases it along all other directions. The simplest approach is to choose X at the mideoint of this interval, Le., 1/2(bl + ‘/2b2). This choice is labeled X in the figure. Figure 1 shows a typical result for a plot of step length vs. X; these calculated values supplement the schematic diagrams presented in earlier paper^.^^^ As shown in the figure, the Cerjan-Miller choice, XCM, gives quite a short step that will be directed more toward the v2 direction (since X = b2). Furthermore, since this choice is greater than Il2b2,eq 6 shows that the new energy (to second order) will increase rather than decrease along the v2 direction; this problem was first pointed out by Simons et a1.6 Choosing X as we do here between bl and 1/2b2yields a step which is still acceptably small and which is directed much more closely along the v1 eigenvector direction. This simple choice of X must be modified in some circumstances. The most obvious problem arises when bl > ‘/*b2,Le., when the desired range of X does not exist. An example of such behavior is shown in Figure 2a. Even if the desired range does exist, problems may arise if b, and ll2b2are too close together, for then X = b, and the step generated may be very large. Both of these problems may be overcome by scaling the coordinate system along the vl direction,6 so that in the new coordinates (denoted by a prime)
+
(x-v,)’ = (x.vl)/n (8) where x is an arbitrary displacement and n is a scaling parameter. This procedure does not change the location of the stationary points on the true energy surface, but does modify the Hessian eigenvalues and the gradient vector to make b,’ = n2bl and (gvl)’ = n(gv,). In practice, we have adopted a suggestion of Simons et aL6 to choose n such that bl’ = b2/4. We perform this trans-
\
Figure 2. Plots of step length (A in A) vs. X for a case where 6, > ‘/*b2, (a) before and (b) after scaling of coordinates.
formation whenever the unscaled eigenvalues are such that bl > b2/3. Figure 2b shows the results of such a transformation: the coordinate scaling leads to modified eigenvalue spectrum that is ”well-behaved” and indeed is very much like that shown in Figure 1. Choice of Alternative Search Directions. In many cases, one may wish to walk uphill in some direction other than that of the lowest eigenvector. For example, in the alanyl dipeptide, the lowest three eigenvector directions correspond to rotations of the methyl groups about their threefold axes; there are indeed transition states associated with these motions, but we are interested in other motions that change the overall conformation of the molecule. Similarly, in cyclmtane, motion along the lowest eigenvalue leads to pseudorotation, and alternative paths must be explored to find other conformational minima. Two approaches can be used to avoid this problem. First, mass weighting of the Hessian matrix will move hydrogen motions to higher frequencies, allowing one to follow heavy-atom conformational changes. This was very helpful for the alanine dipeptide, as discussed below. Since the location of the stationary points on the surface is not affected by mass weighting, arbitrary masses could be used to emphasize desired motions. The calculations reported below all use the “true” atomic masses, but we have found artificial masses to be useful in finding transition states for other systems. l 2 Second, one can choose to “walk” along an eigenvector direction other than the lowest one. This is accomplished by a scaling transformation like that of eq 8, but along a higher mode, so that it becomes the lowest mode in the scaled system. As one approaches a transition state, the “desired” direction generally becomes the lowest mode, and the scaling is no longer necessary. At each step, then, we typically choose to follow that mode with the largest projection on the previous step direction. For example, if we begin to follow the third lowest eigenvector, as the walk continues we generally will switch to the second lowest, and then the lowest, mode until a stationary point is reached. This approach
The Journal of Physical Chemistry, Vol. 89, No. 19, 1985 4023
Large-Molecule Potential Energy Surfaces
Figure 3. Dihedral angles of N-methylalanylacetamide.
4 A
-180
-60
60
180
W
b)
t \
180
-180
-60
60
18 0
CD
Figure 5. Same as Figure 4 for N-methylalanylacetamide. 60
CI
-60
-180
-60
60
180
Figure 4. “Walk” on 4-q energy contour map (in kcal/mol) of Nmethylglycylacetamide (a) along the first eigenvector and (b) second eigenvector directions. Minima are marked by and transition states by A.
was first developed and tested in ref 6, which should be consulted for more details. Our results for cyclooctane show how useful it can be to be able to move away from a minimum in more than one direction. All of the results reported here were obtained by using double-precision arithmetic on a Digital Equipment Corporation VAX 11/750 computer. Reported stationary points have a rootmean-square gradient of less than kcal/(mol A). Results A . Dipeptides. We chose the “dipeptide“ structures Nmethylglycylacetamide and N-methylalanylacetamide for our first test cases. These surfaces are characterized by the presence of two “important” and relatively floppy degrees of freedom, the 4 and $ torsional angles about the central carbon atom (see Figure 3). The 3N-8 remaining internal degrees of freedom are only loosely coupled to these, so that for the most part one can visualize our results on a conventional @-$ map of the type used for many years to study peptide conformations. Nevertheless, the existence
of these extra degrees of freedom makes finding transition states more difficult than for the two-dimensional problem, and we had to try several variations of our basic algorithm before finding a robust procedure capable of finding transition states from a variety of starting points. The potential functions we use are those of the AMBER molecular mechanics programs.’0 For the glycyl dipeptide we used a “united atom” approach for which hydrogens bonded to carbon are not represented explicitly.” This has the simplifying effect of reducing the number of atoms to 11, and of removing the torsional motion of the terminal methyl groups (since CH3 is treated as a single unit). For the alanyl dipeptide and cyclooctane, we used a more recent potential in which all hydrogens are explicitly represented.’* This allows us to study the interference of low-energy but “unimportant” motions (such as methyl torsions) on our ability space. to find transition states in $--I) Although the AMBER functions are simple ones, and are used here primarily as model potentials, they appear to give a good account of the known features of dipeptide energy surfaces. Recently, the AMBER energies of the minima on the alanyl dipeptide surface were compared to ab initio Hartree-Fock results by using 4-31G and 4-21G basis sets.I3 If the Hartree-Fock results are corrected for dispersion effects, the calculations are in good qualitative agreement, with the relative energies for various conformers agreeing within 0.5 kcal/mol. Although this comparison does not indicate that the entire surfaces are nearly the same, it does give us some confidence that the basic features of (10) Weiner, P. K.; Kollman, P. A. J. Comput. Chem. 1981, 2, 287. (11) Weiner, S.J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Alagona, G.; Ghio, C.; Weiner, P. K. J. Am. Chem. SOC.1984, 106, 165. (12) Weiner, S.J.; Kollman, P.A,; Nguyen, D. T.; Case, D. A. manuscript submitted for publication in J . Comput. Chem. (13) Weiner, S.J.; Singh, U. C.; ODonnell, T. J.; Kollman, P. A. J. Am. Chem. SOC.1984, 106, 6243. Schafer, L.; Klimkowski, V. J.; Momany, F. A.; Chuman, H.; Van, Alsenoy, C. Biopolymers 1984,23,2335. For a survey of other recent molecular mechanics studies on the alanyl dipeptide, see: Pettitt, B. M.; Karplus, M. J. Am. Chem. SOC.1985, 107, 1166.
4024
The Journal of Physical Chemistry, Vol. 89, No. 19, 1985
TABLE I: Stationary Points for Dipeptides” re1 energy
0.00 3.21 4.13 4.04* 4.43* 5.20* 7.13*
d A. Glycyl Dipeptide 77 180 66 123 82 115 -1
TABLE 11: Methyl Group Rotations in the Alanyl methyl group re1 xa Xb xc rotated energy 0.00 60.02 61.90 59.93 a 0.62* 120.59 61.90 59.94 b 5.49* 59.97 122.94 60.04 c 0.51* 60.03 61.87 120.41
IL -64 180 35 -126 -4 66 78
-76 69 -161 -62 53 -159 -108 -7 7 62 -120 3 0 141 -95 125 76 125
Dipeptide4
d
IL
284.41 284.33 284.51 284.18
65.86 65.95 68.05 66.25
uEnergies in kcal/mol, angles in degrees. Values marked with an asterisk are transition states. See Figure 3 for the definition of dihedral angles.
B. Alanyl Dipeptide 0.00 0.68 3.22 3.60 4.39 5.22 4.47* 4.49* 5.45* 5.83* 6.74* 7.04* 7.37. 7.46* 7.95’ 8.41* 9.31*
Nguyen and Case
66 -64 169 -40 42 -63 I46 3 2 -64 -79 80 -60 -119 -138 118 72
LI Energies in kcal/mol; angles in degrees; values marked with an asterisk are transition states.
the (gas-phase) dipeptide surfaces are realistic. Figures 4 and 5 show “adiabatic” maps in which the 4 and angles are fixed and all other degrees of freedom are relaxed.” Superimposed on top of these contour maps are the locations of stationary points on the full 3N-6 dimensional surface and the paths taken by our search program to find them. For both glycyl and alanyl dipeptides we began with the molecule in an extended conformation with standard bond lengths and angles, and used ordinary Newton-Raphson steps to locate the nearest local minimum. We then used the procedure outlined above to walk uphill along the lowest eigenvector direction, with starting steps both parallel and antiparallel to this direction. Results are given in Figures 4a and 5a, which show smooth motion connecting (via transition states) the extended conformation with one local minimum characterized by a seven-membered ring closed by a hydrogen bond (C,), and a second minimum near the a-helix region in longer polypeptides. Energies and geometrical parameters for these stationary points are given in Table I. In general, about 10-20 steps were required to move from one stationary point to another. The lowest eigenvector directions for these molecules are primarily directed along $, Le., the paths are mostly vertical in Figures 4a and Sa. By following the second-lowest mode, which is mostly along 4, we can generate “horizontal” paths that search new regions of conformational space. As shown in Figures 4b and 5b, these paths find new stationary points, not previously characterized for dipeptides. However, as Table I indicates, these are of fairly high energy relative to the global C7 minima, and are unlikely to be of great importance. Nevertheless, the results indicate that one can systematically search a surface for new stationary points without needing to know where they are in advance. Table I1 gives results for three transition states found by following normal modes corresponding to rotations of the methyl groups in the alanyl dipeptide. Since the nature of methyl rotations in proteins may be probed by neutron diffraction and by N M R relaxation measurements, calculations of this sort may be of considerable interest in estimating barriers to this motion. As might be expected, rotations of the end methyl groups (angles xa and xo in Figure 3) have low barriers (0.5-0.6 kcal/mol) whereas it takes 5.5 kcal/mol to rotate the “side-chain” methyl when the backbone is near the C7 conformation. The $ angle opens up by
+
2.2’ during this latter motion (see Table 11); if a rigid rotation about Xb is performed with fixed backbone angles, the computed barrier is 6.5 kcal/mol. This difference between ”rigid” and ”adiabatic” barriers is very often important in predicting conformational flexibility in proteins and other large molecules.14 The usual method for estimating the adiabatic barriers involves repeated energy minimization calculations in which an internal degree of freedom (xbin this case) is fixed at various values. The methods described here can take the place of such calculations, and may be especially useful in cases of low symmetry where the location of the barrier is not known in advance, or where a concerted motion of several angles is involved. B . Cyclooctane. The cyclooctane molecule has fewer heavy atoms than the alanyl dipeptide, but is a more complicated conformational problem because of the constraints imposed by the ring system, and because of coupling between ring inversion and pseudorotation processes. There are now five “important” degrees of freedom, p, el, e2, &, and q53, which describe the displacements of the methylene groups from a planar reference ge01netry.I~ (Coordinate p represents an overall deviation from planarity; 8, specifies the relative amount of crown character and 8, the amount of boat-boat and twist-twist character; 42and 43give the contributions of bending and twisting.) Because there are more than two degrees of freedom, the conformations cannot be mapped in the simple way used above for the dipeptides. “Cuts” through the five-dimensional hypersurface can be very useful, however, in depicting parts of the conformational space.I6 The conformations of cyclooctane have been the subject of extensive study,I7 and it is not our intention to give a complete account of this subject here. Most empirical potential energy functions, beginning with the early work of Hendrickson,Is predict the boat-chair form to be the global minimum, and this appears to be in best accord with the vibrational spectrum.I6 There is no general agreement on the relative energies of other conformers, but the twist-chair-chair, twist-boat-chair and C2 boat-boat forms are likely candidates for local minima with energies within 5 kcal/mol of the boat-chair f ~ r m . ’ ~It . is~ possible ~ that the twist-chair-chair form has an energy low enough to have a significant population at room temperature. It is by no means a simple job to generate structures that correspond to the various “allowed” conformers for cyclooctane. Hendrickson has devised a systematic procedure for generating conformations of medium-sized rings containing either a plane or an axis of symmetry.’* This scheme is restricted by its symmetry requirements, and becomes unwieldy to apply if heteroatoms or double bonds are inserted into the structures. Distance geometry methods have also been used in an attempt to generate conformers automatically for cyclooctane and other ring systems.I9 For (14) See, for example: Gelin, B. R.; Karplus, M. J . Am. Chem. Soc. 1975, 97, 6996. Proc. Natl. Acad. Sci. U.S.A. 1975, 72, 200. (15) Pickett, H. M.; Strauss, H. L. J . Chem. Phys. 1971, 55, 324. Offenbach, J. L.; Strauss, H. L.; Graveron-Demilly, D. J . Chem. Phys. 1978, 69, 3441. (16) Pakes, P. W.; Rounds, T. C.; Strauss, H. L. J . Phys. Chem. 1981,85, 2469. (17) For a recent review, see: Strauss, H. L. Annu. Rev. Phys. Chem. 1983, 31,301. See also: Ivanov, P. H.; bsawa, E. J . Comput. Chem. 1984, 4 , 307. (18) Hendrickson, J. B. J. Am. Chem. SOC.1964,86,4854; 1967,89,7036, 7047.
The Journal of Physical Chemistry, Vol. 89, No. 19, 1985 4025
Large-Molecule Potential Energy Surfaces A
w- w - 9 Crown
TFC
TABLE 111: Energies of Cyclooctane (kcal/mol)
cc
J T3
“A” twist-chair-chair chair-chair
crown T3 “B” twist-boat boat-boat boat T2 “C” boat-chair T1 twist-boat-chair chair twist-chair
C
Figure 6. Interconversion pathways of cyclooctane.
cyclooctane this “embedding” procedure was able to find structures near to the “extended” or “noncongested” conformers, but additional input information (in the form of dihedral angle constraints) was required to build structures near to the boat, chair, and twist-chair forms. Furthermore, since the distance geometry approach does not depend upon energy, subsequent energy refinement calculations are needed to find stationary points on the potential energy surface; as discussed above, energy minimization techniques based on Taylor series expansions may often miss transitions states and proceed to nearby local minima.20 Our results for cyclooctane are summarized in Figure 6. We began with a model-built conformation in the boat-chair (BC) form (see ref 18 for a discussion of the nomenclature of cyclooctane conformers). As is commonly the case in attempts to build ring systems from internal coordinates, our initial structure was not very good, but minimized easily to a structure which turned out to be our global minimum. The various pathways we found in moving away from this structure are illustrated in Figure 6, and the geometries and energies of the stationary points are collected in Tables I11 and IV. Following the lowest eigenvector direction away from BC leads to a transition state (“Tl”) which has C1point group symmetry and which is 2.8 kcal/mol above the boat-chair form; proceeding downhill in the same direction yields a twist-boat-chair (TBC) minimum, lying 1.7 kcal/mol above BC. Continuing further leads to a pseudorotated T1, and then back to a pseudorotated BC form. Although it has been recognized for many years that the pseudorotation pathway for cyclooctane alternates between the BC and TBC forms,’618 few estimates of the transition-state geometries or energy barriers to this motion have been made. Although we know of no way (short of a systematic global search) to prove that T1 is the lowest transition state between the BC and TBC forms, the way in which it was found suggests that this is likely to be the case. Because of its Cl symmetry, it could not be found by procedures like Hendrickson’s,’s which require an axis or plane of symmetry. Following higher eigenvector directions away from BC yields other transition states which in turn yield new minima, as shown (19) Weiner, P. K.; Profeta, Jr., S.; Wipff, G.; Havel, T.; Kuntz, I. D.; Langridge, R.; Kollman, P. A. Tetrahedron 1983, 39, 11 13. (20) This is illustrated by the results for cyclooctane reported in ref 19: by specifying four dihedral angles, a chair form with CZhsymmetry was generated, but with considerable internal energy strain. Upon minimization this structure moved to the nearby twist-boat-chair minimum.
sym Hend.“ PRSb MM2‘ AMBERd D2 1.7 0.20 0.97 0.70 C2,
Du D2 S4 D2d
D2d
c1 c, c1
C2 C2, C,
1.9 2.8
0.58 0.86
0.9 1.4 10.3
2.88 4.03
0.0
0.00
2.0 8.3 8.7
1.70 8.27 6.49
1.15 1.16 11.24 3.12 3.55 10.16 13.71 0.00 2.81 1.66 7.50 8.10
0.84* 0.86** 12.16** 3.72 4.85* 10.61* 10.14* 0.00 2.87* 1.78 1.12* 8.44*
“Reference 18. bReference 16. ‘This work, using parameters of ref 21. dThis work * = saddle point, ** = “peak” (stationary point with two negative eigenvalues).
in Figure 6. We have grouped the resulting conformers into three ”families” of structures, which can interconvert relatively easily among themselves and which are separated from other families by higher barriers. At the top, labeled “A”, are the twistchair-chair (TCC), chair-chair (CC), and crown forms. As in previous potential calculations,16 the TCC form has the lowest energy, and is the only local minimum in this family. At the center of Figure 6 is the “B” family, which consists of the boat (B), twist-boat (TB), and boat-boat (BB) forms, with the TB form as the local minimum. Finally, at the bottom of Figure 6 are the BC pseudorotation pathways described above, plus the chair (C) and twist-chair (TC) transition states, which are higher energy transition states along a TBC pseudorotation pathway. This division into families is based on energetic criteria, but is similar to a classification proposed earlier16 based on an analysis of the five “important” coordinates described above. Although a large number of named conformers have been identified for cyclooctane, only four appear to be relative minima, and these are shown in the central column of Figure 6. The barrier between BC and TBC is less than 3 kcal/mol, whereas the barriers between families are 10-12 kcal/mol. Table I11 gives relative energies for the various conformers with four different empirical force fields. The MM2 results were determined here by starting from the AMBER stationary points and refining in the MM2 force field.2’ The largest deviation between MM2 and AMBER is for the “T2” conformation. The agreement between the two sets of energies for the remaining 11 conformers listed in Table 111 is very good, with a root-mean-square difference of 0.6 kcal/mol; This agreement is not overly surprising, since many of the hydrocarbon parameters in the AMBER force field were taken from MM2, but it does serve to document the close relation between the two force fields for hydrocarbons. Although all four force fields give the BC form as the global minimum, the early calculations of Hendrickson are the only ones to place the “C” family at such a low energy relative to BC. The three more recent force fields are all quite similar, although the potential of Pakes, Round and Straws places the TCC form lower than do the MM2 and AMBER force fields. All three of the recent force fields would lead to the conclusion, however, that in thermal equilibrium at room temperature, significant population of the TCC form would be found. Further studies of the vibrational contributions to the thermodynamic properties of the BC and TCC forms are in progress. Discussion
We have described here our experience with a constrained Newton-Raphson approach to finding stationary points for large molecules. Several other methods have been proposed to find transition states on potential energy surfaces.22 The choice of (21) Allinger, N. L. J . Am. Chem. SOC.1977, 99, 8127. Our calculations used program No. 395 from the Quantum Chemistry Program Exchange (Bloomington, Indiana), modified for a Digital Equipment Corporation VAX computer by Profeta, Jr., S . ; Pensak, D. A.
4026 The Journal of Physical Chemistry, Vol. 89, No. 19, 1985
Nguyen and Case
TABLE IV: Dihedral Angles of Cyclooctane“ twist-chair-chair chair-chair crown T3 twist-boat boat-boat boat T2 boat-chair T1 twist-boat-chair chair twist-chair
w1
w2
w3
a4
u5
O6
w7
W8
112.55 92.60 88.62 114.79 67.98 51.51 74.29 56.82 43.22 -1 1.07 -48.45 0.00 38.28
-86.59 -85.58 -89.01 -43.34 -32.48 -51.51 0.00 -70.45 -102.42 -75.29 -47.97 -76.86 -107.65
64.86 85.58 89.66 -12.22 -67.98 -51.51 -74.29 -3.94 69.66 83.09 93.33 119.25 107.65
-86.59 -92.60 -89.27 -43.34 32.48 51.51 0.00 -12.03 -69.65 -77.52 -91.94 -76.86 -38.28
112.55 92.60 88.62 114.79 67.98 51.51 74.29 84.39 102.42 104.82 93.33 0.00 -38.28
-86.59 -85.58 -89.02 -43.34 -32.48 -51.51 0.00 -38.78 -43.22 -58.89 -47.98 76.86 107.65
64.86 85.58 88.06 -12.22 -67.98 -51.51 -74.29 -74.64 -65.96 -50.89 -48.45 -1 19.25 -107.65
-86.59 -92.60 89.27 -43.34 32.48 51.51 0.00 57.15 65.96 98.33 119.14 76.86 38.28
“Angles in degrees.
method depends to some extent on the type of energy surface being studied, and on the ease of computing first and second derivatives. For quantum mechanical calculations, it is still relatively expensive to calculate second derivatives, so that gradient methods (requiring only first derivatives) are attractive. In molecular mechanics calculations, on the other hand, second derivatives are more readily available (although still more expensive to compute than first derivatives), and Newton-Raphson minimization procedures have been used in this context for many years. A choice of method may also depend upon the amount of prior knowledge one has about a potential energy surface. The early work by McIver and KormonickiZ3(using semiempirical potential energy functions) used direct minimization of the norm of the gradient vector, and required a good initial guess for the transition-state structure. Other approaches assume prior knowledge of two minima and attempt to find a stationary point near some path that connects them.24 Such methods will be more difficult to apply to larger molecules, where knowledge of the potential energy surface is meager. We envision cases where only a single minimum is known, probably generated by energy refinement from some postulated or experimental starting structure. It may be expected, however, that a price must be paid for this generality, and that a constrained Newton-Raphson technique may be less efficient than other methods for cases in which reasonable initial guesses to the location of transiton states can be generated. Although the examples presented here contain 10-25 atoms, we have successfully used this technique for systems with as many as 90 atoms. It is worth remembering that for large molecules, most of the computer time in this and other second-derivative calculations is spent doing linear algebra, rather than evaluating derivatives. Consider the computing times for the AMBER mo(22) For a recent review, see: Bell, S.;Crighton, J. S. J . Chem. Phys. 1984, 80, 2464. (23) McIver, J. W.; Kormonicki, A. J . Am. Chem. SOC.1972, 94, 2625. (24) Halgren, T. A.; Lipscomb, W. N. Chem. Phys. Lert. 1977, 49, 225. Bell, S.: Crighton, J. S.; Fletcher, R. Chem. Phys. Lett. 1981. 82, 122; see also ref 22.
lecular mechanics potential for a 90-atom fragment of DNA: evaluation of the energy requires 0.1 s of CPU time in double precision on a VAX 11/780;calculating the gradient (first derivative) takes about l s, and evaluation of the second derivatives about 10 s. In contrast to these, it takes about 120 s to determine the eigenvalues of the Hessian matrix with a Givens-Householder algorithm, and calculating all of the eigenvectors takes an additional 830 s. Solution of a set of 270 simultaneous linear equations (as in eq 2 or 4a, using an LU decomposition) takes about 145 s. Since we are dealing with an iterative method, it is important to make each step as efficient as possible. In the implementation described above, one needs to compute several low-lying eigenvalues and eigenvectors, and then solve a set of linear equations for the search direction; we find that our computational times per step are about twice that of a conventional Newton-Raphson minimization step.25 A basic goal of this work was to find an algorithm which, starting from a single conformer, would be able to locate in a systematic fashion other low-energy stationary points on the surface-and to do this with no input “knowledge” of where these stationary points might be. We feel that we have been successful in this objective, and in addition are able to classify the conformers we find into “families” of structures that can interconvert with fairly low barriers. These techniques should be applicable to a wide variety of conformational analysis problems. Acknowledgment. We are grateful to Agit Banerjee, Hugh Taylor, and Herb Strauss for many useful discussions. This work was supported by the University of California Cancer Research Coordinating Committee. D.A.C. is an Alfred P. Sloan Foundation Fellow. Registry No. Cyclooctane, 292-64-8; N-methylglycylacetamide, 7606-79-3; N-methylalanylacetamide, 227 15-68-0. (25) A related method which may in principle be more efficient for large molecules is given by: Banerjee, A.; Adams, N.; Simons, J.: Shepard, R . J. Phys. Chem. 1985, 89, 52.
