118
J . Phys. Chem. 1989, 93, 118-125
Intrinsic Reaction Coordinate Calculations for Very Flat Potential Energy Surfaces: Application to Singlet Si,H, Isomerization Shiro Koseki and Mark S. Gordon* Department of Chemistry, North Dakota State University, Fargo, North Dakota 58105 (Received: April 4, 1988)
Si2H2is an important species that may appear in the chemical vapor deposition of silicon. This paper reports an intrinsic reaction coordinate (IRC), or a minimum energy path (MEP), of the isomerization from silasilene to bridged disilyne obtained by using the local cubic and quadratic approximations. These new approximationsgenerate a correct IRC for this isomerization, while some conventional methods fail to predict a reasonable IRC because of the very flat potential energy surface. This reaction path bifurcates to two identical IRCs to reach bridged disilyne. The activation energy for this isomerization is predicted to be less than 3 kcal/mol. This energy barrier may disappear at high temperatures. This paper also compares the IRCs on the simple potential energy surfaces of the ammonia inversion and the HCN isomerization generated by using the local cubic and quadratic approximations with those obtained by some conventional methods.
I. Introduction Si2H2is one of the important species that may appear in the mechanism of the chemical vapor deposition (CVD) of silicon.' Snyder et aL2 reported theoretical equilibrium geometries for five prototypical conformations of singlet Si2HZ,namely, linear (D.+ 'E,'), trans-bent (C2h, 'Ag), twisted (C2,IA), and bridged (C,,, 'Al) disilynes, and planar silasilene (CZu,lAI: Si-SiH,). According to their generalized valence bond3 (GVB)/4-3 l G Zresults, linear disilyne and planar silasilene are the least and most stable structures, respectively. However, no vibrational frequencies were reported for these structures. More recently, Binkley4 and Luke et aL5 have reexamined singlet Si2H2using fourth-order perturbation theory (MP4SDTQ): including triple excitations, and the 6-31 ++G(d,p) basis at Hartree-Fock (HF) geometries. These studies demonstrated that bridged disilyne and planar silasilene are the true minima on the singlet potential energy surface (PES) of Si2H2and that bridged disilyne is about 10 kcal/mol lower than planar silasilene. Binkley has also reported that the triplet minima of Si2H2are 3B) and bent silasilene (Cs,3A").4 However, twisted disilyne (C,, 3A2) in a recent notes it was reported that planar silasilene (C,,, is lower in energy than the bent structure when zero-point vibrational energies and/or the effects of electron correlation are included. Furthermore, the unrestricted HF (UHF) wave function for the 3B state of twisted disilyne is unstable to internal pert ~ r b a t i o n ,and ~ trans-bent disilyne is predicted to be a true minimum on the triplet PES. As part of a continuing study of reactions that may be involved in the mechanism of silicon CVD, the goal of the present work is to determine the isomerization path from silasilene to bridged disilyne. This problem is very interesting because of the existence of bifurcation points,1° and it is very difficult to calculate the path because of a very flat PES. (1) (a) Coltin, M. E.; Kee, R. J.; Miller, J. A. J . Electrochem. Soc. 1984, 131, 425. (b) Steinwandel, J.; Hoeschele, J. Chem. Phys. Lett. 1985, 116, 25. (2) Snyder, L. C.; Wasserman, Z. R.; Moskowitz, J. W. In?. J. Quantum Chem. 1982, 21, 565. (3) Bobrowicz, F. W.; Goddard, W. A. In Modern Theoretical Chemistry; Schaefer 111, H. F., Ed.; Vol. 3, Chapter 4. (4) Binkley, J. S . J . Am. Chem. SOC.1984, 106, 603. (5) Luke, B. T.; Pople, J. A.; Krogh-Jespersen, M.-B.; Apeloig, Y.; Karni, M.; Chandrasekhar, J; Schleyer, P. von R. J . Am. Chem. SOC.1986,108,270. (6) (a) Krishnan, R.; Pople, J A. Int. J. Quantum Chem. 1978, 14, 91. (b) Krishnan, R.; Frisch, M. J.; Pople, J. A. J . Chem. Phys. 1980, 72, 4244. (7) (a) Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta 1973,28, 213. (b) Francl, M. M; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. J. Chem. Phys. 1982, 77, 3654. (c) Schleyer, P. v. R., in press. See ref 4. (d) Gordon, M. S. Chem. Phys. Lett. 1980, 76, 163. (8) Koseki, S.; Gordon, M. S . J . Phys. Chem. 1988, 92, 364. (9) (a) Seeger, R.; Pople, J. A. J . Chem. Phys. 1976,65,265. (b) Seeger, R.; Pople, J. A. J. Chem. Phys. 1977, 66, 3045. (10) (a) Hoffman, D. K.; Nord, R. S.; Ruedenberg, K. Theor. Chim. Acta 1986,69, 265. (b) Valtazanos, P.; Ruedenberg, K. Theor. Chim. Acta 1986, 69, 281
0022-3654/89/2093-0118$01.50/0
The concept of an intrinsic reaction coordinate (IRC) or minimum energy path (MEP) on a PES has been proposed by Fukui and others," and methods of generating IRCs have been developed by several groups.12 In addition to providing a useful perspective of reaction paths, the IRC is an important starting point for the prediction of reaction rate constants. However, it will be demonstrated below that conventional methods cannot generate the correct IRC for the isomerization from silasilene to bridged disilyne because the PES is too flat. Recently, Page and McIverI3 have proposed a new method for tracing an IRC. The present paper investigates the accuracy of this new method and compares it with some conventional methods. The PESs for inversion of ammonia and isomerization of H C N are employed as simple examples. Then, the new method is applied to the isomerization of singlet Si2H2. 11. Theory and Computational Methods
The local cubic approximation (LCA) and the local quadratic approximation (LQA) proposed by Page and McIver are briefly reviewed here.13 The LCA determines the first point xlto which a reaction system moves from a transition state %. So that the required curvature K of a reaction coordinate usn can be obtained, derivatives of the force constant matrix at %, F(') = dF/ds, have to be calculated, because the gradient vector is zero. Force constant matrices are calculated analytically by using GAUSSlAN82'4 at structures displaced linearly forward and backward along the vibrational mode Q, corresponding to the imaginary frequency, and these are differentiated numerically to produce the matrix F(l). Then, the curvature K at xo is given by
K = (2Q:FQJ
- F)-l(F(I)Q, - Q,+F(l)Q,Qs) K
= IKI
(1)
(2)
where I is a unit matrix. By using this value, the first point xI along an IRC is predicted to be x1 = x,,
+ Q,6s + K6s2/2
(3)
Note that common methods for obtaining x1 neglect the last term in eq 3. The "stride" 6s is a small, arbitrarily determined step. (11) (a) Truhlar, D. G.; Kupperman, A. J. Am. Chem. SOC.1971, 93, 1840. (b) Fukui, K. J . Phys. Chem. 1970, 74,4161. (c) Fukui, K. Acc. Chem. Res. 1981, 14, 363. (d) Fukui, K. Pure Appl. Chem. 1982, 54, 1825. (12) (a) Ishida, K.; Morokuma, K.; Komornicki, A. J. Chem. Phys. 1977, 66,2153. (b) Schmidt, M. W.; Gordon, M. S.;Dupuis, M. J . Am. Chem. SOC. 1985, 107, 2585. (13) (a) Page, M.; McIver, J. W. J . Chem. Phys., in press. (b) Doubleday, C.; McIver, J. W.; Page, M. submitted for pubiication in J . Phys. Chem. (14) Binkley, J. S.; Frisch, M.; Ragavachari, K.; DeFrees, D. J.; Krishnan, R.; Schlegel, H. B.; Whiteside, R.; Fluder, E.; Pople, S. R.; Pople, J. A. GAUSSIANQ program, Carnegie-Mellon University, Pittsburgh, 1983 (un-
published).
0 1989 American Chemical Society
Very Flat Potential Energy Surfaces
The Journal of Physical Chemistry, Vol. 93, No. 1, 1989 119
Once the first point on the IRC away from the transition state has been determined, successive points on the IRC are obtained by using the LQA: xn+I
= xn
+ Una(t)Un+gn
(n > 1)
(4)
In eq 4, g , is the gradient vector at x,, U , is a unitary matrix constructed from eigenvectors of the force constant matrix F, at x,, and a ( t ) is a diagonal matrix whose elements are given by
ai = (e-hf- l)/f;. = Cfik-l(-t)k/k! k>O
(5)
H B
I I I
f;.is an eigenvalue of the force constant matrix F,. The magnitude of the second term in eq 4 is taken to be the “stride” 6s in this paper. The parameter “t” is determined to obtain an appropriate length of the stride. It is difficult to keep 6s constant, but the deviation is forced to be less than 1% in nearly the entire region of the IRC in this study. With regard to the numerical differentiation of the force constant matrix at xo, it appears that the best criterion is an energy lowering of 106-10-7 hartree at structures displaced linearly from the transition state. Full implementation of the LQA requires the calculation of a force constant matrix at each point on an IRC. Because this is very time-consuming, we have followed the suggestion of Page and c o - ~ o r k e r s : ’A~ ~force constant matrix obtained at a point on an IRC is employed to predict the next several points along the IRC on the assumption of a very slowly changing force constant matrix. For example, the total energy is computed at 60 points on the ammonia inversion IRC, while force constant matrices are calculated at only 10 selected points (described in section 1II.A in detail). In this paper, three “conventional methods” are employed to compare them with the LQA. The first method is the so-called “linear Euler method”-the original approach proposed by Fukui and others.” In this scheme, a new point X,+, along an IRC is predicted by following a gradient vector g , at a current point x,: X,+I
Vn
= x,
+ v,6s
= -gn/lgnl
(6) (7)
Ishida et al. have proposed a stabilization technique.12* This method predicts a new point along an IRC by searching along the bisector vector defined by dn+l
= (vn+l - vn)/Ivn+1 - VnI
(8)
v, is calculated at the nth point x, along an IRC, while v,+~ is obtained at the next point predicted by the linear method. A modified version of this stabilization technique developed in this laboratory is the second conventional method employed in this s t ~ d y . ’ ~ This ~ , ’ ~method is called the “Euler stabilization” (ES2) in ref 16. The third method is the “quadratic” method,16 which generates a new point x,+~ by using X,+I
= x, Vn(l)
+ v,6s + Vn(’)6S2/2 = (v, - v,,)/6s
(9)
(10)
This is called the first-order “fixed-step-size Adams predictor method” (FAP1) in ref 16. The IRCs generated by the methods considered in this paper will be referred to as “LIN” IRC (obtained by using the linear method), “STB” IRC (linear method with stabilization), “QUAD” IRC (quadratic method), and “LQA” IRC (local cubic and quadratic approximations), respectively. The stride 6s has been chosen to be 0.01 for the ammonia inversion and the Si2H2 I
(15) (a) Dupuis, M.; Spangler, D.; Wendoloski, J. J. Program QGOI, National Resource for Computation in Chemistry Software Catalog: Lawrence Berkeley Laboratory, Berkeley, CA, 1980: p 60. (b) Schmidt, M. W.; Boatz, J. A.; Baldridge, K. K.; Koseki, S.; Gordon, M. S.; Elbert, S. T.; Lam, B. QCPE 1987, 7, 115. (16) Garrett, B. C.; Redmon, M. J.; Steckler, R.; Truhlar, D. G.; Baldridge, K. K.; Bartol, D.; Schmidt, M. W.; Gordon, M. S. J . Phys. Chem. 1988.92, 1476.
I I
c3
Figure 1. RHF/6-31G geometrical structures of the transition state (A, D3,J and the equilibrium (B, C3”)in the ammonia inversion. Bond lengths and bond angles are in angstroms and degrees.
isomerization and 0.05 for the H C N isomerization, unless otherwise nokd. We have employed the restricted HartreeFock (RHF) method and the 6-31G basis (RHF/6-31G) to calculate IRCs of the ammonia inversion and the H C N isomerization. The reaction path for the isomerization from silasilene to bridged disilyne, which is our main interest, is investigated on the R H F / 6 - 3 1 G ( d , ~ ) ~ ~ v ~ PES. Boys localized molecular orbitals (LMOs) with the 6-31G(d,p) basis set are used to interpret the electronic structures at the equilibrium geometries of singlet Si2H2. Bent bond lengths of the Si-Si bond in bridged disilyne were determined by using the maximum electron density (MED) pathI7 and the 6-3 1G(2d,p) basis set.l8 The relative energies of the stationary points on this surface are estimated by using two kinds of higher level computational methods. One of these is Maller-Plesset nth-order perturbation calculations6 ( n = 2,4) based on the RHF/6-31G(d,p) wave function. In this calculation, core electrons are frozen, and all single, double, triple, and quadruple excitations from a ground electronic configuration are included (denoted RMP4SDTQ/6-31G(d,p)). The second approach is full optimized reaction space (FORS) MCSCFI9 using the 6-3 lG(d,p) basis set. The active space for the FORS MCSCF calculation includes all valence orbitals and electrons. 111. Results and Discussion
A . Inversion of Ammonia. It is well-known that double-{ and split valence basis sets incorrectly predict ammonia to be planar, or nearly so.2o Because of this, the RHF/6-31G ammonia inversion barrier is only 0.37 kcal/mol, with a predicted equilibrium H N H angle of 116.1’ (Figure 1). This gives rise to a very flat PES, where it is very difficult to generate a reasonable IRC. The frequency of the out-of-plane mode (a;) in the planar structure (D3&is 4241’ crn-’. Since this corresponds to a pure out-of-plane vibration, linear displacements lengthen the N-H bonds. Therefore, a relatively large displacement makes the total energy higher than that of the planar structure. Even if the total energy does not go up at a displaced structure, the force that shortens the N-H bonds is very strong at such a displaced (17) Bader, R. F. W.; Anderson, S. G.; Duke, A. J. J . Am. Chem. Soc. 1979, 101, 1389. (18) The 6-31G(d,p) basis set seems to give reasonable geometries and relative energies. However, it cannot provide a reasonable density map for bridged disilyne. Therefore, the 6-31G(d,p) is employed to obtain the length of the bent bond of bridged disilyne. (19) (a) Ruedenberg, K.; Schmidt, M. W.; Gilbert, M. M.; Elbert, S. T. Chem. Phys. 1982,71,41. (b) Ruedenberg, K.; Schmidt, M. W.; Gilbert, M. M. Chem. Phys. 1982, 71, 51. (c) Ruedenberg, K.; Schmidt, M. W.; Gilbert, M. M.; Elbert, S. T. Chem. Phys. 1982, 71, 6 5 . (20) Rauk, A.; Allen, L. C.; Clementi, E. J . Chem. Phys. 1970,52,4133.
Koseki and Gordon
120 The Journal of Physical Chemistry, Vol. 93, No. 1, 1989 1
I
I 0
0.5
0.5
0
s Cbohr*sqrt(a.rn.u.)l
Figure 2. N-H bond length along the LIN, QUAD, and LQA IRCs of the ammonia inversion: (a) LQA (6s = 0.01), (b) LIN (6s = 0.01), (c) LIN (6s = 0.005), (d) QUAD (6s = 0.01), and (e) QUAD (6s = 0.005).
See text.
1
s Cbohr*sqrt(a.rn.u.)l Figure 4. Potential energy surfaces (PES) along the LIN, STB, QUAD, and LQA IRCs of the ammonia inversion. The labels are the same as in Figures 2 and 3.
a
i
\ ; ,
,
,
,
,
,
,
,
,
,
, a
0.98
0
0
0.5
s Cbahr*sqrt(a.rn.u.)l
Figure 3. N-H bond length along the STB and LQA IRCs of the ammonia inversion: (a) LQA (6s = 0.01), (f) STB (6s = 0.01), and (g) STB (6s = 0.005). Also see Figure 2.
structure. The conventional methods start to generate an IRC at this structure, whose gradient vector has a relatively large contribution from the symmetric N-H stretching vibration. Accordingly, the conventional methods (LIN, STB, and QUAD) give N-H bonds at the next point that are too long. As a result, the N-H bonds oscillate along the entire IRCs obtained by these conventional methods, as shown in Figures 2 and 3. Because of the problem described above, there is an N H bond oscillation with a 0.003-Awidth in the LIN IRC. This oscillation makes the corresponding PES oscillate, and the total energies along this IRC often become higher than that of the planar structure, as demonstrated in Figure 4. A small stride (6s = 0.005) makes the PES smooth, but the PES goes down very slowly and the N-H bonds still apparently oscillate (see Figure 2). An infinitesimal stride might generate a reasonable IRC, but this is impractical. Strange results appear on the QUAD IRCs. The QUAD method is a higher level approximation than LIN, so that QUAD should generate a better IRC. However, the oscillations of the PESs and the N-H bonds along the QUAD IRCs are larger than those along the LIN IRCs, as shown in Figures 2 and 4. The second-order correction works badly, though the QUAD IRC obtained by using a smaller stride (6s = 0.005) improves in the region of s > 0.5. The X-N-H angles (X is a dummy atom on the C,axis.), depicted in Figure 5, decrease very slowly along the QUAD IRCs as well as the LIN IRCs, because of the N-H bond oscillation. This shows that the coupling between the X-N-H bending (reaction coordinate) and N-H stretching modes must be considered explicitly. The stabilization technique provides smooth P E S , as illustrated in Figure 4. However, the N-H bonds still oscillate slightly along
0.5
s Cbohr*sqrt(a.rn.u.)l
Figure 5. X-N-H angles along the LIN, STB, QUAD, and LQA IRCs of the ammonia inversion. "X" is a dummy atom as shown in Figure 1. The labels are the same as in Figures 2 and 3.
I
v, I
4000
E,
Y
A
c 3 U
e 2000
LL
0
0.5
s Cbohr*sqrt(a.rn.u.)l Figure 6. Vibrational frequencies calculated by the diagonalization of the projected force constant matrix along the LIN, STB, QUAD, and LQA IRCs of the ammonia inversion. X shows the point where the vibrational analysis is carried out. Every IRC gives almost pame values for e modes, (YJ and (v4), while an al mode has four different curves, (via) LQA (6s = 0.01), ( ~ l b ) LIN (6s = O.OOS), (vic) QUAD (6'' = 0.005),and ( u l d ) STB (6s = 0.005). u2 corresponds to the reaction mode.
the STB IRCs (Figure 3). The width of this oscillation is about 0.001 A, which can be viewed as acceptable, but the X-N-H angle does not change smoothly in the region 0 < s < 0.1, as shown in Figure 5. When a small stride (6s = 0.005) is employed, the N-H bond oscillation almost disappears after several steps and the
The Journal of Physical Chemistry, Vol. 93, No. I, 1989 121
Very Flat Potential Energy Surfaces 1.20,
1.1ej
1
/
- I
HNC
t-
__*
HCN
I
\\
0 l.16/
1.144. -3
I
I
-2
.
I
-1
1
r
0
1
I
1
I
. 2
I
r
3
.
20004,. -3 -2
. -1
X-N-H angle changes reasonably. Thus, this stabilization technique seems to be very effective, if an appropriate stride is employed. Figure 6 plots the vibrational frequencies obtained by the diagonalization of the projected force constant matrices2' at several points along the IRCs. The a, (N-H symmetric stretching) mode has a bump near the planar structure along the STB IRC (6s = 0.005). This bump is much smaller than those along the LIN and QUAD IRCs (6s = 0.005). This a, mode has the same symmetry as the reaction coordinate, so that the bump appears due to the neglect of the coupling between the reaction coordinate and this mode. This deficiency will have a serious impact on calculated free energy surfaces (FES)13b322and VTST (variational transition-state theory) rate constant^,^^^^' because the number of vibrational states and the density of vibrational states along the IRCs depend very strongly on the vibrational frequencies and the mode-mode interactions. The first displacement calculated by using the LCA is a mixture a,) and N-H symmetric of the reaction (out-of-plane, a? a l ) modes (ratio = 1:0.0024). Therefore, the stretching (a,' first point xIgiven by eq 3 has longer N-H bonds than those of the planar structure by less than A. The subsequent points are generated by using eq 4, and no N-H bond oscillation appears as shown in Figure 3, even though the stride is large (6s = 0.01). The total energy and the X-N-H angle also change very smoothly (Figures 4 and 5 ) , even though force constant matrices are obtained only at selected points, instead of at every point. As a result, there is no bump in the curves of the vibrational frequencies (Figure 6), and reliable reaction dynamics are expected to be obtained by using this IRC. The LIN and STB IRCs starting at the point predicted by the LCA method were also generated in this study. The N-H bonds still oscillate along these IRCs, though the oscillation width becomes smaller. Therefore, the coupling between the N-H bond stretching and reaction modes is very important along the entire path. B . Isomerization HNC HCN. Based on the discussion in the previous section, the stabilization technique appears to be potentially very helpful for generating a qualitatively reasonable IRC (as opposed to one that would be used to calculate a FES and/or VTST rate constant). However, as noted previously, one
-
-
-
(21) (a) Miller, W. H.; Handy, N. C.; Adam, J. E. J . Chem. Phys. 1980, 72, 99. (b) Gray, S. K.; Miller, W. H.; Yarnaguchi, Y.; Schaefer 111, H. F. J . Am. Chem. Soc. 1981, 103, 1900. (22) (a) Garret, B. C.; Truhlar, D. G. J . Am. Chem. Soc. 1979,101,5207. (b) Garret, B. C.;Truhlar, D. G. J . Phys. Chem. 1979, 83, 1079. (c) Doubleday, C.; McIver, J.; Page, M.; Zielinski, T. J . Am. Chem. SOC.1985,107, 5800. (d) Brown, F. B.; Tucker, S. C.; Truhlar, D. G. J . Chem. Phys. 1985, 83, 4451. (e) Garret, B. C.; Truhlar, D. G.; Bowman, J. M.; Wagner, A. F. J . Phys. Chem. 1986, 90, 4305.
1
2
3
s Cbohr*sqrt(a.m.u.)l
s Cbohr*sqrt(a.m.u.)l Figure 7. C-N bond length along the LIN, STB, and LQA IRCs of the HCN isomerization. A large negative value of s leads to HNC, while HCN corresponds to a large positive s. The stride 6s is set to 0.05. (a) LIN and LQA and (b) STB.
. . 0
Figure 8. Vibrational frequencies calculated by the diagonalization of
the projected force constant matrix along the LQA IRC of the HCN isomerization. The stride 6s is set to 0.05. X shows the point where the vibrational analysis is carried out.
-3
-2
-1
0
1
2
3
s Cbohr*sqrt(o.m.u)l
Figure 9. (a) Potential energy surface, (b) potential energy surface corrected by zero-point energy, and (e-h) (Helmholtz) free energy surfaces at several temperatures along the LQA IRC of the HCN isomerization. The temperatures are set to (c) 100,200, 298.15,400, and 500, (d) 1000, (e) 2000, (f) 3000, (8) 4000, and (h) 5000 K.
must be careful not to use too small a value of the stride with this met hod. The PES for the isomerization, HNC HCN, is not flat, but typical and popular. It is easy to generate a reasonable IRC. The PESs along LIN, STB, and LQA IRCs from H N C to HCN are very smooth and very similar to each other, and seem to be very reasonable. However, as illustrated in Figure 7, the stabilization technique forces the C-N bond to change rapidly. The other internal coordinates change reasonably along this IRC, though a small deviation of the STB IRC from the LIN and LQA IRCs can be found in the entire region. As discussed in ref 16, the stride (6s = 0.05) is too small to employ the stabilization technique because of very tiny changing of the direction of the gradient vectors. The LIN IRC is almost the same as the LQA IRC, as illustrated in Figure 7. The reaction coordinate is coupled with the N-H/C-H stretching mode, but the results show that this coupling is not important. It is worthwhile to consider the vibrational frequencies and (Helmholtz) FESs of this isomerization. Figures 8 and 9 plot vibrational frequencies and FESs at several temperatures along the LQA IRC, respectively. The lower frequency, corresponding to the C-N stretching mode, is almost unchanged along the entire IRC. The higher frequency mode is the N-H stretch at the beginning of the isomerization (s < 0). It decreases as the isomerization progresss, and it is replaced by the C-H stretch around the transition state. Figure 9 shows very interesting results. The higher the temperature becomes, the lower the energy barrier
-
Koseki and Gordon
122 The Journal of Physical Chemistry, Vol. 93, No. 1, 1989
B t.C
si
I
-
I
I
I
I
si
I
t
E
k
P l a n a r silaailene
Bridged diailyne
Figure 10. Schematicof the electronic structures of bridged disilyne and planar silasilene, according to localized molecular orbitals. The two pairs of electrons among the silicon and hydrogen atoms in bridged disilyne actually form 7r bonds between the two silicon atoms (see text). a* in planar silasilene shows a pair of a electrons.
is. Indeed, H N C might be bent above 3000 K, while H C N appears always to be linear. C. Isomerizationfrom Planar Silasilene to Bridged Disilyne.
H
Initial Considerations. The L M O pictures of the electronic structures of bridged disilyne and planar silasilene are shown schematically in Figure 10. The two lone pair LMOs of bridged disilyne are almost the same as those of the Si2 molecule (Figure 1la). The Si-Si L M O is largely expanded in the region opposite to the hydrogens (Figure l l b ) . The Si-H-Si LMOs have a maximum density on the hydrogen atoms and look like threecentered two-electron bonds (Figure 1 I C The Si-Si bond of bridged disilyne is relatively short (2.185 ); indeed, it is similar to that in d i ~ i l e n e . ~However, ~ the length obtained by tracing the MED path using the 6-31G(2d,p) basis set is 2.246 A.” As indicated by the dashed line in Figure 12a, this Si-Si bond is bent outward in a direction away from the hydrogens. Figure 12b illustrates the total electron density, plotted in a planar containing both silicons and one hydrogen. It is clear from this picture that the two hydrogen atoms simply attach to the two Si-Si ?r bonds, rather than forming the ordinary three-centered two-electron bonds that one might expect. That there is no H-H bond is illustrated in Figure 12c, in which the total density is plotted in a plane perpendicular to and bisecting the Si-Si axis. While there appears to be no MED path between the Si and H atoms in bridged disilyne, for simplicity, these are referred to as Si-H-Si bridged bonds in subsequent discussion. The H-X-H angle (X is the midpoint between the silicon atoms) is about 102O. This is a little larger than the 90° angle between the 7r orbitals, because of the repulsive interaction between the two hydrogen atoms. On the basis of the relative s:p:d Mulliken populations of 1:0.23:0, the lone-pair L M O centered on the terminal silicon of planar silasilene is virtually an atomic 3s orbital. The contribution of this atom to the Si-Si u bond is nearly pure 3p, (Figure lo), while the 3p, orbitals of both silicon atoms form the ?r bond. The 3p, orbital of the terminal silicon is vacant. So, the terminal silicon has virtually no hybridization. On the other hand, the central silicon seems to be sp2 hybridized, because s:p:d = 1:1.01:0.06 at the central silicon in the Si-Si u bond and the ratio in the Si-H LMOs is s:p:d = 1:1.62:0.15.24 The isomerization could occur via a C,reaction path in which the mirror plane perpendicular to the molecular plane of planar silasilene is retained. Although this least motion path is symmetry allowed, with 12 a’ and 3 a” doubly occupied orbitals, the structure at the energy maximum along this path has a second imaginary
1
(23) Schmidt, M. W.;Truong, P.N.; Gordon, M. S . J . Am. Chem. SOC. 1987, 109, 5217. (24) Mulliken population analysis at the HF/6-3 IG(d,p) level gives an s:p:d ratio of 1:1.99:0.20 for sp’ orbitals of the SiH, molecule, 1:2.05:0.20 for sp’ orbitals of the Si2H6molecule, and 1:I .43:0.13 for sp2 orbitals of the D2* structure of the Si2H4molecule.
C -
Figure 11. (a) Lone-pair localized molecular orbital (LMO), (b) Si-Si LMO, and (c) Si-H-Si LMO in bridged disilyne.
frequency (a”), in addition to the mode that leads to the two isomers. The electronic structure of planar silasilene suggests that, at the beginning of the isomerization from planar silasilene, one of
The Journal of Physical Chemistry, Vol. 93, No. 1, 1989 123
Very Flat Potential Energy Surfaces
a
I
I
A -
I
C
n 1.686
H -H= 1.9 9 7 0, TS’), especially near the VRI point. The vibrational frequency corresponding to the a” mode goes to zero at the VRI point. Therefore, vibrationally excited states have to be considered explicitly near the VRI point even at a low temperature to obtain a reliable FES. Though the FESs around the VRI point may not be reliable, it can be concluded that the transition state shifts to the reactant side as temperature increases, and silasilene might not be observable at high temperatures.2s Reexamination of Relative Energies of Stationary Structures. Table I lists the relative energies of silasilene, TS, TS’, and bridged disilyne at several computational levels. The energy barrier of the isomerization from silasilene to bridged disilyne is 9.8 (RHF/6-31G(d,p)), 1.6 (RMP2/6-31G(d,p)), or 2.7 (RMP4SDTQ/6-3lG(d,p)) kcal/mol, from use of RHF/6-31G(d,p) geometries. The zero-point energy (ZPE) correction obtained at the RHF/6-31G(d,p) level reduces the barrier by 2.7 kcal/mol, so that there might be no barrier at the RMP2+ZPE and RMP4SDTQ+ZPE levels. Accordingly, silasilene easily isomerizes to bridged disilyne, as Ho et al. have The barrier for the reverse reaction is 14.0 (RHF), 13.9 (RMP2), or 14.0 (RMP4SDTQ) kcal/mol (the ZPE correction also reduces this barrier by 2.7 kcal/mol). FORS MCSCF/6-3 lG(d,p) calculations were also carried out at the RHF/6-3 lG(d,p) geometries. The energy difference between bridged disilyne and silasilene is only l .7 kcal/mol. The energy barrier for the isomerization from silasilene to bridged disilyne is 6.0 kcal/mol (7.7 kcal/mol) for the reverse reaction). (25) Ho, P.; Coltrin, M. E.; Binkley, J. S.; Melius, C. F. J . Phys. Chem. 1986, 90, 3399.
This barrier might be reduced by the ZPE correction, as described in the previous paragraph, yielding an activation energy of about 3 kcal/mol, assuming a similar ZPE correction as at the R H F level. The energetic ordering calculated by the MCSCF method is the same as that obtained by the R M P approaches. Thus, it can be concluded that the activation energy for the isomerization of silasilene is no more than 3 kcal/mol. The MCSCF energy difference between silasilene and bridged disilyne is much smaller than the R M P ones. The MCSCF calculations consider correlation effects only in the valence space, while the RMP approaches include electron correlation between the valence and external spaces as well as within the valence space. Even though the amount of configuration mixing increases slightly from bridged disilyne to silasilene (Table 11), the mixing is relatively small along the points of the surface considered here. It is therefore likely that the R M P results are more reliable.
IV. Conclusion The LCA and LQA proposed by Page and McIver can generate reasonable IRCs for the very flat RHF/6-31G PES of the ammonia inversion, the typical RHF/6-31G PES of the H C N isomerization, and the flat RHF/6-31G(d,p) PES of the singlet Si2H2isomerization. These IRCs can be employed to calculate dynamical properties for these reactions. According to the results of the FESs for the H C N isomerization, H N C might be bent at high temperatures. A valley-ridge inflection (VRI) point, where the reaction path bifurcates to two identical IRCs, is found along the IRC for the isomerization from silasilene to bridged disilyne. This point is very close to the transition state TS‘ for the pseudorotation of the two hydrogen atoms in bridged disilyne, as well as to the transition state TS for this isomerization. The energy barrier for this isomerization is about 6 (MCSCF) or 3 kcal/mol (RMP4SDTQ), but it might go to zero at high temperatures. Therefore, only bridged disilyne might be observable in the silicon CVD. Acknowledgment. This work was supported by grants from the National Science Foundation (CHE86-40771) and the Air Force Office of Scientific Research (AFOSR 87-0049). We are indebted to Dr. Michael Page, Dr. Michael W. Schmidt, and Jerry A. Boatz for several helpful discussions. All calculations were performed on an IBM 3081D at the North Dakota State University Computer Center, a VAX 11/750 minicomputer (purchased with the aid of AFOSR Grant 84-0428), a microVAX I1 minicomputer (purchased with the aid of N S F Grant CHE-8511697), and a VAX 8530 minicomputer (purchased with the aid of AFOSR Grant 86-0237). Computer time was also made available on the CRAY XMP48 at the San Diego Supercomputer Center through the N S F grant. Registry No. Si,H,, 87970-40-9.
