Nuclear Instruments and Methods in Physics Research B 115 (I 9%) 501-504
*H &J
NOM B Boam Interactions with Materials 8 Atoms ELSEVJER
Empirical approach for the interatomic potential of carbon
t
E. Pablo Andribet, Javier Domhguez-VZizquez, A. Mari Carmen Perez-Martin, Eduardo V. Alonso ’ , Jose J. Jimhez-Rodriguez * Departamento de Electricidad y Electrhica,
Fact&ad de Fisicas, Universidad Complutense. E-28040 Madrid, Spain
Abstract A new empirical potential for both diamond and graphite structures, inspired by the previously published work by Heggie [J. Phys. Condensed Matter 3 (1991) 3065; Carbon 30 (1992) 711, has been developed. This potential accounts for the weak attractive interaction between planes as well as for a direct interpolation between the sp* and the sp3 orbital configurations. In order to check the validity of the proposed interatomic potential we study by molecular dynamics (MD) carbon structures (diamond and graphite) around the equilibrium and finally, the dynamic behaviour of a vacancy and an interstitial are also studied.
1. Introduction Molecular dynamics is a powerful method to reproduce numerically the time dependent behaviour of a system either isolated or under the influence of external forces. Structural changes and/or the relaxation of defects are deeply related with the interatomic potential. A good knowledge of the interatomic potential allows not only to evaluate but even to predict behaviours and to establish the limits of applicability under different conditions of the system. Carbon structure and in particular the graphite phase is a matter of interest. Its technological applications and its rather anisotropic structure makes of it a tempting system to study. There are, for graphite, some empirical interatomic potentials available in the literature [1,2]. Tersoff [2] does not account for the Van der Waals or weak interaction between layers that makes this potential unusable to study the dynamic behaviour of interstitials between planes. This interatomic potential [2] also does not describe the passage from the sp* to sp3 hybridization, as reported in Refs. [3,4]. An empirical interatomic potential [l] was specifically developed for carbon which accounts for the structural passage from graphite to diamond as described by Fahy et al. [5], unfortunately this empirical
potential is difficult to handle and requires a great computing time that makes this potential unapproachable. We propose, in this paper, a new and simpler empirical approach for the interatomic potential of carbon in which the three-body interaction is taken into account. The gradual transition from graphite to diamond structures is also considered. The next section is devoted to explicit the form of the interatomic potential and to the MD study of carbon structures (diamond and graphite) around the equilibrium. Section 3 is dedicated to discuss the way in which interstitials between planes in graphite structure should be treated following the theoretical [6] and experimental [7,8] indications from the literature. In this section the dynamic behaviour of a vacancy and an interstitial is also studied. Finally, summary and conclusions are exposed in Section 4.
2. The potential After the works of Ferrante et al. [9] and Abel1 [lo] the description of the interaction potential energy, Vij, between an i and a j atom, for a covalent system is widely accepted as Vij = A exp( - c1rij) - B exp( - @,,)
’ This work is dedicated to the memory of Bob Collins, a colleague and a good friend, with whom we had the pleasure to discuss the preliminary stages of this project. ’ Permanent address: DivisiQ Sistemas Integrados, lnvap S. E., 8400-Batiloche, Argentina. * Corresponding author. Fax: + 1.394.5196; E-mail: jJj@ eucmax.sim.ucm.es.
,
(1)
which is the so called Morse potential. rij is the distance between the i and j atoms and (A, a, B, p> is referred to as the set of parameters which will vary from one crystalline structure to another. For practical purposes everybody agrees to limit the range of the potential, since for many application short-range functions permit a consider-
0168-583X/%/$15.00 Copyright 0 1996 Elsevier Science B.V. All rights reserved VI. PROJECTILE INDUCED/ASSISTED PROCESSES SSDI 0168-583X(96)01 527-2
E. Pablo Andribet et al./Nucl.
502
able reduction in computing time. A smoothly vanishes from v, to vI is [2],
which
v < V[
107
f(v, v,, vu) =
function
Instr. and Meth. in Phys. Res. B 115 (1996) 501-504
I[ I
v -
f 1 -cm p-
vu-v,
1,
VI
1 ’
Vl~VslVu, v > v,
(2) where the parameters v, and v, are set as required below. The graphite structure requires beside the interaction potential for atoms within the layer, as described by Eq. (1). an interatomic potential to describe the longer range interaction for atoms in different planes. To describe this long range interaction we shall use the expression [ 111,
0i Fig. 1. 2D-image of the Wigner-Seitz
where A’ = 15.1685 eV A6 and a0 = 3.35 A corresponds to the distance between layers in a graphite structure and coincides with the minimum in the potential energy. The set of parameters (A, (Y, B, p) for diamond and graphite crystalline structures of carbon are, respectively, (2890.45 eV, 4.64 A-‘, 74.56 eV, 1.64 A-‘) and (2710.25 eV, 4.64 A-‘, 130.57 eV, 1.89 A-‘). These values are chosen to fit the corresponding values of the cohesive energy, lattice constant and bulk modulus of each structure. The numerical algorithm followed to obtain these values is known as simulated annealing [12] and the experimental values taken for the fitting are those from Ref. [13]. Eqs. (1) and (3) represent the two-body contribution to the interatomic potential. To add the three-body contribution we shall follow the clever idea proposed by Heggie [ l]: the initial premise is that the Wigner-Seitz cell is the best description of the local environment of an atom. Fig. 1 illustrates a 2D-image of the Wigner-Seitz cell of an i-atom. Let us drive our attention to the interaction between the i and j atoms. Fig. 1 also displays the way to
cell of an i-atom.
build the eijk angle. Note that the closer the k-atom to the i-atom the smaller the Oi, angle and consequently Silk, which is given by, 6;,
2Pijk = tg eijk = rij
(4)
’
where pip is illustrated in Fig. 1. Let us define sjj as sij = min{ sijk} ,
(5)
where the index k run over all the first neighbours of i with the exception of j. Consequently, 6,, is a good measurement of the distance between the i-atom and its neighbours. The Vii interaction will be therefore modified by means of the weight function fls,,, 6,, 6,). Eq. (21, with 6, = 0 and 6, = 1.4 or 1.7, respectively, for diamond or graphite. In others words, the potential energy of the i-atom will be,
Table 1 Energy per atom in the gradual transition from graphite to diamond structures
This work
Ab initio calculation [5] Bond length between
Bond length within
Buckling angle
layers [A]
layers [A]
khl
1.54
1.540
1.80
1.510 1.475 1.438 1.424 I .422 1.42 1.42
109.5 106.4 102.9 98.5 94.6 92.8 90.5 90.0
2.00 2.20 2.45 2.60 2.90 3.35
E [eV/atom]
E [eV/atom]
A
0 0.07 0.19 0.25 0.35 0.32 0.17 0.012
0 0.07 0.20 0.26 0.35 0.3 1 0.22 0.012
0.20 0.52 0.62 0.70 0.84 0.90 0.95 1
E. Pablo Andribet et al./Nucl.
Instr. and Meth. in Phys. Res. B 115 (1996) 501-504
503
.. .
. ?I...
.
.
l
l
l
*lo
. .
-9 .
where to differentiate between the interatomic interactions within the plane, Eq. (l), and the interatomic interactions between planes, Eq. (3), appear multiplied, respectively, by the cut-off functions [I - f(rij, r,,, r,)] and Ar,J, ru. r,), with r, = 2.2 A and r,, = 2.4 A. The explicit dependence of the set of parameters on Xij allows the gradual transition from graphite to diamond structures. Let us term P to any of the four parameters for a given structure, assuming a linear interpolation, one has, P(A)
= (1 - A)P(O)
+ AP( l),
(7)
where A = I corresponds to graphite and A = 0 to diamond structure. The question that emerges is: What is the A value for a given couple of atoms interacting? The answer is in Ref. [5]. In this work Fahy et al. calculated the total energy of the solid using local-density-functional theory with ab initio pseudopotential in the transition, maintaining rhombohedral symmetry, from graphite to diamond structure by hydrostatic pressure. The results obtained in these calculations are summarised in Table 1. This table also displays the A-values to fit the energy per atom in the transition obtained from Eq. (6) and the set of parameters which results from Eq. (7). A reasonable good adjustment of the A-parameter as a function of the bond length between layers, BLBL, is: A = - 0.2789BLBL2 + 1.7757BLBL - 1.8321.
(8)
There are two ways to test the empirical interatomic potential proposed above. First of all around the equilibrium and, secondly far from equilibrium. Concerning the former point we have calculated the energy per atom as a function of the distance to the first neighbours, the characteristic acoustic and optical frequencies of the atoms and the elastic constants. A detailed description of all these calculations will be published elsewhere [ 141. With respect to the latter in the next section we shall study the dynamic behaviour of a vacancy and of an interstitial as a function of time.
3. Vacancy and self-interstitial The dynamic simulations have been carried out by the computer code MD-TOPS (Trajectories Of Particles Simulation) [15]. It uses a third order multi-step algorithm [16] to integrate the classical equations of motion and, due to the high anisotropy of the graphite structure, the force is calculated by evaluating the potential at seven points (with
N .
9-
2 -6
. . -7
. -....
.
00.
I
.
.
4 -2
0
2
.
.
b) -6
-3
0
3
6
x (4 Fig. 2. (a) Relaxation of a vacancy. (b) Relaxation of an interstitial.
three coordinates each); this guarantees a good conservation of the total energy of the system. A neighbour list, first suggested by Verlet [17], is used to make the simulation more efficient in computer time. The crystallite used in the simulation is composed of 212 atoms distributed in five layers. The dimensions of the crystal are, approximately, 13 X. 13 X 13 A. All atoms in the outermost planes are kept in their initial positions. The crystal is large enough to ensure that the deformations of the (0001) planes never reach the outermost planes. When an atom is removed, it is observed that the three first neighbours relax moving closer to their corresponding neighbours remaining in the plane, as it is shown in Fig. 2a. These three atoms remain bound by two u-bonds instead of three, resulting in a shorter bond length. This behaviour agrees with the previously reported by other authors [ 181. The formation energy of the vacancy is found to be 7.1 eV in very close agreement with the value suggested from experiment [ 191 and other calculations [6]. As far as we know from the results reported in the bibliography the behaviour of a self-interstitial is not clear in a graphite structure. Thrower and Mayer [ 191, in a quite complete review article about point defects and self-interstitials in graphite established that the bound interstitial and the interstitial mechanism of diffusion cannot be discussed separately. They reported that the positions of equilibrium for the interstitials are those points between layers placed below/above a hexagonal ring in the upper/lower plane and above/below an atom in the lower/upper layer. This would give rise to a sp3 hybridization which formation energy would be about 7 eV.
VI. PROJECTILE INDUCED/ASSISTED PROCESSES
504
E. Pablo Andribet et al. / Nucl. lnstr. and Meth. in Phys. Res. B I15 (I 996) 501-504
They also suggest the interstitial firstly studied by Wallace [20] would not explain the diffusion through the planes. These interstitials correspond to places in which there is an atom above and another below giving rise to two quasi sp3 hybridization. The formation energy for these interstitials is larger. More recently, Heggie [l] assumed that an interstitial interacts with the closest atoms in a plane as they would in the transition from graphite to diamond of the whole lattice. This is equivalent to admit that the Wallace’s interstitial is more stable. Xu et al. [6] studied the point-defects properties in graphite by a tight-binding-force model and they argue that the sp’ hybridization bond in the planes is too strong to be modified by the presence of an interstitial and consequently it is rather unlikely to find an sp3 hybridization as suggested in Refs. [1,20]. This means that an interstitial will be repelled by the planes in graphite. Indirectly this is also corroborated by experimental results [7,8] in which they have found bumps at the surface induced by single ion bombardment on graphite. Following this argument we have therefore assumed, in the following calculations, that interstitials will not bind to the atoms of the planes. Fig. 2b displays the relaxed structure of graphite containing a self-interstitial with the hexagonal ring below. Notice how the adjacent planes to the interstitials are warped. The effective volume of the interstitials results to 2.41 atomic volumes, in agreement with other predictions [21,22].
4. Summary and conclusions The proposed empirical interatomic potential of carbon is rather a simple potential which makes it a good candidate to be used in molecular dynamics calculations. Besides its simplicity it predicts very well the lattice constant, cohesive energy, bulk modulus, elastic constants and vibrational state distributions. The dynamic study of the lattice relaxation induced by a vacancy seems to have a time dependence behaviour in accordance with other previous calculations [ 181. Finally, the dynamic study of an interstitial in a graphite structure encourages us to proceed with the study of the ion beam modification of a graphite structure.
Acknowledgement
Helpful discussions with Prof. A. Gras-Martf are gratefully acknowledged. This work has been partially supported by DGICYT (project no. PB94-0283).
References [I] M.I. Heggie, J. Phys: Condens. Matter 3 (1991) 3065; Carbon 30 (1992) 71. [2] J. Tersoff, Phys. Rev. B 37 (1988) 6991; Phys. Rev. B 38 ( 1988) 9902. [3] D.W. Brenner, Phys. Rev. B 42 (1990) 9458. [4] K.E. Kbor and S. Das Sarma, Phys. Rev. B 38 (1988) 3318. [5] S. Fahy, S.G. Louie and M.L. Cohen, Phys. Rev. B 34 (1986) 1191; Phys. Rev. B 35 (1987) 7623; S. Fahy and S.G. Lottie, Phys. Rev. B 36 (1987) 3373. [6] C.H. Xu, C.L. Fu and D.F. Pedraza, Phys. Rev. B 48 (1993) 13273. [7] GM. Shedd and P.E. Russell, J. Vat. Sci. Technol. A 9 (1991) 1261. [8] R. Coratger. A. Claverie, A. Chahboum, V. Laudry, F. Ajustron and J. Beauvillain, Surf. Sci. 262 (1992) 208. [9] I. Ferrante, J.R. Smith and J.H. Rose, Phys. Rev. Lett. 50 (1983) 1385; J.H. Rose, J.R. Smith and J. Ferrante, Phys. Rev. B 28 (1983) 1835. [lo] G.C. Abell, Phys. Rev. B 31 (1985) 6184. [l 11 L.A. Girifalco and R.A. Ladd, J. Chem. Phys. 25 (1956) 693. [ 121 S. Kirkpatrick, C.D. Gelatt, Jr. and M.P. Vecchi, Science 220 (1983) 671. [13] M.T. Yin and M.L. Cohen, Phys. Rev. B 29 (1984) 6996. 1141 E.P. Andribet, J. Dominguez-Viizquez, A.M.C. Perez-Martin and J.J. Jimtnez-Rodriguez, to be published. 1151 K. GIrtner et al., Nucl. Instr. and Meth. B 102 (1995) 183. [I61 R. Smith and D.E. Harrison Jr., Comput. Phys. 3 (1989) 68. [17] L. Verlet, Phys. Rev. 159 (1967) 98. [18] E. Kaxiras and K.C. Pandey, Phys. Rev. Lett. 61 (1988) 2693. [19] P.A. Thrower and R.M. Mayer, Phys. Status Solidi A 47 (1978) 11. [20] P.R. Wallace, Solid State Commun. 4 (1966) 521. [21] J.H.W. Simmons, Radiations Damage in Graphite (Pergamon, Oxford, 1965). [22] R. Henson and R.W. Reynolds, Carbon 3 (1965) 277.