- Email: [email protected]
First-principles simulations for structures and optical spectra of carbon cluster C8
Chemical Physics Letters 380 (2003) 330–336 www.elsevier.com/locate/cplett
First-principles simulations for structures and optical spectra of carbon cluster C8 Xiang-Rong Chen b
a,*
, Yu-Lin Bai
a,b
, Xiao-Lin Zhou
a,c
, Xiang-Dong Yang
a,b
a Institute of Atomic and Molecular Physics, Sichuan University, Chengdu, Sichuan 610065, PR China Department of Electronic Information Science and Technology, Yibi University, Yibin 644000, PR China c Department of Physics, Sichuan Normal University, Chengdu 610066, PR China
Received 7 July 2003; in final form 5 September 2003 Published online: 28 September 2003
Abstract We apply a finite-difference pseudopotential density functional theory in real space and the Langevin molecular dynamics annealing technique as well as the adiabatic time-dependent density functional theory within the local density approximation (TDLDA) to the descriptions of structures and optical absorption spectra of carbon clusters C8 . It is shown that the odd-numbered clusters have the linear structures and most of the even-numbered clusters prefer cyclic structures. The calculated spectra exhibit a variety of features that can be used for comparison against future experimental investigations. Ó 2003 Elsevier B.V. All rights reserved.
1. Introduction Carbon is a unique element, which has the ability to form strong chemical bonds with a variety of coordination numbers from two (e.g., linear chains or carbine phase), to three (e.g., graphite) and four (e.g., diamond). Although there were extensive studies before, many interesting problems including the high-temperature and high-pressure phase diagram of carbon and the geometric and electronic structure of various disordered carbon phases remain unresolved. Mean-
*
Corresponding author. Fax: +862885405515. E-mail address: [email protected] (X.-R. Chen).
while, the carbon atomic aggregates (i.e., carbon clusters) have also been the subjects of extensive theoretical and experimental studies, even before the development of fullerene chemistry [1,2]. The interplay between the electronic degrees of freedom and the geometrical structure in small carbon clusters has attracted a lot of interest [3]. Over the past years, the structures and properties of small clusters have been the focus of using a variety of theoretical calculations mainly based on ab initio techniques [4,5], empirical interatomic potentials [6], density functional theory with the local density approximation [7], embedded atom methods (EAM) [8], tight-binding molecular dynamics (TBMD) [9,10]. Most of the larger carbon clusters were found to have either linear or
0009-2614/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2003.09.032
X.-R. Chen et al. / Chemical Physics Letters 380 (2003) 330–336
monocyclic ground states involving multiple bonding, as is clearly due to the presence of strong p-bonding involving carbon. In general, the oddnumbered carbon clusters have linear structures and most of the even-numbered clusters prefer cyclic structures [11,12]. In this Letter, we focus on the carbon cluster C8 , the linear structure of which has hitherto been undetected. Raghavachari and Binkley [11] were the first to find a ring with C4h symmetry using allelectron ab initio molecular orbital techniques [13] as the global minimum on the Hartree–Fock potential surface. The multi-reference configuration €f interaction (MRCI) study by Parasuk and AlmlO [14] has also found the same structure. However, Hutter et al. [2] have predicted the C8 ring to have only C2h symmetry. We here have applied a finite-difference pseudopotential density functional theory method in real space and the Langevin molecular dynamics annealing technique [15] to determine the structures of cluster C8 in the ground state, and applied the adiabatic time-dependent density functional theory (TDDFT) within the local density approximation (LDA) [16] to calculate its optical spectrum. This work has not been performed before. In the next section, we make a brief description of the computational methods. The results and discussion are presented in Section 3.
2. Computational method 2.1. Pseudopotential density functional theory method In the ab initio pseudopotentials constructed within the LDA, the effect of the core electrons and nuclei is replaced with an effective ionic potential. The resulting potential and wave functions obtained using this procedure are smoothly varying. Therefore, simple bases, such as plane waves or real space grid methods, can be directly applied. The pseudopotential method has been applied to a variety of semiconductor clusters, resulting in good agreement with experimental results [17–19]. While pseudopotentials provide an accurate and efficient method for computing structural
331
energies, the issue of selecting energetically viable structures remains. For relatively small clusters, one can make an inventory of all possible structures and investigate each one. Thus, the ÔtrueÕ ground state can be obtained. Unfortunately, as the cluster size exceeds a few atoms, the number of competing structures becomes quite large. The computational burden imposed by the construction of such an inventory and the subsequent calculations excludes this explicit method. For the larger clusters, it is useful to employ computer simulations that exploit a simulated annealing strategy. In our work, we optimize the cluster structures using the Langevin molecular dynamics (LMD) coupled with the simulated annealing procedure. The cluster is heated to a high temperature and then gradually cooled. Initially, the cluster samples numerous configurations without regard to the corresponding structural energy. As the cluster is cooled, only energetically favorable structures are sampled. We control the temperature of the system without rescaling the velocities, as is often done in Newtonian molecular dynamics. Energy can flow into and out of the system as required by the temperature of the heat bath. In this approach, the atoms are collectively moved in order to sample the configuration space. The minima in a potential energy surface are located based on the interatomic forces. Conversely, the presence of random forces helps the system escape from metastable states. The forces can be determined from classical interatomic potentials or quantum mechanically from ab initio pseudopotentials. In this work, the exchange-correlation potential term is approximated by the Ceperley–Alder functional [20,21]. The Troullier–Martins nonlocal pseudopotentials [22] in Kleinman–Bylander form [23] are used. The cutoff core radii adopted are 1.50, 1.54, and 1.54 a.u. for C 2s, 2p, and 3d electronic states, respectively [24–26], in constructing the Troullier–Martins scheme pseudopotentials. With respect to computing the quantum forces, the cluster in question is placed in a spherical domain. Outside of this domain, the wave function is required to vanish. The radius of the sphere is such that the outmost atom is 7 a.u. from the boundary. The grid spacing was 0.5 a.u.,
332
X.-R. Chen et al. / Chemical Physics Letters 380 (2003) 330–336
which is evaluated roughly with a plane wave 2 cutoff of ðp=hÞ or about 30 Ry. For our study of the structure of the C8 cluster, the initial simulation temperature is taken to be 4100 K, and the final temperature 100 K. The annealing schedule lowers the temperature 500 K each 50 time steps. The time step is taken to be 80 a.u. (1 a.u. ¼ 0.0241888428 fs). A friction coefficient is taken to be 0.0006 a.u. After reaching the temperature of 100 K, the cluster is quenched to 0 K. The ground state structure of this cluster is found through a direct minimization by a steepest descent procedure. It is very important to choose the initial atomic configuration. If the atoms are too far apart, they may not form a stable cluster as the simulations proceed. If the atoms are too close together, they may form a metastable cluster from which the ground state may be kinetically inaccessible. Typically, a random placement of the atoms must reside within 1.05 and 1.30 times the dimer bond length. If the initial atomic configuration cannot form a stable cluster, there are two methods to cope with it. One method is that, you may check the energy per atom of the first step to check that it nearly equals (no more or less than 4 eV) that of the dimer. When the energy per atom is much greater than that of the dimer, you must reform the initial configuration. The other method is that, the time step may be smaller within the first fifty steps. The detailed account of this method has been reported elsewhere [15].
approximate owing to the use of the LDA instead of the exact spatial dependence, and an adiabatic approximation instead of the exact time dependence of the exchange-correlation functional. The time-dependent LDA (TDLDA) technique requires a substantially smaller computational effort and is able to handle much larger clusters. This technique has recently been applied to Na, Si, GaAs, and CdSe clusters and was found to yield results in very good agreement with experiment [27–29]. We here focus on applying the TDLDA technique to the calculation of the optical absorption spectra of the C8 cluster. We apply it within a frequency-space approach, in which the electronic transition energies Xn are obtained from the solution of the following eigenvalue problem: h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii pffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2ijr dik djl drs þ 2 fijr xijr Kijr;kls fkls xklr Fn
2.2. Time-dependent density functional theory within the local density approximation
where /ðrÞ is a Kohn–Sham wave function, i; j; r; k; l; s are the occupied states, unoccupied state, and spin indices, respectively, qðrÞ is the Kohn– Sham charge density, and mxc ðrÞ is the LDA exchange-correlation potential. This computational technique is based on the higher-order finite difference pseudopotential method [30]. The solution proceeds in three parts: First, the Kohn–Sham energies and wave functions are found by means of a conventional ground state calculation. Second, the wave functions are used to construct the coupling matrix of Eq. (2). Third, the TDLDA eigenvalue equation, Eq. (1), is solved and the absorption spectrum constructed. All three steps were performed using a real-space grid confined in a spherical domain [16].
It is known that the density functional theory states are uniquely determined by the ground state charge density. However, the classic formulation of the density functional formalism is restricted to the time-independent case only. A proper treatment of electronic excitations is not possible within the time-independent framework and requires a generalization of DFT to the time-dependent phenomena. This limitation has led to the development of TDDFT. Within the TDDFT, the main theorem of the density functional formalism is extended to time-dependent systems. Formally, this is an exact approach, whereas in practice it is
¼ X2n Fn ;
ð1Þ
where xijr ¼ ejr eir are the Kohn–Sham transition energies, and fijr ¼ nir njr are the differences between the occupation numbers of the ith and jth states. The eigenvectors Fn are related to the oscillator strength, and Kijr;kls is a coupling matrix given by Z Z 1 dmxc r ðrÞ K¼ /ir ðrÞ /jr ðrÞ þ jr r0 j dqs ðr0 Þ /ks ðr0 Þ /ls ðr0 Þ dr dr0 ;
ð2Þ
X.-R. Chen et al. / Chemical Physics Letters 380 (2003) 330–336
3. Results and discussion We first determine the ground state structure of the C8 cluster by simulated annealing via the Langevin molecular dynamics technique. We illustrate the simulated annealing procedure in Figs. 1 and 2. The initial geometry is a pair of parallel
333
rhombus containing several bonds. After 169 steps (327.03 fs), it forms an octagon structure (see Fig. 2(r)), but the structure is still somewhat removed from the ground state. After 368 steps (711.15 fs), the ground state structure of cyclic D2h (see Fig. 2(t)) is essentially formed. The result is different from that obtained using all-electron ab
-6.0
0.6 Ek (eV)
Binding Energy (eV)
0.8 -6.5 -7.0
0.4
-7.5 0.2 -8.0
0.0 0
100
200
300
400
500
0
100
200
300
400
500
Fig. 1. Binding energy and kinetic energy of the C8 cluster during a Langevin simulation. The initial temperature is 4100 K and the final temperature is 100 K. The time step is 80 a.u.
Fig. 2. The geometries of the C8 cluster at various time steps during a Langevin simulation.
334
X.-R. Chen et al. / Chemical Physics Letters 380 (2003) 330–336
initio molecular orbital techniques [13] by Raghavachari and Binkley [11] and different from the multi-reference configuration interaction (MRCI) € f [14], and also difstudy by Parasuk and AlmlO ferent from that predicted by Hutter et al. [2] We also calculate the ground state structures of carbon clusters Cn ðn ¼ 2–7Þ. Our results reveal several interesting aspects. The bond length of the (the experimental value is 1.2429 C2 dimer is 1.27 A [3]) and the corresponding binding energy is 5.44 A eV (the experimental value is about 6.2 eV [11]). The C3 cluster is an isosceles triangle with C2v symmetry. The corresponding bond length is 1.24 (the experimental value is about 1.30 A [3]), apex A angle is h ¼ 137:97°, and the binding energy is 14.68 eV (the experimental value is 13.9 eV [11]). For the C4 cluster, the structure is a rhombus with , the minor diagonal length the side length 1.44 A and the binding energy 22.65 eV (the ex1.42 A perimental value is 19.4 eV [31]). The configuration of the C5 cluster is a linear structure. The four bond , respeclengths are 1.23, 1.25, 1.28, and 1.32 A tively. This structure is a little difference from that of the C5 cluster in [11]. The ground state structure of the C6 cluster is a monocyclic structure, as is different from the linear structure, but agrees well with that obtained through ab initio calculations of Raghavachari and Binkley [11]. The monocyclic structure is found to be about 0.2 eV lower than the monocyclic structure in energy. As regards the structure of the C7 cluster, it is also a linear structure, just like the structure of the C5 cluster [11]. In summary, our calculated results agree well with those obtained by Raghavachari and Binkley [11]. Xu et al. [32] have applied the tight-binding model to study these small carbon clusters. They have obtained very regular structures for these small carbon clusters. We do not agree with them.
Therefore, we believe that our scheme has made a successful prediction of the small carbon clusters. Moreover, all the results show that there is an odd– even alternation in the small carbon cluster geometries: the odd-numbered clusters have linear structures and many of the even-numbered clusters prefer cyclic structures. In Table 1, we list the average bond lengths together with our calculated total energies of these carbon clusters. It is obvious from our calculated , of the C4 results that the bond length, 1.44 A cluster structure is the longest among these small carbon clusters; the even-numbered small carbon clusters have longer bond lengths and the oddnumbered clusters prefer slightly shorter bond lengths. However, the absolute value of the total energies per atom increases roughly with the number of atoms of the cluster. For the TDLDA calculations of the C8 cluster with D2h symmetry as illustrated in Fig. 2(t), we carefully have tested for convergence of the computed excitation energies and associated absorption spectra with respect to the size of the spherical domain, the grid spacing, and the total number of electronic states included in the calculations. Convergence is obtained with a domain radius of 33 a.u., a grid spacing of 0.6 a.u., and a number of unoccupied Kohn–Sham orbitals that is at least three times greater than the number of occupied ones. The calculated TDLDA absorption spectrum (solid line) together with the absorption spectrum from time-independent LDA (dashed line) are shown in Fig. 3. All spectra are broadened using a Gaussian in order to simulate the effect of a finite temperature as well as the finite resolution of an experiment. We use a Gaussian width of 0.1 eV, which is a typical accuracy for optical absorption measurements.
Table 1 ) and calculated total energies ET (eV/atom) of Cn ðn ¼ 1–8Þ The average bond lengths L (A n a
L Lb ET a
2
3
1.27 1.2429 )150.90
1.24 1.2968 )152.80
Our work. b Average experimental data from [3].
4 1.44 )153.02
5
6
1.27 1.2855 )153.65
1.30 1.2780 )153.58
7
8 1.28
1.31
)153.94
)154.21
Absorption cross section
X.-R. Chen et al. / Chemical Physics Letters 380 (2003) 330–336
335
ties. The TDLDA computed spectrum, which exhibit a variety of features that could be used for cluster identification, differ significantly from those computed using a simple LDA approach. It should allow for a comparison with future experimental investigations.
2.0 1.5 1.0 0.5 0.0
0
5
10
15
20
Fig. 3. The calculated optical absorption spectra (arbitrary units) for the C8 cluster with D2h symmetry, where the solid and dashed lines correspond to the case of TDLDA and LDA, respectively.
It is noted that the calculated TDLDA absorption spectrum are characterized by three extremely intense absorption peaks near 7.6, 11.3, and 17.5 eV, respectively, and display a significantly blue shift with respect to the Kohn–Sham eigenvalue spectrum, which are consistent with the well-known band gap underestimate of LDA. The first TDLDA absorption peak display a blue shift by about 3.0 eV with respect to the Kohn–Sham eigenvalue spectrum, the second display a blue shift by about 1.8 eV and the third by about 0.6 eV. This implies that TDLDA works well for systems with a significant degree of ionicity. There is a qualitative difference between the current results and those of Gan Asn [33] and Si2n [27] clusters. Similarity to III–V and IV–IV clusters is that the long absorption tails have been found in the calculated absorption spectra in the region of lower transition energies. As we know that, these tails, which extend deeply into the region of lower transition energies, are likely due to the existence of free surfaces in the clusters. Such behavior is absent in quantum dots that are passivated at the boundaries [34]. In summary, the ground state structure for the C8 cluster is calculated using the finite-difference pseudopotential density functional theory method in real space and the Langevin molecular dynamics annealing technique. We have found that the C8 ring has D2h symmetry and the TDLDA formalism can provide an efficient alternative to more complex theoretical methods for excited state proper-
Acknowledgements The authors thank Prof. J.R. Chelikowsky for providing us PARSEC_CODE, a serial version of real space, time-dependent density functional code. We also acknowledge support for this work by the National Natural Science Foundation of PR China under Grant No.10274055 and by the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry.
References [1] W. Weltner Jr., R.J. Van Zee, Chem. Rev. 89 (1989) 1713. [2] J. Hutter, H.P. Luthi, F. Diederich, J. Am. Chem. Soc. 116 (1994) 750. [3] J.M.L. Martin, P.R. Taylor, J. Phys. Chem. 100 (1996) 6047. [4] G. Galli, R.M. Martin, R. Car, M. Parrinello, Phys. Rev. Lett. 62 (1989) 555, ibid, 63 (1989) 988. [5] Z.Y. Lu, C.Z. Wang, K.M. Ho, Phys. Rev. B 61 (2000) 2329. [6] C.Z. Wang, C.T. Chan, K.M. Ho, Phys. Rev. B 42 (1990) 11276. [7] S. Fahy, S.G. Louie, Phys. Rev. B 36 (1987) 3373. [8] V.G. Grigoryan, M. Springborg, Phys. Chem. Chem. Phys. 3 (2001) 5135. [9] J.L. Wang et al., Solid State Commun. 117 (2001) 593. [10] J.J. Zhao, J.L. Wang, G.H. Wang, Phys. Lett. A 275 (1993) 281. [11] K. Raghavachari, J.S. Binkley, J. Chem. Phys. 87 (1987) 2191. [12] Y.L. Bai, X.R. Chen, X.D. Yang, P.F. Lu, Commun. Theor. Phys., to be published. [13] W.J. Hehre, L. Radom, P.V.R. Schleyer, J.A. Pople, Ab initio Molecular Orbital Theory, Wiley, New York, 1979. € f, Theor. Chim. Acta 83 (1992) 227. [14] V. Parasuk, J. AlmlO [15] C. Troparevsky, J.R. Chelikowsky, J. Chem. Phys. 114 (2001) 943, and references therein. € gu €t, J.R. Chelikowsky, Phys. Rev. B 65 [16] I. Vasiliev, S. O (2002) 115416, and references therein. [17] N. Binggeli, J.R. Chelikowsky, Phys. Rev. B 50 (1994) 11764. [18] N. Binggeli, J.R. Chelikowsky, Appl. Phys. Lett. 66 (1995) 1316. [19] J.R. Chelikowsky, J. Phys. D 33 (2000) R33.
336 [20] [21] [22] [23] [24]
[25] [26] [27] [28]
X.-R. Chen et al. / Chemical Physics Letters 380 (2003) 330–336 D.M. Ceperley, B.J. Alder, Phys. Rev. Lett. 45 (1980) 566. J.P. Perdew, A. Zunger, Phys. Rev. B 23 (1981) 5048. N. Troullier, J.L. Martins, Phys. Rev. B 43 (1991) 1993. L. Kleinman, D.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. X.R. Chen, A. Oshiyama, S. Okada, Phys. Rev. B 67 (2003) 033408; Chem. Phys. Lett. 371 (2003) 528. X.R. Chen, A. Oshiyama, S. Okada, Q.Q. Gou, Chin. Phys. Lett. 20 (2003) 404. X.R. Chen, X.L. Zhou, J. Zhu, Q.Q. Gou, Phys. Lett. A 315 (2003) 403. € €t, J.R. Chelikowsky, Phys. Rev. Lett. 82 I. Vasiliev, S. O gu (1999) 1919, and references therein. € € t, J.R. Chelikowsky, in: P. Jena, S.N. I. Vasiliev, S. O gu Khanna, B.K. Rao (Eds.), Cluster and Nanostructure Interfaces, World Scientific, Singapore, 2000, p. 259.
[29] M. Claudia Troparevsky, Leeor Kronik, James R. Chelikowsky, Phys. Rev. B 65 (2001) 033311. [30] (a) J.R. Chelikowsky, N. Troullier, K. Wu, Y. Saad, Phys. Rev. B 50 (1994) 11355; (b) J.R. Chelikowsky, N. Troullier, X. Jing, D. Dean, N. Binggeli, K. Wu, Y. Saad, Comput. Phys. Commun. 85 (1995) 325; (c) Y. Saad, A. Stathoupolos, J.R. Chelikowsky, K. Wu, S. € gu €t, BIT 36 (1996) 563. O [31] K. Raghavachari, R.A. Whiteside, J.A. Pople, J. Chem. Phys. 85 (1986) 6623. [32] C.H. Xu, C.Z. Wang, C.T. Chan, K.M. Ho, J. Phys.: Condens. Matter 4 (1992) 6047. € gu €t, J.R. Chelikowsky, Phys. Rev. B 60 [33] I. Vasiliev, S. O (1999) 8477. € gu € t, J.R. Chelikowsky, Phys. Rev. Lett. 86 [34] I. Vasiliev, S. O (2001) 1813.
First-principles simulations for structures and optical spectra of carbon cluster C8 Xiang-Rong Chen b
a,*
, Yu-Lin Bai
a,b
, Xiao-Lin Zhou
a,c
, Xiang-Dong Yang
a,b
a Institute of Atomic and Molecular Physics, Sichuan University, Chengdu, Sichuan 610065, PR China Department of Electronic Information Science and Technology, Yibi University, Yibin 644000, PR China c Department of Physics, Sichuan Normal University, Chengdu 610066, PR China
Received 7 July 2003; in final form 5 September 2003 Published online: 28 September 2003
Abstract We apply a finite-difference pseudopotential density functional theory in real space and the Langevin molecular dynamics annealing technique as well as the adiabatic time-dependent density functional theory within the local density approximation (TDLDA) to the descriptions of structures and optical absorption spectra of carbon clusters C8 . It is shown that the odd-numbered clusters have the linear structures and most of the even-numbered clusters prefer cyclic structures. The calculated spectra exhibit a variety of features that can be used for comparison against future experimental investigations. Ó 2003 Elsevier B.V. All rights reserved.
1. Introduction Carbon is a unique element, which has the ability to form strong chemical bonds with a variety of coordination numbers from two (e.g., linear chains or carbine phase), to three (e.g., graphite) and four (e.g., diamond). Although there were extensive studies before, many interesting problems including the high-temperature and high-pressure phase diagram of carbon and the geometric and electronic structure of various disordered carbon phases remain unresolved. Mean-
*
Corresponding author. Fax: +862885405515. E-mail address: [email protected] (X.-R. Chen).
while, the carbon atomic aggregates (i.e., carbon clusters) have also been the subjects of extensive theoretical and experimental studies, even before the development of fullerene chemistry [1,2]. The interplay between the electronic degrees of freedom and the geometrical structure in small carbon clusters has attracted a lot of interest [3]. Over the past years, the structures and properties of small clusters have been the focus of using a variety of theoretical calculations mainly based on ab initio techniques [4,5], empirical interatomic potentials [6], density functional theory with the local density approximation [7], embedded atom methods (EAM) [8], tight-binding molecular dynamics (TBMD) [9,10]. Most of the larger carbon clusters were found to have either linear or
0009-2614/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2003.09.032
X.-R. Chen et al. / Chemical Physics Letters 380 (2003) 330–336
monocyclic ground states involving multiple bonding, as is clearly due to the presence of strong p-bonding involving carbon. In general, the oddnumbered carbon clusters have linear structures and most of the even-numbered clusters prefer cyclic structures [11,12]. In this Letter, we focus on the carbon cluster C8 , the linear structure of which has hitherto been undetected. Raghavachari and Binkley [11] were the first to find a ring with C4h symmetry using allelectron ab initio molecular orbital techniques [13] as the global minimum on the Hartree–Fock potential surface. The multi-reference configuration €f interaction (MRCI) study by Parasuk and AlmlO [14] has also found the same structure. However, Hutter et al. [2] have predicted the C8 ring to have only C2h symmetry. We here have applied a finite-difference pseudopotential density functional theory method in real space and the Langevin molecular dynamics annealing technique [15] to determine the structures of cluster C8 in the ground state, and applied the adiabatic time-dependent density functional theory (TDDFT) within the local density approximation (LDA) [16] to calculate its optical spectrum. This work has not been performed before. In the next section, we make a brief description of the computational methods. The results and discussion are presented in Section 3.
2. Computational method 2.1. Pseudopotential density functional theory method In the ab initio pseudopotentials constructed within the LDA, the effect of the core electrons and nuclei is replaced with an effective ionic potential. The resulting potential and wave functions obtained using this procedure are smoothly varying. Therefore, simple bases, such as plane waves or real space grid methods, can be directly applied. The pseudopotential method has been applied to a variety of semiconductor clusters, resulting in good agreement with experimental results [17–19]. While pseudopotentials provide an accurate and efficient method for computing structural
331
energies, the issue of selecting energetically viable structures remains. For relatively small clusters, one can make an inventory of all possible structures and investigate each one. Thus, the ÔtrueÕ ground state can be obtained. Unfortunately, as the cluster size exceeds a few atoms, the number of competing structures becomes quite large. The computational burden imposed by the construction of such an inventory and the subsequent calculations excludes this explicit method. For the larger clusters, it is useful to employ computer simulations that exploit a simulated annealing strategy. In our work, we optimize the cluster structures using the Langevin molecular dynamics (LMD) coupled with the simulated annealing procedure. The cluster is heated to a high temperature and then gradually cooled. Initially, the cluster samples numerous configurations without regard to the corresponding structural energy. As the cluster is cooled, only energetically favorable structures are sampled. We control the temperature of the system without rescaling the velocities, as is often done in Newtonian molecular dynamics. Energy can flow into and out of the system as required by the temperature of the heat bath. In this approach, the atoms are collectively moved in order to sample the configuration space. The minima in a potential energy surface are located based on the interatomic forces. Conversely, the presence of random forces helps the system escape from metastable states. The forces can be determined from classical interatomic potentials or quantum mechanically from ab initio pseudopotentials. In this work, the exchange-correlation potential term is approximated by the Ceperley–Alder functional [20,21]. The Troullier–Martins nonlocal pseudopotentials [22] in Kleinman–Bylander form [23] are used. The cutoff core radii adopted are 1.50, 1.54, and 1.54 a.u. for C 2s, 2p, and 3d electronic states, respectively [24–26], in constructing the Troullier–Martins scheme pseudopotentials. With respect to computing the quantum forces, the cluster in question is placed in a spherical domain. Outside of this domain, the wave function is required to vanish. The radius of the sphere is such that the outmost atom is 7 a.u. from the boundary. The grid spacing was 0.5 a.u.,
332
X.-R. Chen et al. / Chemical Physics Letters 380 (2003) 330–336
which is evaluated roughly with a plane wave 2 cutoff of ðp=hÞ or about 30 Ry. For our study of the structure of the C8 cluster, the initial simulation temperature is taken to be 4100 K, and the final temperature 100 K. The annealing schedule lowers the temperature 500 K each 50 time steps. The time step is taken to be 80 a.u. (1 a.u. ¼ 0.0241888428 fs). A friction coefficient is taken to be 0.0006 a.u. After reaching the temperature of 100 K, the cluster is quenched to 0 K. The ground state structure of this cluster is found through a direct minimization by a steepest descent procedure. It is very important to choose the initial atomic configuration. If the atoms are too far apart, they may not form a stable cluster as the simulations proceed. If the atoms are too close together, they may form a metastable cluster from which the ground state may be kinetically inaccessible. Typically, a random placement of the atoms must reside within 1.05 and 1.30 times the dimer bond length. If the initial atomic configuration cannot form a stable cluster, there are two methods to cope with it. One method is that, you may check the energy per atom of the first step to check that it nearly equals (no more or less than 4 eV) that of the dimer. When the energy per atom is much greater than that of the dimer, you must reform the initial configuration. The other method is that, the time step may be smaller within the first fifty steps. The detailed account of this method has been reported elsewhere [15].
approximate owing to the use of the LDA instead of the exact spatial dependence, and an adiabatic approximation instead of the exact time dependence of the exchange-correlation functional. The time-dependent LDA (TDLDA) technique requires a substantially smaller computational effort and is able to handle much larger clusters. This technique has recently been applied to Na, Si, GaAs, and CdSe clusters and was found to yield results in very good agreement with experiment [27–29]. We here focus on applying the TDLDA technique to the calculation of the optical absorption spectra of the C8 cluster. We apply it within a frequency-space approach, in which the electronic transition energies Xn are obtained from the solution of the following eigenvalue problem: h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii pffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2ijr dik djl drs þ 2 fijr xijr Kijr;kls fkls xklr Fn
2.2. Time-dependent density functional theory within the local density approximation
where /ðrÞ is a Kohn–Sham wave function, i; j; r; k; l; s are the occupied states, unoccupied state, and spin indices, respectively, qðrÞ is the Kohn– Sham charge density, and mxc ðrÞ is the LDA exchange-correlation potential. This computational technique is based on the higher-order finite difference pseudopotential method [30]. The solution proceeds in three parts: First, the Kohn–Sham energies and wave functions are found by means of a conventional ground state calculation. Second, the wave functions are used to construct the coupling matrix of Eq. (2). Third, the TDLDA eigenvalue equation, Eq. (1), is solved and the absorption spectrum constructed. All three steps were performed using a real-space grid confined in a spherical domain [16].
It is known that the density functional theory states are uniquely determined by the ground state charge density. However, the classic formulation of the density functional formalism is restricted to the time-independent case only. A proper treatment of electronic excitations is not possible within the time-independent framework and requires a generalization of DFT to the time-dependent phenomena. This limitation has led to the development of TDDFT. Within the TDDFT, the main theorem of the density functional formalism is extended to time-dependent systems. Formally, this is an exact approach, whereas in practice it is
¼ X2n Fn ;
ð1Þ
where xijr ¼ ejr eir are the Kohn–Sham transition energies, and fijr ¼ nir njr are the differences between the occupation numbers of the ith and jth states. The eigenvectors Fn are related to the oscillator strength, and Kijr;kls is a coupling matrix given by Z Z 1 dmxc r ðrÞ K¼ /ir ðrÞ /jr ðrÞ þ jr r0 j dqs ðr0 Þ /ks ðr0 Þ /ls ðr0 Þ dr dr0 ;
ð2Þ
X.-R. Chen et al. / Chemical Physics Letters 380 (2003) 330–336
3. Results and discussion We first determine the ground state structure of the C8 cluster by simulated annealing via the Langevin molecular dynamics technique. We illustrate the simulated annealing procedure in Figs. 1 and 2. The initial geometry is a pair of parallel
333
rhombus containing several bonds. After 169 steps (327.03 fs), it forms an octagon structure (see Fig. 2(r)), but the structure is still somewhat removed from the ground state. After 368 steps (711.15 fs), the ground state structure of cyclic D2h (see Fig. 2(t)) is essentially formed. The result is different from that obtained using all-electron ab
-6.0
0.6 Ek (eV)
Binding Energy (eV)
0.8 -6.5 -7.0
0.4
-7.5 0.2 -8.0
0.0 0
100
200
300
400
500
0
100
200
300
400
500
Fig. 1. Binding energy and kinetic energy of the C8 cluster during a Langevin simulation. The initial temperature is 4100 K and the final temperature is 100 K. The time step is 80 a.u.
Fig. 2. The geometries of the C8 cluster at various time steps during a Langevin simulation.
334
X.-R. Chen et al. / Chemical Physics Letters 380 (2003) 330–336
initio molecular orbital techniques [13] by Raghavachari and Binkley [11] and different from the multi-reference configuration interaction (MRCI) € f [14], and also difstudy by Parasuk and AlmlO ferent from that predicted by Hutter et al. [2] We also calculate the ground state structures of carbon clusters Cn ðn ¼ 2–7Þ. Our results reveal several interesting aspects. The bond length of the (the experimental value is 1.2429 C2 dimer is 1.27 A [3]) and the corresponding binding energy is 5.44 A eV (the experimental value is about 6.2 eV [11]). The C3 cluster is an isosceles triangle with C2v symmetry. The corresponding bond length is 1.24 (the experimental value is about 1.30 A [3]), apex A angle is h ¼ 137:97°, and the binding energy is 14.68 eV (the experimental value is 13.9 eV [11]). For the C4 cluster, the structure is a rhombus with , the minor diagonal length the side length 1.44 A and the binding energy 22.65 eV (the ex1.42 A perimental value is 19.4 eV [31]). The configuration of the C5 cluster is a linear structure. The four bond , respeclengths are 1.23, 1.25, 1.28, and 1.32 A tively. This structure is a little difference from that of the C5 cluster in [11]. The ground state structure of the C6 cluster is a monocyclic structure, as is different from the linear structure, but agrees well with that obtained through ab initio calculations of Raghavachari and Binkley [11]. The monocyclic structure is found to be about 0.2 eV lower than the monocyclic structure in energy. As regards the structure of the C7 cluster, it is also a linear structure, just like the structure of the C5 cluster [11]. In summary, our calculated results agree well with those obtained by Raghavachari and Binkley [11]. Xu et al. [32] have applied the tight-binding model to study these small carbon clusters. They have obtained very regular structures for these small carbon clusters. We do not agree with them.
Therefore, we believe that our scheme has made a successful prediction of the small carbon clusters. Moreover, all the results show that there is an odd– even alternation in the small carbon cluster geometries: the odd-numbered clusters have linear structures and many of the even-numbered clusters prefer cyclic structures. In Table 1, we list the average bond lengths together with our calculated total energies of these carbon clusters. It is obvious from our calculated , of the C4 results that the bond length, 1.44 A cluster structure is the longest among these small carbon clusters; the even-numbered small carbon clusters have longer bond lengths and the oddnumbered clusters prefer slightly shorter bond lengths. However, the absolute value of the total energies per atom increases roughly with the number of atoms of the cluster. For the TDLDA calculations of the C8 cluster with D2h symmetry as illustrated in Fig. 2(t), we carefully have tested for convergence of the computed excitation energies and associated absorption spectra with respect to the size of the spherical domain, the grid spacing, and the total number of electronic states included in the calculations. Convergence is obtained with a domain radius of 33 a.u., a grid spacing of 0.6 a.u., and a number of unoccupied Kohn–Sham orbitals that is at least three times greater than the number of occupied ones. The calculated TDLDA absorption spectrum (solid line) together with the absorption spectrum from time-independent LDA (dashed line) are shown in Fig. 3. All spectra are broadened using a Gaussian in order to simulate the effect of a finite temperature as well as the finite resolution of an experiment. We use a Gaussian width of 0.1 eV, which is a typical accuracy for optical absorption measurements.
Table 1 ) and calculated total energies ET (eV/atom) of Cn ðn ¼ 1–8Þ The average bond lengths L (A n a
L Lb ET a
2
3
1.27 1.2429 )150.90
1.24 1.2968 )152.80
Our work. b Average experimental data from [3].
4 1.44 )153.02
5
6
1.27 1.2855 )153.65
1.30 1.2780 )153.58
7
8 1.28
1.31
)153.94
)154.21
Absorption cross section
X.-R. Chen et al. / Chemical Physics Letters 380 (2003) 330–336
335
ties. The TDLDA computed spectrum, which exhibit a variety of features that could be used for cluster identification, differ significantly from those computed using a simple LDA approach. It should allow for a comparison with future experimental investigations.
2.0 1.5 1.0 0.5 0.0
0
5
10
15
20
Fig. 3. The calculated optical absorption spectra (arbitrary units) for the C8 cluster with D2h symmetry, where the solid and dashed lines correspond to the case of TDLDA and LDA, respectively.
It is noted that the calculated TDLDA absorption spectrum are characterized by three extremely intense absorption peaks near 7.6, 11.3, and 17.5 eV, respectively, and display a significantly blue shift with respect to the Kohn–Sham eigenvalue spectrum, which are consistent with the well-known band gap underestimate of LDA. The first TDLDA absorption peak display a blue shift by about 3.0 eV with respect to the Kohn–Sham eigenvalue spectrum, the second display a blue shift by about 1.8 eV and the third by about 0.6 eV. This implies that TDLDA works well for systems with a significant degree of ionicity. There is a qualitative difference between the current results and those of Gan Asn [33] and Si2n [27] clusters. Similarity to III–V and IV–IV clusters is that the long absorption tails have been found in the calculated absorption spectra in the region of lower transition energies. As we know that, these tails, which extend deeply into the region of lower transition energies, are likely due to the existence of free surfaces in the clusters. Such behavior is absent in quantum dots that are passivated at the boundaries [34]. In summary, the ground state structure for the C8 cluster is calculated using the finite-difference pseudopotential density functional theory method in real space and the Langevin molecular dynamics annealing technique. We have found that the C8 ring has D2h symmetry and the TDLDA formalism can provide an efficient alternative to more complex theoretical methods for excited state proper-
Acknowledgements The authors thank Prof. J.R. Chelikowsky for providing us PARSEC_CODE, a serial version of real space, time-dependent density functional code. We also acknowledge support for this work by the National Natural Science Foundation of PR China under Grant No.10274055 and by the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry.
References [1] W. Weltner Jr., R.J. Van Zee, Chem. Rev. 89 (1989) 1713. [2] J. Hutter, H.P. Luthi, F. Diederich, J. Am. Chem. Soc. 116 (1994) 750. [3] J.M.L. Martin, P.R. Taylor, J. Phys. Chem. 100 (1996) 6047. [4] G. Galli, R.M. Martin, R. Car, M. Parrinello, Phys. Rev. Lett. 62 (1989) 555, ibid, 63 (1989) 988. [5] Z.Y. Lu, C.Z. Wang, K.M. Ho, Phys. Rev. B 61 (2000) 2329. [6] C.Z. Wang, C.T. Chan, K.M. Ho, Phys. Rev. B 42 (1990) 11276. [7] S. Fahy, S.G. Louie, Phys. Rev. B 36 (1987) 3373. [8] V.G. Grigoryan, M. Springborg, Phys. Chem. Chem. Phys. 3 (2001) 5135. [9] J.L. Wang et al., Solid State Commun. 117 (2001) 593. [10] J.J. Zhao, J.L. Wang, G.H. Wang, Phys. Lett. A 275 (1993) 281. [11] K. Raghavachari, J.S. Binkley, J. Chem. Phys. 87 (1987) 2191. [12] Y.L. Bai, X.R. Chen, X.D. Yang, P.F. Lu, Commun. Theor. Phys., to be published. [13] W.J. Hehre, L. Radom, P.V.R. Schleyer, J.A. Pople, Ab initio Molecular Orbital Theory, Wiley, New York, 1979. € f, Theor. Chim. Acta 83 (1992) 227. [14] V. Parasuk, J. AlmlO [15] C. Troparevsky, J.R. Chelikowsky, J. Chem. Phys. 114 (2001) 943, and references therein. € gu €t, J.R. Chelikowsky, Phys. Rev. B 65 [16] I. Vasiliev, S. O (2002) 115416, and references therein. [17] N. Binggeli, J.R. Chelikowsky, Phys. Rev. B 50 (1994) 11764. [18] N. Binggeli, J.R. Chelikowsky, Appl. Phys. Lett. 66 (1995) 1316. [19] J.R. Chelikowsky, J. Phys. D 33 (2000) R33.
336 [20] [21] [22] [23] [24]
[25] [26] [27] [28]
X.-R. Chen et al. / Chemical Physics Letters 380 (2003) 330–336 D.M. Ceperley, B.J. Alder, Phys. Rev. Lett. 45 (1980) 566. J.P. Perdew, A. Zunger, Phys. Rev. B 23 (1981) 5048. N. Troullier, J.L. Martins, Phys. Rev. B 43 (1991) 1993. L. Kleinman, D.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. X.R. Chen, A. Oshiyama, S. Okada, Phys. Rev. B 67 (2003) 033408; Chem. Phys. Lett. 371 (2003) 528. X.R. Chen, A. Oshiyama, S. Okada, Q.Q. Gou, Chin. Phys. Lett. 20 (2003) 404. X.R. Chen, X.L. Zhou, J. Zhu, Q.Q. Gou, Phys. Lett. A 315 (2003) 403. € €t, J.R. Chelikowsky, Phys. Rev. Lett. 82 I. Vasiliev, S. O gu (1999) 1919, and references therein. € € t, J.R. Chelikowsky, in: P. Jena, S.N. I. Vasiliev, S. O gu Khanna, B.K. Rao (Eds.), Cluster and Nanostructure Interfaces, World Scientific, Singapore, 2000, p. 259.
[29] M. Claudia Troparevsky, Leeor Kronik, James R. Chelikowsky, Phys. Rev. B 65 (2001) 033311. [30] (a) J.R. Chelikowsky, N. Troullier, K. Wu, Y. Saad, Phys. Rev. B 50 (1994) 11355; (b) J.R. Chelikowsky, N. Troullier, X. Jing, D. Dean, N. Binggeli, K. Wu, Y. Saad, Comput. Phys. Commun. 85 (1995) 325; (c) Y. Saad, A. Stathoupolos, J.R. Chelikowsky, K. Wu, S. € gu €t, BIT 36 (1996) 563. O [31] K. Raghavachari, R.A. Whiteside, J.A. Pople, J. Chem. Phys. 85 (1986) 6623. [32] C.H. Xu, C.Z. Wang, C.T. Chan, K.M. Ho, J. Phys.: Condens. Matter 4 (1992) 6047. € gu €t, J.R. Chelikowsky, Phys. Rev. B 60 [33] I. Vasiliev, S. O (1999) 8477. € gu € t, J.R. Chelikowsky, Phys. Rev. Lett. 86 [34] I. Vasiliev, S. O (2001) 1813.