REVIEW PHYSICAL B
VOLUME 61, NUMBER 2
1 JANUARY 2000-II
Locally activated Monte Carlo method for long-time-scale simulations
M. Kaukonen, J. Perajoki, and R. M. Nieminen
Laboratory of Physics, Helsinki University of Technology, P.O. Box 1100, FIN-02015, Helsinki, Finland
G. Jungnickel and Th. Frauenheim
Laboratory of Physics, University of Paderborn, Warburger Strasse 100, 33098 Paderborn, Germany Received 11 August 1998; revised manuscript received 2 July 1999 We present a technique for the structural optimization of atom models to study long time relaxation processes involving different time scales. The method takes advantage of the benets of both the kinetic Monte Carlo KMC and the molecular dynamics simulation techniques. In contrast to ordinary KMC, our method allows for an estimation of a true lower limit for the time scale of a relaxation process. The scheme is fairly general in that neither the typical pathways nor the typical metastable states need to be known prior to the simulation. It is independent of the lattice type and the potential which describes the atomic interactions. It is adopted to study systems with structural and/or chemical inhomogeneity which makes it particularly useful for studying growth and diffusion processes in a variety of physical systems, including crystalline bulk, amorphous systems, surfaces with adsorbates, uids, and interfaces. As a simple illustration we apply the locally activated Monte Carlo to study hydrogen diffusion in diamond.
I. INTRODUCTION
The molecular-dynamics MD Ref. 1 simulation method is an extremely powerful tool to study microscopic motion based on the Newtonian dynamics of atoms interacting through a model potential. The equations of motion are solved using a mesh of discrete time steps. The time step must be short compared the phonon frequencies, usually a fraction of a femtosecond. Hence, the total simulation time, which is of the order of 1000 to 100 000 steps dependent on the model potential, is only in the picosecond region. At maximum, MD methods can be used to simulate atomic processes occurring on a time scale of nanoseconds. For a variety of physical situations, however, this is far too short to study the true dynamics of a system, in particular for diffusive processes such as atom migration on surfaces, certain formation processes during growth, and defect migration in bulk material. A recent extension of standard MD schemes due to Voter2 focuses on the simulation of such processes. By adding an articial boosting potential to the true local energy landscape of the atoms, the energy barriers to congurations which are normally not accessible by ordinary MD can be overcome. Thus, the atoms are forced to do movements which are related to barrier heights incompatible with normal thermal activation energies, yielding an extended time scale. Another class of methods that has been proposed to relax a system over a large period of time is based on the knowledge of the local energy barriers. In a method by Barkema and Mousseau3,4 the local energy barriers are explicitly searched for with an inverse conjugate gradient method using a modied force vector. Atoms are allowed to make jumps over the actual saddle points found according to a standard Metropolis Monte Carlo MC algorithm with a ctitious temperature 2500 K . This temperature is the parameter which controls the acceptance rate of certain relaxation processes and hence the time scale under consideration. This
0163-1829/2000/61 2 /980 8 /$15.00 PRB 61
parameter is, however, not known a priori and difcult to determine. Furthermore, the method suffers from the fairly general problem whether all saddle points relevant for the evolution of the system can be found. This scheme which is known as the activation-relaxation technique ART has been used in identifying relevant local relaxation processes in amorphous silicon at low temperatures.6 If the dynamics of a system can be described as a sequence of rather independent infrequent events, long time scales can be modeled using transition-state theory TST . Since its development in the thirties, TST has been applied to a wide range of phenomena.79 Voter developed a TST based kinetic Monte-Carlo KMC method for describing the dynamics of such infrequent events in a regular lattice and applied it to the study of rhodium clusters on Rh 100 .10 Here, we develop a method based on similar ideas which, however, is more general and can be used to investigate diffusion reactions without assuming a regular lattice. Since this requires some knowledge of the local energy landscape in the vicinity of a moving atom it also contains features common to the ART described above. In TST based methods, rate constants for infrequent events usually depend on the predetermination of reactants and products, e.g., on the knowledge of the local energy minima prior to and after a chemical reaction. Then, various schemes11,12 may be applied to reach transition states which are characterized by the saddle points of the energy landscape. From this information one can extract the transition probability for an event and the related time scale. However, if signicant reactions are missed initially or if the potential changes remarkably during the evolution of the system, the information gained from such simulations is quite restricted. Our method, therefore, starts from ideas similar to ART in that we focus on the determination of the most relevant if not all energy barriers that an atom sees in its immediate neighborhood at any time the atom is going to make a move. In the locally activated Monte-Carlo LAMC technique,13
980 2000 The American Physical Society
PRB 61
LOCALLY ACTIVATED MONTE CARLO METHOD FOR . . .
981
we imagine a model structure as a system for which the short range and the short medium range order14 are the prime factors responsible for the actual form of the local potential. Therefore, we concentrate on nding the energy barriers in the vicinity of an atom and the related smallest energy paths. This is done by efciently mapping the energy landscape in a few directions around each atom and dening a local distribution function for the probability of the atoms to escape from their current positions. Within LAMC an event is the instantaneous jump of selected atoms called the movers over one of their nearby barriers. For this, the local distribution of escape rates determines the probability in which direction the selected atom will move. Once the escape direction is chosen due to this distribution function the transition probability is assumed to be unity. Hence, movers are forced to jump even if the energy barrier in the selected direction is rather high. Since the typical escape rate for the event is mapped there is a welldened control over the real time scale in which this event takes place. Using the LAMC scheme together with a MD simulation introduces small random disturbances but allows to advance the clock after a jump according to the average escape rate of the mover and, therefore, extends the time scale enormously. The paper is organized as follows. In Sec. II we present the ideas behind LAMC in more detail. As an example we discuss in Sec. III the diffusion of a hydrogen atom in diamond, which is an important process for understanding the chemical vapor deposition CVD frequently used to deposit diamond thin lms.15 Note that this application involves modeling of heteropolar interactions although the problem is simplied due to the homopolar symmetric host in which the hydrogen atom is allowed to move. A discussion of the power of LAMC follows in Sec. IV.
II. METHOD
of any given structure provided the dynamical behavior becomes largely determined by infrequent events. Classically, the fundamental assumption in TST is that there exists a dividing surface in phase space with two properties: i it separates reactants from products and ii any trajectory crossing this surface will not recross it. The related rate constants which describe the equilibrium ux of particles through the dividing surface can be approximated to a good extent by simple transition state theory STST 10:
k STST n p
0 exp
E saddle E min /k b T ,
1
In traditional TST based MC studies, the possible reactions that may occur in the system are assumed to be known a priori. The global evolution of the system is separated into single atom events generic moves for which the typical barriers are predetermined and assumed to remain unchanged as the system relaxes. Usually the barriers are calculated for the generic moves of an atom in an otherwise ideal host matrix by the most accurate methods available16,17 before actually doing any structural optimization of the model. A Monte Carlo step consists of randomly selecting one of the atoms and one of the generic moves and of evolving the system according to the transition rate for this move. The rate is given as an exponential of its predetermined barrier height. The barrier controls the acceptance rate in much the same way as the total energy difference between the initial and a nal state of the model in an ordinary Metropolis MC method. In practical situations such as surface growth or relaxation of amorphous materials it is rather unlikely that the typical barriers remain constant over a longer period of time. Also, it is extremely difcult if not impossible to predict the most relevant generic events in particular in systems with many different types of atoms. Therefore, we present a method which within the limits of TST is suitable for the relaxation
where n p is the number of possible exit directions, 0 is the harmonic frequency, E saddle and E min are the energies at the transition state and at the minimum, respectively, T is the temperature and k b is the Boltzmann constant. The second basic TST assumption in practice is violated to a certain extent, since each crossing of the dividing surface does not necessarily correspond to a reactive event. Thus Eq. 1 gives an upper bound to the true rate constant. MD methods have been used to calculate dynamical corrections to STST rate by determining the fraction of TST surface crossings that lead to a true reactive event.18 These studies show that the STST rates are very close to the dynamically exact rate constants. Provided there exists a systematic way of nding all or at least all the lowest and signicant saddle points of the energy landscape in the immediate neighborhood of an atom one would be able to evolve the system in accord with the STST expression for the rate constant. We wish to implement this idea into a Monte Carlo type algorithm which can be easily combined with MD in order to study a structure dynamically under the inuence of long-term processes. Generally, there is no explicit restriction for the model potential which is used to evaluate the atomic interactions and, hence, the saddle points. The actual choice for this potential may be critical though and the most accurate quantum-chemical potentials should be given the preference. This is particularly essential in sensitive bonding situations such as in carbon when classical potentials1923 usually applied for large models may frequently fail to reproduce the true barriers in the structure. For example, Tersoffs classical potential which has been applied very successfully in a number of carbon studies21,2426 was shown to result in additional local minima and associated barriers in the energy landscape when studying the relaxation of an hypothetical icosahedral carbon cluster.27 In contrast, true density-functional or Hartree-Fock based self-consistent potentials require computer resources that would restrict the size of the models under consideration enormously and are almost impractical for real diffusion or growth modeling. For the total-energy calculations in this study we therefore utilize as a reasonable compromise the recently developed density-functional based tight-binding approach DFTB .28,29 This method derives its name from the use of self-consistent density-functional calculations for pseudoatoms in order to construct transferable tight-binding TB potentials for the non-self-consistent solution of the Kohn-Sham equations of the many body system. The main idea of the scheme is to superpose local atomiclike orbitals to make up the molecular
982
M. KAUKONEN et al.
PRB 61
basin in the energy landscape. To avoid extensive calculai tions, however, we consider this number n dd to be nite and i small. The n dd search directions are obtained in the followi ing way: take a random diffusion direction; the rest n dd 1 directions are chosen by uniformly sampling an imaginary sphere around the moving atom. The total number of search directions determines the success rate with which one will be able to nd all or at least the i essential local barriers in the next step of LAMC. The n dd for each mover can be made to depend on the local geometry i around the mover, for example, n dd may depend on the number of the nearest neighbors of the mover.
C. Search for the saddle points FIG. 1. The owchart of the LAMC procedure described in Chap. II. 1. Finding true diffusion barriers by the projected conjugate gradient (PCG) method
states where the set of local orbitals is chosen in such a way as to predict the total charge density of the full structure as well as possible. Then, the Kohn-Sham equations for the complete model need not be solved self-consistently but still give a reasonable total energy. The DF-TB method has been successfully applied to various carbon systems, ranging from small clusters28,30 to buckminster fullerenes and related oligomers,31 amorphous carbon systems,32 and carbon surfaces.3335 Moreover, heteropolar interactions of carbon with hydrogen, boron, and nitrogen have been modelled accurately with this technique.3537 A ow chart of the LAMC method that we introduce is presented in Fig. 1 and will be discussed below.
A. Choice of the diffusing atoms
In general all atoms in a given model structure may be involved in infrequent jumps diffusive events over barriers in their neighborhood. The number of atoms that are explicitly considered may be restricted in order to study only such processes that are associated with a certain atom type or a subregion of the full model such as a surface . In the following those atoms where events are initiated are called the movers which does not mean that the remaining atoms in the system are kept xed at their positions. The number movers will be designated by N diff . The reason for these denitions is to prevent other events such as the self-diffusion in the bulk from interfering with the particularly interesting cases such as the migration on a surface.
B. Choice of the global search directions
i The denition of general search directions r direction for which events may take place is a central feature of our method. We will call such a direction a diffusion direction hereafter although we do not necessarily restrict the scheme to diffusion in the classical sense. Diffusion directions may be assigned to any of the N diff atoms. The maximum number i of the possible diffusion directions per atom n dd basically determines whether we are able to nd all the local barriers that an atom sees when being activated from its harmonic
The global search directions dened so far in general do not contain the true barriers of the system although the saddle points are expected to be close, in particular when the number of search directions is large. To nd the relevant diffusion barriers we utilize a method recently introduced to specify changes in barrier heights on top of diamond surfaces due to the presence of dopants in subsurface layers.35 In this particular scheme single energy barriers are found by a series of conjugate gradient CG relaxation processes with modied forces for the diffusing atom. The principle is similar to what has been introduced to ART.3,4 In our method, however, the force on the diffusing atom is projected onto equidistant planes which are always perpendicular to the global search direction. Therefore, the diffusing atom can only relax within these virtual planes, whereas all other atoms can fully relax due to the true interatomic forces acting on them. Between two consecutive CG steps the diffusing atom is pushed from the current to the next plane on a straight line connecting the position in the current plane with the nal point on the global search direction. Then the whole system is CG relaxed under the restrictions described above. Note that all atoms can react to the changes in the position of a single mover and that the particles can even get around huge barriers that may exist along the global search direction. Figure 2 shows a snapshot of a typical situation when searching for a single barrier. The solid circles indicate the initial and nal position of the moving atom, dotted circles mark its positions along the migration path at various steps of the calculation. The start and end points along the global search direction which is marked by the straight dashed line i i for the ith single mover are indicated by r initial ,r nal, respectively. The relaxed position of the mover in one of the planes i is r current . As discussed above, in any diffusion step the mover is initially set to the crossing point of the next plane along the i i global search path and the vector r nal r current . Therefore, the mover will be always focused back onto the global search direction so that the search path cannot diverge which is the major difculty in another study5 that attempts to nd all the local barriers. The position of atom i in the next virtual plane is given by
PRB 61
LOCALLY ACTIVATED MONTE CARLO METHOD FOR . . .
983
in the structure to change their positions while constraining the mover onto the plane perpendicular to r initial r nal. The projected force for such a diffusing atom is calculated as
k F projected F k k k r nal r initial
F r r
2
,
3
where r r nal r initial , and k denotes the x,y, or z coordinates. The projected force is substituted for the true total force acting on the diffusing atom. After this procedure the CG total energy minimization operates in the usual way. Note also that projecting forces before actually calling the CG algorithm is sufcient for all of the subsequent CG steps since the conjugated directions constructed during the minimization of the energy functional depend linearly on the initial data. The method enables us to nd a minimum-energy path related to any of the given global search directions independently. Therefore, the scheme can take full advantage of parallel computer architectures. If the initial position of a mover is a local minimum this path must either contain at least one saddle point to overcome or the total energy increases steadily.
2. Finding the saddle point congurations
FIG. 2. The projected conjugate gradient method. Here the diffusing atom is moved in four steps from the initial position r in to nal position r in . Forces F projected acting on the diffusing atom are restricted to the planes perpendicular to the diffusion direction during the CG minimization.
i i r nal r current
i i r next r current
n total n current
,
2
where n total is the total number of steps between the initial and the nal point along the global search direction. This is the parameter which determines the accuracy of the barrier height nally found. It is chosen such that the step length i.e., the distance between the adjacent planes in Fig. 2 i i r nal r initial/n total becomes so small that the desired accuracy in the height of the energy barrier can be obtained. There is a certain degree of freedom how to pick the nal position of a mover along a selected global search direction. In order to make the algorithm very much independent of this choice the position is set quite far away ( 100 ) from the initial local minimum. This is possible since the search for saddle points in the vicinity of the global search direction is terminated as soon as the rst relevant saddle is found. The maximum number of search steps determines the distance between the virtual planes that restrict the motion of i the mover. For the example in Sec. III we used r nal i r current 100 , n total 2000 yielding a spacing of 0.05 . The number of search steps already made (n current) ranges from zero to n total 1. The single search step is nished by relaxing the system with the CG method allowing all atoms
The saddle points are searched by moving diffusing atoms or movers along global search directions as described above. The movement of a single mover is continued until i a saddle point is reached with an energy at least E min higher than the total energy at the starting point or, ii the total energy becomes much larger ( E max) than the initial energy. For this work, E min and E max were chosen to be 0.2 and 7.0 eV, respectively. i The saddle point position of the ith mover r saddle is dened as the position vector of atom i where the total energy decreases for the rst time after leaving its initial position. The lower limit E min is used to suppress very frequent events. The detailed study of the short term behavior of the system is the standard task of ordinary MD simulations. Here, this limitation was set to 0.2 eV and we investigated migration processes with exceeding barriers. Therefore, large escape rates associated with the smallest energy barriers do effectively not contribute to the time scale. For case ii , any k esc attributed to associated directions is set to zero. The (i, j) (i, N saddle saddle point positions r saddle and k escj) are saved.
D. Choice of a diffusion event
After determination of the
N diff
N saddle
i 1
n dd i
4
saddle points related to the events i.e., all the barriers which surround all atoms declared to be movers , the escape rates i k esc are calculated using Eq. 1 assuming that the attempt frequency is the same for all events. In principle it is possible to estimate the true attempt frequency for each process using the harmonic approximation.10 This calculation of 0 , however, is demanding and beyond the current computer re-
984
M. KAUKONEN et al.
PRB 61
i sources. The number of diffusion directions n dd which are taken into account for an atom i is only restricted by the computer power. It may depend on the local atomic geometry. The saddle points found for each global search direction of every mover dene a set of local escape rates of the atoms from their harmonic basins. This set of escape rates can be regarded as a distribution function from which the most popular event, i.e., a mover plus its escape direction be can found by drawing with a probability given by the distribution function. In the simplest variant, the selection of the diffusing atom, the mover, and its diffusion direction, is a single step process. The particular diffusion event is selected from the pool of all possible ones by weighting the selection with the probability of the jump. In this way the physical reactions which can occur within the predened energy window of the saddle points E min Esaddle E max are explicitly taken into account in the STST-way according to Voter.10 Here, lower lying saddle points are more frequently chosen than others which may be a deciency in particular in heteropolar or inhomogeneous systems. In the worst situation the time increments are always almost the same and only a subset of events is considered. The evolution of the system may then become articially conned to a small part of the full system. In contrast, to increase the importance of slower events for the evolution of the system the procedure may be altered to a random selection of the mover followed by drawing the event in accord with the pool of the escape rates of the selected mover. So the sum in Eq. 5 is restricted to the diffusion directions belonging to the selected single mover only. All the diffusing atoms have the same probability to be chosen for the mover and the uctuation of time increments after completed diffusion jumps becomes remarkable Eq. 5 . Here, the long term evolution is governed by the diffusion events with large barriers and rapid diffusion events are partially suppressed. However, when choosing this superlong time scale, the detailed balance condition is lost. Thus the ensemble averages cannot be accumulated in a reliable way. We note that the energy window has the role of a time lter in a way similar to previous KMC studies.38,39
where the index j includes the diffusion directions of the atom which made the most recent jump. k esc is the same as k STST in Eq. 1 , except that the number of possible exit directions n p is set to one because all the saddle points are searched individually. This time increment t hop yields a lower limit of the time duration while it assumes that any of the diffusion jumps has taken place. An alternative way to dene the time increment is given by Battaile et al.17 However, in the LAMC, the pool of the escape rates may change after every diffusion jump, so it is not clear that the time increment suggested by Battaile et al. converges to the correct value.
F. Update of the distribution of escape rates
For those atoms whose local environment changes considerably due to the diffusion event the local escape rates need to be updated before drawing from the pool of escape rates again in order to dene the next successful event. This is done in the same way as described above in Secs. II A, II B, and II C. The update is only necessary for atoms that changed their coordination numbers in the recent event remarkably or that have moved during the subsequent CG steps more than a critical distance. Note that atoms may change their status from being dened as movers to nonmovers or vice versa during the simulation according to the denition of a mover.
III. SIMULATION DEMONSTRATIONS A. Results for hydrogen diffusion in diamond
E. Completion of the event
After the event has been chosen the chosen atom is set to (i, j) the precalculated Sec. II C saddle point position r saddle i.e. slightly behind the saddle point, as dened in Sec. II C . Thereafter the whole structure is CG relaxed without any further restriction or modications of the forces in order to nd the nal state after the successful diffusion event. Placing the mover at the actual saddle point is the activation step of the event in the ART terminology. The time dimension can be included into the simulations, if the attempt frequency 0 in Eq. 1 is known from either experiments or calculations of the vibrational spectra at the local minima and the saddle point congurations. The clock is incremented after an event by
n dd 1 j k esc
t hop
,
5
j
1
To illustrate our LAMC method we study hydrogen diffusion in diamond. We consider this process as a sequence of uncorrelated jumps from one interstitial site to another. A single hydrogen atom is moving between 64 carbon atoms enclosed in a cubic cell with periodic boundary conditions. Thus, the number of diffusing atoms N diff Sec. II A is just one. The number of global diffusion directions per atom n dd Sec. II B is set to six. The rst initial diffusion directions Sec. II B is selected randomly and the rest ve are sampled uniformly on an imaginary sphere. The six virtual nal points r nal Sec. II C are selected to be points towards the six initial diffusion directions and 100 away from the diffusing atom. The total number of diffusion n total steps is set to 2000, yielding the spacing of 0.05 between the PCG planes Sec. II C and Fig. 2 . Six carbon atoms on the faces of the supercell are xed to prevent a center of mass motion. This leads to higher energy barriers at the boundary of the supercell and to an articial reection of the diffusing atom. These effects do not change our conclusions signicantly. An alternative way to prevent the motion of the supercell is to add a constant force component to all atoms in the supercell after the true forces are calculated. The properties of interstitial hydrogen in diamond have been studied earlier by Estle et al.40 at the approximate ab initio Hartree-Fock level with the method of partial retention of diatomic differential overlap PRDDO , representing the bulk host by small saturated cluster models. Their calculations indicate that the lowest-energy site for hydrogen is an interstitial in the relaxed bond-centered BC site which appears to be 2.7 eV below the tetrahedral T site. By linear
PRB 61
LOCALLY ACTIVATED MONTE CARLO METHOD FOR . . .
985
TABLE I. Energies of various sites for H relative to the energy at a BC site and the energy barriers. All energies are in eV. This Work ET E AB
barr E T-to-BC barr E BC-to-BC barr E T-to-AB barr E AB-to-BC
Ref. 40 2.7 5.1
Ref. 42 5.3a
Ref. 41 2.7 2.5
Ref. 43 1.9
1.6 1.7 0.4 2.6 0.4 0.2
0.1 0.1 0.1 0.1 0.1 0.1
1.9
a
T site is not stable.
interpolation between the atomic positions of the relaxed BC and T site models, they could estimate the barrier inbetween to be of the order of 5 eV. The BC site has also been found more stable compared to the T site by Chu et al.41 on the same level of theory using saturated cluster models containing up to 44 carbon atoms. There, however, the barrier between the BC and the T site has been determined to be only about 2.5 eV above the BC site. A barrier of about 1.9 eV has been determined for the migration from one BC site to a neighboring BC site in the study by Mehandru et al.42 using the semiempirical atom superposition and electron-delocalization molecular orbital technique with a 46 C-atom cluster model. The densityfunctional pseudopotential self-consistent eld calculation using the local-density approximation SCF-LDA calculation by Briddon et al.43 utilizing a C26H3 O cluster model focuses on the stability of molecular hydrogen inside the diamond crystal. For monatomic hydrogen they nd the BC site to be more stable compared to the T site by 1.9 eV. The site energies and related barrier heights are summarized in Table I and compared to values obtained in the previous studies. We started our investigation with the determination of the equilibrium congurations for hydrogen in either the interstitial BC or the T site by CG relaxing idealized geometries, allowing all the atoms in the system to relax. When hydrogen is in the BC site relaxation forces the C-C bond containing the hydrogen atom in the middle to stretch by 0.80 giving rise to C-H distances of 1.17 . The C-C distance is 52% greater than the normal bond which is slightly larger than previously reported by Estle et al.40 42% , Mehandru et al.42 43.5% , and Briddon et al.43 43% . Within DF-TB the BC site is more stable than the T site by 1.6 eV. This energy difference is about 1 eV lower than found in the Hartree-Fock calculations40 but in very good agreement with the SCF-LDA results.43 We then applied the PCG algorithm to the calculation of the characteristic barriers separately. Our calculations do not fully support the results by Mehandru et al.42 with respect to the transition from a BC to a neighboring BC site. The barrier appears to be about 30% higher than previously reported and seems to be the most unlikely transition compared to the other pathways between T and antibonding AB , AB and BC, or T and BC sites. The latter three barriers are almost isoenergetic and about 2 eV lower than the former saddle point. This strongly suggests that the diffusion of hydrogen in diamond does not occur between neighboring BC sites.
The BC sites self-trap hydrogen and can additionally bind other Hs in a nearby AB conguration.43 This causes hydrogen to stay at BC sites for longer periods than at T or AB sites. After activation into one of those congurations H has almost equal chances to rapidly diffuse between them or back to a BC site where it is trapped again. This appears to be consistent with the picture for silicon.41,44 On the diffusion path from the BC site to the T site the hydrogen atom maintains its bonding to one of the two neighboring C atoms with an increasing bond length from 1.1 to 1.3 . The energy barrier is reached close to within 0.25 the T site. This saddle point conguration is characterized as a stretched tetrahedron with a H atom in the middle. The nearest neighbor H-C distances are about 1.3, 1.5, 1.7, and 1.8 . Note that within our method the system is continuously relaxed during the search of the saddle points and that it may be viewed as a search with least constraints. Typically, the relaxations are of the order of 0.4 for the C atoms neighboring the BC hydrogen and 0.1 for C atoms closest to the hydrogen at a T site. Therefore, it is not surprising that the barriers are noticeably lower than in the previous studies by Estle et al. and Chu et al., where a linear interpolation technique to constrain the geometries was applied to nd upper limits for the transition state energies. When studying the evolution of the system with the LAMC algorithm, diffusion takes place mostly between BC and T sites as expected. During a simulation of 20 diffusion jumps, only once a different state was reached. This new conguration is of particular interest and will be called the antibonding site which has not been observed for monatomic hydrogen in diamond before. For silicon the same conguration has been determined to be a local minimum44 in the energy landscape, too. The AB state occurs when the hydrogen atom is neither in an exact bond-centered nor a tetrahedral position. The AB state is reached by moving the H atom 0.15 from the tetrahedral position to the 111 direction. The neighboring C atom in this direction moves towards the H atom yielding a true C-H bond of 1.08 . The AB site is only marginally 0.1 eV higher in energy than the T site, and the barrier from a T site to the AB site is found to be 0.4 0.1 eV. We calculated the vibrational spectra for the BC and T sites, as well as for the transition-state between them, in order to get the attempt frequencies for BC to T and T to BC reactions. Equation 1 yields k BC-to-T 4 7.76 1012 exp( 2.0 eV/kT)(1/s) and k T-to-BC 6 1.16 1013 exp( 0.4 eV/kT)(1/s). At 1100 K these result to k BC-to-T 21360(1/s) and k T-to-BC 7.02 1011(1/s). With this information the diffusion constant can be estimated as44,45 D 1 6 R i R j 2n ik i j , 6
i, j
where n i is the probability for the hydrogen atom to be located at a given site in the lattice. 16 8 Here n BC 24 , n T 24 exp( 1.6 eV/kT) 1.56 10 8 at 1100 K, and R i R j R T R BC 1.47 . We thus arrive at an estimate of D 9.0 10 13 cm2 /s. The diffusion constant D can be evaluated another way by calculating it directly from
986
M. KAUKONEN et al.
PRB 61
0.60 and to the reverse direction 0.12 eV Munro et al. 0.03 eV . However, the supercell in the calculations by Munro et al. consists of 216 atoms compared to ours of 64. These results demonstrate that the current method can also handle diffusion processes involving bond rotation.
IV. DISCUSSION
FIG. 3. The square of the displacement of the hydrogen atom in diamond; only BC sites are included.
D
1 lim 6 t
d dt
x t
x 0
2
,
7
where the average is an ensemble average. The time evolution of the square of the displacements of the diffusing H atom is shown in Fig. 3. In this evaluation only the BC positions are taken into account. By calculating the diffusion constant from Eq. 7 we obtain again D 9.0 10 13cm2 /s. This is expected since the Eqs. 6 and 7 should give the same result when the simulation time approaches innity. The Eq. 7 is of course more convenient in a general diffusion case when neither the diffusion paths nor the metastable states are known a priori. Our estimated diffusion constant for the ideal crystalline host is about three times larger than experimental values.46 However, in experiments diamond contains point defects, such as vacancies, impurities, and grain boundaries. The BC interstitial itself has a partially lled level close to the conduction band which becomes lled when a second hydrogen in an AB type of conguration or other impurities such as nitrogen are present nearby.43 Hence, lattice defects can trap a diffusing H atom and thus reduce its diffusion rate.
B. Self-interstitial diffusion in silicon
Additionally we ha...