Research: Science and Education
Energy Contour Plots: Slices through the Potential Energy Surface That Simplify Quantum Mechanical Studies of Reacting Systems
W
Andrew G. Leach*,† Department of Chemistry and Biochemistry, University of California, Los Angeles, CA 90095-1569;
[email protected] E. Goldstein Department of Chemistry, California State Polytechnic University at Pomona, Pomona, CA 91768
†
Current address: AstraZeneca, Mereside, Alderley Park, Macclesfield, Cheshire, SK10 4TG, United Kingdom.
www.JCE.DivCHED.org
•
dictions of a particular level of theory and a “route map” such that the key structures (minima and transition states) can be optimized readily, leading to quantitative predictions of, for example, the relative rates of competing reactions. Throughout this article, points that were highlighted by the contour plots as interesting minima or transition states are marked with a + sign. While the reactions highlighted in this article are good starting points, the techniques that are highlighted are more important as they provide entrances into the study of many reactions in a generic sense and therefore should allow students to work independently but with a clear sense of how to proceed. Description of the Method When beginning an ab initio study of a chemical reaction, one is in effect presented with a potential energy surface and required to locate the stationary points that are relevant to the reaction of interest. Any given set of N atoms (with N > 2) define a potential energy surface of 3N − 6 dimensions. For any systems of experimental interest to organic chemists, this will be sufficiently complex to be beyond comprehension. To study the surface, it is therefore essential to be able to project it into fewer dimensions in a way that captures the key points. When studying a reaction, the optimal projection would be in to the one dimension that corresponds to the reaction coordinate to obtain the kind of plots that are familiar to all chemists such as that shown in Figure 1.
transition state
E
Computational chemistry lab classes form an expanding part of many undergraduate chemistry courses.1 Similarly, a familiarity with the tools provided by theory is an essential part of thorough graduate studies in all areas of chemistry. This reflects the increasing use of the techniques of theoretical chemistry as accompaniments to experimental research (1). This expansion has been fueled by the decreasing cost and increasing speed of the requisite computational resources. Authors who previously may not have considered incorporating calculations into their research now regularly use calculations to provide a more complete picture (1). Companies marketing quantum mechanics packages have designed programs explicitly for these undergraduate classes and for graduate-level applications for those not specializing in computations (2). The focus of undergraduate computational lab classes and even graduate-level courses is frequently on simple examples that are somewhat distant from material in other classes. While organic chemistry classes focus on understanding reactivity, which depends on the relative energy of transition states, intermediates, and products, computational and theoretical chemistry lab classes tend to involve relatively simple calculations, predominantly concerning themselves with minima on the potential energy surface: reactants, products, and conformational spaces or with simple reactions such as those between atoms and diatomics. The theory behind such calculations can usually be presented in fairly rigorous detail. Molecular mechanics and even self consistent field calculations have fairly straightforward theoretical bases that could form the material of accompanying lectures. However, calculations that will provide even a qualitative picture of reactivity usually rely on higher levels of theory that are not so easily or transparently communicated. Nonetheless, studies employing computations can form the basis of highly productive undergraduate-level research projects as well as PhDlevel work. It would therefore be attractive to introduce these techniques and provide a toolkit for students that simplifies the location of the key species for a reaction of interest. This article introduces a particularly effective method for discovering the qualitative predictions that a particular level of theory makes: energy contour plots as simplified versions of the potential energy surface. Articles advocating effective means to visualize the potential energy surface have appeared in this Journal before but for chemically simpler systems (3, 4). These surfaces provide an overview of the qualitative pre-
reactant
product
Reaction Coordinate Figure 1. A projection into one dimension of the potential energy surface for a reaction. Each point on the curve is a minimum in every dimension other than the one shown.
Vol. 83 No. 3 March 2006
•
Journal of Chemical Education
451
Research: Science and Education
In this ideal case, the minimum energy path linking a reactant and product is mapped: each point on the line is also a minimum in every other dimension other than the one shown, such that the maximum corresponding to the transition state is a maximum in one dimension only and the minima are true local minima. The plot can only be obtained once the reaction pathway has been identified. A number of approaches have generally been taken to work around this: • A long handed plot of energy variation with internal coordinates can be used but historically this was limited to the study of simple reacting systems (3). • Most computational studies nowadays use an optimization algorithm to locate the nearest stationary point to an initial guess (5). This relies on a good guess at the geometry, which can usually be achieved relatively easily for minima but can be challenging even for an expert. When a transition state is located, it is often not certain that it is the correct structure and a recent example showed that this method can misidentify transition states (1b). IRC (intrinsic reaction coordinate) calculations or related variants (6) are often used to verify the nature of transition states but we have found such calculations to be problematic to perform and often inconclusive. • Other methods involve providing reactant and product geometries but these are difficult to implement and variable in performance (7).
The vastly increased computer power now available even to undergraduate students leads us to propose turning to systematic searching of the potential energy surface using energy contour plots as a complement to searching for minima and transition states (8). In this way, good starting geometries can be provided for a subsequent optimization step. The role of the optimized points will be well understood and they will be located more rapidly. This is particularly impor-
tant for those new to computational chemistry and those without sufficient chemical intuition to make appropriate guesses of a starting geometry. Once an archetypal reaction has been studied, more complex and highly substituted examples may be studied by analogy and without the need to plot further surfaces. To arrive at an adequate projection, an appropriate internal coordinate (or coordinates) must be chosen. These may be interatomic distances, angles, or dihedrals. If studying a bond dissociation for instance, a plot of how the energy varies with the distance corresponding to the bond length (X⫺Y) can map the relevant features (Figure 2). This distance is not strictly the reaction coordinate. Although the reaction coordinate is a linear combination of all the internal coordinates, it would be dominated by the X⫺Y distance. In most reactions, the reaction coordinate will be dominated by more than one internal coordinate. For instance, the Diels–Alder reaction (Figure 3), discussed in more detail below, will be described well by the two interatomic distances corresponding to the bonds that form during the cycloaddition (C1⫺C6 and C4⫺C5). Interatomic distances are easier to visualize than angles or torsions but these other coordinates can be used also (4e, 9). A further constraint must usually be placed on the geometry to maintain broadly the same structure to avoid discontinuities. For instance, the diene fragment must maintain an approximately s-cis geometry in the case of the Diels–Alder reaction. The choice of intramolecular coordinates for projection will often be clear. A number of suggestions in Table 1 should provide a guide to likely choices. The detail of the number of steps and the stepsize for the scan will be determined by the quantity of time available. For instance, the Diels–Alder reaction between butadiene and ethylene was studied with the C1–C6 and C4–C5 distances being scanned from 1.5 Å to 3 Å in 0.1 Å steps and took a few days on a single processor Linux node. This range covers fully formed bonds going to completely rup-
Diels–Alder reaction
ⴙ
E
retro Diels–Alder reaction
diene fragment C2
Y
X
C1
X−Y
C4
C6
R X−Y
Journal of Chemical Education
C5
1
dienophile fragment
Figure 2. A one-dimensional cut through the potential energy surface for a bond dissociation.
452
C3
•
Figure 3. The Diels–Alder reaction between butadiene and ethylene. The numbering adopted throughout is indicated as well as the atomic arrangement that is maintained throughout the energy contour plot in which the section labeled diene fragment must maintain an s-cis arrangement. A possible diradical intermediate, 1, for the stepwise reaction is shown.
Vol. 83 No. 3 March 2006
•
www.JCE.DivCHED.org
Research: Science and Education
tured bonds. Although 1.5 Å to 2.5 Å is often a good first guess for the required range, one must be prepared to adapt to the output from the scan to ensure that regions of interest are fully mapped—if relevant features are still being revealed or are only partly revealed the plot should be extended. For instance, in the plot shown in Figure 9 for a different reaction, a scan of C1⫺N down to 1.5 Å was not sufficient to see one edge of the valley along the right hand side just as the initial scan of C2⫺N out to only 2.35 Å was not sufficient to see the valley edge along the top. For interatomic distances, 0.1 Å steps are often ideal but larger steps can be used so long as the resulting plots are scrutinized more closely. Torsions can frequently be scanned in 15⬚ steps. The selection of these details depends on the balance between thoroughness and speed that will be determined by the resources available. A description of how we prepared the plot in Figure 4 is provided in the Supplemental Material.W Although it is not the subject of this contribution, it is critical to choose an appropriate level of theory to use in the study (14). This requires a consideration of the balance between accuracy and computer time. No method will be perfect and the shortcomings should be investigated before undertaking a study. Energy contour plots are not limited to just stationary points and so the shortcomings of a method in other regions may also be important.2 For instance, whereas B3LYP/6-31G* performs well at reasonable computational cost for hydrocarbon pericyclic reactions, it does not deal adequately with weak intermolecular interactions such that areas of the plot corresponding to large separations may not be chemically accurate (15). Reports concerning the application of theory to other classes of reaction suggest that other
Figure 4. An RB3LYP/6-31G* energy contour plot generated by scanning the C1⫺C6 and C4⫺C5 distances of cyclohexene. The + symbol indicates a point of interest for subsequent optimization.
Table 1. Reaction Classes Likely To Be Suitable for Energy Contour Plots Suggested Internal Coordinates
Reaction Class
r1
Atom + diatomic, e.g., H + H2
H
SN2 reaction
r2 H
r1
Sigmatropic shifts, e.g., Cope reaction
ⴚ
r2
Require nucleophile (Nu) and leaving group (LG) on opposite sides of reaction center. Colinearity need not be imposed as the reaction center will impose preference. Usually require inclusion of solvation to achieve adequate theoretical treatment. Intramolecular reactions may be best target (e.g., Payne rearrangement) (11).
LG
r1
r1
Constrain to be co-linear in order to maintain two dimensions (3). Theoretically challenging despite simplicity of reaction (10).
H
Nu
Cycloadditions, e.g., Diels–Alder reaction
Notes
Excellent target as accurate theory at reasonable cost available (12). Reaction coordinate dominated by two internal coordinates.
r2
r2
r1
r2
Electrocyclic reactions
Excellent target as accurate theory at reasonable cost available (12). Reaction coordinate dominated by two internal coordinates (13).
Excellent target as accurate theory at reasonable cost available (12). φ1 φ2 φ
Conformational changes
Only one internal coordinate per rotor being studied.
NOTE: This table should act as a guide and not a definitive list.
www.JCE.DivCHED.org
•
Vol. 83 No. 3 March 2006
•
Journal of Chemical Education
453
Research: Science and Education
In this article we will use some examples from cycloaddition chemistry to show how plotting appropriate energy contour plots provides an easy way to begin to study a particular reaction. The authors have studied a number of cycloaddition reactions and it is the lessons learned during the course of these studies concerning good ways to work methodically through a quantum mechanical study that we wish to communicate in this section. Cycloaddition reactions are an ideal case for the application of such an approach because the two interatomic distances corresponding to the bonds formed in the cycloaddition form two natural coordinates to describe the reaction and dominate the real reaction coordinate. The examples below all employ B3LYP/6-31G* that should provide a good description of cycloaddition reactions (12, 16).3 The first example is the Diels–Alder reaction, which will be described in detail and will demonstrate the utility of energy contour plots as enabling tools for students to perform computational studies of a reaction of interest. The subsequent two examples arose during the course of research investigations and will highlight the application of energy contour plots in research.
Example 1: The Diels–Alder Reaction The best known of the pericyclic reactions is the Diels– Alder reaction, which is an allowed [4+2] cycloaddition. The reaction was discovered in 1928 and was recognized to be of such importance that its inventors were awarded the Nobel Prize for Chemistry in 1950 (17), consequently, there have been many theoretical studies of the reaction (18). Experimental kinetic isotope effects for the reaction between isoprene and maleic anhydride were measured and agree well with those computed by B3LYP/6-31G* suggesting that this level of theory provides a good description of the transition state and the reactants (19), as well as predicting the reaction energetics and activation energy (12). As with most pericyclic reactions, one could propose that stepwise reactions might also operate alongside the concerted reactions that are treated by the Woodward–Hoffmann rules (20). How would one go about studying the Diels–Alder reaction with theoretical methods in order to learn more about this reaction? One good way is to start with the simplest Diels–Alder reaction, between butadiene and ethene (Figure 3). The product of this reaction is cyclohexene, which has a C2 symmetric lowest energy conformation that can be easily optimized at 454
Journal of Chemical Education
•
23.4
6.7 0.0
3.3
/ (kcal/mol)
Examples: Cycloadditions
an appropriate level of theory and is a good starting point for an energy contour plot. For this method, it is not essential to have the structure optimized at the same level of theory that will be used for the scan, as constraints will be imposed during the scan. It is key that the general atomic arrangement is correct, as this will determine the structures obtained. An energy contour plot of how the energy varies as the C1⫺C6 and C4⫺C5 bonds are stretched where the atomic arrangement remains approximately cyclohexene-like will reveal what is predicted by the chosen level of theory for this reaction. We have calculated such a surface using B3LYP/631G*. Figure 4 shows the contour plot using entirely restricted calculations in which no open-shell (radical and diradical like) structures are treated. In this plot, as the 3D pictures reveal, the atoms maintain a roughly cyclohexene like arrangement. It can be seen that as the two bonds are stretched the energy increases and that this increase is minimized by stretching the two bonds simultaneously. In the center of the plot, the energy reaches a maximum, which is symmetrical. This is the Diels–Alder transition state and is concerted and synchronous: the two bonds form to the same extent. Continuing to stretch the bonds towards the top left hand corner moves towards a complex of s-cis-butadiene and ethene.4 With this plot in hand, the key transition state is identified. The geometry corresponding to the center of the plot can be extracted from the output file and optimized fully at the same level of theory. This optimization will be quick and reliable because the starting guess at the geometry is good and, most importantly from the perspective of students, the
E
levels may be more appropriate (14d, e). Some have even modified density functional theory methods to better predict kinetic parameters (14f ). Reactions involving charged species or species with significant charge separation may require a treatment that includes solvation. Reactions involving charge transfer will be challenging and may require inclusion not only of solvation but of open-shell species also. Energy contour plots may help deduce that there is a problem with the chosen level of theory by highlighting the absence of barriers or discontinuities in ground-state electronic structure that usually characterize these problem cases. These problems may also manifest themselves in a crashed job or never-ending optimization step.
−38.6
Reaction Coordinate Figure 5. The B3LYP/6-31G* energetics of the Diels–Alder reaction between butadiene and ethene. The energies are enthalpies in kcal/mol relative to s-trans-butadiene and ethene. The enthalpies of reactants are the sum of the enthalpies for butadiene and ethene calculated separately.
Vol. 83 No. 3 March 2006
•
www.JCE.DivCHED.org
Research: Science and Education
optimized structure will be well understood; the plot in Figure 4 has revealed that it is the transition state linking cyclohexene with s-cis-butadiene and ethene. A similar process using the C1⫺C2⫺C3⫺C4 dihedral as coordinate could be used to plot a 1D surface for the interconversion of s-cis
Figure 6. A UB3LYP/6-31G* energy contour plot obtained by stretching the C1⫺C6 and C4⫺C5 bonds in cyclohexene and with the guess=(mix,always) command to generate open-shell structures when they are lower in energy than closed-shell ones. The + symbol indicates a point of interest for subsequent optimization.
ⴙ
ⴙ
2
3
endo
exo
Scheme I. The [6+4] cycloaddition reaction between cyclopentadiene and cycloheptatriene.
www.JCE.DivCHED.org
•
and s-trans-butadiene in order to locate the transition state linking these two isomers. In this way, the complete picture of the concerted reaction can be obtained and is shown in schematic form in Figure 5. s-cis-Butadiene is higher in energy than s-trans-butadiene so the observed activation energy should correspond to the difference between s-trans-butadiene and the Diels–Alder transition state. This corresponds to ∆H‡0K = 23.4 kcal兾mol in excellent agreement with experimental data, where ∆H‡0K is found to be 24.5 kcal兾mol (12). The study of the Diels–Alder reaction can be extended by using a method that permits open-shell electronic structures as well as closed shell. We have used unrestricted B3LYP and the guess=(mix,always) command in Gaussian98.3 The modified surface shown in Figure 6 is obtained. Instead of the one reaction channel seen in the RB3LYP surface, three channels are observed, the concerted transition state is still present in the center but now two alternative, higher energy channels are open. In these areas, the lowest energy structure is open shell; these are channels leading towards diradicals like 1 (see Figure 3). The diradicals themselves will not be observed in a plot like this because there is no barrier to their reclosure to cyclohexene without a conformational change, which is without the scope of this plot. However, it can be seen that these channels are higher in energy than the concerted transition state by about 4 kcal兾mol. We can discern more information from the shape of these channels. The concerted channel is rather broad and the diradical ones narrow. A narrow channel corresponds to a vibrational degree of freedom (orthogonal to the path along the “valley” floor) that has a high force constant that in turn corresponds to a high zero point energy, 0.5(k兾µ)1/2, and widely spaced vibrational energy levels. This means that the effective barrier is even higher to the diradical channels than the valley bottoms suggest. This is in relatively good agreement with Zewail’s observations from femtosecond spectroscopy (21). This simple example shows how the beginnings of a mechanistic study of a reaction may be built up without knowledge of the structure of the transition states of interest. With a well understood reaction mapped out with optimized structures, a great deal may be learned. The following two examples illustrate the approach as used by ourselves in the context of actual research projects. For both we give a brief description of our motivation with references provided to the published details. This is followed by a description of the energy contour plots along with what was deduced from them. Some of the conclusions ultimately drawn from these studies over and above what was learned from the contour plots are given for each example by way of a summary.
Example 2: The [6+4] Cycloaddition between Cyclopentadiene and Cycloheptatriene A long-standing interest in the [6+4] cycloaddition prompted a detailed study of the reaction between cyclopentadiene and cycloheptatriene (Scheme I). Woodward and Hoffmann had predicted that the [6+4] cycloaddition should be allowed according to the rules of orbital symmetry control (20, 22). Soon after they made this prediction, a number of [6+4] cycloadditions were observed (23). The [6+4] cycloaddition leading to 2 and 3 is in competition with four alternative allowed [4+2] cycloadditions,
Vol. 83 No. 3 March 2006
•
Journal of Chemical Education
455
Research: Science and Education
three forbidden [2+2] cycloadditions, a forbidden [4+4] cycloaddition, and a forbidden [6+2] cycloaddition. Concerted and stepwise processes involved in all of these reactions were studied (24). Many reactions are computed to be close in energy. In line with Woodward and Hoffmann’s prediction, the [6+4] cycloaddition is computed to be one of the two most favored cycloadditions. Woodward and Hoffmann had further predicted, based on secondary orbital interactions, that the [6+4] cycloaddition would be strongly exo selective (25). This was borne out by the calculations and will be discussed in more detail. The transition state for the [6+4] exo cycloaddition, leading to 3, was easily located and corresponded to an activation energy of 26.0 kcal兾mol. Many prospecting searches for the transition state for the endo [6+4] cycloaddition proved
fruitless. Constraining the search to CS symmetry, that is, enforcing a synchronous bond formation, yielded a second-order saddle point 4 at an energy of 37.0 kcal兾mol relative to reactants. It was unclear whether a transition state existed but was being missed or whether the transition state did not exist. An energy contour plot was able to clarify the situation and revealed that no endo [6+4] transition state exists in this case. A contour plot in which the C1⫺C8 and C6⫺C11 distances were constrained to a range of values between 1.5 Å and 3.1 Å starting from, and with the geometry maintaining an atomic arrangement similar to, the [6+4] endo adduct, 2, was generated and is shown in Figure 7. In this figure the top left hand corner, in which both distances are long, represents the reactants having been brought into one another’s proximity. The bottom right hand corner, in which both distances are short, represents the [6+4] endo adduct 2. Top right and bottom left corners represent the case where only one of the two bonds (C1⫺C8 or C6⫺C11) of the [6+4] cycloadduct are formed. Unlike the situation for the retroDiels–Alder reaction of cyclohexene, this does not correspond to a diradical like 1. Inspection of the geometry in the output file showed that it is the [4+2] cycloadduct 5. The diradical is thwarted by the ready availability of other bonds to form. Finally, at the center, the ridge indicated corresponds to the second-order saddle point 4 that was located when optimizing in CS symmetry. The lowest energy path, marked with solid arrows, causes one of the two bonds of the [6+4] cycloadduct to form. It is unfavorable to form the second bond compared to forming a different C⫺C bond (either C4⫺C9 or C3⫺C10) corresponding to a [4+2] cycloaddition leading to 5. The [6+4] adduct is only reached through a sigmatropic shift. This occurs through a second transition state that corresponds to the
[6+4] endo cycloaddition C11
[4+2] exo cycloaddition
C6
C10 C9
C8
C1
C4 C3
attractive interaction corresponding to forming bond repulsive 2° interaction attractive 2° interaction
Figure 7. The B3LYP/6-31G* energy contour plot generated by scanning the C1⫺C8 and C6⫺C11 distances corresponding to the two bonds that would be formed in the [6+4] cycloaddition between cyclopentadiene and cycloheptatriene leading to 2. The + symbol indicates a point of interest for subsequent optimization.
456
Journal of Chemical Education
•
Figure 8. The HOMO of cyclopentadiene and LUMO of cycloheptatriene interaction corresponding to the [6+4] endo cycloaddition leading to 2 (left hand side) and to the alternative [4+2] exo cycloaddition leading to 5 (right hand side).
Vol. 83 No. 3 March 2006
•
www.JCE.DivCHED.org
Research: Science and Education
ridge part way along the bottom edge (or along the right hand edge), shown in Figure 7 by the break in the arrows. Figure 7 demonstrates that there is no [6+4] endo transition state. Subsequently we were able to rationalize this based on frontier orbital interactions and steric clashes. Woodward and Hoffmann had predicted that repulsive secondary orbital interactions would lead to a preference for exo cycloaddition (25). As indicated in the left hand box of Figure 8, there is a repulsive interaction between C3 (or C4) of the LUMO of cycloheptatriene and C9 (or C10) of the HOMO of cyclopentadiene. This is the HOMO–LUMO pair that are closest in energy. This interaction is a factor in the exo selectivity of the [6+4] cycloaddition reaction and is a key factor responsible for the high energy ridge in the center of Figure 7. With a slight geometric distortion, the repulsive secondary HOMO–LUMO interactions of the [6+4] endo cycloaddition can be swapped for attractive primary and secondary orbital interactions in the [4+2] cycloaddition (right hand side box of Figure 8). A second factor in the contour plot is the steric repulsion between the two CH2 groups that is maximized in the synchronous [6+4] endo cycloaddition, 4.
Example 3: Aziridine N-Oxides in the Ene Reaction of Nitroso Compounds A theoretical study of the ene reactions of nitroso compounds was recently performed by one of the authors (26). A key aspect of this was to establish the role of aziridine Noxide intermediates like 6 (Scheme II). To calculate how aziridine N-oxide intermediates are formed, an energy contour plot was generated using restricted calculations. The two C⫺N bonds of the aziridine N-oxide that would be formed between HNO and propene were stretched. This energy map is shown in Figure 9. On this energy contour plot, the aziridine N-oxide is in the bottom right hand corner of the plot where both C⫺N distances are short. The reactants (propene and HNO) are in the top left hand corner where both C⫺N distances are long. Structures in which only one of the two C⫺N bonds have formed are located in the top right and bottom left hand corners of the plot. This energy contour plot established that the lowest energy path linking propene and HNO to the corresponding aziridine N-oxide is not concerted; the concerted path would pass through the center of the plot and corresponds to the arrow that clearly passes through the high energy area. The lowest energy path is stepwise and involves formation of a bond at the least substituted end of the alkene first, leading to a local minimum (an intermediate) in the top right hand corner of the plot. This is followed by a second bond formation to yield the aziridine N-oxide in the bottom right hand corner of the plot. This plot permitted the key stationary points to be optimized rapidly and a clear view of their role was known from the outset. Subsequent calculations revealed that the open chain intermediate could be optimized with either restricted or unrestricted calculations and a detailed set of studies led us to style this intermediate a “polarized diradical”. This has been described in depth elsewhere (26). By contrast, Adam et al. studied the same reaction type but with p-NO2C6H4NO as the nitroso enophile and obtained only one of the two transition states for this process www.JCE.DivCHED.org
•
(1b). It was subsequently shown that they had missed the second transition state and the open-chain intermediate for this process. It is unlikely that this would have occurred if energy contour plots for the reaction had been calculated.
H
O R
OH
ⴙ
N
R
R
N
−
+
O
N
6 intermediate
Scheme II. The ene reaction of nitroso compounds and a possible aziridine N-oxide intermediate, 6.
Figure 9. The B3LYP/6-31G* energy contour plot generated by scanning the C1⫺N and C2⫺N bonds corresponding to the aziridine N-oxide intermediate formed in the reaction between HNO and propene. The + symbol indicates a point of interest for subsequent optimization.
Vol. 83 No. 3 March 2006
•
Journal of Chemical Education
457
Research: Science and Education
Conclusions This article has highlighted the utility of energy contour plots for a number of purposes principally by way of illustrative examples. They permit a theoretical study of a reacting system to be started in a systematic way; a prerequisite for the uninitiated. A more detailed picture of a reaction according to a particular level of theory can be obtained than a set of optimized minima and transition structures provides. Energy contour plots can establish whether a particular transition state and therefore reaction pathway exists at a given level of theory and can provide useful pointers for understanding the absence of expected transition states. Energy contour plots can establish a mechanism and the likely structure of intermediates and transition states before the search for the optimized species themselves is undertaken. In this way, they permit optimizations to be undertaken more easily, and with a more detailed understanding than would be obtained with a simple hunt for stationary points. The figures plotted in the course of the studies subsequently provide eloquent illustrations when presenting these studies in publication or lecture form. Acknowledgments The authors thank K. N. Houk and his research group for providing the stimulating environment in which the research was undertaken. The National Computational Science Alliance and the UCLA Academic Technology Service are thanked for allocations of computer time and AGL thanks AstraZeneca and the U.K. Fulbright Commission for a fellowship. W
Supplemental Material
Instructions to compute and plot the energy contour plots are available in this issue of JCE Online. Notes 1. For instance, the ACS guidelines for Undergraduate Professional Education in Chemistry (spring 2003 version) lists chemical computation as an essential function of a viable chemistry department and states that “an early introduction to instrumental and computational techniques can give students an immediate sense of how chemistry is done today.” 2. Agreement between dynamics calculations and experimental results provides support for the validity of the non-stationary parts of the PES. Such calculations are extremely demanding even for the relatively efficient density functional methods. A recent study of the reaction of atomic oxygen with methane found that B3LYP provides excellent dynamics results (8). A study of the ene reaction of singlet oxygen (another pericyclic reaction) found that B3LYP provides a surface that predicts kinetic isotope effects in excellent agreement with experiment. In this case, the kinetic isotope effects are not determined by stationary points but by the shape of the surface around a reaction path bifurcation. Kraka and co-workers recently studied the Diels-Alder reaction and concluded that B3LYP provides an excellent description of the whole of the reaction pathway (18).
458
Journal of Chemical Education
•
3. All calculations in this article were performed in Gaussian 98 (Revision A.7), Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian, Inc., Pittsburgh PA, 1998. 4. Cyclohexene itself is C2 symmetrical while the Diels–Alder transition state and this complex are Cs symmetrical. As the bonds are stretched this change in symmetry occurs smoothly. In the forward sense of the Diels–Alder reaction and not revealed by this plot because of the choice of coordinates, after passing through the Cs symmetrical transition state, the reaction must choose to proceed to one or other of the enantiomeric C2 symmetric conformations of cyclohexene (a Cs symmetrical cyclohexene would actually be the transition state linking these two conformations). The two bond lengths are identical in this product structure and so they would appear identical in this plot (27).
Literature Cited 1. See for instance (a) Aggarwal, V. K.; Harvey, J. N.; Richardson, J. J. Am. Chem. Soc. 2002, 124, 5747–5756. (b) Adam, W.; Bottke, N.; Engels, B.; Krebs, O. J. Am. Chem. Soc. 2001, 123, 5542–5548. 2. See for instance (a) Spartan Student Edition Home Page. http:/ /www.wavefun.com/news/spartanstudented.html. (b) Student HyperChem. http://www.hyper.com/products/Student/ default.htm. (c) Gaussian Home Page. http://www.gaussian.com (all accessed Nov 2005). 3. (a) Hulse, J. E.; Jackson, R. A.; Wright, J. S. J. Chem. Educ. 1974, 51, 78–82. (b) Moss, S. J.; Coady, C. J. J. Chem. Educ. 1983, 60, 455–461. (c) Fernández, G. M.; Sordo, J. A.; Sordo, T. L. J. Chem. Educ. 1988, 65, 665–667. 4. (a) Bauer, S. H. J. Chem. Educ. 1999, 76, 440–443. (b) Macomber, R. S. J. Chem. Educ. 1998, 75, 1346–1350. (c) Leventis, N.; Hanna, S. B.; Sotirior-Leventis, C. J. Chem. Educ. 1997, 74, 813–814. (d) Lehman, J. L.; Goldstein, E. J. Chem. Educ. 1996, 73, 1096–1098. (e) Tonge, K. H. J. Chem. Educ. 1988, 65, 65–68. 5. (a) Pulay, P. Mol. Phys. 1969, 17, 197–204. (b) Schlegel, H. B. J. Comp. Chem. 1982, 3, 214. 6. Silva, M. A.; Goodman, J. M. Tetrahedron Lett. 2003, 44, 8233–8236. 7. Foresman, J. B.; Frisch, A. Exploring Chemistry with Electronic Structure Methods, 2nd ed.; Gaussian Inc.: Pittsburgh, PA, 1993; p 46. 8. A recent example is found in Ensing, B.; Buda, F.; Gribnau, M. C. M.; Baerends, E. J. J. Am. Chem. Soc. 2004, 126, 4355– 4365. 9. Suhrada, C. P.; Houk, K. N. J. Am. Chem. Soc. 2002, 124, 8796–8797.
Vol. 83 No. 3 March 2006
•
www.JCE.DivCHED.org
Research: Science and Education 10. Grant, D. M.; Wilson, P. J.; Tozer, D. J.; Althorpe, S. C. Chem. Phys. Lett. 2003, 375, 162–166. 11. Jung, M. E.; van den Heuvel, A.; Leach, A. G.; Houk, K. N. Org. Lett. 2003, 5, 3375–3378. 12. Guner, V.; Khuong, K. S.; Leach, A. G.; Lee, P. S.; Bartberger, M. D.; Houk, K. N. J. Phys. Chem. A 2003, 107, 11445– 11459 and references therein. 13. Black, K. A.; Leach, A. G.; Kalani, M. Y. S.; Houk, K. N. J. Am. Chem. Soc. 2004, 126, 9695–9708. 14. For properties of ground states, see (a) Curtiss, L. A.; Rhagavachari, K. Theor. Chem. Acc, 2002, 108, 61–70. (b) Henry, D. J.; Parkinson, C. J.; Mayer, P. M.; Radom, L. J. Phys. Chem. A 2001, 105, 6750–6756. (c) Foresman, J. B.; Frisch, A. Exploring Chemistry with Electronic Structure Methods, 2nd ed.; Gaussian Inc.: Pittsburgh, PA, 1993. For activation energies see ref 12 and references therein. (d) Kang, J. K.; Musgrave, C. B. J. Chem. Phys. 2001, 115, 11040–11051. (e) Parthiban, S.; Oliviera, G.; Martin, J. M. L. J. Phys. Chem. A 2001, 105, 895–904. (f ) Lynch, B. J.; Truhlar, D. G. J. Phys. Chem. A 2003, 107, 3898–3906 and references therein. 15. (a) Koch, W.; Holthausen, M. C. Hydrogen Bonds and Weakly Bound Systems. In A Chemist’s Guide to Density Functional Theory, 2nd ed.; Wiley-VCH: Weinheim, Germany, 2001; pp 217–238. (b) Wesolowski, T. A.; Parisel, O.; Ellinger, Y.; Weber, J. J. Phys. Chem. A 1997, 101, 7818. (c) Lundell, J.; Latajka, Z. J. Phys. Chem. A 1997, 101, 5004–5009. 16. (a) Becke, A. D. J. Chem. Phys. 1993, 98, 1372–1377. (b) Becke, A. D. J. Chem. Phys. 1993, 98, 5648–5652. (c) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623–11627. (d) Lee, C.; Yang,
www.JCE.DivCHED.org
•
W.; Parr, R. G. Phys. Rev. B 1988, 37, 785–789. 17. For more information see the Web site of the Nobel foundation http://www.nobel.se (accessed Nov 2005). 18. Kraka, E.; Wu, A.; Cremer, D. J. Phys. Chem. A 2003, 107, 9008–9021 and references therein. 19. (a) Diau, E. W.-G.; De Feyter, S.; Zewail, A. H. Chem. Phys. Lett. 1999, 304, 134–144. (b) Beno, B. R.; Houk, K. N.; Singleton, D. A. J. Am. Chem. Soc. 1996, 118, 9984–9985. 20. Woodward, R. B.; Hoffmann, R. Angew. Chem., Int. Ed. Engl. 1969, 8, 781. Woodward, R. B.; Hoffmann, R. The Conservation of Orbital Symmetry; Verlag Chemie: Weinheim, Germany, 1970. 21. Horn, B. A.; Herek, J. L.; Zewail, A. H. J. Am. Chem. Soc. 1996, 118, 8755–8756. 22. Hoffmann, R.; Woodward, R. B. J. Am. Chem. Soc. 1965, 87, 2046–2048. 23. (a) Houk, K. N.; Woodward, R. B. J. Am. Chem. Soc. 1970, 92, 4143–4145. (b) Houk, K. N.; Woodward, R. B. J. Am. Chem. Soc. 1970, 92, 4145–4147. (c) Bhacca, N. S.; Luskus, L. J.; Houk, K. N. Chem. Commun. 1971, 109–111. 24. Leach, A. G.; Goldstein, E.; Houk, K. N. J. Am. Chem. Soc. 2003, 125, 8330–8339. 25. Hoffmann, R.; Woodward, R. B. J. Am. Chem. Soc. 1965, 87, 4388–4389. 26. (a) Leach, A. G.; Houk, K. N. Org. Biomol. Chem. 2003, 1, 1389–1403. (b) Leach, A. G.; Houk, K. N. J. Am. Chem. Soc. 2002, 124, 14820–14821. (c) Leach, A. G.; Houk, K. N. Chem. Commun. 2002, 1243–1255. 27. (a) Basilevsky, M. V. Chem. Phys. 1982, 67, 337–346. (b) Basilevsky, M. V. Theo. Chim. Acta 1987, 72, 63–67.
Vol. 83 No. 3 March 2006
•
Journal of Chemical Education
459