Surface Diffusion of a Zn Adatom on a Zn(001) Surface: A DFT Study

Aug 22, 2007 - Kei IOKIBE , Takeshi SAKAMOTO , Hiroto TACHIKAWA , Kazuhisa AZUMI. Journal of The Surface Finishing Society of Japan 2009 60 (9), 592- ...
0 downloads 0 Views 937KB Size
Download PDF
13510

J. Phys. Chem. C 2007, 111, 13510-13516

Surface Diffusion of a Zn Adatom on a Zn(001) Surface: A DFT Study Kei Iokibe, Kazuhisa Azumi, and Hiroto Tachikawa* DiVision of Materials Chemistry, Graduate School of Engineering, Hokkaido UniVersity, Sapporo 060-8628, Japan ReceiVed: April 22, 2007; In Final Form: June 27, 2007

The adsorption and diffusion of a Zn atom on a Zn(001) surface has been investigated theoretically by using first-principles periodic boundary condition calculations to elucidate the mechanism of crystal growth of zinc in a Zn surface. Three surface models, terrace, step, and kink, of the Zn(001) surface were considered as adsorption sites in the present study. The most stable site in the terrace model was the on-top site, where the zinc adatom is bound to only one zinc atom of the surface. The activation barrier between the on-top site and the next on-top site was negligibly low (the activation energy is 44 meV at the PW91/LANL2DZ level), suggesting that the Zn adatom diffuses easily on the Zn(001) surface. The Zn adatom was more stabilized at the step site (752 meV) and kink site (935 meV) with respect to the on-top site. It was found that the magnitude of 4s-4p orbital mixing of the Zn adatom (hybridization) is strongly related to the binding energy. The mechanism of the adsorption of Zn on Zn(001) surfaces was discussed on the basis of theoretical results.

1. Introduction In recent decades there has been considerable interest in the electrodeposition of zinc, since zinc plating efficiently protects steel from corrosion.1 More recently, zinc alloys and zinc oxide have also attracted considerable attention due to their possible applications in wide-band-gap semiconductor and UV devices.1-3 For example, zinc telluride (ZnTe), one of the zinc metal alloys, is a semiconductor of the II-VI family which is utilized in a variety of optoelectronic devices such as pure-green-lightemitting devices, terahertz detectors, solar cells, waveguides, and modulators.4 Also, zinc oxide (ZnO) is used in the window layers of solar cells and ultraviolet-light-emitting devices because of its wide direct band gap of 3.37 eV at room temperature.5 These materials are synthesized mainly by the epitaxial crystal growth technique and have been extensively studied over the past few decades. However, the mechanism of the epitaxial crystal growth is still in controversy. This is due to the fact that the mechanism of epitaxial growth of pure zinc metal is not clearly understood. Also, the electronic states of the zinc surface and adsorbed atom are unclear. In the present study, density functional theory (DFT) calculations under the periodic boundary condition (PBC) are applied to the surface diffusion of a Zn adatom on a Zn(001) surface to shed light on the mechanism of zinc crystal growth on the Zn(001) surface. Three surface models, terrace, step, and kink, are considered as diffusion regions of the Zn adatom. We focus our attention mainly on the electronic states of the Zn adatom in the surface models. In previous studies,6 we investigated the clustering mechanism of zinc clusters using ab initio and DFT methods and found that the binding energy of the Zn adatom is strongly related to the magnitude of 4s-4p hybridization: namely, the contribution of a 4p electron to the ground-state wave function dominates the binding nature of the Zn atom. In this work, we examine the validity of the 4s-4p hybridization * To whom correspondence should be addressed. E-mail: hiroto@ eng.hokudai.ac.jp. Fax: +81 11706-7897.

model to the Zn adatom on the Zn(001) surfaces as well as in the atomic cluster. To our knowledge, the adsorption of Zn on the Zn(001) surface has not been studied with the first-principles method. For the other systems, several theoretical works have been carried out mainly using the DFT calculation. Kim and Chung investigated the surface diffusion of a Co atom on an Al(001) surface using the DFT method.7 The energy barrier for the surface diffusion (migration of the Co adatom to an adjacent hollow site passing the bridge site) was calculated as 1.01 eV. Large displacement of neighboring Al atoms was accompanied by the surface diffusion of the Co adatom. Yarovsky and coworkers investigated the interactions of atomic sulfur with Fe(100) and Fe(110) surfaces.8 Vibrational frequency calculations were applied to determine the nature of stationary points at the adsorption sites. The binding energies of the S atom to atop, bridge, and hollow adsorption sites of the Fe(100) surface were calculated to be 4.11, 4.90, and -6.10 eV, respectively. From these works, it was shown that the DFT method is effective to determine the electronic states of an adatom (metal) adsorbed on a metal surface. In the present work, the DFT calculation together with the PBC is applied to the Zn adatom on the Zn(001) surface. 2. Method of Calculations Three surface models, terrace, step, and kink, were considered as diffusion regions of the Zn adatom on Zn(001) surfaces. The structure of the surfaces was constructed from the experimental value of the unit cell (a ) 2.659 Å and b ) 4.935 Å). Two and three layers were considered as the zinc surface in the present calculation. In the DFT calculations, the exchange component of Perdew and Wang’s 1991 functional theory (PW91) was used with the Los Alamos effective-core potential (ECP) basis set LANL2DZ9 under the PBC. First, the PW91/LANL2MB calculation was preliminarily carried out to obtain the potential energy surface of diffusion of the Zn adatom, and the height of the Zn adatom from the surface was roughly estimated. After that, the PW91/LANL2DZ

10.1021/jp073096t CCC: $37.00 © 2007 American Chemical Society Published on Web 08/22/2007

Diffusion of a Zn Adatom on a Zn(001) Surface

J. Phys. Chem. C, Vol. 111, No. 36, 2007 13511

Figure 1. Terrace model of the Zn(001) surface and potential energy surface. (a) Structure of nine units of the supercell composed of two layers (12 atoms in the first layer and 12 atoms in the second layer). (b) Top view of the terrace model. Bold and normal squares indicate the supercell and plot area used in the calculation, respectively. (c) Contour map of the potential energy surface relevant for the diffusion of the Zn atom on the terrace calculated at the PW91/LANL2DZ level with the PBC. The contour is drawn every 0.02 eV.

calculation was applied to obtain a more accurate potential energy surface. It should be noted that the DFT calculation at the PW91/LANL2DZ level of theory gives reasonable structures and electronic states for zinc clusters ZnN and ZnN-H2O (N ) 1-32), as shown in our previous paper.6 All calculations were carried out using the Gaussian 03 program package10 at the Information Initiative Center at Hokkaido University. 3. Results A. Diffusion of the Zn Adatom on the Zn(001) Surface (Terrace Model). We consider the adsorption of a single Zn atom on a pure Zn(001) surface (terrace model) as the first step. The potential energy surface relevant for diffusion of the Zn atom on the Zn(001) surface is calculated at the PW91/ LANL2DZ level under the PBC. The result is plotted as a

contour map in Figure 1. Figure 1a shows nine units of the supercell used in the PBC calculation. The supercell is composed of 12 atoms in the first layer and 12 atoms in the second layer. The height of the Zn atom is fixed to 2.64 Å in the calculation. The plot area is also given by a square where the contour map consists of 948 calculated points. From the DFT-PBC calculation, it is found that two stable sites of the Zn adatom are located on the terrace of the surface. These are the on-top site and 3-fold site (denoted by “a” and “b”, respectively). In the on-top site, the Zn atom is bound to one Zn atom on the Zn(001) surface, while the Zn atom in the 3-fold site is bound equivalently to three Zn atoms in the first layer. In addition to these sites, the Zn(001) surface has a hollow site (denoted by “c”) where the adatom is connected to three Zn atoms in the first layer and one Zn atom in the second layer.

13512 J. Phys. Chem. C, Vol. 111, No. 36, 2007

Iokibe et al.

TABLE 1: Binding and Relative Energies (meV) of the Zn Adatom to Several Points on the Zn(001) Model Surfacesa binding energy (meV)

relative energy (meV)

Zn(001) + Zn a b c a′ b′ a′′

Nonoptimized Structureb 0.0 a′′′ 493 (2.64) 0.0 TS1 485 (2.64) 8 c′ 455 (2.64) 38 TS2 493 (2.64) 0.0 TS3 485 (2.64) 8 d 493 (2.64) 0.0 e

489 (2.84) 479 (2.84) 494 (2.84) 484 (2.84) 731 (2.84) 751 (2.84) 928 (2.74)

0.0 10 -5 5 -242 -262 -439

Zn(001) + Zn a b c a′

Optimized Structurec 0.0 b′ 529 (2.84) 0.0 a′′ 485 (2.64) 44 TS3 463 (2.78) 66 d 529 (2.84) 0.0 e

485 (2.64) 529 (2.84) 754 (2.59) 752 (2.82) 935 (2.76)

44 0.0

point

binding energy (meV)

relative energy (meV)

point

a The values are calculated at the PW91/LANL2DZ level under the PBC. b Values in parenthesis are distances (Å) of the Zn atom from the first layer used in the energy calculation. b Optimized distances (Å) are given in parentheses.

The binding and relative energies of the sites are summarized in Table 1. The zero level of the binding energy corresponds to a dissociation limit of the Zn adatom: namely, Zn(001) + Zn. All sites (a, b, and c) have negative values, indicating that the Zn adatom adsorbed on the Zn surface is more stable in energy than the dissociation limit. The binding energies of the Zn adatom to the on-top, 3-fold, and hollow sites are 493, 485, and 455 meV, respectively, indicating that the on-top site is the most stable point on the Zn(001) surface, although the energy differences are negligibly small. Also, it is found that the hollow site is a transition state (TS) between the on-top and 3-fold sites. To obtain more realistic energetics, the heights of the Zn adatom are optimized in the points. The binding and relative energies corrected by the geometry optimization are given in Table 1 (lower columns). The binding energies of the Zn adatom adsorbed on the on-top, 3-fold, and hollow sites are slightly changed: 529, 485, and 463 meV, respectively. There are two different diffusion paths of the Zn atom from one on-top site (a) to another on-top site (a′′). In path I, the trajectory moves along the path expressed by a f b f a′ f b′ f a′′, while the trajectory in path II moves as a f c f b′ f a′′. At the PW91/LANL2DZ level, the energies (in meV) are changed as 0.0 (a) f 44.0 (b) f 0.0 (a′) f 44.0 (b′) f 0.0 (a′′) in path I and as 0.0 (a) f 66.0 (c) f 44.0 (b′) f 0.0 (a′′) in path II. Thus, both paths have significantly low activation barriers. This result strongly indicates that the Zn adatom at room temperature can move freely on the terrace region on the Zn(001) surface with a low activation energy. B. Diffusion of the Zn Atom in the Step Site on the Zn(001) Surface (Step Model). Figure 2 shows the step model used in the DFT-PBC calculation. Nine units of the supercell are shown in Figure 2a. The supercell is composed of 12 atoms in the first layer and 9 atoms in the second layer, while the step site consists of 3 zinc atoms. The height of the Zn atom is fixed to 2.84 Å in the calculation. The potential energy surface relevant for the diffusion of Zn on the step site of the Zn(001) surface is given in Figure 2. The point a′′′ indicates the on-top site of the Zn(001) surface corresponding to point a in Figure 1. Point c′ is also the on-top site, but the stable position is slightly deviated from the standard position of the on-top site. The 3-fold site near the step (denoted by “×” in Figure 2) is not stable, but it is located on the slope connected between the

terrace and step sites. The most stable site in the step model is point d, whose stabilization energy is 262 meV relative to point a′′′ at the PW91/LANL2DZ level. The minimum energy path from the terrace site to the step site (a′′′ f d) is illustrated by solid lines in Figure 2. This path is expressed by point a′′′ (0.0) f TS1 (10) f c′ (-5) f TS2 (5) f TS3 (-242) f point d (-262), where each value in parentheses indicates the relative energy at each point (meV). After the geometry optimization, the binding energy at point d is changed to 263 meV. The energetics indicates that the Zn atom can move easily from the terrace region to the step region and the Zn atom is stabilized in the step region. The stabilization energy is 263 meV relative to the on-top site (a′′′). The activation energy between site d and the next site d via TS3 is calculated to be only 2 meV, which is negligibly low. These results indicate that the diffusion of the Zn atom takes place without an activation barrier along the edge of the step (the diffusion path is expressed as point d f TS3 f point d). C. Diffusion of the Zn Atom in the Kink Site on the Zn(001) Surface (Kink Model). A model of the kink site used in the present study is illustrated in Figure 3. Nine units of the supercell used in the PBC calculation are illustrated in Figure 3a. The supercell is composed of 13 atoms in the first layer, 10 atoms in the second layer, and 4 atoms in the kink site: namely, a model of the kink site is composed of three layers. The potential energy surface for the diffusion of the Zn adatom around the kink site of the Zn(001) surface is given in Figure 3c. In the calculation, the height of the Zn adatom is fixed to 2.74 Å from the first layer. Point d indicates the stable point of the step corresponding to point d in Figure 2. Point e is the most stable point on the kink region. The stabilization energy of point e is 439 meV relative to point a in the terrace. This energy is modified to be 446 meV by the optimization of the height of the Zn adatom. The optimized height is obtained to be 2.76 Å at point e. The potential energy curve from point d to point e decreases monotonically without an activation barrier. The dashed curve indicates the minimum energy path from the kink site to the terrace site. The large energy needs to escape from the kink site to the terrace region because the kink site is strongly bound in energy. These results suggest that the diffusion of the Zn atom from the step site to the kink site easily occurs without an activation barrier. Once the Zn adatom binds to the kink site, the escape of the Zn from the kink site is difficult by thermal activation. D. Nature of Binding of the Zn Atom to the Zn(001) Surface. In a previous study,6 we investigated the clustering mechanism of the zinc cluster ZnN (N ) 2-32) using DFT methods and showed that the binding energy is strongly dependent on the magnitude of 4s-4p orbital mixing of the zinc atom (hybridization). This value has a close relationship to the magnitude of the contribution of the 4s f 4p excited state to the ground-state wave function, which is expressed by (4s) + δ(4p), where δ means the magnitude of the contribution of the 4p electron. The increase of the population of the 4p electron in the adatom (N4p) causes an increase of the binding energy. In the case of the isolated Zn atom, the contribution of the 4p orbital is zero because the ground-state configuration of the Zn atom is 4s24p0. In the larger clusters, the Zn atom in the surface region has N4p ) 0.1, while the internal Zn atom has a large value (N4p ) 0.6-0.8).6 In this section, the investigation of the nature of the binding of the Zn atom to the Zn(001) surface using the Zn cluster model of the Zn(001) surface is described. The cluster model used for the step model is illustrated in Figure 4A. This model is

Diffusion of a Zn Adatom on a Zn(001) Surface

J. Phys. Chem. C, Vol. 111, No. 36, 2007 13513

Figure 2. Step model of the Zn(001) surface and potential energy surface. (a) Structure of nine units of the supercell which is composed of three layers (12 atoms in the first layer, 9 atoms in the second layer, and 3 atoms in step region). (b) Top view of the terrace model. (c) Contour map of the potential energy surface relevant for the diffusion of the Zn atom around the step site. The contour is drawn every 0.02 eV.

composed of 20 zinc atoms. The height of the Zn adatom is fixed to 2.4675 Å during the calculation. The potential energy curve shows that the energy decreases gradually as the step is approached. The energy is stabilized at the step region (denoted by c) as well as the PBC calculation. The natural electron population in the 4p orbital of the adatom (N4p) is calculated to be 0.1 at the terrace region (point a), and it increases gradually as the step is approached (a f b f c). The value is N4p ) 0.2 at the step site (point c). These results indicate that the contribution of the 4p electron makes the binding energy of the Zn adatom to the Zn(001) surface. The magnitudes of the contribution of the 4p electron at the terrace and step sites are calculated to be N4p ) 0.1 and 0.2, respectively. The variation of the natural electron population (N4p) along the edge of the step site is given in Figure 4B. The cluster model

of the kink site is composed of 33 zinc atoms. The height of the Zn adatom is fixed to 2.4675 Å during the calculation. The values at the step site and kink site are calculated to be N4p ) 0.21 and 0.30, respectively. Thus, N4p increases drastically as N4p ) 0.0 (gas phase), 0.1 (surface terrace region), 0.2 (step site), and 0.3 (kink site). It can therefore be concluded that the binding nature of the adatom is strongly dependent on the magnitudes of the contribution of the 4p electron. 4. Discussion A. Summary of the Present Calculations. In the present study, the potential energy surfaces relevant for diffusion of the Zn adatom on the terrace, step, and kink regions of the Zn(001) surface were calculated at the PW91/LANL2DZ level under the PBC to elucidate the diffusion mechanism of Zn on

13514 J. Phys. Chem. C, Vol. 111, No. 36, 2007

Iokibe et al.

Figure 3. Kink model of the Zn(001) surface and potential energy surface. (a) Structure of nine units of the supercell which is composed of three layers (13 atoms in the first layer, 10 atoms in the second layer, and 4 atoms in step region). (b) Top view of the terrace model. (c) Contour map of the potential energy surface relevant for the diffusion of the Zn atom around the step site. The contour is drawn every 0.02 eV.

Zn(001). Two stable sites were found on the terrace region. In the on-top site, the Zn atom is bound to one zinc atom of the Zn(001) surface. The Zn adatom in the 3-fold site is equivalently bound to three zinc atoms. The Zn adatom can diffuse on the Zn(001) surface with a very low activation energy. The Zn adatom is more stabilized in step and kink sites. B. Mechanism of the Epitaxial Crystal Growth of the Zinc Surface. On the basis of the present calculations, the model of epitaxial crystal growth of the zinc surface is proposed in this section. A schematic illustration of the model is given in Figure 5. The zinc atom has a 4s2 electronic configuration, and the contribution of the 4p orbital is zero (N4p ) 0.0). When the zinc atom binds to the Zn(001) surface (point a), the 4p state is mixed with the 4s2 state (N4p ) 0.1). The adatom diffuses freely on the surface with a very low activation barrier and reaches the step region (point b). The contribution of the 4p electron in the Zn adatom is calculated to be N4p ) 0.2 at the step site. The energy of the system becomes lower (the increase of the binding energy from the terrace region (529 meV) to the step site (752 meV) is calculated to be 223 meV at the PW91/ LANL2DZ level). The increase of the binding energy is caused

by the increase of the contribution of the 4p orbital. The adatom can diffuse easily along the edge of the step and enters the kink site (point c). The energy is further stabilized (935 meV) and relaxed in this region (N4p ) 0.3). Once the zinc atom is bound to the Zn atom in the kink site, the dissociation of the Zn atom does not occur due to the large stabilization energy. Thus, the change of the contribution of the 4p electron is important in the crystal growth of the zinc atom. C. Comparison with Previous Studies. In a previous study,6 we investigated the clustering mechanism of the zinc cluster ZnN (N ) 2-32) using DFT methods to elucidate the change from van der Waals to metallic properties in the Zn system. We showed that the binding energy is strongly dependent on the magnitude of 4s-4p orbital mixing in the adatom (4s f 4p electronic excitation). This value is strongly related to the magnitude of the contribution of the 4s f 4p excited state to the ground-state wave function, which is expressed by (4s) + δ(4p). In the case of the isolated Zn atom, the contribution of the 4p orbital is zero because the ground-state configuration of the Zn atom is 4s24p0. In the larger clusters, the Zn atom in the

Diffusion of a Zn Adatom on a Zn(001) Surface

J. Phys. Chem. C, Vol. 111, No. 36, 2007 13515

Figure 4. Potential energy curves and natural populations in the 4p orbital (N4p) of the adatom plotted along the diffusion paths of the Zn adatom (A, left) from the terrace site to the step site and (B, right) from the step site to the kink site. The values are calculated at the PW91/LANL2DZ level. The height of Zn was fixed to 2.4675 Å in the calculations. The calculation was carried out at the PW91/LANL2DZ level without the PBC. The lower panels indicate the step model (A) and kink model (B) using the calculations.

Figure 5. Model of the diffusion of the Zn adatom on the Zn(001) surface and epitaxial growth of the Zn surface.

surface region has N4p ) 0.1, while the internal Zn atom has a large value (N4p ) 0.6-0.8). The present calculations showed that the values of the Zn adatom on the terrace, step, and kink regions are drastically changed to 0.1, 0.2, and 0.3, respectively. Also, it was found that the binding energy increases with increasing magnitude of the contribution of the 4p orbital. This feature is the same in both the cluster and the surface. Therefore, it can be concluded that the electronic properties of zinc metal are dependent on the magnitude of orbital mixing from 4p to 4s orbitals (i.e., magnitude of hybridization). D. Final Remarks. In the present study, we examined only two and three layers of the zinc surface to the diffuse region of the Zn atom. Although this model is the smallest of the Zn(001) surface, a reasonable feature could be obtained for the diffusion process of the Zn adatom. To discuss the diffusion process

quantitatively, larger sized surface models and more flexible basis sets would be required as a future work. Such calculation will be possible after the development of a high-speed CPU machine. The dynamics effects on the diffusion of the Zn adatom were not considered in this study. In the near future, a direct ab initio or semiempirical MO-molecular dynamics (MD) technique will be applied to the diffusion process of Zn on the Zn(001) surface as well as in our previous works.11-14 Despite the several assumptions introduced here, the results enable us to obtain valuable information on the mechanism of the diffusion of the zinc atom on Zn(001) surfaces. Acknowledgment. We are indebted to the Computer Center at the Institute for Molecular Science (IMS) for the use of the computing facilities. H.T. also acknowledges partial support

13516 J. Phys. Chem. C, Vol. 111, No. 36, 2007 from a Grant-in-Aid for Scientific Research (C) from the Japan Society for the Promotion of Science (JSPS). References and Notes (1) Pereira, M. S.; Barbosa, L. L.; Souza, C. A. C.; De Moraes, A. C. M.; Carlos, I. A. J. Appl. Electrochem. 2006, 36, 727. (2) Look, D. C. J. Electron. Mater. 2006, 35, 1299. (3) Ankiewicz, A. O.; Carmo, M. C.; Sobolev, N. A.; Gehlhoff, W.; Kaidashev, E. M.; Rahm, A.; Lorenz, M.; Grundmann, M. J. Appl. Phys. 2007, 101, Article No. 024324. (4) Kume, Y.; Guo, Q.; Tanaka, T.; Nishio, M.; Ogawa, H.; Shen, W. J. Cryst. Growth 2007, 298, 441. (5) Terasako, T.; Ishiko, Y.; Saeki, K.; Yudate, S.; Shirakata, S. 2007, 298, 481. (6) (a) Iokibe, K.; Tachikawa, H.; Azumi, K. J. Phys. B 2007, 40, 427. (b) Tachikawa, H.; Iokibe, K.; Azumi, K.; Kawabata, H. Phys. Chem. Chem. Phys. 2007, 9, 3978. (7) Kim, C.; Chung, Y.-C. J. Solid State Chem. 2005, 178, 47. (8) Todorova, N.; Spencer, M. J. S.; Yarovsky, I. Surf. Sci. 2007, 601, 665. (9) (a) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 270. (b) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 284. (c) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 299.

Iokibe et al. (10) Ab initio MO program Gaussian 03, revision B.04: Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A., Gaussian, Inc., Pittsburgh, PA, 2003. (11) Tachikawa, H.; Shimizu, A. J. Phys. Chem. B 2006, 110, 20445. (12) Tachikawa, H. J. Chem. Phys. 2006, 125, 133119. (13) Tachikawa, H. J. Chem. Phys. 2006, 125, 144307. (14) Tachikawa, H.; Abe, S. J. Chem. Phys. 2007, 126, 294310.