In the Classroom
Computational Studies of Chemical Reactions: The HNC–HCN and CH3NC–CH3CN Isomerizations Arthur M. Halpern Department of Chemistry, Indiana State University, Terre Haute, IN 47809;
[email protected] Advances in computational techniques and the increasing accessibility of these methods to nonspecialists have opened new vistas in undergraduate education. Thus it is now possible for students to explore the application of computational chemistry to molecular structure and spectroscopy and chemical reactivity and dynamics. Although these activities are optimally coupled with formal coursework in computational chemistry or quantum chemistry, it is nevertheless possible for students to become engaged in special projects in order for them to gain new insights into both the capability of computational chemistry techniques, as well as the nature of chemical reactivity and dynamics. For example, many students presented with a simple chemical reaction might wonder, “How does that happen?” Exposing them to the capabilities of computational chemistry methods to explore reaction mechanisms helps them reach satisfying answers to this question. The projects described here can be readily used as the context for problem-based learning experiences for individual students or for those working in groups. In the latter case, students (with faculty guidance, of course) can teach and learn from each other. At the completion of the project the student (or team) should be encouraged to organize and present the results. The experiences gained can readily act as a springboard from which some students can jump into other topics or projects in computational chemistry. Students at various stages in their education can benefit from these learning experiences. They could be juniors taking physical chemistry, seniors in a special topics course, or beginning graduate students. However, particularly motivated beginning students could also profitably participate in such a study group. The advent of computational chemistry packages and related visualization applications that run on stand-alone computers provide the instructor with new, powerful tools to help satisfy students’ curiosity and will perhaps inspire them to ask more sophisticated and detailed questions about molecular properties and reactivity. This article first describes the application of computational methods to the isomerization of hydrogen isocyanide to hydrogen cyanide,
HNC HCN Although this reaction appears very straightforward, it has many of the attributes associated with more complex atom migration processes. Because this reaction involves triatomics (and only two heavy atoms), it is not only easy for students to visualize the process, once elucidated, but it is also very amenable to ab initio calculations at a reasonably high level without demanding impractical computing time or resources. In addition, this isomerization contains all of the fundamental issues associated with chemical reactivity, such as the identification of the stable reactant and product, the transition state, the construction of a potential energy surface, and the calculation of rate constants. Students can then use this information to obtain basic information about the reaction, inwww.JCE.DivCHED.org
•
cluding its thermochemistry, activation energies of forward and back reactions, and estimates of the reaction rates. In other words, this reaction provides students with a comprehensive learning experience in chemical physics that spans and augments much of the standard topics in physical chemistry. This article also describes the logical extension of this exercise to the isomerization of the methyl-substituted compounds, methylisocyanide and methylcyanide,
CH3NC
CH3CN
Students will encounter the challenges and complexities associated with this process, which has been experimentally studied in more detail than the HNC–HCN isomerization. Computational Methods We recommend using a density functional theory (DFT) method (1) because it is a good compromise between rigor and efficiency. The DFT approach attempts to account for exchange energy and electron correlation effects. Various DFT applications use different parameterized functional forms to account for these interaction energies. A widely used DFT method is B3LYP (2), which we will employ in these calculations. A reasonable choice of basis set for this problem is the correlation-consistent polarized valence set of functions at the double-zeta level, cc-pVDZ (3) , which is built up from Gaussian-type functions with two different exponents containing s through f functions. Thus the computational method is B3LYP/cc-pVDZ. We used these methods as implemented in GAUSSIAN 98 (G98) (4). This application can also be run on a Windows platform (G98W) (5). The description of the input files pertains to G98W. Users of G98 or other applications can alter the files accordingly. Setting Up the Problem Students can be provided with an overview of the project that establishes its overall scope, objectives, and methods. This approach should help them avoid getting lost in the details and thus keep up a desired level of excitement, challenge, and (eventual) accomplishment. The following protocol can be used to frame the problem and identify the tactics of its study: • Characterize the reactant and product, optimize their structures, obtain vibrational frequencies, and compare with experiment • Obtain an approximate transition-state structure • Optimize and confirm the transition-state structure • Perform intrinsic reaction coordinate (IRC) (6) calculations on the forward and reverse reactions from the transition state • Visualize the reaction potential energy surface • Use the results to calculate reaction enthalpies, equilibrium constants, and rate constants
Vol. 83 No. 1 January 2006
•
Journal of Chemical Education
69
In the Classroom
The HNC–HCN Isomerization
Reactants and Products The reactant and product are the triatomic molecules, HNC and HCN. By examining the Lewis structures of these compounds, students should be able to predict linear structures. The input information for G98W contains the following components:1 Routing Section: b3lyp/cc-pvdz opt freq Title Section: HCN Optimization (or other entry) Charge and multiplicity: 0,1 Molecule Specification: The molecule specification section consists of one line per atom containing its atomic symbol (or atomic number) and, in this case, its Cartesian coordinates (in Å). Here we have chosen to orient the molecule along the z axis with the C atom at the origin. The z coordinates represent initial guess values, based, for example, on average bond lengths (7) corresponding to H⫺C and C⬅N bond lengths of 1.1 and 1.2 Å, respectively.2 H 0.0 0.0 -1.1 C 0.0 0.0 0.0 N 0.0 0.0 1.2
This calculation first obtains the optimized geometry and then calculates the vibrational frequencies3 and other fundamental information such as the rotational constant and thermodynamic properties. In the case of HNC, the symbols “C” and “N” can simply be switched and the Title Section changed accordingly. These jobs should take less than one minute to complete. The combined results for the HCN and HNC stationary (ground electronic) states are listed in Table 1 along with experimental data. It should be noted that while HCN is a well characterized molecule, HNC is somewhat elusive, even though it is predicted to be stable (with respect to isomerization to HCN) under isolated conditions at room temperature. It is evidently reactive, perhaps rapidly forming HCN in a bimolecular process, although it has been observed in interstellar space, in
Table 1. Optimized Properties of HNC and HCN and Comparisons with Experiment HNC
Property
Calc
a
HCN a
Exp
Calc
rN(C)H/Å
1.0044
0.9938b
1.0770
1.0655g
rCN/Å
1.1780
1.1690b
1.1579
1.1532g
772.1
713.46h
∼v /cm᎑1 1 ∼v /cm᎑1 2 ∼v /cm᎑1
467.5
c
477
2100.3
2029.2
2199.9
2096.7h
3780.2
3620c
3466.0
3311.5h
E /hartree
᎑93.406222
---
᎑93.430176
---
He/hartree
᎑93.386892
---
᎑93.410271
---
Ge/hartree
᎑93.410219
---
᎑93.433126
---
ZPEf/hartree
0.0155271
---
0.0164258
3 d
c
Exp
a
--b
This work. All calculations are at the B3LYP/cc-pVDZ level. Ref 8. Ref 9. dElectronic energy only. eIncludes zero-point energy and f g h thermal correction at 298 K. Zero-point energy. Ref 10. Ref 11. c
70
Journal of Chemical Education
•
comet Hale–Bopp (12) and in other conditions at high temperature at which the equilibrium composition of HNC is high enough to permit its detection. From the calculated standard enthalpies, students can easily determine that the heat of reaction is predicted to be ᎑61.5 kJ兾mol. This result is in remarkably good agreement with an experimental value of ᎑61.9 ± 2 kJ兾mol that has been deduced from ion cyclotron double resonance spectroscopy (13). Students should be encouraged by the good overall agreement with experiment, especially the thermochemistry of the reaction, and be already convinced of the usefulness of computational methods.
Transition State In locating the transition state for the isomerization, students will have to visualize what happens during the process. They will realize that the H atom essentially migrates from one end of the molecule to the other along the pi-system of the N⫺C bond, and there will be a high-energy structure, the transition state, in which the H atom lies at the vertex of a triangle, not necessarily equidistant from the N and C atoms. One straightforward way to estimate the transition-state structure is to perform a scan of the potential energy surface describing the migration of the hydrogen atom. This can be readily accomplished in a calculation that utilizes internal, rather then Cartesian coordinates. Using this approach, the H⫺C⫺N bond angle can be systematically varied, allowing the H atom to pass through the transition-state area. The structure of the internal coordinates, commonly referred to as a z-matrix (14), that is used in the molecule specification section of the calculation is: N C 1 rcn H 2 rch 1 theta rcn 1.2 rch 1.1 theta 175.0 s 15 -10.
The calculation will obtain partially optimized structures and their energies as the H⫺C⫺N bond angle is varied from 175⬚ to 35⬚ in steps of 10⬚. The output file will contain a section with the header “Summary of the Optimized Potential Surface Scan”. By perusing the energies listed, the student can readily determine that the maximum energy occurs when rcn = 1.1922 Å, rch = 1.1833 Å, and theta = 75.0⬚. The actual energy is unimportant because it does not represent the fully optimized transition state. This information is now used to obtain an optimized transition-state structure and its vibrational frequencies. For this compound job, the routing section contains: b3lyp/cc-pvdz opt=(ts,calcfc,z-matrix) freq
The charge and multiplicity and geometry specification are the same as for the scan, but the initial guesses for the transition-state coordinates are updated by the values cited earlier. The keywords direct the program to locate and optimize the transition state in internal coordinates and to calculate the force constants. As the student should expect, the calculation produces a triangular structure, which is depicted in Figure 1. Students should verify that the calculation produces a single imaginary frequency, which arises from a vibrational
Vol. 83 No. 1 January 2006
•
www.JCE.DivCHED.org
In the Classroom
mode having a negative force constant. This result is a requirement for a first-order saddle point, or transition state. The mode involved corresponds, as anticipated, to the inplane rocking motion of the H atom in the triangular transition-state structure. At this point it is very instructive for students to use a visualization application that animates the vibrational modes obtained in the calculation. The properties of this transition state are summarized in Table 2. This species has a smaller H⫺N⫺C bond angle (54.0⬚) as compared with the H⫺C⫺N bond angle (72.1⬚), and, accordingly, a smaller C⫺H bond length (1.197 Å) relative to the N⫺H bond length (1.406 Å), and thus may be said to lie on the potential energy surface closer to the product (HCN) than the reactant (HNC), in an apparent exception to the Hammond postulate (16), which suggests that in an exothermic process the transition state more closely resembles the reactant.
ICR Calculations At this point, with the transition state located and characterized, the student can now perform intrinsic reaction coordinate (IRC) calculations of the HNC–HCN isomerization. The IRC potential, which was introduced by Fukui (6), represents here the minimum energy pathway (MEP) that connects a transition state with a stable entity. The “forward” and “reverse” directions along the IRC potential lead either to HNC or HCN depending on the nature of the algorithm. The IRC, moreover, accounts quantitatively for the displacement of all atoms along this path. The value of the ith IRC coordinate, qi, along the MEP is expressed as qi =
∑ m j ( x j,i
− x j ,0
j
( yj,i + ( z j, i +
)
N
C
Figure 1. The B3LYP/cc-pVDZ optimized transition state for the HNC–HCN isomerization. The properties of this structure are contained in Table 2.
Table 2. Optimized Properties of the HNC–HCN Transition State and Comparison with Literature Values Calca
Ref 15
rCH/Å
1.1968
1.1835
Property rCN/Å
1.1931
1.1867
rNH/Å ∠HCN
1.4062
1.3875
72.084
71.733
∠HNC
54.045
---
∼v /cm᎑1 1 ∼v /cm᎑1 2 ∼v /cm᎑1
1119.3i
---
2058.2
---
2628.3
---
᎑93.355436
---
H /hartree
᎑93.340983
---
Gc/hartree
᎑93.365779
---
ZPE/hartree
0.0106765
---
3
Eb/hartree c
a
This work. All calculations are at the B3LYP/cc-pVDZ level. bElectronic energy only. cIncludes zero-point energy and thermal corrections at 298 K.
2 2
) 2 z j, 0 )
− y j, 0 −
H
1
2
(1)
where the sum is over all atoms. We denote mj as the mass (in amu) of the j-th atom, and x0 and xi are the x-coordinate values of the initial and ith IRC point, respectively, and so on, for y and z. The units of the IRC used in this work are Å amu1兾2. A significant advantage to using a potential expressed in mass-weighted coordinates is that the Hamiltonian matrix can be diagonalized to obtain the eigenvalues of the system without having to make assumptions about the reduced mass of the system, that is, whether it is independent of, or some function of, the IRC. The IRC calculations utilize the optimized transitionstate structure obtained in the previous run. The routing section is as follows: b3lyp/cc-pvdz irc(calcfc,forward,maxpoints=100,stepsize=5)
These keywords instruct the calculation to calculate the force constants of the input structure, allowing the calculation to jump off the saddle point in the “forward” direction, and to advance along the IRC (in this case) in steps of 0.05 bohr amu1/2. The molecular specification section should contain
www.JCE.DivCHED.org
•
the Cartesian coordinates of the optimized transition-state structure, which can be cut and pasted from the output file of the transition-state optimization run (the geometry following the header “Standard Orientation” nearest the end of the file). The geometry in this case is: N 0.081339 -0.570256 0.000000 C 0.081339 0.622850 0.000000 H -1.057413 0.254691 0.000000
This calculation should produce 71 points that advance to 3.44587 bohr amu1兾2, and an energy that corresponds to the near-equilibrium structure of HCN.4 An identical calculation is performed in which “reverse” replaces “forward” in the IRC argument. In this case, 82 points are generated that run from ᎑4.02563 bohr amu1兾2 (at the equilibrium structure of HNC) back to an IRC value of 0, which corresponds to the transition state.
Potential Energy Surface Using these two IRC arrays, the student should be able to assemble the IRC potential that spans a full revolution of the H atom: HNC → HNC‡ → HCN → HCN‡ → HNC. It is helpful to conceive of the isomerization process as the motion of the H atom from one end of the C⫺N moiety to the other under the influence of a double minimum poten-
Vol. 83 No. 1 January 2006
•
Journal of Chemical Education
71
In the Classroom TS 50
40
E / (kcal molⴚ1)
tial. This IRC potential is illustrated in Figure 2 in which the IRC origin is set to correspond to the metastable HNC molecule, and zero energy is assigned to the more stable HCN isomer. Students might find it enlightening to visualize the minimum energy trajectory that the H atom follows in isomerizing from HNC to HCN. This cyclic path is shown in Figure 3, in which the distance between the H atom and the center of mass of the C⫺N diatom, xcm, is plotted radially as a function of the angular variable, θ, defined as the angle between the line connecting the H atom to xcm and the line connecting xcm to the N atom. Thus in one-half cycle, θ spans from 0 (the reactant, HNC) to π radians (the product, HCN). The diamond symbols indicate the locations of the transition state. It is also instructive to transform the linear IRC potential shown in Figure 2 to its radial form. This graph is illustrated in Figure 4. This representation helps students visualize the energy as the H atom moves from the HNC species, at 0 radians, through the transition state, indicated by the triangles, and on to HCN at π radians. The means of constructing the radial plots is described in the Appendix.
30
20
10
0 -5
-4
-3
-2
-1
0
1
IRC / (Å amu
1/2
2
3
4
5
)
Figure 2. The IRC potential for the HNC–HCN isomerization obtained at the B3LYP/cc-pVDZ level of theory. The origin corresponds to the reactant, HNC. A step size of 0.04 bohr amu1/2 was used in these calculations using Gaussian 98.
The CH3NC–CH3CN Isomerization Now that students have explored the details of the HNC–HCN reaction, they might be inspired to apply this knowledge to the somewhat more complex, but experimentally more well-studied system, the methylisocyanide– methylcyanide isomerization,
In this case, it is the methyl group that undergoes the 1,2 shift. The basic methodology described above is followed in this study. Students will realize, in analogy with HNC or
HCN, that the three heavy atoms are in a linear arrangement (both molecules have C3v symmetry). However, for these molecules they should be advised to create molecule specification files in z-matrix format rather than in Cartesian coordinates. This approach is not only more convenient in setting up the molecule specification file, but it also will allow one to exert control of the methyl rotational orientation in the scan and in the transition-state structure calculations. However, for the equilibrium structures, the linear C⫺N⫺C (or C⫺C⫺N) structure requires the use of a “dummy atom” (rep-
π/2
π/2
CH3NC
CH3CN
TS TS
H
2Å 1Å π
C
xcm
0
N
π
0
TS
TS 3π/2
3π/2
Figure 3. Radial plot of the minimum energy path followed by the H atom in the isomerization process. The circles represent distances between the H atom and the center of mass, xcm, of the CN diatom of 1 Å and 2 Å. The diamonds labeled TS indicate the locations of the transition state.
Figure 4. Radial plot of the IRC potential illustrated in Figure 2. The circle represents an energy of 80 kcal/mol and the marks labeled TS indicate the transition state.
72
Journal of Chemical Education
•
Vol. 83 No. 1 January 2006
•
www.JCE.DivCHED.org
In the Classroom
resented as X) in the z-matrix to avoid computational problems (17). The appropriate molecule specification for CH3NC is X N 1 1.0 C 2 rcn C 2 rcc H 4 rch H 4 rch H 4 rch rcn 1.1 rcc 1.5 rch 1.1 a 109.
1 1 2 2 2
90. 90. a 1 a 1 a 1
with respect to the C⫺N⫺C plane. This scan reveals a maximum energy with the geometrical parameters rcn 1.19950 rcc 1.76545 rcha 1.09642 rchb 1.09384 a1 75.0 a2 122.91697 a3 99.36112 d 122.09221
3 180. 180. 60. -60.
A transition-state optimization is now performed on the updated geometry specification with the route section containing b3lyp/cc-pvdz opt=(ts,calcfc,z-matrix) freq.
This input structure is shown in Figure 5. The X atom (which does not enter into the energy calculation) is fixed at a 1 Å distance from the N atom, and the two C atoms are held at 90⬚ with respect to the N⫺X “bond.” The dihedral angles of the three H atoms, defined with respect to the C, N, and X species, are set to the values indicated in order to produce the three-fold rotational symmetry of the CH3 moiety. The three C⫺H bond lengths and H⫺C⫺N angles are made equivalent, as required by the symmetry of the molecule. The molecule specification for CH3CN is identical, except that the N and (first) C atoms are switched. The results of these calculations are summarized in Table 3. The calculated and experimental enthalpies of reaction, ᎑99.54 and ᎑99.2 (24) kJ兾mol, respectively, are in excellent agreement (see Table 5). The location and optimization of the transition state proceed as described above for HNC–HCN reaction. First one performs a relaxed scan of the potential energy surface on CH3CN in which the N⫺C⫺CH3 bond angle is brought from 165⬚ to 25⬚ in 10⬚ steps. A z-matrix containing a dummy atom is useful here because the job allows for the proper movement and rotational control of the CH3 group. The Route Section is identical to one used for the previous optimizations. The molecule specification for this scan is: X N 1 1.0 C 2 rcn 1 90. C 3 rcc 2 a1 1 180. H 4 rcha 3 a2 1 0.0 H 4 rchb 3 a3 1 d H 4 rchb 3 a3 1 -d rcn 1.1 rcc 1.7 rcha 1.0 rchb 1.0 a1 165. S 14 -10. a2 110. a3 110. d 120.
www.JCE.DivCHED.org
H6 C4
C3
Figure 5. The structure of the CH3NC molecule, and the identities of the atoms, corresponding to the z-matrix used in the optimization calculation. The “dummy atom” is represented by X1.
Table 3. Optimized Properties of CH3NC and CH3CN and Comparison with Experiment CH3NC
Property
CH3CN
Calca
Exp
Calca
Exp
rMe(N)C/Å
1.4211
1.419b
1.4615
1.4585h
rCN/Å
1.1786
1.171b
1.1611
1.1564h
1.1005
b
1.105
1.0905h
b
rCH/Å
7
c
᎑1
1.091
110.24
109.94h
269.1
d
263
383.9
364.7i
968.1
945d
110.04
109.6
938.3
915.4i
d
1043.9
1040.9i
d
1388.3
1390i
1459.4
d
1467
1443.3
1448.0i
2240.3
2166d
2371.7
2270.6i
3042.5
d
3049.4
2922.7i
d
1128.6 1433.4
1129 1429
2966
3124.0
3014
3130.9
3009.2i
E /hartree
᎑132.723599
---
᎑132.716509
---
Hf/hartree
᎑132.674040
---
᎑132.711952
---
᎑132.702027
---
᎑132.739502
---
0.0447581
---
0.0449976
---
v8 /cm e
•
N2
H5
∼a HC(N)C/deg ∼v c/cm᎑1 1 ∼v /cm᎑1 2 ∼v c/cm᎑1 3 ∼v /cm᎑1 4 ∼v c/cm᎑1 5 ∼v /cm᎑1 6 ∼v /cm᎑1
In this calculation we anticipate that the transition state will possess a symmetry in which two of the three C⫺H bonds are equivalent. We also express the H⫺C⫺N bond angles in terms of two different variables, a2 and a3. In addition, we have set the conformation of the CH3 group so that one of the C⫺H bonds forms a dihedral angle of 0⬚
X1
H7
f
G /hartree ZPEg/hartree a
This work. All calculations are at the B3LYP/cc-pVDZ level. bRef 18. c Doubly degenerate mode. dRef 19. eElectronic energy only. f Including zero-point energy and thermal corrections at 298 K. gZeropoint energy. hRef 20. iRef 21.
Vol. 83 No. 1 January 2006
•
Journal of Chemical Education
73
In the Classroom
C
N
C
Figure 6. The B3LYP/cc-pVDZ optimized transition state for the CH3NC–CH3CN isomerization. The properties of this structure are contained in Table 4.
Table 4. Optimized Properties of the CH3NC–CH3CN Transition State and Comparison with Literature Values Property
Calca
Literaturee
rCH/Å
1.200
1.178
rMeC/Å
1.787
1.740
rMeCN/Å
1.836
1.898
rCHa/Å
1.096
1.078
b
rCH /Å
1.094
1.076
aMeCN/deg
72.85
78.54
aMeNC/deg
68.49
72.74
aCMeN/deg
38.66
37.47
435.1i
458i
206.2
255
∼v /cm᎑1 1 ∼v /cm᎑1 2 ∼v /cm᎑1 3 ∼v /cm᎑1 4 ∼v /cm᎑1 5 b
613.2
677
968.8
1063
971.4
1083
᎑132.658959
---
H /hartree
᎑132.612336
---
Gc/hartree
᎑132.641362
---
0.0420351
---
E /hartree c
ZPEd/hartree
a This work. All calculations are at the B3LYP/cc-pVDZ level. bElectronic energy only. cIncluding zero-point energy and thermal corrections at 298 K. dZero-point energy. eRef 22.
The optimization and frequency run produces a result with the structure shown in Figure 6. This transition state is confirmed by the presence of a single imaginary frequency (i.e., 435.1i cm᎑1). Students should examine this transition-state structure in the context of the expectation provided by the Hammond postulate (16). It is interesting to note that the methyl C⫺H bond lying most directly over the C⫺N bond is in an eclipsed conformation, as constrained in the Molecular Specification. If the CH3 group is rotated by 180o, thus producing the staggered conformation, a second-order saddle point, not a transition state, is obtained in which there are two imaginary frequencies, corresponding to the rotation of the CH3 group, as well as its movement along the reaction coordinate. The energy of this structure lies about 8 kJ兾mol above that of the transition state. The properties of this transition state are summarized in Table 4 and compared with earlier calculations (23). 74
Journal of Chemical Education
•
The IRC potential can be obtained in the same fashion as described for the HNC–HCN reaction. In this case one can use a larger step size, 0.2 bohr amu1兾2, in view of the larger mass of the migrating CH3 group. At this point the student will have carried out reasonably complete calculations of the HNC–HCN and CH3NC– CH3CN isomerization reactions. With the information obtained it is possible to estimate the reaction enthalpies, equilibrium constants, and the rate constants for these processes. Table 5 lists the reaction enthalpies and equilibrium constants and their comparisons with experiment, where appropriate. To obtain the rate constants, we will first apply RRK theory, assuming that energy is randomly distributed among the vibrational degrees of freedom. Using the expression for the high-pressure limit (24), kRRK K = νi exp
− E0 kB T
(2)
where νi is the frequency (in Hz) of the critical oscillator (i.e., the magnitude of the imaginary frequency) and E0 is the difference in zero-point energies of the transition state and reactant. Students obtain the necessary information from the E0 values of the reactants and transition states and the frequency calculations of the transition states (see Tables 1–4).5 From eq 2 one obtains (for HNC) kRRK = 2.4 × 10᎑8 s᎑1 and 8.4 s᎑1 at 298 and 500 K, respectively. For CH3NC the results are 4.2 × 10᎑16 s᎑1 and 1.35 × 10᎑4 s᎑1 at 298 and 500 K, respectively. These results are listed in Table 5. One can also estimate the rate constant from transitionstate theory (TST) formulation (28) using kTST =
− ∆G ‡ ° kB T exp h RT
(3)
where kB and h are Boltzmann’s and Planck’s constants, respectively, ∆G‡⬚ is the change in the standard Gibbs free energy between the transition state and HNC species (including the respective zero-point energies). The G98W frequency calculations provide these values (see Tables 1, 2, and 3), which led to ∆G‡ = 116.7 kJ兾mol and kTST = 2.19 × 10᎑8 s᎑1 at 298 K and 11.3 s᎑1 at 500 K, in close agreement with the RRK results at the respective temperatures. Students can readily conclude that HNC is a kinetically stable species at ambient temperatures, although it is thermodynamically unstable with respect to the unimolecular isomerization to HCN. The H and G values at 500 K are obtained from runs using the equilibrium HNC and transition-state structures in which the keyword frequency = readisotope is used. An additional set of statements is included (following a blank line) after the Molecule Specification as follows: 500. 1.0 1.0 14 12 1
The first line corresponds to the temperature (in K), pressure (in atm), and the vibrational frequency scaling factor (taken here as 1.0),3 and the next three lines indicate the atomic masses numbers (in the same order as the molecule specification).
Vol. 83 No. 1 January 2006
•
www.JCE.DivCHED.org
In the Classroom Table 5. Summary of Thermodynamic and Kinetic Properties of the HNC–HCN and CH3NC–CH3CN Isomerization Reactions Property
HNC–HCN
CH3NC–CH3CN
298 K
500 K
298 K
500 K
∆Eo/(kJ/mol)
᎑61.38a
᎑62.13
᎑99.58
᎑99.58
∆Ho/(kJ/mol)
᎑62.13 ᎑58.91
᎑99.54 (᎑99.2)d ᎑98.41
᎑99.58
∆Go/(kJ/mol)
᎑61.38 (᎑61.9)b ᎑60.12
∆G‡/(kJ/mol)
116.7
114.5
159.3 (160)e
157.6
4.13
1.77 x 1017
Ko
3.49 x 1010
᎑16
1.44 x1010
᎑1
2.42 x 10 (1.1 x 10᎑8)c
8.4 (16.1)c
4.2 x 10 (2.7 x 10᎑16)f
1.35 x 10᎑4 (4.48 x 10᎑4)g
kTST/s᎑1
2.19 x 10᎑8
11.3
5.22 x 10᎑16
3.58 x 10᎑4
kRRK/s
᎑8
᎑97.24
a Energies obtained at the B3LYP/cc-pVDZ level. bRef 13. cData obtained from ref 23. dRef 24. eSee ref 26. This is reported Ea value. Using the relationship Ea = RT + ∆E‡o, one obtains 164 kJ/mol. fExtrapolated value based on A∞ and E∞ values of ref 26. gCalculated from A∞ and E∞values of ref 27.
As an additional component of this project, students can repeat the rate constant calculations for the isotopomer reaction, DNC–DCN. Although the electronic energy difference between the reactant and transition state is unaltered by the heavy hydrogen atom substitution, the rotational constants, vibrational frequencies and partition functions will, of course, differ. These values can be obtained from a frequency calculation in which the keyword frequency = readisotope is used as indicated above, and the appropriate temperature and atomic masses entered. For this reaction at 500 K, kTST and kRRK are 5.31 and 3.13 s᎑1, respectively, and compare with values of 11.3 and HCN reaction (see Table 5). These 8.4 s᎑1 for the HNC calculated rate constants are in excellent agreement with values obtained from calculations reported by Peri´c et al. (23). As in the case of the HNC–HCN reaction, we can estimate the forward rate constants using RRK or transition-state theory. However, for the CH3NC–CH3CN reaction there is a considerable information about the kinetics and thermodynamics from experimental studies. The calculated and experimental rate constants are displayed in Table 5. The experimental value at 500 K, which is more reliable, is larger by a factor of ca. 2 than kRRK, which is quite reasonable. Summary The isocyanide–cyanide isomerization reaction provides an extensive set of challenging and rewarding learning experiences for students. A variety of projects can designed that allow them to more fully utilize the fundamental tools of computational chemistry. The linkages to experimental studies in thermodynamics and kinetics nicely round out this project. Acknowledgment The author acknowledges with pleasure many helpful conversations with his colleagues, E. D. Glendening and B. R. Ramachandran. www.JCE.DivCHED.org
•
Notes 1. Input text is denoted in Courier font. 2. The value for C⬅C is used for C⬅N. 3. Frequencies reported here are unscaled. B3LYP frequencies are sometimes scaled by 0.96 in an attempt to account for anharmonic effects. 4. The number of points and the range of the IRC depend on the step size used. 5. Thus E0 = E + ZPE, where E and ZPE are obtained from Tables 1–4.
Literature Cited 1. Levine, I. N. Quantum Chemistry, 4th ed.; Prentice Hall: Englewood Cliffs, NJ, 1991; pp 520–525. 2. (a) Lee, C.; Yang, W.; Parr, R. G. Phy. Rev. B 1988, 37, 785– 789. (b) Mielich, B.; Savin, A.; Stoll, H.; Preuss, H. Chem. Phys. Lett. 1989, 157, 200. (c) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. 3. Woon, D. E.; Dunning, T. H., Jr. J. Chem. Phys. 1993, 98, 1358–1371. 4. Gaussian 98; Revision A.7; Frisch, M. J. et al.; Gaussian, Inc.: Pittsburgh PA, 1998. 5. Gaussian 98W; Revision A.9; Frisch. M. J. et al.; Gaussian, Inc.: Pittsburgh PA, 1998. 6. (a) Fukui, K. J. Phys. Chem. 1970, 74, 4161–4163. (b) Fukui, K. Acc. Chem. Res. 1981, 14, 363–368. 7. Atkins, P.; de Paula, J. Physical Chemistry, 7th ed.; W. H. Freeman and Company: New York, 2002; p 423. 8. Creswell, R. A.; Robiette, A. G. Mol. Phys. 1978, 36, 869– 876. 9. Milligan, D. E.; Jacox, M. E. J. Chem. Phys. 1967, 47, 278– 285. 10. Winnewisser, G.; Maki, A. G.; Johnson, D. R. J. Mol. Spec. 1971, 39, 149–158. 11. Herzberg, G. Molecular Spectra and Molecular Structure. III. Electronic Spectra and Electronic Structure of Polyatomic Molecules; D. van Nostrand Company, Inc.: Princeton, NJ, 1966; p 588.
Vol. 83 No. 1 January 2006
•
Journal of Chemical Education
75
In the Classroom 12. Ziurys, L. M.; Savage, C.; Brewster, M. A.; Apponi, A. J.; Pesch, T. C.; Wyckoff, S. Astrophys. J. 1999, 527, L67–L71. 13. Pau, C.-F.; Hehre, W. J. J. Phys. Chem. 1981, 86, 321–322. 14. Foresman, J. B.; Frisch, Æ. Exploring Chemistry with Electronic Structure Methods; Gaussian, Inc.: Pittsburgh, PA, 1993; pp 33–36. 15. van Mourik, T.; Harris, G. J.; Polyansky, O. L.; Tennyson, J.; Császár, A. G; Knowles, P. J. J. Chem. Phys. 2001, 115, 3706– 3718. 16. Loudon, G. M. Organic Chemistry, 4th ed.; Oxford University Press: New York, 2002; pp 147–149. 17. Foresman, J. B.; Frisch, Æ. Exploring Chemistry with Electronic Structure Methods; Gaussian, Inc.: Pittsburgh, PA, 1993; p 146 ff. 18. Halonen, L.; Mills, I. M. J. Mol. Spectrosc. 1978, 73, 494– 502. 19. Shimanouchi, T. Tables of Molecular Vibrational Frequencies; National Bureaux of Standards: Washington, DC, 1972; p 1. 20. LeGuennec, M.; Wladarczak, G.; Burie, J.; Demaison, J. J. Mol. Spectrosc. 1992, 154, 305–323. 21. Duncan, J. L.; McKean, D. C.; Tullini, F.; Nivellini, G. D.; Pena, J. P. J. Mol. Spectrosc. 1978, 69, 123–140. 22. Saxe, P.; Yamaguchi, Y.; Pulay, P.; Schaefer, H. F., III. J. Am. Chem. Soc. 1980, 102, 3718–3723. 23. Peri´c, M.; Mladenovi´c, M.; Peyerimhoff, S. D.; Buenker, R. J. Chem. Phys. 1984, 86, 85–103. 24. Baghal-Vayjooee, M. H.; Callister, J. L.; Prichard, H. O. Can. J. Chem. 1977, 55, 2634–1636. 25. Steinfeld, J. I.; Francisco, J. S.; Hase, W. L. Chemical Kinetics and Dynamics, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 1989; pp 340–341. 26. Schneider, F. W.; Rabinovitch J. Am. Chem. Soc. 1962, 84, 4215–4230.
76
Journal of Chemical Education
•
27. Collister, J. L.; Prichard, H. O. Can. J. Chem. 1977, 54, 2380– 2384. 28. Steinfeld, J. I.; Francisco, J. S.; Hase, W. L. Chemical Kinetics and Dynamics, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 1989; pp 300–301.
Appendix The x, y coordinates of the three atoms along the IRC are cut and pasted into a spreadsheet. For the ith step along the IRC, the distance from the H atom to the center of mass of the NC diatom, ri, is obtained from the following expression
ri =
x H, i
+
14..00 xN,i + 12.00 xC,i − 26.00
y H, i −
2
14.00 yN, i + 12.00 yC,i 26.00
2
1
2
where xH,i and yH,i denote the x, y coordinates of the H atom, and so forth. One then obtains an array, r(N ), where N corresponds to the total number of points in a full-cycle IRC, such as shown in Figure 2. Next the IRC values are converted to radians, θ, by multiplying by 2π兾N. Two new arrays are now constructed, Xi = ri sinθ and Yi = ri cosθ. The angular plot is rendered by plotting Yi versus Xi. To avoid the distortion that is caused by a default rectangular format, the y and x axes should be adjusted to have the same values per unit length. The angular plot of the IRC potential is created in the analogous fashion.
Vol. 83 No. 1 January 2006
•
www.JCE.DivCHED.org
