J. Phys. Chem. 1990, 94, 4480-4482
4480
A Model Potential Function for Carbon Systems: Clusters T. Takai,+ C. Lee,' T. Halicioglu,* and W. A. Tiller Department of Materials Science and Engineering, Stanford University, Stanford, California 94305 (Received: November 10, 1989; In Final Form: January 2, 1990)
I n this investigation a new model potential energy function for carbon systems was developed. The function has a total of nine parameters and it is based on a two-body and an angle-dependent three-body interaction. The parameters for the function were evaluated from fits to various properties of C2, graphite, and the diamond crystal. The potential function reproduces exactly the experimental values for the mininum point properties of C,. Energy- and structure-related properties for bulk graphite and diamond were reproduced with less than 4%deviations. The function was employed to calculate stable structures of carbon clusters with varying shapes and sizes. For clusters of C, it was found that for small n linear clusters are energetically more stable. For increasing n, first cyclic then hexagonal honeycomb structures become the ground-state configuration. Among several three-dimensional polyhedral structures, only C6, was found to have the lowest energy. Calculated results are in a general agreement with the available data
In atomic level simulation studies, the potential energy function describing interactions among the particles of a system plays a central role in the accuracy of the calculated properties. Because high-level first principles calculations are applicable to systems with only a very small number of atoms, presently, in calculating energy- and structure-related properties of systems containing a large number of atoms, semiempirical or model potentials are employed. Despite the fact that energy values calculated by the semiempirical or model functions are less accurate than those obtained by ab initio calculations, today, they are used quite widely due to their simplicity and ability to produce meaningful results for larger systems (almost approaching macroscopic dimensions). A number of model potential functions have been developed for silicon, for instance, and used successfully in various simulation calculation^.'-^ Recently, Tersoff developed a new model function for carbon systems with I 1 adjustable parameters, describing energy- and structure-related properties of the diamond crystal structure quite accurately.* He also showed that the function can describe properties of the individual basal planes of graphite. A variant of the Tersoff function was also used more recently by Brenner.9 He evaluated a different set of parameters for carbon species and successfully calculated various properties of the diamond crystal and the individual graphite planes. In both cases, however, this function fails to represent interactions between the planes of the graphite structure, because of its short-range nature. In the present study, a different model potential function for carbon systems is developed. The function has a range somewhat longer than the Tersoff potential, and it is used here for simulating the structural properties of carbon clusters. In general, it is assumed that, in the absence of external forces, a function @(Fl,...,FN) exists describing the total potential energy of a system of N atoms as a function of their positions which may be expanded aslo
j> I
j>ik>j
where u(Ti,Fj) represents the two-body interaction between two particles at positions Fi and Fj, and U ( F j , i j , & ) represents the three-body interaction among particles at positions and Fk, respectively. In this so-called many-body expansion of a, it is usually believed that the series has a rapid convergence; therefore, the higher moments may be negligibIe.'J0 In the present work, the total potential energy of the system is assumed to be given by sums of two- and three-body functions; accordingly, four and higher body terms are neglected. 'Department of Industrial Arts, Faculty of Education, Kagawa University, Takamatsu, Kagawa, 760 Japan. Department of Metallurgical Engineering, Inha University, Inchon, 160 Korea.
0022-3654/90/2094-4480$02.50/0
TABLE I: Comparison of Simulated and Measured Property Data for the Carbon System sDecies DroDertv exptl [ref] calcd cz bond length, A 1.2425 1291 I .2425 -6.21 [ i s ] - -6.214 binding energy, eV force const, eV/A2 75.98 [29] 75.98 -7.427 [30] -7.597 graphite crystal cohesive energ unit cell dimensions, a 2.4508 [31] 2.3592 C 6.7076 [31] 6.7137 diamond crystal cohesive energ , e V -7.408 [30] -7.126 unit cell dimensions, a 3.5668 [32] 3.6172
1'
1
The specific representation of the two- and three-body interactions is always a central issue in simulation calculations. In this study, two-body interactions are assumed to be given by
where rij represents the distance between the particles i andj, and the parameters are denoted by ql through q5. While an exponential function describes the repulsive arm of the potential in this formalism, an inverse tangent function is employed for the attractive arm. The convergence of the attractive arm is very important in a realistic simulation. The use of the inverse tangent function, in that respect, provides an added flexibility for adjusting the converging features of the potential for the long-range regions. The three-body part, on the other hand, is represented by an angle-dependent function u(?i,?j,?k) = & (COS oi h)(COS o j -f- h)(COS 6 k + h))e-~(ri14r,t2+rJi2) (3)
+
+
where, 8,,0j,0k and rij,rik,rjkrepresent the angles and the sides of the triangle formed by the three particles. The three-body adjustable parameters are z, p , h, and b. The angle dependence of the three-body part is introduced via cosine functions and the exponential function provides the convergence. The parameters in eqs 2 and 3 were evaluated by a least-squares fit to the experimental data given in Table I. In the fitting (1) Tersoff, J. Phys. Reu. Lett. 1986, 56, 632.
(2) Dodson, B. W. Phys. Reu. 1987, B 35, 2795. (3) Stillinger, F. H.; Weber, T. A. Phys. Reo. 1985, B 31, 5262. (4) Biswas, R.; Hamann, D. R. Phys. Reu. Lett. 1985, 55, 2001. (5) Dodson, B. W. Dodson, Phys. Reu. 1986, B 33, 7361. (6) Blaisten-Barojas, E.; Levesque, D . Phys. Reu. 1986, B 34, 3910. ( 7 ) Pearson, E.; Takai, T.; Halicioglu, T.; Tiller, W. A. J . Crysr. Growth 1984, 70, 33. (8) Tesrsoff, J. Phys. Reu. Lett. 1988, 61, 2879. (9) Brenner, D. W. Mater. Res. SOC.Symp. Proc. 1989, 141, 59. (IO) Murrell, J. N.; Carter, S.;Farantos, S. C.; Huxley, P.; Varandas, A . J. C . Molecular Potential Energy Functions; Wiley: New York, 1984.
0 1990 American Chemical Society
Model Potential Function for Carbon Systems
The Journal of Physical Chemistry, Vol. 94, No. 1 1 , 1990 4481 -4.0 I
TABLE 11: Calculated Values of the Parameters for the Two-and Three-Body Interactions
two-body part 41 = 10.149804 q2 = 7.936986 (I/A) q3 = q4 = qs =
261.527033 (eV) 0.527263 ( I / & 3.071221 (A)
three-body part 2 = 20.0 eV h = 0.205 p = 1.340 b = 0.588 118,
procedure additionally, we considered the following constraints: (i) For the diamond crystal, the stability condition given by
-4.5
-
-5.0
-
-1.5
-
-6.0
-
-6.5
-
-1.0
-
t! I
e
($)=O was imposed at V = Vowhere Vodenotes the equilibrium volume. (ii) For the graphite case, the stability condition was established via the pressure tensors along the a and c axes (corresponding to the parallel and perpendicular directions to the basal plane) which were both forced to vanish. (iii) Minimized total energies for common crystalline structures (such as fcc, hcp, bcc,and sc) were constrained to remain higher than the minimized energy of the diamond crystal. (iv) Finally, the minimum energy of the linear C3 molecule was forced to be the lowest among all triangular C3 structures. Numerical values of the determined parameters are given in Table 11. This parameter set now satisfies all the conditions imposed above and, at the same time, it can reproduce all the input data within acceptable limits (see Table I). For Cz, the bond length as well as the dissociation energy and the force constant were reproduced correctly. For graphite, the cohesive energy was estimated within 2.2%. Calculations for the structural data (Le., the unit cell corresponding to the minimum energy configuration) produced primary cell dimensions a and c within 3.7% and 0.1% of the experimental values. For the diamond structure, the calculated cohesive energy was found to be within 3.8% and the lattice constant was reproduced within a 1.4%. The calculations carried out in this investigation include four different structural families of carbon clusters: (i) one-dimensional linear chain, (ii) planar cyclic, (iii) two-dimensional hexagonal honeycomb, and (iv) three-dimensional spherical shell (polyhedral) structures. In each case, the total energy of the cluster was minimized with respect to relaxation of the atomic coordinates. The total energy of a cluster calculated here, using eqs 1-3, is the binding energy of the cluster which is equivalent to the energy released in the formation of the cluster from isolated atoms. Throughout this paper, therefore, we used the basic thermodynamic definition for the binding energy which is a negative quantity for a bound molecule (Le., a negative energy of dissociation). As anticipated, the energy per atom decreases as the number of atoms, n, in the cluster increases. The variation in the minimized energy as a function of n is shown in Figure 1 for the four different types of clusters. In the smaller n region (for n = 3-5) the linear form is the energetically most favorable structure. Table 111 compares the present results for optimized energy and structure of clusters (up to five carbon atoms) with the values in the literature. In general, the calculated energies compare well with the published results. Even though the binding energies obtained in this work are slightly lower than those calculated by Raghavachari and Binkley" based on high-level first principles methods, they were found to be in good agreement with the availabe experimental r e s ~ l t s . ~ ~Bond ~ ' ~ lengths for cyclic structures were found to be about 9% longer than the corresponding linear forms. This is consistent with results obtained from first principles calculations." By the present potential function, however, the optimized bond lengths (for cyclic as well
K.;Binkley, J. S. J . Chem. Phys. 1987, 87, 2191. (12) Drowart, J.; Burns, R. P.;DeMaria, G.;Inghram, M. G. J . Chem. Phys. 1959, 31, 1131. ( I 3) Stull, D. R.;Prophet, H. JANAF Thermochemicul Tables; US.Natl. Bur. Stand.; Natl. Stand, Ref. Data Ser. 37; U.S.GPO: Washington, DC, (11) Raghavachari,
1971.
I
..- I
-?
0.0
I
I
10.0
I
10.0
30.0
I
I
I
50.0
40.0
60.0
70.0
n Figure 1. Calculated energies for clusters of carbon atoms as a function of the cluster size. Open circles are for linear chains, squares represent
cyclic forms, triangles are for two-dimensional hexagonal honeycomb structures, and three-dimensionalpolyhedral shell structures are represented by unconnected solid cirles. Energies are in eV per atom and n denotes the number of atoms in the cluster. TABLE 111: Calculated and Literature Values of Energies and Interatomic Distances for Small Clusters of Carbon"
cluster type trimer (linear)
calculated energy, bond length, eV
A
-12.85
1.225
trimer -12.56 (equilateral triangle) tetramer -19.56
1.335
(linear)
tetramer (rhombus)
-18.89
tetramer -17.04 (tetrahedra) -26.06 pentamer (linear) pentamer (cyclic) a
-25.68
rl = 1.225
r2 = 1.209 r = 1.360 e = 68.0
energy, eV
-12.7 -13.9 -13.7 -12.3
literature bond length,
A
[ I l l 1.278 [ I l l (121 1.277 [34] [I31 (331
-19.0 [I21 r , = 1.330 [I41 -19.4 [I31 r2 = 1.305 [I41
-17.3 (111 r = 1.425 [ I l l e = 61.5 [ i l l r = 1.46 (141 e = 63.0 [MI
1.474 rl = 1.211
-23.9 [ I l l rl = 1.271 [ I l l r2 = 1.226 -26.5 [I21 r2 = 1.275 [ I l l -26.7 ( 131 r = 1.335 r = 1.396 [35] r = 1.417 (351
Numbers in brackets are reference numbers.
as for linear forms) were calculated to be approximately 3.5% shorter than the reported ab initio results. Energetically, the linear form of C3 is about 0.3 eV more strongly bound than the equilateral triangular structure. For C4 also, the linear form was found to be the energetically most stable structure followed by the rhombus form with an energy difference of 0.67 eV. The three-dimensional tetrahedral structure, on the other hand, is in a much higher configurational energy state, and therefore, it is well separated from the linear and rhombus forms. Which one of these two low-energy forms represents the ground-state configuration is still a controversial issue.'"ls Even though exper(14) Bernholdt, D. E.; Magers, D. H.; Bartlett, R. J. J . Chem. Phys. 1988, 89, 3612.
4482
The Journal of Physical Chemistry, Vol. 94, No. 1 I, 1990
imental results seem to favor the linear structure to be the ground state, the issue is not yet completely resolved. In the case of C5, also, the present calculations indicate that the linear structure is energetically more favorable than the cyclic form. This outcome is consistent with the results obtained from various ab initio methods.' '-I9 In the present investigation, planar cyclic structures, from n = 6 to approximately n = 22, were found to have the lowest energy value. Within this size interval, the energy values (electron volts per carbon atom) varies from -5.61 1 for n = 6 to -6.571 for n = 22. These values are about 0.17-0.3 eV per atom lower than the corresponding energies for the linear forms (see Figure I ) . Earlier theoretical calculations2*22 on C, suggest that the linear structures should be favorable as long as n does not become too large (n IIO). Other c a c u l a t i ~ n s ,however, ~ ~ * ~ decrease this limit to n I 4. Calculations for 6 In II O by Raghavachari and Binkley," on the other hand, produced cyclic c6, C8, Clo and linear C,, C9 forms to be the ground-state configurations. In another theoretical report by Ray," who carried out ab initio calculations up to n = 20, it has been shown that cyclic ring structures are preferred over linear chains for C, clusters with n > 9. For n > 23, clusters with the hexagonal structure become energetically more favorable. This type of carbon cluster is constructed simply as multiples of Kekult structures and is geometrically similar to fragments of the graphite basal plane. Calculations were also carried out for the IO In I62 region. For Cl0a naphthalene-like structure was considered and, for larger members, only close-packed hexagonal forms (pyrene-like structures) were included. As shown in Figure 1, the energy values per carbon atom, in this case, varies from -6.060 eV for n = 10 to -6.992 eV for n = 62. For the case of the three-dimensional polyhedral structures, calculations included clusters with n = 12, 20, 30, and 60 corresponding to the shells of a regular icosahedron. While, CI2,C20, and Cj0 polyhedra are energetically less stable than the other three clusters forms, c 6 0 in the polyhedral form was found to be the energetically most stable configuration. The energy per atom for the optimized structure of c 6 0 (in the soccerball geometry) was calculated as -7.13 1 eV, corresponding to a polyhedral radius of 3.448 A. This calculated radius is in a good agreement with the experimental value of 3.5 A reported by Kroto et al.25 Its low binding energy, also, is quite consistent with the unusual stability (15) Ritchie, J . P.; King, H . F.; Young, W . S. J . Chem. Phys. 1986, 85, 5175. (16) Magers, D. H.; Harrison, R. J.; Bartlett, R. J. J . Chem. Phys. 1986, 84, 3284. (17) Ray, A . K. J . Phys. E: At. Mol. Phys. 1987, 20, 5233. ( I 8) Koutecky, J.; Fantucci, P. Chem. Reo. 1986, 86, 539. (19) Botschwina, P.; Sebald, P. Chem. Phys. Lett. 1989, 160, 485. (20) Kraetschmer, W.; Sorg, N.; Huffman, D. R . Surf. Sci. 1985, 156, 814. (21) Clementi, E.; Clementi, A. J . Chem. Phys. 1962, 36, 2824. (22) Clementi, E. J . Am. Chem. Soc. 1961, 83, 4501. (23) Koutecky, J.; Pacchioni, G. Ber. Bunsen-Ges. Phys. Chem. 1984, 88, 233. (24) Whiteside, R. A.; Krishnan, R.; DeFrees, D. J.; Pople, J. A , ; von Rague-Schleyer, P. Chem. Phys. Lett. 1981, 78, 538. (25) Kroto, H . W.; Heath, J. R.; OBrien, S. C.; Curl, R. F.; Smalley, R. E. Nature 1985, 318, 162.
Takai et al. of Cm as reported by many researchers.26-28 A model potential energy function based on two- and three-body interactions was developed in this study for systems containing carbon atoms. In the description of the total energy, as anticipated, the functional forms of the two- and three-body interactions are very important. The choice of functions of this purpose is, in general, a cumbersome job. In the present investigation a number of functions with different analytical forms and with a varying number of parameters were taken into consideration. Each of these functions was fitted to experimental data in a similar way and the one with the best performance is reported here. The total potential energy function reported in this work has nine parameters. They were estimated from fits to various properties of carbon dimer, graphite, and diamond. The potential function with the calculated parameters is able to reproduce exactly the experimental properties for C2 at its minimum point and it provides the graphite and diamond structures as stable forms with graphite being the energetically most stable form. The potential function was used in this study primarily to calculate energy- and structure-related properties for various clusters of carbon. Despite existing controversies in the literature about the ground-state configurations of carbon clusters, the overall results obtained in this study are in accord with the available data from both experiment and theory. The present results clearly indicate that, for smaller n of C, clusters, the linear forms are favorable. However, for increasing n, cyclic and eventually hexagonal honeycomb structures become energetically more stable. While the three-dimensional polyhedral structures for smaller n are not very stable, the energy for the special C a polyhedral cluster was found to be the lowest, in agreement with various experimental results. Although qualitatively correct, the optimized cohesive energy values for graphite and diamond structures as well as their respective unit cell dimensions produced by the present function exhibit some deviations (less than 4%) with respect to their experimental values. Some inaccuracies may also be expected in calculating other bulk properties such as elastic constants, phonon frequencies, etc. Thus, some "fine tuning" is probably still needed for this potential function to reproduce different properties to a high degree of accuracy.
Acknowledgment. This work was supported by N A S A Ames Research Center via Grant NCC2-297. (26) Curl, R. F.; Smalley, R. E. Science 1988, 242, 1017. (27) Kroto, H . W. Science 1988, 242, 1 1 39. (28) Iijima, S. J . Phys. Chem. 1988, 91, 3466. (29) Huber, K. P.; Herzberg, G . Molecular Spectra and Molecular Structure IV. Constants of Diatomic Molecules; Van Nostrand Reinhold: New York, 1979. (30) Hultgren, R.; Desai, P. D.; Hawkins, D. T.; Gleiser, M.; Kelly, K. K; Wagman, D. D. Selected Values of the Thermodynamic Properties of the Elements; American Society for Metals: Metals Park, OH, 1973. (31) Brandt, N . B.; Chudinov, S. M.; Ponomarev, Y . G. Semimetals: Graphite and irs Compounds; North-Holland: New York, 1988. (32) Wyckoff, R. W . Crystal Structures; Wiley: New York, 1963; Vol. 1.
(33) Chelikowsky, J. R.; Chou, M . Y. Phys. Reo. 1988, B 37, 6504. (34) Gausset, L.; Herzberg, G.; Lagervist, A,; Rosen, A. Discuss. Faraday Soc. 1963, 35, 113. (35) Ewing, D. W.; Pfeiffer, G. V. Physics and Chemistry of Small Clusters; Jena, P., Rao, B. K., Khanna, S. N., Eds.; N A T O AS1 Series B: Plenum Press: 1987; Vol. 158, p 67.
