3-1 Theoretical results for OH+H
2H+H
2O reaction
3-1.1 Theoretical calculation for MEP
There are several theoretical groups devoted on the research of OH+H2H+H2O reaction; Smith and Zenner1 use two different semi-empirical potential energy surface; LEPS and BEBO and they neglect tunneling with barrier same as experiment activation energy, Schatz used large scale POL-CI3 wave function with basis 3s3p1d/3s1p and 4s3p2d/3s2p and Truong4 used the PMP4/6-311++G(2d f ,2pd)//QCISD/6-311+G(d,p) method. And in this report, Gaussian 09 is used to calculate the reaction path with CCSD/aug-cc-pvdz, CCSD(T)/ aug-cc-pvdz//CCSD/
aug-cc-pvdz and CCSD(T)/aug-cc-pvtz//CCSD/aug-cc-pvdz. All the results are shown in table 3.1. From table 3.1, we could see that the energy correction by CCSD(T) is necessary, which match pretty well with Schatz’s larger basis and Troung, no matter aug-cc-pvdz or aug-cc-pvtz, but seems shown difference with Zener and the activation energy (Ea) at 300K, which Zener took the activation energy as barrier height in the semi-empirical potential energy surface. It could be contributed to the tunneling effect and temperature effect. Schatz4 and Truong5 have calculated the rate constant with the effect of tunneling, then they get pretty good results with experiments, that is to say, the results of CCSD(T) energy correction is ok. And the comparison for MEP of my work between before and after energy correction by CCSD(T) shown in figure 3.1 indicates that the saddle point does not change to much compared to the one without energy correction and almost no difference between the two different basis in energy correction, CCSD(T)/aug-cc-pvdz//CCSD/aug-cc-pvdz and CCSD(T)/aug-cc-pvtz//
CCSD /aug-cc-pvdz. So we use the MEP calculated by CCSD(T)/ aug-cc-pvdz //
CCSD aug-cc-pvdz.
3-1.2 Geometry of transition state
Figure 3.2 displays the geometries of transition state with CCSD/aug-cc-pvdz method, which predicted that the O-H length of HOH’H’’ is 0.976 A, which almost same as the length of OH radical. This result match the concept that several studies bring out, which OH can be seen as a spectator group during the reaction. The geometry is almost same as Schatz2 and Truong4’s result, but different from Zenner1, which is no longer planar.
3-1.3 Normal modes along reaction path
Table 3.2 presents the frequencies along reaction path calculated by G09 with CCSD/aug-cc-pvdz. In general, there should be 7 zero or very small frequencies and 3N-7 non-zero frequencies, here is 5. From table 3.2, you could see that it’s not the case. For points near transition state (s = 0) are still ok, but points away from s = 0 show strong defects. Since there are two general problems exist in the G09 calculation for points on reaction coordinate:
1. Linear transformation of Cartesian coordinate to internal coordinate:
For non local maximum or minimum point, the transformation will be no longer linear, since the first order gradient will also contribute to it like Truhlar metioned6. Bellow is the energy based on Cartesian (equation (3.1)) and internal coordinates (equation (3.2))
+L
∆
∆ +
∆ +
=
∑ ∑
Nj i
j i ij i
N
i Gi R F R R
V V
3
, 3
0 2
1 (3.1)
And
,in which ∆Ri and ∆qi are, respectively, the Cartesian and internal coordinate.
The internal coordinate can be expressed in terms of power serious of the Cartesian coordinate as bellow:
∑
From equation (3.1) to (3.3) yields the expressions for Cartesian gradient and force const matrix in terms of their internal coordinate counterpart as
g
In above equation, it could be seen easily that it’s not linear transformation, which is the only first term of equation (3.5)
2. The process of projecting out the unit component along reaction path, rotation and translation is not conducted.
Because of problem 1, we choose to use Molpro to calculate the frequencies by the same method and basis, CCSD/aug-cc-pvdz, with the geometries from G09. Since Molpro directly diagonalize the transformation matrix in Cartesian coordinate without changing to internal coordinate. The results of Molpro are indicated in table 3.3. It seems better but there is still some defects in it, because of problem 2 mentioned above. This program for this is not that easy to construct. In order to find the 3N-7 nonzero frequencies, we’ve done the following steps:
1. Transform the first order gradient in Cartesian coordinate from Molpro to normal
mode coordinate with the transformation matrix from Molpro
2. Do mode-scanning for each coordinate with the transformation matrix obtained in Molpro for each mode at each point on reaction path.
From the gradient, reduced mass and the mode scan in normal mode, we could distinguish which is frequency of reaction path (gradient max one), rotation and translation. Then we obtained the organized frequencies in table 3.4 and. Figure 3.3.
Those frequencies belong to the modes of reactants (OH, H2) and products (H2O);
W3~W4 in table 3.4, are reasonable. But W1 and W2 have unphysical imaginary frequencies during the reaction path. It’s due to the problem 2 mentioned above. The contamination of reaction coordinate, rotation and translation is stronger for smaller frequencies (those contribute to the rotation of reactants and products, that is, W1 and W2). But the program for dealing problem2 is no easy to write. Even if we have done this, there will be still a big problem, which is shown by Truhlar6. You still may get unphysical imaginary frequencies not until you use curvilinear internal coordinate.
This is even harder than dealing with problem 2, since it’s a non linear process. So we modified the W1 and W2 by following the trends of frequencies around s = 0 and goes to frequencies of rotation in the asymptotic region (Originally, it should be zero, but we know it’s contributed by the rotation in the asymptotic region, which won’t be defined as frequency). Then, we fit the curves in the reaction region (s > 0) and product region( s < 0). This process is shown in figure3.4. (Only show W1 and W2, since W3~W4 don’t need to modification.) The fitting process will also be done for W3, W4, W5 and MEP. The fitting result for W1, W2 and W3 of OH+H2 part (all change to rotation at asymptotic region) is shown as bellow:
W1:
OH+H2 part ( The rotation of OH, s > 0) :
) , where 18.601cm-1 is the rotation constant B of OH in asymptotic region.
H+H2O part (The total rotational J of H2O, here we approximate it as prolate with
, where 11.92cm-1 is the rotation constant B’ of H2O in asymptotic region.
W2: , where 0cm-1 is the rotation constant B’ of OH/H2 in asymptotic region.
H+H2O part (The k component of H2O ):
, where 14.875cm-1 is the rotation constant A’ of H2O in asymptotic region.
W3:
, where 57.66cm-1 is the rotation constant B of H2 in asymptotic region.
We fit these frequencies directly to certain rotation, which in fact is not the case. Since the motion of those rotations in asymptotic region are mixed together in the complex HOHH. In reality, it’s not a big problem for two reasons: (1) The couplings between these rotations could be considered during the process of non-adiabatic transition. (2) We could take the rotation as an average field, that is, we only need to know this rotation belong to which vibration without knowing what kind of rotation it is, which is more general for large system. So from the result we can gain the fitted frequencies
along the reaction path as table 3.5 and figure 3.5
3-2.Contruct the adiabatic curves of J = 0 on MEP
By the fitted frequencies along the MEP, we could construct the adiabatic curves of J =0 on it with second and third term in the right side of the equation (2.9).
But there are some modifications, which is necessary. Those frequencies that become rotation in the asymptotic region don’t always preserved the form of Eni = wi(ni+1/2).
It evolves into the form of rotation in regions that are away from transition state.
Besides that, for W2 of OH+H2, the rotational quantum number itself is restricted to the rotational quantum numbers of OH and H2 rotation, and for W1 of H+H2O, it’s restricted by rotational total J of H2O. The restriction seems unreasonable, but in fact ok as the two points I mentioned above. So the modification will be done as bellow:
w J s a w n w J J s a w
J s a
EJ = ( )⋅ ⋅ +(1− ( ))⋅( 2+ )⋅ = ⋅ − ( )⋅ 2⋅ (3.11) Or
w K s a w
K s a
EK = ( )⋅ ⋅ +(1− ( ))⋅( 2)⋅ (3.12) , in which J is the rotational total quantum number and K is the quantum number of z component of J in molecular axis. a(s) is percentage for the portion of vibration and (1-a(s)) then is the portion of rotation. ZPE is directly considered with MEP, so no 1/2 for vibration in equation (3.11) and (3.12). The way of how to obtain a(s) is pretty easy, just follow the trends of how frequencies drop in equation (3.6)~(3.10). Take W1 for example, you could get a(s) for :
W1 of OH+H2:
) 0.56019 /
exp(
)
(s s
a = − (3.13) W1 of H+H2O:
[
( 0.24983)/0.48136]
) So En for adiabatic curves would be:OH+H2 side : momentum equals to zero.
H+H2O side :
< 0). Thus from (3.15) and (3.16) we can gain the adiabatic curves as figure 3.6
3-3. Cumulative reaction probability
3-3.1 Cumulative reaction probability for OH+H
2HCl+H
2O
Before talking about the state to state dynamics, we would like to show the cumulative reaction probability first. Cumulative reaction probability is the property that indicates how much reactants become products at certain total energy, which could be shown as below:
∑∑
→, in which N(E)(J) is cumulative reaction probability at total energy E, pif(J) stands
for the reaction probability from certain state i of reactants to certain state f of products and both of them are the results of total angular momentum equals to J. From equation (3.18), it is easy to understand that N(E) is macroscopic properties, which means that the result of N(E) decides whether the potential energy curves we use is right or not, since it’s a total effect, that is , it’s easier to achieve than pi
f(J), which is a microscopic property. Here, we only show the result of N(E)(J=0), since other N(E)(J) could be obtain by J shift approximation, that is, we only need to make sure N(E)(J=0) is correct. Below we would like to display several results of N(E)(J=0) from different types of model :
Type I :
We only consider the energy levels at the two sides of adiabatic curves, s=-2.5bohr and s=2.0bohr with the cumulative reaction probability as bellow
f r
f r
N N
N J N
E
N +
= ⋅
= )0 )(
( (3.19)
, where Nr, Nf are, respectively, the number of states available at E for s=2.5bohr and s=2.0bohr. This N(E)(J=0) is the result of barrierless and every states are equal partition if the states are available at that E . We do the same thing for reactants and products, where Nr, Nf would then, respectively, be the number of states available at E for reactants and products. Then we obtain the result as figure 3.7. The outcome is larger than Miller’s result7, which is reasonable, since no barrier in this case. Th S=-2.5bohr/s=2.0bohr, these two positions are pretty close to the asymptotic region, and from figure 3.7, we could find that it goes well with the trend of reactant and product. That is, the linkage between the two sides of the energy curves and asymptotic region is pretty good.
Type II :
We consider the entire curves as shown in figure 3.6, but we didn’t evaluate the total S matrix as equations (2.26) to (2.35) in reality. We only consider the effect of tunneling of each curve. Because of the macroscopic properties of cumulative reaction probability, it is easy to assume that the N(E) is contributed mainly by the tunneling property of each adiabatic curve not by the non-adiabatic transition between them, which is a microscopic property that would easily vanish in the macroscopic properties. The strategy we applied is pretty easy, if there is only one barrier in each adiabatic curve, we just evaluate the tunneling probability of each curves and then do summation. But if any curve inside the adiabatic curves has more than one barrier, evaluating the S matrix becomes necessary. Since OH+H2 side, reactant side, is higher than H2O+H side in each adiabatic curve, that is, the possibilities of tunneling would be decided only by the energy levels of OH+H2 in the asymptotic region. From figure 3.8 which indicated the adiabatic curves that connects the energy levels of reactants lower than 15kcal/mol. It is obvious that the assumption of one barrier is sufficient, even if there is some pretty small well in higher levels, which won’t give big effect.
The transmittance probability is evaluated as bellow:
∫
−= 2
1
) ( 1 T 2
T µE V s ds
δ h (3.20)
,where T1 and T2 are the turning points in each adiabatic curve. If the adiabatic curves are in parabolic shape, then the transmittance probability are in the form as equation (3.21) and (3.22):
2 max 2
tan
1 for E E
e
ptransmit ce e total >
= + −
− δ δ
(3.21)
2 max 2
tan
1 forE E
e
ptransmit ce e total <
= + δ δ (3.22) , in which Etotal is the total energy and Emax is the max energy in each adiabatic
curve. In type II, we obtain the result as figure 3.9. The deviation from Miller’s result7 is not small. There are three possible reasons:
1. The linkage between reactants and products is not appropriate, since there two rotational related vibratioanl modes in product side, but three for reactant side.
Possible solution would be the way of dealing with the mode, which shows vibrational behavior in one side but rotational behavior in other side. We should consider the evolution of this kind of mode from reactant to product, not just take one side as vibration and the evolution only in another side. Besides, the way of dealing with the evolution from vibration to rotation should be modified also, and this would be done in the future.
2. The non-adiabatic transition may cause some effect to the cumulative reaction probability.
3. The accuracy of the frequencies along the reaction coordinate would influence the density of states. For those big frequencies, the defects may be small ,but not the case for those small frequencies, especially those rotational related modes in the asymptotic regions.
Type III:
Since the adiabatic curves are not appropriate yet, in order to obtain the cumulative reaction probability, we’ve done following three assumptions:
1. The shape of upper levels is same as the ground state, and take all of curves as parabolic model
2. Take the energy level of transition state as the max energy, Emax, of each adiabatic curve.
3. From assumption 1 and 2, we could make a subroutine of delta (δ ) vs deltaE for ground energy curve, where deltaE is the difference between Etotal and Emax.
Then take this subroutine as the reference for upper energy levels. So for each adiabatic curves we only need to know two things, one is Etotal > Emax or Etotal< Emax, another is difference between Emax and Etotal. Then from the subroutine, we could get its correspond delta(δ ). For Etotal > Emax, equation (3.21) is applied to get Ptransmitance. For Etotal < Emax, equation (3.22) is applied to get Ptransmitance.
We call this model as adiabatic model. The energy diagram of this model is shown in figure 3.10. The delta (δ ) vs deltaE for ground energy curve is displayed in figure 3.11 and the fitted delta (δ ) vs deltaE is shown in equation (3.23):
4 The cumulative reaction probability of J = 0 is in figure 3.12. We could found that it’s much closer to Miller’s result7 than type I, but still has not small difference. There would three possible reasons:
1. From figure 3.10, it is easy to see that the energy density of reactant side is larger than transition state, that is, the delta would be smaller than ground energy curve at the same delta E for the upper curves, as the parabolic model still holds. So we’ve done some modification for upper curves as bellow:
max
of adiabatic curve considered and ∆ v equals to Emax-v0. The result of this
modification is shown in figure 3.13. No big difference between before and after
modifications. But after directly multiple adiabatic models result by a factor 3.6 or
the modification order change from 0.5 to 2.95 as bellow:
2.95 max
we found that the result match Miller’s result7 better. These are displayed in figure
3.14. The possible explanations would be the parabolic model is not appropriate or
the density of state in transition state is not correct which would be explained in
reason 3. We fit the ground energy curve in polynomial equation, and obtained the
equation (3.28)
From equation (3.28), we could find that s2 is not the only dominated term, that is,
parabolic model no longer good, but whether this is the main reason, it is still in
request.
2. Non-adiabatic transition between each adiabatic curve may cause some effects, which would be shown in type IV.
3. The accuracy of the frequencies along the reaction coordinate would influence the density of states. For those big frequencies, the defects may be small ,but not the
case for those small frequencies, especially those rotational related modes in the asymptotic regions.
Type IV:
In order to consider the influence of non-adiabatic transition, we use diabatic model to evaluate the energy curves for the transition state. In diabatic model, all modes in transition state preserve the properties of the motion in reactants. Since diabatic model keep all the properties from reactant to transition state, that is, the non-adiabatic coupling is considered entirely in each curves. So the five vibrational frequencies 546.6cm-1, 609.94cm-1, 1059.12cm-1, 2479.16cm-1 and 3736.43 cm-1 in transition state would be changed as bellow:
1. 546.6cm-1: correlates to the rotation between OH and H2 in asymptotic region, which is zero wave number in asymptotic region, but 4.89cm-1 in transition state, which the radius for moment inertia is calculated by the mass center of OH and H2
2. 609.94cm-1: correlates to the rotation of OH in asymptotic region, which is 18.601cm-1 in asymptotic region, but 18.6cm-1 in transition state. No big differences between asymptotic region and transition state, since OH can be seen as a spectator during the reaction.
3. 1059.12cm-1: correlates to the rotation of H2 in asymptotic region, which is 57.66cm-1 in asymptotic region, but 50cm-1 in transition state.
4. 2479.16cm-1: correlates to the vibration of H2 in asymptotic region, which is 4344.54cm-1 in asymptotic region, but 2479.16cm-1 in transition state.
5. 3736.43cm-1: correlates to the rotation of OH in asymptotic region, which is 3714.33cm-1 in asymptotic region, but 3736.43cm-1 in transition state.
No big differences between asymptotic region and transition state, the reason is same as the concept mentioned in 609.94cm-1
So, for the diabatic model, in order to preserve the properties, that is, the quantum numbers of those five motions in transition state should follow the quantum numbers of reactants. This is pretty different from adiabatic model. In adiabatic model, we count the number of possible tunneling amounts from reactants, and then the Emax of each adiabatic curve are taken from low level to high level, that is, the relation between each Emax and each state of reactants are decided by the sequence of the energy level of transition state and reactants. But for diabatic model, for every state of reactants, we should remember the quantum numbers of those 5 motions, and its corresponded Emax should have the same quantum numbers in order to preserve the properties. So the energy level of transition state would be like equation (3.29)：
50 diabatic model as figure 3.15. We could find out that the density of state in transition state is similar to the density in reactants, since it is diabatic model. Then the process to obtain cumulative reaction probability is entirely same as type III except the
linkage between Emax and state of reactants. So we use the same equation (3.23) for delta (δ ) vs deltaE, then we get the cumulative reaction probability of diabatic model as figure 3.16. In this figure, it also shows the adiabatic model. It’s obvious to see that miller’s result is between adiabatic model and diabatic model, that is to say, non-adiabatic coupling may give some contribution to the cumulative reaction probability, since diabatic model consider the non-adiabatic coupling entirely in each curve, but entirely no non-adiabatic coupling in adiabatic model.
linkage between Emax and state of reactants. So we use the same equation (3.23) for delta (δ ) vs deltaE, then we get the cumulative reaction probability of diabatic model as figure 3.16. In this figure, it also shows the adiabatic model. It’s obvious to see that miller’s result is between adiabatic model and diabatic model, that is to say, non-adiabatic coupling may give some contribution to the cumulative reaction probability, since diabatic model consider the non-adiabatic coupling entirely in each curve, but entirely no non-adiabatic coupling in adiabatic model.