- Email: [email protected]
A new global reaction route map on the potential energy surface of H2CO with unrestricted level
Chemical Physics Letters 460 (2008) 55–58
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
A new global reaction route map on the potential energy surface of H2CO with unrestricted level Satoshi Maeda, Koichi Ohno * Department of Chemistry, Graduate School of Science, Tohoku University, Aramaki, Aoba-ku, Sendai 980-8578, Japan
a r t i c l e
i n f o
Article history: Received 26 February 2008 In final form 2 June 2008 Available online 6 June 2008
a b s t r a c t The potential energy surface (PES) of H2CO has been revisited by the anharmonic downward distortion following algorithm, where an unrestricted ab initio method was employed during the automated exploration. Some structures previously reported on the PES of restricted theories were found to be open-shell species, and new local minima and first-order saddle points were discovered in high energy regions related to radical fragments which are important in high energy processes such as combustion. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction
2. Method
Equilibrium structures (EQ) and transition state structures (TS) on the potential energy surface (PES) of H2CO were studied extensively [1–11]. Especially, the TS for the formaldehyde decomposition into H2 and CO was studied by many groups at various computation levels, and the best theoretical estimate [7] of its barrier height of 81.9 (± 0.3) kcal/mol shows excellent agreement with a recent experimental value [12] of 81.7 (± 0.5) kcal/mol. Moreover, this PES was used as a benchmark of automated global exploration techniques on a given ab initio PES [4,6,8–10]. Such global explorations have been made on PESs of restricted theories: the RHF surface was explored by Bondensgård and Jensen by using the gradient extremal following method [4], the RHF surface was revisited by Quapp et al. by means of their reduced gradient following method [6], the RMP2 and the RB3LYP surfaces were examined by the anharmonic downward distortion following (ADD-following) algorithm [8–10]. In this study, we revisited the PES of H2CO by using the ADD-following method [8–10]. The present purpose is its application to open-shell parts of the ground state PES including high energy regions, and the use of an unrestricted method coupled with a series of the MO stability check during the exploration was tested by using the PES of H2CO as a benchmark system. Some structures previously reported on the PES of restricted theories were found to be open-shell species, and new EQs and TSs have been discovered in high energy regions of the PES.
2.1. Anharmonic downward distortion following algorithm The ADD-following was proposed as an uphill walking approach along reaction pathways from an EQ to TSs and dissociation channels (DC) [8–10]. Since typical reaction paths always change their curvatures from a concave to a convex on going to a TS or a DC, slopes should always decline from respective harmonic curve because of energy lowering interactions leading to a TS or a DC. Therefore, one can assume that a reaction pathway can be found as an ADD around an EQ. To detect an ADD maximal direction in a multidimensional coordinate space, we developed the scaled hypersphere search (SHS) method [8], by which an ADD is searched by a comparison of the reference harmonic PES with the actual anharmonic PES on the iso-energy hypersurface of the harmonic PES. The iso-energy hypersurface constructs a hyperellipsoid centered at an EQ. In the SHS method, the scaled normal coordinates qi (=ki1/2Qi) are employed instead of the usual normal coordinates Qi, where ki is the eigenvalue of Qi. The use of qi converts a hyperellipsoidal iso-energy surface of the harmonic PES into a simple hypersphere. Harmonic energy on this scaled hypersphere is constant, and therefore energy minima on it correspond to ADD maxima. It follows that an ADD maximal direction can be found by a simple constrained energy minimization algorithm. By following ADD maximal directions with expanding hypersphere radii, reaction pathways can be followed from an EQ to TSs and DCs. Further details of the ADD-following by the SHS method and automated procedure of the global reaction route mapping (GRRM) can be found in previous papers [9,10]. 2.2. Exploring on open-shell surface
* Corresponding author. Fax: +81 22 795 6580. E-mail address: [email protected] (K. Ohno). 0009-2614/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2008.06.003
We have applied the ADD-following method to unimolecular systems [13–16], H-bond cluster systems [17–19], and a large
56
S. Maeda, K. Ohno / Chemical Physics Letters 460 (2008) 55–58
Here, we consider the use of the UB3LYP method coupled with the MO stability check, where the MO stability check is performed in every prediction (sphere expansion) steps of the SHS procedure. The use of the MO stability check is necessary, since a search starting from a closed-shell part of the PES tend to remain on the closed-shell surface. The radical channel region is further explored in this study. Although a certain threshold is used in the GRRM program for recognition of DC, 1.5 times larger threshold value with respect to the default for unimolecular systems is used in the application to the UB3LYP surface. Since the radical channel region of PES is fairly flat, a half step size for the sphere expansion procedure is used not to overlook loose TSs in such regions. Totally, 32 376 times force calculations and 2222 times Hessian calcula-
organometallic system [20], within the restricted single-reference theory framework. In cases of H-bond clusters and an organometallic system, we explored only low energy parts of PESs connected via weak bond rearrangement pathways (i.e., H-bond and coordination bond reorganization reactions) by using the large-ADDfollowing algorithm [17], and the restricted theory framework is sufficient for such reactions. On the other hand, high energy species with unsaturated electronic structures may have a significant biradical character, and therefore one can expect that our global reaction route maps for unimolecular systems may overlook some biradical EQs and TSs. In this study, to explore on open-shell surfaces by the GRRM program [8–10,17,20,21], we introduced some modifications which can be used by options.
O
O
O
O C
C
C
C TS5/6 (Cs) 382.5 kJ/mol
EQ5 (Cs) 379.0 kJ/mol
EQ6 (Cs) 381.9 kJ/mol
EQ8 (Cs) 581.9 kJ/mol
O
C O TS4/5 (Cs) 379.1 kJ/mol
TS3/6 (Cs) 402.6 kJ/mol
TS3/8 (C1) 627.5 kJ/mol
O C
C
O
TS3/4 (C1) 329.1 kJ/mol
O C
O
C
O
C
TS7/8 (C1) 654.8 kJ/mol
C
EQ4 (Cs) 238.7 kJ/mol
TS3/7 (C1) 555.4 kJ/mol
EQ3 (Cs) 220.0 kJ/mol
O C TS2/5 (Cs) 382.4 kJ/mol
O C O O C TS0/3 (Cs) 360.7 kJ/mol
C TS2/4 (Cs) 450.7 kJ/mol
O O C O
EQ0 (C2v) 0.0 kJ/mol
C O
TS0/6 (Cs) 382.7 kJ/mol
C
TS0/2 (Cs) 346.9 kJ/mol
C
O O
O
O
C EQ2 (Coov) 32.0 kJ/mol
TS1/2 (Cs) 33.1 kJ/mol
C EQ1 (Cs) 31.1 kJ/mol
C TS1/7 (Cs) 555.4 kJ/mol
C EQ7 (Cs) 554.5 kJ/mol
Fig. 1. The global reaction route map of the potential energy surface of H2CO at a computation level of UB3LYP/cc-pVDZ. EQn are equilibrium structures (local minima), and TSn/m are transition structures (first-order saddle points) connecting EQn and EQm. Symmetry of each structure is indicated in parentheses. A chemical bond, a weak bond, and a rearranging bond at each TS are shown by a double solid line, a dotted line, and a thin line, respectively.
57
S. Maeda, K. Ohno / Chemical Physics Letters 460 (2008) 55–58
2.3. Computations
tions were made for this calculation. These numbers are larger than triples of the case of the previous application to the RB3LYP surface [10]. This is due to the decrease of the sphere expansion step size and the increase of the DC recognition threshold. In this study, we further consider the CASSCF surface for comparisons. In this case, we restrict all geometrical displacements to be smaller than 0.1 Å, because the SCF convergence in multi-reference methods is very sensitive to a quality of the MO guess. Since the SCF convergence was sometimes very slow in regions of biradical species, the CASSCF method is applied only to the chemical bonding region of the PES, where 11 495 times force calculations and 859 times Hessian calculations were required.
The cc-pVDZ basis set was employed in both UB3LYP and CASSCF calculations. Six orbitals and six electrons were included in the active space of the CASSCF calculation, i.e., [6,6]-CASSCF. Geometries and intrinsic reaction coordinate (IRC) [22] connections were refined at the UCCSD/aug-cc-pVDZ level, and then single point energies at the UCCSD geometries were computed at the UCCSD(T)/aug-cc-pVTZ level. The automated exploration was performed by the GRRM1.20 program [8–10,17,20,21], in which energy, gradient, and Hessian are obtained from the GAUSSIAN 03 programs [23].
O
O
O C C
TS0/5 (Cs) 394.8 kJ/mol (383.5 kJ/mol)
EQ5 (Cs) 394.5 kJ/mol (383.3 kJ/mol)
C O
EQ8 (Cs) 574.5 kJ/mol (547.0 kJ/mol)
C O C
O TS0/3 (Cs) 435.7 kJ/mol (434.1 kJ/mol) C
TS4/5 (Cs) 404.9 kJ/mol (404.4 kJ/mol)
O
O
C
C
EQ4 (Cs) 235.9 kJ/mol (231.1 kJ/mol)
O
TS1/8 (C1) O 689.6 kJ/mol (683.3 kJ/mol) C
O C
TS3/4 (C1) 339.9 kJ/mol (330.3 kJ/mol)
EQ3 (Cs) 215.8 kJ/mol (210.5 kJ/mol)
O
C TS2/5 (Cs) 394.6 kJ/mol (383.4 kJ/mol)
TS3/8 (C1) 642.7 kJ/mol (633.1 kJ/mol)
C O
O
C TS0/3 (Cs) 359.2 kJ/mol (369.3 kJ/mol)
C TS2/4 (Cs) 454.0 kJ/mol (460.4 kJ/mol)
TS1/3 (C1) 564.8 kJ/mol (543.8 kJ/mol)
O C O EQ0 (C2v) 0.0 kJ/mol (0.0 kJ/mol)
C O
TS0/2 (Cs) 364.9 kJ/mol (373.6 kJ/mol)
C
O O C
EQ2 (Coov) 19.9 kJ/mol (16.0 kJ/mol)
TS1/2 (Cs) 20.6 kJ/mol (16.6 kJ/mol)
C EQ1 (Coov) 20.1 kJ/mol (16.1 kJ/mol)
Fig. 2. The global reaction route map of the potential energy surface of H2CO at a computation level of UCCSD(T)/aug-cc-pVTZ//UCCSD/aug-cc-pVDZ. Values in parentheses are UCCSD/aug-cc-pVDZ energies. Other notations are the same as Fig. 1.
58
S. Maeda, K. Ohno / Chemical Physics Letters 460 (2008) 55–58
3. Results and discussions Fig. 1 shows a new global reaction route map for H2CO based on the UB3LYP method, where EQn is a local minimum and TSn/m is a first-order saddle point for an IRC connection between EQn and EQm. In Fig. 1, EQ5, EQ6, EQ7, EQ8, TS0/6, TS1/7, TS2/5, TS3/4, TS3/6, TS3/7, TS3/8, TS4/5, TS5/6, and TS7/8 have non-zero spin expectation values of 0.5–1.0, and they can be regarded to be open-shell species. In the case of the [6,6]-CASSCF surface, EQ0, EQ3, EQ4, EQ8, TS0/2, TS0/3, TS2/4, TS3/4, TS3/6, TS3/8, TS4/5, and TS7/8 were obtained in the chemical bonding region, and others were recognized as DC during the exploration because of the smaller threshold for the DC recognition. Concerning with the chemical bonding region, the [6,6]-CASSCF surface and the UB3LYP surface are qualitatively equivalent to each other. Around stable isomers of EQ0, EQ3, and EQ4, four TSs, i.e., TS0/6, TS3/6, TS3/7, and TS4/5, were newly added to the map for the PES of H2CO. All these four TSs are open-shell species, and all of previous automated explorations based on restricted theories overlooked them [4,6,8–10]. Although EQ3 and EQ4 seem to be the carbene with an unsaturated carbon atom, these two structures are closed-shell species. Since EQ8 contains a highly unsaturated atomic carbon fragment, EQ8 has a significant open-shell character, and the CO distance in EQ8 of 1.87 Å is enlarged with respect to the previous geometry [10] of 1.46 Å based on the RB3LYP method. The longer OH distance in EQ7 of 2.67 Å is also enlarged with respect to the previous geometry [10] of 1.76 Å. Another remarkable change from previous global maps is an appearance of the path via TS0/6, EQ6, TS5/6, EQ5, and TS2/5 into the molecular products of H2 and CO. A similar reaction route is known as the roaming pathway for the photo-dissociation of formaldehyde [24]. Another TS where an H atom go around C=O bond axis, i.e., an out-of-plane route, was discovered for the roaming pathway recently [11], and the present one is an alternative in-plane type. These open-shell parts are mostly related to radical fragments, where recombination reactions between such fragments are very important in some reaction fields such as flames, the interstellar space, the atmosphere, and so on. The unrestricted hybrid-DFT method can be applied to problems involving biradical species such as drawing radical dissociation curves, although it may not be very accurate because of the spin contamination [25]. Another well-known drawback of the DFT method [26] is that it cannot precisely describe the dispersion force which is important to compute an accurate binding energy between a pair of fragments. Hence, geometries and IRC connections were re-computed at the UCCSD level using the UB3LYP structures as initial guesses. Fig 2 gives the GRRM of H2CO at the UCCSD(T)/aug-cc-pVTZ//UCCSD/aug-cc-pVDZ level. Cartesian geometries, relative potential energy values, harmonic zero-point energy values, and harmonic frequencies for all stationary points on the UCCSD surface are listed in Supplementary data. Here, two weak complexes (EQ6 and EQ7) are not seen on the UCCSD surface, and TSs related to them are directly connected with next EQs. However, overall atom movements along IRC paths via these disappeared intermediates are very similar to that on the UCCSD surface. For an example, the IRC starting from TS0/5 on the UCCSD surface passes structures like EQ6 and TS0/6. In other words, the global IRC path network on the UB3LYP surface well reproduce that on the UCCSD surface. Although geometries of weak complexes such as EQ6 and EQ7 are very sensitive to levels of ab initio calcu-
lations, we were able to reveal qualitative overview of energy landscapes on the UCCSD surface by the initial search at the UB3LYP/cc-pVDZ level. The present treatment, i.e., the ADD-following via an unrestricted hybrid-DFT with a series of the MO stability check, will be used for an initial automated exploration of whole area on a PES including such radical fragment regions. 4. Conclusions The PES of H2CO were explored by the ADD-following method coupled with the UB3LYP method and the CASSCF method, and then geometries and IRC connections were refined at the UCCSD level. This is the first automated systematic global reaction route mapping on the PES of H2CO via unrestricted theories, where we introduced some modifications to explore on open-shell surfaces by the ADD-following method. Global topography on these three surfaces are qualitatively equivalent to each other. Although, in low energy regions of the PES, there are few changes in the topography on the PES with respect to those of restricted theories, new structures have been discovered in high energy regions involving radical fragments. It follows that the ADD-following via an unrestricted hybrid-DFT with a series of the MO stability check may provides new structures in such regions even for systems well studied via restricted theories. Acknowledgements S. Maeda is supported by Research Fellowship of the Japan Society for Promotion of Science for Young Scientists. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.cplett.2008.06.003. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26]
J.D. Goddard, H.F. Schaefer III, J. Chem. Phys. 70 (1979) 5117. G.E. Scuseria, H.F. Schaefer III, J. Chem. Phys. 90 (1989) 3629. L. Deng, T. Ziegler, L. Fan, J. Chem. Phys. 99 (1993) 3823. K. Bondensgård, F. Jensen, J. Chem. Phys. 104 (1996) 8025. H. Nakano, K. Nakayama, K. Hirao, M. Dupuis, J. Chem. Phys. 106 (1997) 4912. W. Quapp, M. Hirsch, O. Imig, D. Heidrich, J. Comput. Chem. 19 (1998) 1087. D. Feller, M. Dupuis, B.C. Garrett, J. Chem. Phys. 113 (2000) 218. K. Ohno, S. Maeda, Chem. Phys. Lett. 384 (2004) 277. S. Maeda, K. Ohno, J. Phys. Chem. A 109 (2005) 5742. K. Ohno, S. Maeda, J. Phys. Chem. A 110 (2006) 8933. L.B. Harding, S.J. Klippenstein, A.W. Jasper, Phys. Chem. Chem. Phys. 9 (2007) 4055. J. Troe, J. Phys. Chem. A 111 (2007) 3862. X. Yang, S. Maeda, K. Ohno, J. Phys. Chem. A 109 (2005) 7319. X. Yang, S. Maeda, K. Ohno, Chem. Phys. Lett. 418 (2006) 208. X. Yang, S. Maeda, K. Ohno, J. Phys. Chem. A 111 (2007) 5099. Y. Watanabe, S. Maeda, K. Ohno, Chem. Phys. Lett. 447 (2007) 21. S. Maeda, K. Ohno, J. Phys. Chem. A 111 (2007) 4527. Y. Luo, S. Maeda, K. Ohno, J. Phys. Chem. A 111 (2007) 10732. S. Maeda, K. Ohno, J. Phys. Chem. A 112 (2008) 2962. S. Maeda, K. Ohno, J. Phys. Chem. A 111 (2007) 13168. The execution module of the GRRM1.20 program is available by contact with the corresponding author (see http://qpcrkk.chem.tohoku.ac.jp/ for details). K. Fukui, Acc. Chem. Res. 14 (1981) 63. M.J. Frisch et al., GAUSSIAN 03, Revision , C.02. Gaussian, Inc., Wallingford CT, 2004. D. Townsend et al., Science, 306 (2004), p. 1158. F. Jensen, Introduction to Computational Chemistry, Wiley, Chichester, 1998. E.R. Johnson, G.A. DiLabio, Chem. Phys. Lett. 419 (2006) 333.
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
A new global reaction route map on the potential energy surface of H2CO with unrestricted level Satoshi Maeda, Koichi Ohno * Department of Chemistry, Graduate School of Science, Tohoku University, Aramaki, Aoba-ku, Sendai 980-8578, Japan
a r t i c l e
i n f o
Article history: Received 26 February 2008 In final form 2 June 2008 Available online 6 June 2008
a b s t r a c t The potential energy surface (PES) of H2CO has been revisited by the anharmonic downward distortion following algorithm, where an unrestricted ab initio method was employed during the automated exploration. Some structures previously reported on the PES of restricted theories were found to be open-shell species, and new local minima and first-order saddle points were discovered in high energy regions related to radical fragments which are important in high energy processes such as combustion. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction
2. Method
Equilibrium structures (EQ) and transition state structures (TS) on the potential energy surface (PES) of H2CO were studied extensively [1–11]. Especially, the TS for the formaldehyde decomposition into H2 and CO was studied by many groups at various computation levels, and the best theoretical estimate [7] of its barrier height of 81.9 (± 0.3) kcal/mol shows excellent agreement with a recent experimental value [12] of 81.7 (± 0.5) kcal/mol. Moreover, this PES was used as a benchmark of automated global exploration techniques on a given ab initio PES [4,6,8–10]. Such global explorations have been made on PESs of restricted theories: the RHF surface was explored by Bondensgård and Jensen by using the gradient extremal following method [4], the RHF surface was revisited by Quapp et al. by means of their reduced gradient following method [6], the RMP2 and the RB3LYP surfaces were examined by the anharmonic downward distortion following (ADD-following) algorithm [8–10]. In this study, we revisited the PES of H2CO by using the ADD-following method [8–10]. The present purpose is its application to open-shell parts of the ground state PES including high energy regions, and the use of an unrestricted method coupled with a series of the MO stability check during the exploration was tested by using the PES of H2CO as a benchmark system. Some structures previously reported on the PES of restricted theories were found to be open-shell species, and new EQs and TSs have been discovered in high energy regions of the PES.
2.1. Anharmonic downward distortion following algorithm The ADD-following was proposed as an uphill walking approach along reaction pathways from an EQ to TSs and dissociation channels (DC) [8–10]. Since typical reaction paths always change their curvatures from a concave to a convex on going to a TS or a DC, slopes should always decline from respective harmonic curve because of energy lowering interactions leading to a TS or a DC. Therefore, one can assume that a reaction pathway can be found as an ADD around an EQ. To detect an ADD maximal direction in a multidimensional coordinate space, we developed the scaled hypersphere search (SHS) method [8], by which an ADD is searched by a comparison of the reference harmonic PES with the actual anharmonic PES on the iso-energy hypersurface of the harmonic PES. The iso-energy hypersurface constructs a hyperellipsoid centered at an EQ. In the SHS method, the scaled normal coordinates qi (=ki1/2Qi) are employed instead of the usual normal coordinates Qi, where ki is the eigenvalue of Qi. The use of qi converts a hyperellipsoidal iso-energy surface of the harmonic PES into a simple hypersphere. Harmonic energy on this scaled hypersphere is constant, and therefore energy minima on it correspond to ADD maxima. It follows that an ADD maximal direction can be found by a simple constrained energy minimization algorithm. By following ADD maximal directions with expanding hypersphere radii, reaction pathways can be followed from an EQ to TSs and DCs. Further details of the ADD-following by the SHS method and automated procedure of the global reaction route mapping (GRRM) can be found in previous papers [9,10]. 2.2. Exploring on open-shell surface
* Corresponding author. Fax: +81 22 795 6580. E-mail address: [email protected] (K. Ohno). 0009-2614/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2008.06.003
We have applied the ADD-following method to unimolecular systems [13–16], H-bond cluster systems [17–19], and a large
56
S. Maeda, K. Ohno / Chemical Physics Letters 460 (2008) 55–58
Here, we consider the use of the UB3LYP method coupled with the MO stability check, where the MO stability check is performed in every prediction (sphere expansion) steps of the SHS procedure. The use of the MO stability check is necessary, since a search starting from a closed-shell part of the PES tend to remain on the closed-shell surface. The radical channel region is further explored in this study. Although a certain threshold is used in the GRRM program for recognition of DC, 1.5 times larger threshold value with respect to the default for unimolecular systems is used in the application to the UB3LYP surface. Since the radical channel region of PES is fairly flat, a half step size for the sphere expansion procedure is used not to overlook loose TSs in such regions. Totally, 32 376 times force calculations and 2222 times Hessian calcula-
organometallic system [20], within the restricted single-reference theory framework. In cases of H-bond clusters and an organometallic system, we explored only low energy parts of PESs connected via weak bond rearrangement pathways (i.e., H-bond and coordination bond reorganization reactions) by using the large-ADDfollowing algorithm [17], and the restricted theory framework is sufficient for such reactions. On the other hand, high energy species with unsaturated electronic structures may have a significant biradical character, and therefore one can expect that our global reaction route maps for unimolecular systems may overlook some biradical EQs and TSs. In this study, to explore on open-shell surfaces by the GRRM program [8–10,17,20,21], we introduced some modifications which can be used by options.
O
O
O
O C
C
C
C TS5/6 (Cs) 382.5 kJ/mol
EQ5 (Cs) 379.0 kJ/mol
EQ6 (Cs) 381.9 kJ/mol
EQ8 (Cs) 581.9 kJ/mol
O
C O TS4/5 (Cs) 379.1 kJ/mol
TS3/6 (Cs) 402.6 kJ/mol
TS3/8 (C1) 627.5 kJ/mol
O C
C
O
TS3/4 (C1) 329.1 kJ/mol
O C
O
C
O
C
TS7/8 (C1) 654.8 kJ/mol
C
EQ4 (Cs) 238.7 kJ/mol
TS3/7 (C1) 555.4 kJ/mol
EQ3 (Cs) 220.0 kJ/mol
O C TS2/5 (Cs) 382.4 kJ/mol
O C O O C TS0/3 (Cs) 360.7 kJ/mol
C TS2/4 (Cs) 450.7 kJ/mol
O O C O
EQ0 (C2v) 0.0 kJ/mol
C O
TS0/6 (Cs) 382.7 kJ/mol
C
TS0/2 (Cs) 346.9 kJ/mol
C
O O
O
O
C EQ2 (Coov) 32.0 kJ/mol
TS1/2 (Cs) 33.1 kJ/mol
C EQ1 (Cs) 31.1 kJ/mol
C TS1/7 (Cs) 555.4 kJ/mol
C EQ7 (Cs) 554.5 kJ/mol
Fig. 1. The global reaction route map of the potential energy surface of H2CO at a computation level of UB3LYP/cc-pVDZ. EQn are equilibrium structures (local minima), and TSn/m are transition structures (first-order saddle points) connecting EQn and EQm. Symmetry of each structure is indicated in parentheses. A chemical bond, a weak bond, and a rearranging bond at each TS are shown by a double solid line, a dotted line, and a thin line, respectively.
57
S. Maeda, K. Ohno / Chemical Physics Letters 460 (2008) 55–58
2.3. Computations
tions were made for this calculation. These numbers are larger than triples of the case of the previous application to the RB3LYP surface [10]. This is due to the decrease of the sphere expansion step size and the increase of the DC recognition threshold. In this study, we further consider the CASSCF surface for comparisons. In this case, we restrict all geometrical displacements to be smaller than 0.1 Å, because the SCF convergence in multi-reference methods is very sensitive to a quality of the MO guess. Since the SCF convergence was sometimes very slow in regions of biradical species, the CASSCF method is applied only to the chemical bonding region of the PES, where 11 495 times force calculations and 859 times Hessian calculations were required.
The cc-pVDZ basis set was employed in both UB3LYP and CASSCF calculations. Six orbitals and six electrons were included in the active space of the CASSCF calculation, i.e., [6,6]-CASSCF. Geometries and intrinsic reaction coordinate (IRC) [22] connections were refined at the UCCSD/aug-cc-pVDZ level, and then single point energies at the UCCSD geometries were computed at the UCCSD(T)/aug-cc-pVTZ level. The automated exploration was performed by the GRRM1.20 program [8–10,17,20,21], in which energy, gradient, and Hessian are obtained from the GAUSSIAN 03 programs [23].
O
O
O C C
TS0/5 (Cs) 394.8 kJ/mol (383.5 kJ/mol)
EQ5 (Cs) 394.5 kJ/mol (383.3 kJ/mol)
C O
EQ8 (Cs) 574.5 kJ/mol (547.0 kJ/mol)
C O C
O TS0/3 (Cs) 435.7 kJ/mol (434.1 kJ/mol) C
TS4/5 (Cs) 404.9 kJ/mol (404.4 kJ/mol)
O
O
C
C
EQ4 (Cs) 235.9 kJ/mol (231.1 kJ/mol)
O
TS1/8 (C1) O 689.6 kJ/mol (683.3 kJ/mol) C
O C
TS3/4 (C1) 339.9 kJ/mol (330.3 kJ/mol)
EQ3 (Cs) 215.8 kJ/mol (210.5 kJ/mol)
O
C TS2/5 (Cs) 394.6 kJ/mol (383.4 kJ/mol)
TS3/8 (C1) 642.7 kJ/mol (633.1 kJ/mol)
C O
O
C TS0/3 (Cs) 359.2 kJ/mol (369.3 kJ/mol)
C TS2/4 (Cs) 454.0 kJ/mol (460.4 kJ/mol)
TS1/3 (C1) 564.8 kJ/mol (543.8 kJ/mol)
O C O EQ0 (C2v) 0.0 kJ/mol (0.0 kJ/mol)
C O
TS0/2 (Cs) 364.9 kJ/mol (373.6 kJ/mol)
C
O O C
EQ2 (Coov) 19.9 kJ/mol (16.0 kJ/mol)
TS1/2 (Cs) 20.6 kJ/mol (16.6 kJ/mol)
C EQ1 (Coov) 20.1 kJ/mol (16.1 kJ/mol)
Fig. 2. The global reaction route map of the potential energy surface of H2CO at a computation level of UCCSD(T)/aug-cc-pVTZ//UCCSD/aug-cc-pVDZ. Values in parentheses are UCCSD/aug-cc-pVDZ energies. Other notations are the same as Fig. 1.
58
S. Maeda, K. Ohno / Chemical Physics Letters 460 (2008) 55–58
3. Results and discussions Fig. 1 shows a new global reaction route map for H2CO based on the UB3LYP method, where EQn is a local minimum and TSn/m is a first-order saddle point for an IRC connection between EQn and EQm. In Fig. 1, EQ5, EQ6, EQ7, EQ8, TS0/6, TS1/7, TS2/5, TS3/4, TS3/6, TS3/7, TS3/8, TS4/5, TS5/6, and TS7/8 have non-zero spin expectation values of 0.5–1.0, and they can be regarded to be open-shell species. In the case of the [6,6]-CASSCF surface, EQ0, EQ3, EQ4, EQ8, TS0/2, TS0/3, TS2/4, TS3/4, TS3/6, TS3/8, TS4/5, and TS7/8 were obtained in the chemical bonding region, and others were recognized as DC during the exploration because of the smaller threshold for the DC recognition. Concerning with the chemical bonding region, the [6,6]-CASSCF surface and the UB3LYP surface are qualitatively equivalent to each other. Around stable isomers of EQ0, EQ3, and EQ4, four TSs, i.e., TS0/6, TS3/6, TS3/7, and TS4/5, were newly added to the map for the PES of H2CO. All these four TSs are open-shell species, and all of previous automated explorations based on restricted theories overlooked them [4,6,8–10]. Although EQ3 and EQ4 seem to be the carbene with an unsaturated carbon atom, these two structures are closed-shell species. Since EQ8 contains a highly unsaturated atomic carbon fragment, EQ8 has a significant open-shell character, and the CO distance in EQ8 of 1.87 Å is enlarged with respect to the previous geometry [10] of 1.46 Å based on the RB3LYP method. The longer OH distance in EQ7 of 2.67 Å is also enlarged with respect to the previous geometry [10] of 1.76 Å. Another remarkable change from previous global maps is an appearance of the path via TS0/6, EQ6, TS5/6, EQ5, and TS2/5 into the molecular products of H2 and CO. A similar reaction route is known as the roaming pathway for the photo-dissociation of formaldehyde [24]. Another TS where an H atom go around C=O bond axis, i.e., an out-of-plane route, was discovered for the roaming pathway recently [11], and the present one is an alternative in-plane type. These open-shell parts are mostly related to radical fragments, where recombination reactions between such fragments are very important in some reaction fields such as flames, the interstellar space, the atmosphere, and so on. The unrestricted hybrid-DFT method can be applied to problems involving biradical species such as drawing radical dissociation curves, although it may not be very accurate because of the spin contamination [25]. Another well-known drawback of the DFT method [26] is that it cannot precisely describe the dispersion force which is important to compute an accurate binding energy between a pair of fragments. Hence, geometries and IRC connections were re-computed at the UCCSD level using the UB3LYP structures as initial guesses. Fig 2 gives the GRRM of H2CO at the UCCSD(T)/aug-cc-pVTZ//UCCSD/aug-cc-pVDZ level. Cartesian geometries, relative potential energy values, harmonic zero-point energy values, and harmonic frequencies for all stationary points on the UCCSD surface are listed in Supplementary data. Here, two weak complexes (EQ6 and EQ7) are not seen on the UCCSD surface, and TSs related to them are directly connected with next EQs. However, overall atom movements along IRC paths via these disappeared intermediates are very similar to that on the UCCSD surface. For an example, the IRC starting from TS0/5 on the UCCSD surface passes structures like EQ6 and TS0/6. In other words, the global IRC path network on the UB3LYP surface well reproduce that on the UCCSD surface. Although geometries of weak complexes such as EQ6 and EQ7 are very sensitive to levels of ab initio calcu-
lations, we were able to reveal qualitative overview of energy landscapes on the UCCSD surface by the initial search at the UB3LYP/cc-pVDZ level. The present treatment, i.e., the ADD-following via an unrestricted hybrid-DFT with a series of the MO stability check, will be used for an initial automated exploration of whole area on a PES including such radical fragment regions. 4. Conclusions The PES of H2CO were explored by the ADD-following method coupled with the UB3LYP method and the CASSCF method, and then geometries and IRC connections were refined at the UCCSD level. This is the first automated systematic global reaction route mapping on the PES of H2CO via unrestricted theories, where we introduced some modifications to explore on open-shell surfaces by the ADD-following method. Global topography on these three surfaces are qualitatively equivalent to each other. Although, in low energy regions of the PES, there are few changes in the topography on the PES with respect to those of restricted theories, new structures have been discovered in high energy regions involving radical fragments. It follows that the ADD-following via an unrestricted hybrid-DFT with a series of the MO stability check may provides new structures in such regions even for systems well studied via restricted theories. Acknowledgements S. Maeda is supported by Research Fellowship of the Japan Society for Promotion of Science for Young Scientists. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.cplett.2008.06.003. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26]
J.D. Goddard, H.F. Schaefer III, J. Chem. Phys. 70 (1979) 5117. G.E. Scuseria, H.F. Schaefer III, J. Chem. Phys. 90 (1989) 3629. L. Deng, T. Ziegler, L. Fan, J. Chem. Phys. 99 (1993) 3823. K. Bondensgård, F. Jensen, J. Chem. Phys. 104 (1996) 8025. H. Nakano, K. Nakayama, K. Hirao, M. Dupuis, J. Chem. Phys. 106 (1997) 4912. W. Quapp, M. Hirsch, O. Imig, D. Heidrich, J. Comput. Chem. 19 (1998) 1087. D. Feller, M. Dupuis, B.C. Garrett, J. Chem. Phys. 113 (2000) 218. K. Ohno, S. Maeda, Chem. Phys. Lett. 384 (2004) 277. S. Maeda, K. Ohno, J. Phys. Chem. A 109 (2005) 5742. K. Ohno, S. Maeda, J. Phys. Chem. A 110 (2006) 8933. L.B. Harding, S.J. Klippenstein, A.W. Jasper, Phys. Chem. Chem. Phys. 9 (2007) 4055. J. Troe, J. Phys. Chem. A 111 (2007) 3862. X. Yang, S. Maeda, K. Ohno, J. Phys. Chem. A 109 (2005) 7319. X. Yang, S. Maeda, K. Ohno, Chem. Phys. Lett. 418 (2006) 208. X. Yang, S. Maeda, K. Ohno, J. Phys. Chem. A 111 (2007) 5099. Y. Watanabe, S. Maeda, K. Ohno, Chem. Phys. Lett. 447 (2007) 21. S. Maeda, K. Ohno, J. Phys. Chem. A 111 (2007) 4527. Y. Luo, S. Maeda, K. Ohno, J. Phys. Chem. A 111 (2007) 10732. S. Maeda, K. Ohno, J. Phys. Chem. A 112 (2008) 2962. S. Maeda, K. Ohno, J. Phys. Chem. A 111 (2007) 13168. The execution module of the GRRM1.20 program is available by contact with the corresponding author (see http://qpcrkk.chem.tohoku.ac.jp/ for details). K. Fukui, Acc. Chem. Res. 14 (1981) 63. M.J. Frisch et al., GAUSSIAN 03, Revision , C.02. Gaussian, Inc., Wallingford CT, 2004. D. Townsend et al., Science, 306 (2004), p. 1158. F. Jensen, Introduction to Computational Chemistry, Wiley, Chichester, 1998. E.R. Johnson, G.A. DiLabio, Chem. Phys. Lett. 419 (2006) 333.