- Email: [email protected]
Grand canonical Monte Carlo simulations of intergranular glassy films in β silicon nitride
Materials Science and Engineering A 422 (2006) 123–135
Grand canonical Monte Carlo simulations of intergranular glassy films in  silicon nitride T.S. Hudson a,∗ , D. Nguyen-Manh b , A.C.T. van Duin c , A.P. Sutton a,d,1 a Department of Materials, Parks Road, Oxford OX1 3PH, UK EURATOM/UKAEA Fusion Association, Culham Science Centre, Abingdon OX14 3DB, UK c Beckman Institute (139-74), California Institute of Technology, Pasadena, CA 91125, USA Helsinki University of Technology, Laboratory of Computational Engineering, PO Box 9203, FIN-02015 HUT, Finland b
d
Received 14 February 2005; received in revised form 17 November 2005
Abstract Grand canonical Monte Carlo simulations have been performed on intergranular glassy films in  silicon nitride to determine their equilibrium structure, width and composition at 2000 K. A reactive force field for Si O N systems has been developed and tested, and chemical potentials of oxygen and nitrogen were obtained by requiring that they reproduce, in conjunction with the reactive force field, the experimentally determined compositions of the glass pockets at triple junctions and the interior of the adjoining  silicon nitride grains. Two silicon oxynitride films developed during the course of very long simulations, at the end of which their composition and structure were analysed. The thermodynamic equilibration of the film thickness was explored to give insight into the experimental uniformity of these films. Incorporation of a significant concentration of oxygen into the adjoining Si3 N4 grains was observed, producing somewhat diffuse interfaces on either side of the film. We believe this is the first thermodynamically consistent, atomic-scale, simulation of intergranular glassy films capable of exchanging matter and heat with reservoirs. © 2006 Elsevier B.V. All rights reserved. Keywords: Amorphous; Silicon nitride; Thin film; Interatomic potential; Monte Carlo; Grand canonical ensemble
1. Introduction Polycrystalline  silicon nitride (Si3 N4 ) is a high temperature structural ceramic with high strength and fracture toughness. The relatively high toughness of this material is derived from the fact that fracture occurs primarily intergranularly with a long crack path. The three factors that cause this are the toughness of individual grains, the large aspect ratios of those grains, and the weaker intergranular interfaces through which cracks can propagate [1–3]. The weak interfaces consist of silicon oxynitride glasses which form an intergranular glassy films (IGF) from the oxide surface films during the liquid-sintering process, and are often also modified using sintering additives such as rare-earth oxides [4]. One striking feature of these interfaces is their unifor-
∗ Corresponding author at: School of Chemistry, Bldg. F11, University of Sydney, NSW 2006, Australia. Tel.: +61 2 93514813; fax: +61 2 93513329. E-mail addresses: [email protected] (T.S. Hudson), [email protected] (D. Nguyen-Manh), [email protected] (A.C.T. van Duin), [email protected] (A.P. Sutton). 1 Present address: Department of Physics, Imperial College, London SW7 2AZ, UK.
0921-5093/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.msea.2006.01.014
mity, forming with a constant thickness throughout the sample. In the case of pure Si3 N4 –SiO2 , this thickness is observed to be 1.0 ± 0.1 nm [5], although another study has reported IGF thicknesses ranging from 0.6 to 3.0 nm for a material containing very low quantities of SiO2 [6]. The film thickness is almost independent of grain orientation, but statistical studies have shown that the most common planes for grain termination are those near the {1 0 1¯ 0} prism planes [7]. Great progress has been made in probing the composition and structure of the IGF experimentally. Using electron energy loss spectroscopy (EELS) the film has been determined to contain an anionic ratio for nitrogen of [N]/([N] + [O]) of about 30% [8]. High-resolution transmission electron microscopy has directly imaged the region immediately adjacent to the grains, showing the structural influence of these crystals into the film [9]. Similarly in materials prepared using rare-earth elements in the sintering glasses, high-angle annular dark field (HAADF) experiments have shown the sites these elements occupy in the IGF are strongly influenced by the presence of the crystal, often out to a few atomic layers into the glass [10]. Work has also been done using coherent convergent electron microscopy [11,12] to understand the short and medium
124
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
range order within the amorphous layer via radial distribution functions [12,13]. Modelling, in particular molecular dynamics (MD), has also contributed to our understanding of the detailed structure within these IGFs. One example demonstrated the role of N within the IGF, and the competition between the detrimental energy of mixing and the advantageous role this N plays in reducing the number of dangling bonds at the edges of the IGF [14]. Some studies have investigated segregation of metal ions within the IGF, and the decay of the crystallinity going into the film [15,16]. Other studies investigated the radial distribution functions and composition profiles of the film [13]. A common problem in these studies however is the need for both the composition and thickness of the IGF to be predetermined, because MD is a closed-system method using a fixed set of atoms. In this study, we use the grand canonical Monte Carlo simulation method [17,18] which is able to overcome this problem of predetermined IGF composition and thickness. It allows composition to vary in a thermodynamically correct fashion during the simulation. Chemical potentials are determined, which govern the exchange of material with reservoirs, and so the IGF is able, at least in principle, to establish the equilibrium thickness and composition. In reality, the IGFs are able to exchange matter with real reservoirs, such as the intergranular pockets and the adjoining crystalline regions, in order to establish equilibrium composition and thickness at a given temperature and pressure. 2. Reactive force field for Si O N systems Reactive force fields (ReaxFF) are many-body empirical potentials that seek to extend the force fields widely used for simulations of large molecules, to allow the breaking and forming of bonds [19,20]. They are then appropriate for simulations of condensed network materials, with the capacity to treat differently coordinated polymorphs correctly [19]. They do this by including a large number of physically motivated terms in the interaction energy, and by careful interrelation of these terms. Because these terms cover a wide range of physical effects, a single well-parametrised ReaxFF should be applicable to simulations of a wide range of chemical and material systems. The functional form of the ReaxFF is described in Appendix A. We have followed the notation of the previous descriptions and implementations of the ReaxFF [19–21]. In cases where they are conflicting we have followed the most recently used definitions, unless otherwise noted. 2.1. Parametrising the reactive force field ReaxFF parameters for Si O were previously developed and tested [19]. An attempt has also been made at Si O N parametrisation [13], but it was fitted almost exclusively to data from single molecules and so the applicability to condensed phase systems was unknown. Errors in the earlier Si O N parametrisation [13] were identified in an MD study of silicon nitride IGFs. We constructed a model IGF configuration using this parametrisation by annealing and then cooling a thin slab of silicon dioxide between 
silicon nitride crystals. A second model IGF was constructed by further annealing this IGF configuration [22] using a simple empirical potential (using potentials similar to those of Su and Garofalini [23]). The energies and stabilities of both configurations were compared [24] using density functional theory (DFT). The ReaxFF potential produced a configuration that had a much higher energy according to DFT as compared with the structure predicted by the simpler empirical potential. Further, DFT relaxations of both structures showed that the ReaxFF structure underwent a much more extensive relaxation than that predicted by the simpler potential, and despite this, the relaxed energy was still significantly higher. One reason for this poor performance was the preference of the ReaxFF parametrisation for an excessive concentration of N2 dimers to form within the -Si3 N4 crystal, resulting in under-coordinated Si atoms. We now discuss further inadequacies of the parametrisation developed in [13]. The ReaxFF parameters had been fitted [13] to a training set of Si O N configurations obtained by quantum chemical relaxations which included -Si3 N4 at and around the equilibrium volume under atmospheric pressure, and had very successfully reproduced the bulk modulus (see Fig. 1a). However, upon further investigation, we observed an instability to a fictitious new structure at higher volumes. This fictitious new configuration was calculated according to that ReaxFF parametrisation to have a lower energy than -Si3 N4 , even at the original volume. Thus according to the parametrisation of [13], the  phase is metastable, which is physically not the case. It was concluded that the instability was an indication that the ReaxFF parametrisation of [13] predicts unphysical results. The cause of the instability was elucidated by close examination of the key ReaxFF terms at the instability. The bond bending terms predicted an optimum N Si N bond angle of 155◦ , rather than the 109◦ of the tetrahedral configuration. This, in turn, was caused by a large overestimation of the character of Si N bonds (which are in reality primarily simple single bonds in this four-coordinated-Si network solid). The fitting of the earlier parametrisation [13] had predicted significant bonding at bond lengths that were as large as the onset of significant bonding. Nevertheless, our earlier parametrisation provided a good starting point for a new parameter fitting. The parameters governing the bond lengths at which bonding became significant were manually altered. Additionally, the N N and double bonding were prohibited because they caused excessive N N dimerisation in these solid-state systems when compared with experiment [25]. Although in consequence our new parametrisation is not transferable to molecular systems involving double or triple-bonded nitrogen, it provides greater accuracy for the current system under study. Torsional terms in the ReaxFF were also turned off, because although they contributed only a small energy relative to most other terms, they occasionally gave spuriously low energies during MC simulation, invalidating the run. Removing these terms significantly reduced the number of free parameters in the fit. The earlier training set [13] was augmented by the DFT energies of the silicon nitride IGFs, and those of crystalline ␣-Si3 N4 , -Si3 N4 , ␣-Si2 N2 O, and -Si2 N2 O. Fitting was performed using a combination of a genetic algorithm and a parameter-by-
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
125
Fig. 1. (a) Energy–volume curves for -Si3 N4 predicted by the original parametrisation: red — rigid expansion of equilibrium cell; black — expansion with atomic coordinates relaxed, showing instability at high volumes to a new fictitious structure; green — fictitious structure expanded with atomic relaxations, note that this fictitious structure has a lower energy than the  structure at all volumes. (b) Energy–volume curves predicted by the new parametrisation: red — rigid expansion of the equilibrium cell; black — expansion with atomic coordinates relaxed, showing no instability. The small difference between red and black curves at the ground state indicates that the equilibrium structure predicted by this ReaxFF parametrisation is slightly different to the DFT structure included in the training set.
parameter gradient relaxation. Only parameters specifically relating to Si N and N O bonding (and valence angle bond bending with N involved) were varied. The original general (speciesindependent) ReaxFF parameters (Table B.1) were kept fixed, as were the Si O parameters. This meant that a total of 33 parameters were varied. The final values of all parameters are shown in Appendix in Tables B.1–B.6, with the parameters allowed to vary indicated in bold. Once the fit converged with reasonable accuracy, the same volume expansion test as in Fig. 1a was applied to the new parametrisation. As shown by Fig. 1b, the  structure predicted by DFT is reproduced at all volumes except for a very minor relaxation. Furthermore, these parameters reproduce the ordering of the four IGF structures added to the training set as predicted by DFT calculation [24].
energies with the difference between the initial total energy and the final total energy. This was confirmed to a high degree of precision. The calculations of atomic charges and hence electrostatic energy could not be streamlined in this way. Electrostatic interactions are much longer range, and for every accepted MC move a recalculation involving the entire system is required. An approximate solution to this problem was implemented, whereby after each accepted MC move a region out to the electrostatic cutoff radius was recalculated, saving significant computational time. After every 50 MC attempts, charges on the entire system were recalculated to remove any residual error that had built up. The frequency of this reset was chosen so that average errors in the total energy after this number of steps were still significantly less than the thermal energy kT.
2.2. Testing the implementation of the reaction force field
3. Grand canonical Monte Carlo (GCMC) simulations
To test that the potential was implemented correctly in our Monte Carlo simulation program, the total energies of large Si N O IGF systems were compared to the energies calculated for the same systems by a previous molecular dynamics implementation of the ReaxFF [21,19]. For the same ReaxFF parametrisation, these two codes gave identical results apart from tiny numerical errors on the scale of computational floating-point accuracy. To ensure that the energy of the entire system did not have to be recalculated after each individual Monte Carlo (MC) move, “incremental” energy functions were introduced, which calculated only local changes in each term of the energy following each type of MC move. Because of the many-body nature of many terms of the ReaxFF, a change to a single atom often affects the energy of a number of nearby atoms, sometimes as far as the fourth nearest neighbour. The only exception to this is the electrostatic energy, which will be discussed later. Recalculating the changes to this limited region is computationally much more efficient than recalculating the energy of the entire system. The incremental local energy functions were tested by performing many MC moves on a system and comparing the sum of the incremental changes of local
3.1. Implementation Six types of MC move were included to alter the structure. Three are standard MC moves: translation of a single atom, swapping random pairs of unlike atoms, and swapping atoms with unlike neighbours. The other three are only applicable to the grand canonical ensemble, because they alter the composition of the system: insertions, deletions, and transmutations. The standard MC moves all have a straightforward acceptance criterion. The Boltzmann factor exp(− E kT ) is calculated, where E is the change in potential energy as a result of the move, k the Boltzmann constant, and T is the temperature of the system. If this quantity is greater than 1, the move is immediately accepted. Otherwise the quantity is treated as a probability of the move being accepted, and a pseudo-random number is generated to determine whether or not this occurs. The moves to simulate the grand canonical ensemble have a similar acceptance procedure, but care must be taken to satisfy detailed balance in this ensemble. The probability NA →NA +1 of inserting an atom of type A into a configuration with NA A atoms, NB B atoms, NC C atoms, . . . (NA + 1, NB , NC , . . .) multiplied by the occupation probability NA of this configura-
126
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
tion, must be equal to the probability NA +1→NA of removing an atom of type A from the configuration (NA + 1, NB , NC , . . .) multiplied by the occupation probability NA +1 of this configuration: NA +1 NA →NA +1 = (1) NA NA +1→NA The occupation probabilities are known from the grand canonical partition function [26] to have the ratio: NA +1 1 µA − E V = × 3 × exp (2) NA NA + 1 ΛA kT where V is the volume of the system, µA the chemical potential of A atoms, and ΛA denotes the thermal de Broglie wavelength of each atom of type A, given by the expression: ΛA = √
h 2πmA kT
(3)
where h is Planck’s constant, and mA is the mass of an atom of type A. The probability NA →NA +1 of inserting an atom of type A is equal to the probability of choosing to attempt an insertion, Pins , multiplied by the probability of choosing type A, 1τ , where τ is the number of elements involved in the simulation, multiplied by the acceptance probability min(1, χ). χ is the equivalent of the Boltzmann acceptance factor exp(− E kT ) for standard MC moves, and we derive it for the multicomponent GCMC moves here. The probability NA +1→NA of deleting an atom of type A from the (NA + 1, NB , NC , . . .) configuration is equal to the probability of choosing to attempt a deletion Pdel , multiplied by NA +1 the probability of choosing an atom of type A to delete, N+1 where N = tτ Nt (since we select the atom at random), multiplied by the probability of acceptance, min(1, χ1 ). Therefore, detailed balance (Eq. (1)) requires: 1 V µA − E exp NA + 1 Λ3A kT =
N +1 1 Pins (1/τ) min(1, χ) =χ Pdel (NA + 1/N + 1) min(1, (1/χ)) NA + 1 τ
(4)
Therefore,
τV µA − E exp kT (N + 1)Λ3A τ µ ˜ A − E = exp N +1 kT
χ=
(5)
where we have chosen to set Pins = Pdel . In the final step we have incorporated the factor ΛVA into the chemical potential to make it an effective chemical potential µ ˜ A = µA + kT ln( V3 ). ΛA
This is appropriate as long as our simulations are performed at a fixed volume and temperature. Eq. (5) implies that when considering an insertion, the accepτ ins tance probability is calculated as N+1 exp( µ˜ A −E ). For a delekT tion in a configuration, in which the total number of atoms is initially N + 1, the acceptance probability χ1 can be used. However, we must replace the factor of N + 1 with N if there are N atoms
in the initial configuration from which one is to be deleted. Also note that Edel (N → N − 1) = −Eins (N − 1 → N). So we have an acceptance probability for the deletion step N → N − 1 del of Nτ exp( −µ˜ A −E ). kT A similar argument for the transmutation move type leads to the conclusion that the probability of accepting an atom attempting to change from type A to type B should simply be −EA→B exp ( µ˜ B −µ˜ AkT ). These grand canonical steps can change the system’s composition, so there is a natural question about charge conservation in partially ionic systems such as silicon oxynitrides. A system with a net charge would cause major problems if the periodic images of the simulation cell were allowed to interact correctly. Because the ReaxFF uses a variable charge method to calculate its electrostatics, we are entitled to consider the overall system as always neutral, as long as we only ever insert neutral atoms. The charge will then distribute itself to minimise the total energy. Of course the most favourable energies will be found at stoichiometric compositions, so the system should remain near neutrality. 3.2. Chemical potentials We now turn to the task of determining the effective chemical potentials to use in the IGF simulation. When the IGF is at equilibrium with the glassy pockets at the triple junctions, and with the grains themselves, the effective chemical potentials of silicon, oxygen and nitrogen are constant throughout the entire systems. Because it is easier to extract accurate experimental information about the bulk phases, we may tune our chemical potentials to match experimental data for those phases, and apply these same chemical potentials to the IGF. The target concentrations of ([Si] = 18.48 nm−3 , [N] = 1.44 nm−3 , [O] = 34.78 nm−3 ) in the glassy pockets were based on EELS measurements [25] which also determined that the nitrogen was incorporated into the silicate network rather than being present as molecular N2 . In the  silicon nitride crystal, the corresponding concentrations of ([Si] = 40.6 nm−3 , [N] = 53.6 nm−3 , [O] = 0.56 nm−3 ) were based on the crystal lattice parameters, together with experimental measurements of the oxygen concentration using hot-gas extraction method on liquid-phase sintered -Si3 N4 [27]. Simulations were carried out of both the ␣-SiO2 phase containing nitrogen, and the -Si3 N4 phase containing oxygen. The initial ␣-SiO2 configuration was prepared by Burlakov [28] and contained only silicon and oxygen atoms. An algorithm was adopted whereby if the concentration of any one element was drifting away from the target composition, µ ˜ was adjusted slightly to counteract the drift. The amount by which it was adjusted diminished as the simulation went on, so that eventually convergence was achieved. The values established are shown in Table 1. The value of µ ˜ Si could not be evaluated in the -Si3 N4 simulation, owing to the extremely high energy penalty for creating any Si defect. Thus only our value of µ ˜ Si for the glass is of significance. Assuming the experimental target compositions
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
127
Table 1 Effective chemical potentials at 2000 K ␣-SiO2 phase -Si3 N4 phase Used for IGF
µ ˜ Si (eV)
µ ˜ N (eV)
µ ˜ O (eV)
4.5 N/A 4.5
2.9 1.6 1.9
1.5 2.2 1.9
Fig. 4. (a) Starting configuration of run 2 with a 1.325 nm gap between prism faces of a single crystal of -Si3 N4 . (b) Snapshot of final configuration after 1.44 × 107 GCMC moves.
3.3. Interfacial simulations
are in fact at equilibrium with one another, the differences in the values of µ ˜ N and µ ˜ O for the ␣-SiO2 and -Si3 N4 phases indicate the extent to which our model potential is not thermodynamically consistent. Using the effective chemical potentials obtained from the ␣-SiO2 simulations in a simulation of the Si3 N4 phase produces an excess of O atoms. Similarly using the effective chemical potentials obtained from the -Si3 N4 phase in the ␣-SiO2 phase produces an excess of N atoms. This implies that our potential fails to penalise sufficiently defects in one or both phases, indicating that there is further improvement possible for our potential. For our further simulations we used the intermediate values shown in Table 1, which after equilibration in the -Si3 N4 phase produce an anionic ratio for oxygen of 3%, and after equilibration in the ␣-SiO2 phase produce an anionic ratio for nitrogen of 12%.
In order to simulate a grain boundary, a -Si3 N4 crystal was cleaved at the plane most commonly observed experimentally to terminate Si3 N4 crystals. A gap of 0.3 nm was introduced by separating the two halves of the cleaved single crystal, as shown in Fig. 2a. Because the crystal period in this direction ˚ it is not possible to restore this system to a peris 13.25 A, fect crystal by filling this gap with additional material. Thus the system is forced to accommodate a planar interface. Periodic boundary conditions were imposed both parallel and perpendicular to the interface, and the sides of the periodic supercell (5.59820 nm × 1.52946 nm × 0.87513 nm) were not allowed to change. Our GCMC simulation was then run to determine the equilibrium composition and thickness of this interface. Fig. 2b shows a snapshot of the simulation after 2.64 × 107 MC moves. The system quickly filled the gap, but was unable to generate a well connected bond topology, and so gradually expanded the amorphous oxide layer to relieve this misfit. By tracking the energy and composition of the system (as shown in Fig. 3), we established that this system had not reached an equilibrium by the end of 2100 h of CPU time on a Pentium 4 single-processor machine. A range of interfaces can be generated in this manner, each with a different initial gap between 0 and the crystal period. To
Fig. 3. Cell composition throughout run 1.
Fig. 5. Cell composition throughout run 2.
Fig. 2. (a) Starting configuration of run 1 with a 0.3 nm gap between prism faces of a single crystal of -Si3 N4 . (b) Snapshot of final configuration after 2.64 × 107 GCMC moves.
128
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
Despite the failure of run 2 to close the gap and return to perfect crystal, when compared to run 1 this simulation did show significant differences. We note that the IGF did not penetrate any further into the crystal than the width of the initial gap. The concentration of O defects in the surrounding Si3 N4 is also noticeably lower than in the previous simulation, and does not continue increasing, possibly suggesting a weaker elastic field. Thus in contrast to the run 1, we conclude that the film thick-
Fig. 6. (a) Average composition profile of the IGF during the final period of run 1. (b) Average cumulative frequency of oxygen through the IGF during the final period of run 1.
explore the behaviour of a system which began with a wider gap, in run 2 we chose to study a gap of exactly one crystal period (as shown in Fig. 4a). This has the advantage that the global free energy minimum structure is known — it fills the gap with perfect crystal, leaving almost no defects, and the bulk concentration of oxygen. Thus if our GCMC could explore configuration space very efficiently, it should converge to this structure. However, this did not happen. Just as metastable amorphous materials are unable to crystallise on a physical timescale at low temperatures, our GCMC was unable to find the global minimum within the simulation time of 1200 CPU hours. However, the simulation did converge to a metastable configuration within this time (see Fig. 5). Fig. 4b shows a final snapshot of the simulation. The periodic cell had dimensions of 6.62275 nm × 1.52946 nm × 0.87513 nm.
Fig. 7. (a) Average composition profile of the IGF during the final period of run 2. (b) Average cumulative frequency of oxygen through the IGF during the final period of run 2.
Fig. 8. Average coordination of each type of atom as a function of coordination radius during the final period of run 1. A point f on the line m at radius r means that a fraction f of that type of atom have at least m neighbours within the radius r. (a) silicon; (b) nitrogen; (c) oxygen.
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
129
ness is sufficient to accommodate a well connected oxynitride network, and that any further gains that might be obtained by making it thicker are outweighed by the energetic cost of disrupting crystal. The final period of each run was analysed to extract important physical information about the IGF. Fig. 6a shows the average composition profile through the IGF in run 1. Fig. 6b shows the
Fig. 10. Partial radial distribution functions in the IGF only during the final period of run 1.
average cumulative frequency of oxygen. The film thickness is considered to be the distance between the 10th percentile and the 90th percentile of this distribution (after correcting for the oxygen that is clearly in the crystalline phase), and at the end of ˚ We note that there is substantional this run it is calculated at 16 A. oxygen penetration into the -Si3 N4 phase. Fig. 7 is the equivalent data extracted from run 2. This data gives a metastable film thickness of 1.4 nm. The topology of the system was also studied, specifically with regard to the types of defects present. Fig. 8 shows the fraction of atoms with a particular number of neighbours within a given radius in run 1. Fig. 9 shows the equivalent for run 2. The radius at which first nearest neighbours are typically situated is clearly seen by the first rapid increase in the number of neighbours. At a radius slightly larger than this onset, nearly all silicon atoms have four neighbours, nitrogen atoms have three neighbours, and oxygen atoms have two neighbours, consistent with the valency of each element. A small fraction of the silicon atoms have three neighbours at this radius, and another small fraction have five, indicating the presence of a small quantity of over- and under-coordinated silicon defects. Nitrogen does not show any
Fig. 9. Average coordination of each type of atom as a function of coordination radius during the final period of run 2. A point f on the line m at radius r means that a fraction f of that type of atom have at least m neighbours within the radius r. (a) silicon; (b) nitrogen; (c) oxygen.
Fig. 11. Partial radial distribution functions in the IGF only during the final period of run 2.
130
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
under-coordinated defects, but there are a small fraction of these atoms which have four neighbours and are thus overcoordinated. Oxygen is an interesting case, showing a large fraction of the oxygen atoms to have three neighbours. The distance to the third neighbour varies. The large quantity of over-coordinated oxygen atoms is consistant with the observation that many of them substitute at nitrogen sites in the -Si3 N4 lattice (Figs. 2b and 4b). We have compiled average radial distribution functions (RDFs) for each pair of elements. Because RDFs within the Si3 N4 crystal are not of interest, only atoms within the IGF have been included in the analysis. Normalisation based on average density is performed as usual, however, because we analyse only a thin film, the large-radius asymptote is not to 1 as per normal. Figs. 10 and 11 show the results of these calculations for run 1 and 2, respectively. As expected, these results show a network primarily consisting of bonds between silicon and nitrogen, or silicon and oxygen. Some bonds between pairs of anions are also apparent, which are the primary defects within the IGF layer (bonds between silicon and silicon are not observed).
defects which would show little difference when directly imaged in an electron microscope. Thus unless elementally sensitive techniques such as EELS are used, care should be taken when making conclusions about interfaces that appear to be sharp. Acknowledgements The authors wish to thank Victor Burlakov for providing a simulated amorphous SiO2 cell. Thanks also to Stephen Garofalini and Wai-Yim Ching for their demonstration that our earlier ReaxFF parametrisation [13] needed improvement. Financial support from the European Community Growth Program (Contract No. G5RD-CT-2001-00586) is gratefully acknowledged. There are no competing interests in this work. Appendix A. Reaction force field functional form The terms used in the ReaxFF developed and used in this study are given by Esystem = Ebond + EvdW + ECoulb + Eqdist + Eval + Eover/under + Elp
4. Conclusions We have provided an improved parametrisation for the Si N O ReaxFF. Although not transferrable to systems including molecular N2 , this potential has been successfully applied to condensed Si N O systems. There is room for improvement however, as demonstrated by the thermodynamic inconsistency shown in Section 3.2 when determining chemical potentials in phases which appear experimentally to be in equilibrium. We suggest that further improvements could be obtained by including O ↔ N substitutional defect energies in the fit. We note that the requirement of thermodynamic consistency, in the sense that the chemical potential is constant for each element in phases that are in equilibrium with each other at a given temperature and pressure, is an extremely demanding test of any description of atomic interactions. It is a test that is rarely applied. We have developed and applied a GCMC scheme to simulate thermodynamic equilibrium in IGFs in the Si N O system. However, after many hours of simulation there remains some doubt as to whether ground-state equilibria have been reached, and in the case of run 2, despite reaching a metastable configuration, we are certain that the thermodynamic ground state had not been reached. This suggests that the attainment of true thermodynamic equilibrium in these network ceramics may still be out of reach of modern computation. Nevertheless, we have observed an interface at which an intergranular film formed and thickened into the adjacent crystals. We have also observed an initially wide interface at which an IGF formed and showed no tendency to thicken. The resulting thicknesses of the two simulated IGFs are comparable, and are consistent with experimentally observed IGF widths. Both simulations showed somewhat diffuse interfaces between the IGF and the prism faces of the crystals. The diffuseness was associated with the substitution of O atoms at N crystal sites up to a couple of layers into the crystal. These are substitutional
(A.1)
This is based, with minor alterations, on the form of the ReaxFF potential outlined by van Duin et al. [19], and we have followed the original notation as often as possible. Ebond , EvdW , and ECoulb are two-body1 energy terms representing the bond, van der Waal’s, and Coulomb interaction energies, respectively. Eqdist represents the energy required to ionize neutral atoms to the charges calculated by the electronegativity equalization method (EEM) [29]. Eval is a three-body bond bending term. Eov/und is a many-body term designed to penalise over or under coordination. Finally there are small contributions to the energy, Elp , due to the presence or absense of lone pairs. The differences when compared to Eq. (1) of the overview of ReaxFFs by van Duin et al. are as follows: Epen has been removed from the notation because it was not used in those potentials anyway; Etorsion has been removed to reduce the number of parameters in our potential, chosen because the magnitude of this four-body term is significantly smaller in these systems than other terms; Econj , energy due to conjugation, is not physically relevant to the Si N O system; Eqdist was already included in the practical implementation of the ReaxFF potentials, but was not explicitly mentioned in that article. A.1. Bond energy Two-body (“raw”) bond orders are calculated as follows: pbo,2 pbo,4 rij rij BOij = exp pbo,1 + exp p bo,3 r0 r0 rij pbo,6 + exp pbo,5 (A.2) r0 1
The bond energy is based on bond order. Although it is initially a two-body parameter, it is rescaled to an extent depending on the under or over coordination of both atoms (see Section A.1 for details).
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
where the three terms correspond to the bond, first bond, and second bond (referred to as ), bond orders, respectively, rij the distance between the bonded atoms, r0x the onset radii of each type of bond, and pbo,x are physically related to the onset shape and width of each type of bond. There is a subtlety in this calculation concerning the minimum fraction of a full bond order that is considered to constitute a bond. Many of the other terms in the ReaxFF rely on knowing which pairs of atoms are bonded, and it would be inefficient to consider atoms with a tiny bond order between them to be bonded at all for the purposes of those other terms. A cutoff in the bond order is thus introduced, with a value of 0.01λ30 . However, to ensure continuous potentials for molecular dynamics (MD), bonds with a total bond order of more than the cutoff are stretched such that the bond order becomes 0 as they approach the cutoff. Any bond order less than the cutoff is considered to be 0. These two-body bond orders are then adjusted (Eq. (A.3)) to account for either of the participating atoms being overcoordinated or undercoordinated: BOij = BOij f1 (i , j )f4 (i , BOij )f4 (j , BOij )
(A.3)
where i is the deviation of the sum of the uncorrected bond orders on atom i from its valency Vi (silicon 4, nitrogen 3, oxygen 2), and Vi + f2 (i , j ) 1 f1 (i , j ) = 2 Vi + f2 (i , j ) + f3 (i , j ) Vj + f2 (i , j ) (A.4) + Vj + f2 (i , j ) + f3 (i , j ) f2 (i , j ) = exp(−λ1 i ) + exp(−λ1 j ) f3 (i , j )
1 = − ln λ2
(A.5)
1 [exp(−λ2 i ) + exp(−λ2 j )] 2
where f13 (rij ) =
1 1 + exp(−λ3 (λ4 (BOij )2 − n ) + λ5 )
+
1 γw
λ29 1/λ29 (A.10)
A.3. Electrostatic energy The electrostatic contributions to the total energy of the system consist of two terms, Eqdist , the energy required to place charge on the atoms, and ECoulb , the Coulomb energy of interaction between these charges. The charges themselves are determined by the EEM [29], which requires a large matrix inversion, and is the most computationally expensive part of the energy calculations. Once the charge qi on each atom i is determined, the charge distribution energy can be calculated as follows. Eqdist,i = −χi qi + ηi qi2
(A.11)
where χi is the electronegativity, and ηi is hardness of atom i. The Coulomb interaction is given by ECoulb = T (r) where γcco =
(r 3
C q i qj + γcco )1/3
1 3/2 3/2 γEEM,i γEEM,j
(A.12)
is a screening parameter and may be
different for each pair of interacting elements. A 7th order taper function T (r) is also used to smoothly cut this interaction off in the shell between an inner cutoff radius of λ12 , and an outer cutoff radius of λ13 : T (r) = a0 + a1 r + a2 r 2 + a3 r 3 + a4 r 4
(A.7)
Note some corrections to small typographic errors in Eqs. (3a), (3c), and (3f) of van Duin et al. [20]. Once the bond orders have been corrected in this many-body manner, bond energy is evaluated as: p
+ a5 r 5 + a6 r6 + a7 r 7
(A.8)
a0 = a1 =
a3 =
A.2. Van der Waal’s energy The van der Waal’s energy between each pair of atoms is calculated using: f13 (rij ) EvdW = Dij exp αij 1 − rvdW 1 f13 (rij ) − 2 exp αij 1 − (A.9) 2 rvdW
(A.13)
where the coefficients are given by
a2 =
Ebond = −De BOij exp[pbe,1 (1 − BOijbe,2 )] − De BO ij − De BO ij
λ rij29
and the van der Waal’s interaction radius, rvdW , and γw may be different for each pair of interacting elements. The van der Waal’s energy includes short-range core repulsion as well as long range dispersion. Again note corrections to typographic errors in [20].
(A.6) f4 (n , BOij ) =
131
a4 = a5 = a6 = a7 =
−35λ312 λ413 + 21λ212 λ513 + 7λ12 λ613 + λ713 (λ13 − λ12 )7 140λ312 λ313 (λ13 − λ12 )7 −210(λ312 λ213 + λ212 λ313 ) (λ13 − λ12 )7 140(λ312 λ13 + 3λ212 λ213 + λ12 λ313 ) (λ13 − λ12 )7 3 −35(λ12 + 9λ212 λ13 + λ313 ) (λ13 − λ12 )7 2 84(λ12 + 3λ12 λ13 + λ213 ) (λ13 − λ12 )7 −70(λ12 + λ13 ) (λ13 − λ12 )7 20 (λ13 − λ12 )7
(A.14)
132
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
Eval = f7 (BOij )f7 (BOjk )f8 (j ){ka
A.4. Lone pair energy The energy associated with lone pairs such as those on oxygen and nitrogen atoms is calculated using the expression: Elp =
plp lp,i 1 + exp(−75lp,i )
(A.15)
lp,i is the difference between the actual number of lone pairs on atom i from the optimal number of lone pairs for that element, where the actual number of lone pairs is given by nlp,i = int
ei 2
+ exp −λ16
2 − ei
+ 2 int
ei 2
2
(A.16) ei is the difference between the total number of outer shell electrons, ne,i , and the sum of the bond orders around atom i. Note the corrected signs in comparison to Eq. (4) of van Duin et al. [19].
− ka exp[−kb (Θ0 − Θijk )]}
(A.20)
where p
f7 (BOij ) = 1 − exp(−Aj,26 BOijv,7 )
(A.21)
f8 (j ) = Aj,29 − (Aj,29 − 1) ×
2 + exp(λ15 j ) 1 + exp(λ15 j ) + exp(−pv,2 j )
(A.22)
Note that since previous publications of the ReaxFF, some parameters that were once global have been split into a parameter for each type of bond angle between particular elements (ka , kb , pv,2 , pv,7 ). There is also a subtlety about which pairs of bonds should be counted in this energy. If both bonds are very weak, just above the bond order cutoff, then it would be a waste of computational time to compute bond bending energies. Thus another cutoff and rescaling is included in the bond orders for valence angle energy calculations, similar to that done earlier for bond energies. This only makes a small change to the energy calculations.
A.5. Bond bending energy The bond bending energy calculations must first determine the optimum bond angle, Θ0 , based on bond orders, because atoms with different hybridisations prefer different geometries. This is achieved through the equations: Θ0 = π − Θ0,0 {1 − exp[−λ18 (2 − SBO2)]}
A.6. Over and under-coordination energies The final major term in the ReaxFF energies for this study is due to over or under coordination. The under-coordination energy is expressed as:
(A.17) Eunder,i = −punder
where
SBO2 =
⎧ ⎪ ⎪ ⎪ ⎨
0,
(1 − exp(λ7 i )) (1 + exp(−Ai,25 i ))
f6 (SOVi ) (A.23)
if SBO < 0
SBOλ17 ,
if SBO > 0 and SBO < 1
⎪ 2 − (2 − SBO)λ17 , ⎪ ⎪ ⎩ 2,
where
if SBO > 1 and SBO < 2 if SBO > 2
f6 (SOVi ) =
1 1 + λ9 exp(λ10 SOVi )
(A.24)
(A.18) and
⎡
SBO = ⎣1 −
neighb(j) n=1
+
neighb(j)
and
⎤ exp(−BO8jn )⎦ (j − λ34 lp,j )
BOjn,π
SOVi =
neighb(i) n=1
(A.19)
(A.25)
The over-coordination energy is expressed as: neighb(i)
n=1
Note the corrected bracketing in Eq. (A.19) when compared to Eq. (7c) of van Duin et al. [19]. j is the sum of the bond orders from atom j minus the valency of the atom, minus the optimal number of lone pairs. Once Θ0 has been found, the bond bending or “valence angle” energy can be calculated. For the angle Θijk formed between atoms i, j and k, centred on j, this is expressed as:
(BOn − Vn − lp,n )(BO n + BOn )
Eover,i =
n=1 ov,i
ov,i pbe,3 Dbe,1 BOij + Vi + 10−8 1 + exp(Ai,25 ov,i ) (A.26)
where ov,i =
neighb(i) −lp,i BOni − Vi + 1 + λ33 exp(λ32 SOVi ) n=1
(A.27)
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
133
Appendix B. Reaction force field parametrisation Tables B.1–B.6 Table B.1 General parameters for the reaction force field Parameter
Value
Description
λ1 λ2 λ3 λ4 λ5 λ7 λ9 λ10 λ12 λ13 λ15 λ16 λ17 λ18 λ29 λ30 λ32 λ33 λ34
50.0000 5.5469 21.2839 3.0000 6.5000 1.0159 8.0878 13.0238 0.0000 10.0000 33.8667 25.6125 1.0563 2.0384 1.5591 0.1000 1.6356 5.6937 2.5067
Over-coordination parameter Over-coordination parameter Valency angle conjugation parameter Triple bond stabilisation parameter Triple bond stabilisation parameter Under-coordination parameter Under-coordination parameter Under-coordination parameter Lower Taper-radius Upper Taper-radius Valency undercoordination Valency angle/lone pair parameter Valency angle Valency angle parameter van der Waal’s shielding Cutoff for bond order (×100) Over-coordination parameter Over-coordination parameter Valency/lone pair parameter
These apply to the whole system, and are not element specific. All were previously determined [19,20], and have not been fitted in this study.
Table B.2 Atom parameters Index
Symbol
Si
N
O
Description
Ai,1 Ai,2 Ai,4 Ai,5 Ai,6 Ai,7 Ai,8 Ai,9 Ai,10 Ai,12 Ai,14 Ai,15 Ai,17 Ai,18 Ai,25 Ai,26 Ai,29
r0 V rvdW D γEEM r0 ne ␣ γw −punder χ η r0 plp Ai,25 Ai,26 Ai,29
2.0132 4.0000 2.2203 0.1306 0.8218 1.5629 4.0000 12.1305 1.8838 11.4056 1.8038 7.3852 −1.0000 0.0000 −3.6886 2.9867 2.5791
1.2360 3.0000 1.9333 0.1373 0.8596 1.1602 5.0000 10.0725 8.0467 34.5504 6.8418 6.3404 1.0365 14.7406 −15.0000 2.6491 2.8793
1.1693 2.0000 1.9580 0.0875 1.0804 1.0200 6.0000 10.1761 7.2516 39.2737 8.5000 7.8386 0.9809 18.4195 −7.3000 2.6656 2.9225
cov.r Valency Rvdw Evdw gammaEEM cov.r2 # alfa gammavdW Eunder chiEEM etaEEM cov.r3 Elp ov/un val1 vval4
The index is from the position this variable takes in the structured input to the ReaxFF program by van Duin [21]. Some parameters take this index as a symbol because they were instead denoted as global parameters in previous publications [19,20].
Table B.3 Bond parameters for Si Si, Si N, and Si O Index
Symbol
Si Si
Si N
Si O
Description
Bi,1 Bi,2 Bi,3 Bi,4 Bi,5 Bi,6
De De De pbe,1 pbo,5 ovc13
113.7903 53.9894 30.0000 0.2535 −0.3000 1.0000
144.4831 20.4624 10.0000 − 0.3794 − 0.3000 1.0000
193.1177 41.1424 43.3991 −0.2085 −0.3000 1.0000
Edis1 Dbe1 LPpen Dbe2 n.u. Dbe3 pbe1 pbo5 13corr
134
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
Table B.3 (Continued ) Index
Symbol
Si Si
Si N
Si O
Description
Bi,7 Bi,8 Bi,9 Bi,10 Bi,11 Bi,13 Bi,14 Bi,15
pbo,6 pbe,3 pbe,2 pbo,3 pbo,4 pbo,1 pbo,2 ovc
16.0000 0.0742 0.2586 −0.1963 7.5407 −0.0693 7.9365 0.0000
36.0000 1.5282 3.2297 −0.1012 12.2585 −0.1039 8.1682 1.0000
36.0000 0.7695 0.9220 −0.3683 4.2644 −0.5191 4.4529 1.0000
pbo6 pbe3 pbe2 pbo3 pbo4 pbo1 pbo2 ovcorr
The index is from the position this variable takes in the structured input to the ReaxFF program by van Duin [21].
Table B.4 Bond parameters for N O, N N, and O O Index
Symbol
N O
N N
O O
Description
Bi,1 Bi,2 Bi,3 Bi,4 Bi,5 Bi,6 Bi,7 Bi,8 Bi,9 Bi,10 Bi,11 Bi,13 Bi,14 Bi,15
De De De
147.4287 109.4565 61.7111 0.9507 −0.1616 1.0000 21.7802 0.7500 0.7441 −0.3665 7.2396 −0.2439 4.8277 1.0000
155.9636 0.0000 0.0000 0.2592 −0.1418 1.0000 13.1260 0.3000 0.3989 −0.1336 10.7257 −0.0760 5.2866 1.0000
118.8570 42.8420 167.6132 0.9092 −0.1587 1.0000 10.5117 0.8589 0.8634 −0.2501 6.5155 −0.1613 5.6591 1.0000
Edis1 Dbe1 LPpen Dbe2 n.u. Dbe3 pbe1 pbo5 13corr pbo6 pbe3 pbe2 pbo3 pbo4 pbo1 pbo2 ovcorr
pbe,1 pbo,5 ovc13 pbo,6 pbe,3 pbe,2 pbo,3 pbo,4 pbo,1 pbo,2 ovc
The index is from the position this variable takes in the structured input to the ReaxFF program by van Duin [21]. Table B.5 Off-diagonal parameters Bond
D
rvdW
γw
r0
r0
r0
O Si N Si
0.1233 0.1536
1.7552 1.7871
10.6931 11.0924
1.6096 1.7610
1.2935 1.2495
N/A N/A
If available for a bond type, these override values based on elemental properties.
Table B.6 Bond bending parameters Type O O N O O N Si O O Si O N N O Si N O N
O O O N N N Si Si Si O O Si Si Si N N N O
O N N O N N Si Si O Si Si Si N N Si Si Si Si
θ0,0
ka
kb
pv,2
pv,7
78.4963 79.0601 74.5183 74.8323 77.6183 75.2419 69.3456 70.3016 90.0344 22.1715 73.4663 1.4881 83.1400 77.4641 34.4338 125.3711 69.8728 69.8728
61.0192 38.7432 54.0833 35.2858 29.5041 20.4666 21.7361 15.4081 7.7656 3.6615 25.0761 4.0912 0.0562 4.5724 40.7839 12.9852 32.7155 27.1273
1.0093 2.1471 1.9745 1.4738 1.2928 3.1280 1.4283 1.3267 1.7264 0.3160 0.9143 13.6686 2.2347 1.0849 0.4583 179.7577 1.5875 1.5875
0.4897 0.4897 0.4897 2.7993 2.7993 2.7993 −0.2101 2.1459 0.7689 4.1125 2.2466 12.0668 0.0079 0.7689 6.5339 18.4450 2.2466 2.2466
1.0400 1.0400 1.0400 1.0400 1.0400 1.0400 1.3241 1.0400 1.0400 1.0400 1.0400 1.0400 1.0400 1.0400 1.0400 1.0400 1.0400 1.0400
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
References [1] P.F. Becher, G.S. Painter, E.Y. Sun, C.-H. Hsueh, M.J. Lance, Acta. Mater. 48 (2000) 4493. [2] G.S. Painter, P.F. Becher, H.-J. Kleebe, G. Pezzotti, Phys. Rev. B 65 (2002) 064113. [3] S. Ii, C. Iwamoto, K. Matsunaga, T. Yamamoto, M. Yoshiya, Y. Ikuhara, Phil. Mag. 84 (2004) 2767. [4] C.-M. Wang, X. Pan, M.J. Hoffman, R.M. Cannon, M. R¨uhle, J. Am. Ceram. Soc. 79 (1996) 788. [5] I. Tanaka, H.-J. Kleebe, M.K. Cinibulk, J. Bruley, D.R. Clarke, M. R¨uhle, J. Am. Ceram. Soc. 77 (1994) 911. [6] H. Gu, J. Am. Ceram. Soc. 85 (2002) 33. [7] M. Doeblinger, private communications. [8] H. Gu, in: Interface Science and Materials Interconnection, Silicon Nitride Based Ceramic Materials, Proceedings of JIMIS08. Japan Institute of Metals, Sendai, Japan, (1996), pp. 217–220. [9] A. Ziegler, C. Kisielowski, M.J. Hoffmann, R.O. Ritchie, J. Am. Ceram. Soc. 86 (2003) 1777. [10] G.B. Winkelman, C. Dwyer, T.S. Hudson, D. Nguyen-Manh, M. Doblinger, R.L. Satet, M.J. Hoffmann, D.J.H. Cockayne, Phil. Mag. Lett. 84 (2004) 755. [11] W. McBride, D.J.H. Cockayne, D. Nguyen-Manh, Ultramicroscopy 96 (2003) 191. [12] M. Doblinger, C.D. Marsh, D. Nguyen-Manh, D. Ozkaya, D.J.H. Cockayne, Inst. Phys. Conf. Ser. 179 (2004) 401.
135
[13] D. Nguyen-Manh, D.J.H. Cockayne, M. Doeblinger, C. Marsh, A.P. Sutton, A.C.T. van Duin, Inst. Phys. Conf. Ser. 179 (2004) 139. [14] M. Yoshiya, K. Tatsumi, I. Tanaka, H. Adachi, J. Am. Ceram. Soc. 85 (2002) 109. [15] S.H. Garofalini, W.W. Luo, J. Am. Ceram. Soc. 86 (2003) 1741. [16] X.T. Su, S.H. Garofalini, J. Mater. Res. 19 (2004) 752. [17] D.J. Adams, Mol. Phys. 28 (1974) 1241. [18] D.J. Adams, Mol. Phys. 29 (1975) 307. [19] A.C.T. van Duin, A. Strachan, S. Stewman, Q. Zhang, X. Xu, W.A. Goddard III, J. Phys. Chem. A 107 (2003) 3803. [20] A.C.T. van Duin, S. Dasgupta, F. Lorant, W.A. Goddard III, J. Phys. Chem. A 105 (2001) 9396. [21] A. C. T. van Duin, ReaxFF molecular dynamics program. [22] S. H. Garofalini, private communications. [23] X. Su, S.H. Garofalini, J. Mater. Res. 19 (2004) 752. [24] W.-Y. Ching, private communications. [25] H. Gu, R.M. Cannon, H.J. Seifert, M.J. Hoffmann, I. Tanaka, J. Am. Ceram. Soc. 85 (2002) 25. [26] M.P. Allen, D.J. Tildesley, Computer Simulations of Liquids, Clarendon Press, Oxford, UK, 1987, pp. 39–42. [27] M. Kitayama, K. Hirao, A. Tsuge, M. Toriyama, S. Kanzaki, J. Am. Ceram. Soc. 82 (1999) 3263. [28] V. Burlakov, private communications. [29] W.J. Mortier, S.P. Ghosh, S. Shankar, J. Am. Chem. Soc. 108 (1986) 4315.
Grand canonical Monte Carlo simulations of intergranular glassy films in  silicon nitride T.S. Hudson a,∗ , D. Nguyen-Manh b , A.C.T. van Duin c , A.P. Sutton a,d,1 a Department of Materials, Parks Road, Oxford OX1 3PH, UK EURATOM/UKAEA Fusion Association, Culham Science Centre, Abingdon OX14 3DB, UK c Beckman Institute (139-74), California Institute of Technology, Pasadena, CA 91125, USA Helsinki University of Technology, Laboratory of Computational Engineering, PO Box 9203, FIN-02015 HUT, Finland b
d
Received 14 February 2005; received in revised form 17 November 2005
Abstract Grand canonical Monte Carlo simulations have been performed on intergranular glassy films in  silicon nitride to determine their equilibrium structure, width and composition at 2000 K. A reactive force field for Si O N systems has been developed and tested, and chemical potentials of oxygen and nitrogen were obtained by requiring that they reproduce, in conjunction with the reactive force field, the experimentally determined compositions of the glass pockets at triple junctions and the interior of the adjoining  silicon nitride grains. Two silicon oxynitride films developed during the course of very long simulations, at the end of which their composition and structure were analysed. The thermodynamic equilibration of the film thickness was explored to give insight into the experimental uniformity of these films. Incorporation of a significant concentration of oxygen into the adjoining Si3 N4 grains was observed, producing somewhat diffuse interfaces on either side of the film. We believe this is the first thermodynamically consistent, atomic-scale, simulation of intergranular glassy films capable of exchanging matter and heat with reservoirs. © 2006 Elsevier B.V. All rights reserved. Keywords: Amorphous; Silicon nitride; Thin film; Interatomic potential; Monte Carlo; Grand canonical ensemble
1. Introduction Polycrystalline  silicon nitride (Si3 N4 ) is a high temperature structural ceramic with high strength and fracture toughness. The relatively high toughness of this material is derived from the fact that fracture occurs primarily intergranularly with a long crack path. The three factors that cause this are the toughness of individual grains, the large aspect ratios of those grains, and the weaker intergranular interfaces through which cracks can propagate [1–3]. The weak interfaces consist of silicon oxynitride glasses which form an intergranular glassy films (IGF) from the oxide surface films during the liquid-sintering process, and are often also modified using sintering additives such as rare-earth oxides [4]. One striking feature of these interfaces is their unifor-
∗ Corresponding author at: School of Chemistry, Bldg. F11, University of Sydney, NSW 2006, Australia. Tel.: +61 2 93514813; fax: +61 2 93513329. E-mail addresses: [email protected] (T.S. Hudson), [email protected] (D. Nguyen-Manh), [email protected] (A.C.T. van Duin), [email protected] (A.P. Sutton). 1 Present address: Department of Physics, Imperial College, London SW7 2AZ, UK.
0921-5093/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.msea.2006.01.014
mity, forming with a constant thickness throughout the sample. In the case of pure Si3 N4 –SiO2 , this thickness is observed to be 1.0 ± 0.1 nm [5], although another study has reported IGF thicknesses ranging from 0.6 to 3.0 nm for a material containing very low quantities of SiO2 [6]. The film thickness is almost independent of grain orientation, but statistical studies have shown that the most common planes for grain termination are those near the {1 0 1¯ 0} prism planes [7]. Great progress has been made in probing the composition and structure of the IGF experimentally. Using electron energy loss spectroscopy (EELS) the film has been determined to contain an anionic ratio for nitrogen of [N]/([N] + [O]) of about 30% [8]. High-resolution transmission electron microscopy has directly imaged the region immediately adjacent to the grains, showing the structural influence of these crystals into the film [9]. Similarly in materials prepared using rare-earth elements in the sintering glasses, high-angle annular dark field (HAADF) experiments have shown the sites these elements occupy in the IGF are strongly influenced by the presence of the crystal, often out to a few atomic layers into the glass [10]. Work has also been done using coherent convergent electron microscopy [11,12] to understand the short and medium
124
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
range order within the amorphous layer via radial distribution functions [12,13]. Modelling, in particular molecular dynamics (MD), has also contributed to our understanding of the detailed structure within these IGFs. One example demonstrated the role of N within the IGF, and the competition between the detrimental energy of mixing and the advantageous role this N plays in reducing the number of dangling bonds at the edges of the IGF [14]. Some studies have investigated segregation of metal ions within the IGF, and the decay of the crystallinity going into the film [15,16]. Other studies investigated the radial distribution functions and composition profiles of the film [13]. A common problem in these studies however is the need for both the composition and thickness of the IGF to be predetermined, because MD is a closed-system method using a fixed set of atoms. In this study, we use the grand canonical Monte Carlo simulation method [17,18] which is able to overcome this problem of predetermined IGF composition and thickness. It allows composition to vary in a thermodynamically correct fashion during the simulation. Chemical potentials are determined, which govern the exchange of material with reservoirs, and so the IGF is able, at least in principle, to establish the equilibrium thickness and composition. In reality, the IGFs are able to exchange matter with real reservoirs, such as the intergranular pockets and the adjoining crystalline regions, in order to establish equilibrium composition and thickness at a given temperature and pressure. 2. Reactive force field for Si O N systems Reactive force fields (ReaxFF) are many-body empirical potentials that seek to extend the force fields widely used for simulations of large molecules, to allow the breaking and forming of bonds [19,20]. They are then appropriate for simulations of condensed network materials, with the capacity to treat differently coordinated polymorphs correctly [19]. They do this by including a large number of physically motivated terms in the interaction energy, and by careful interrelation of these terms. Because these terms cover a wide range of physical effects, a single well-parametrised ReaxFF should be applicable to simulations of a wide range of chemical and material systems. The functional form of the ReaxFF is described in Appendix A. We have followed the notation of the previous descriptions and implementations of the ReaxFF [19–21]. In cases where they are conflicting we have followed the most recently used definitions, unless otherwise noted. 2.1. Parametrising the reactive force field ReaxFF parameters for Si O were previously developed and tested [19]. An attempt has also been made at Si O N parametrisation [13], but it was fitted almost exclusively to data from single molecules and so the applicability to condensed phase systems was unknown. Errors in the earlier Si O N parametrisation [13] were identified in an MD study of silicon nitride IGFs. We constructed a model IGF configuration using this parametrisation by annealing and then cooling a thin slab of silicon dioxide between 
silicon nitride crystals. A second model IGF was constructed by further annealing this IGF configuration [22] using a simple empirical potential (using potentials similar to those of Su and Garofalini [23]). The energies and stabilities of both configurations were compared [24] using density functional theory (DFT). The ReaxFF potential produced a configuration that had a much higher energy according to DFT as compared with the structure predicted by the simpler empirical potential. Further, DFT relaxations of both structures showed that the ReaxFF structure underwent a much more extensive relaxation than that predicted by the simpler potential, and despite this, the relaxed energy was still significantly higher. One reason for this poor performance was the preference of the ReaxFF parametrisation for an excessive concentration of N2 dimers to form within the -Si3 N4 crystal, resulting in under-coordinated Si atoms. We now discuss further inadequacies of the parametrisation developed in [13]. The ReaxFF parameters had been fitted [13] to a training set of Si O N configurations obtained by quantum chemical relaxations which included -Si3 N4 at and around the equilibrium volume under atmospheric pressure, and had very successfully reproduced the bulk modulus (see Fig. 1a). However, upon further investigation, we observed an instability to a fictitious new structure at higher volumes. This fictitious new configuration was calculated according to that ReaxFF parametrisation to have a lower energy than -Si3 N4 , even at the original volume. Thus according to the parametrisation of [13], the  phase is metastable, which is physically not the case. It was concluded that the instability was an indication that the ReaxFF parametrisation of [13] predicts unphysical results. The cause of the instability was elucidated by close examination of the key ReaxFF terms at the instability. The bond bending terms predicted an optimum N Si N bond angle of 155◦ , rather than the 109◦ of the tetrahedral configuration. This, in turn, was caused by a large overestimation of the character of Si N bonds (which are in reality primarily simple single bonds in this four-coordinated-Si network solid). The fitting of the earlier parametrisation [13] had predicted significant bonding at bond lengths that were as large as the onset of significant bonding. Nevertheless, our earlier parametrisation provided a good starting point for a new parameter fitting. The parameters governing the bond lengths at which bonding became significant were manually altered. Additionally, the N N and double bonding were prohibited because they caused excessive N N dimerisation in these solid-state systems when compared with experiment [25]. Although in consequence our new parametrisation is not transferable to molecular systems involving double or triple-bonded nitrogen, it provides greater accuracy for the current system under study. Torsional terms in the ReaxFF were also turned off, because although they contributed only a small energy relative to most other terms, they occasionally gave spuriously low energies during MC simulation, invalidating the run. Removing these terms significantly reduced the number of free parameters in the fit. The earlier training set [13] was augmented by the DFT energies of the silicon nitride IGFs, and those of crystalline ␣-Si3 N4 , -Si3 N4 , ␣-Si2 N2 O, and -Si2 N2 O. Fitting was performed using a combination of a genetic algorithm and a parameter-by-
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
125
Fig. 1. (a) Energy–volume curves for -Si3 N4 predicted by the original parametrisation: red — rigid expansion of equilibrium cell; black — expansion with atomic coordinates relaxed, showing instability at high volumes to a new fictitious structure; green — fictitious structure expanded with atomic relaxations, note that this fictitious structure has a lower energy than the  structure at all volumes. (b) Energy–volume curves predicted by the new parametrisation: red — rigid expansion of the equilibrium cell; black — expansion with atomic coordinates relaxed, showing no instability. The small difference between red and black curves at the ground state indicates that the equilibrium structure predicted by this ReaxFF parametrisation is slightly different to the DFT structure included in the training set.
parameter gradient relaxation. Only parameters specifically relating to Si N and N O bonding (and valence angle bond bending with N involved) were varied. The original general (speciesindependent) ReaxFF parameters (Table B.1) were kept fixed, as were the Si O parameters. This meant that a total of 33 parameters were varied. The final values of all parameters are shown in Appendix in Tables B.1–B.6, with the parameters allowed to vary indicated in bold. Once the fit converged with reasonable accuracy, the same volume expansion test as in Fig. 1a was applied to the new parametrisation. As shown by Fig. 1b, the  structure predicted by DFT is reproduced at all volumes except for a very minor relaxation. Furthermore, these parameters reproduce the ordering of the four IGF structures added to the training set as predicted by DFT calculation [24].
energies with the difference between the initial total energy and the final total energy. This was confirmed to a high degree of precision. The calculations of atomic charges and hence electrostatic energy could not be streamlined in this way. Electrostatic interactions are much longer range, and for every accepted MC move a recalculation involving the entire system is required. An approximate solution to this problem was implemented, whereby after each accepted MC move a region out to the electrostatic cutoff radius was recalculated, saving significant computational time. After every 50 MC attempts, charges on the entire system were recalculated to remove any residual error that had built up. The frequency of this reset was chosen so that average errors in the total energy after this number of steps were still significantly less than the thermal energy kT.
2.2. Testing the implementation of the reaction force field
3. Grand canonical Monte Carlo (GCMC) simulations
To test that the potential was implemented correctly in our Monte Carlo simulation program, the total energies of large Si N O IGF systems were compared to the energies calculated for the same systems by a previous molecular dynamics implementation of the ReaxFF [21,19]. For the same ReaxFF parametrisation, these two codes gave identical results apart from tiny numerical errors on the scale of computational floating-point accuracy. To ensure that the energy of the entire system did not have to be recalculated after each individual Monte Carlo (MC) move, “incremental” energy functions were introduced, which calculated only local changes in each term of the energy following each type of MC move. Because of the many-body nature of many terms of the ReaxFF, a change to a single atom often affects the energy of a number of nearby atoms, sometimes as far as the fourth nearest neighbour. The only exception to this is the electrostatic energy, which will be discussed later. Recalculating the changes to this limited region is computationally much more efficient than recalculating the energy of the entire system. The incremental local energy functions were tested by performing many MC moves on a system and comparing the sum of the incremental changes of local
3.1. Implementation Six types of MC move were included to alter the structure. Three are standard MC moves: translation of a single atom, swapping random pairs of unlike atoms, and swapping atoms with unlike neighbours. The other three are only applicable to the grand canonical ensemble, because they alter the composition of the system: insertions, deletions, and transmutations. The standard MC moves all have a straightforward acceptance criterion. The Boltzmann factor exp(− E kT ) is calculated, where E is the change in potential energy as a result of the move, k the Boltzmann constant, and T is the temperature of the system. If this quantity is greater than 1, the move is immediately accepted. Otherwise the quantity is treated as a probability of the move being accepted, and a pseudo-random number is generated to determine whether or not this occurs. The moves to simulate the grand canonical ensemble have a similar acceptance procedure, but care must be taken to satisfy detailed balance in this ensemble. The probability NA →NA +1 of inserting an atom of type A into a configuration with NA A atoms, NB B atoms, NC C atoms, . . . (NA + 1, NB , NC , . . .) multiplied by the occupation probability NA of this configura-
126
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
tion, must be equal to the probability NA +1→NA of removing an atom of type A from the configuration (NA + 1, NB , NC , . . .) multiplied by the occupation probability NA +1 of this configuration: NA +1 NA →NA +1 = (1) NA NA +1→NA The occupation probabilities are known from the grand canonical partition function [26] to have the ratio: NA +1 1 µA − E V = × 3 × exp (2) NA NA + 1 ΛA kT where V is the volume of the system, µA the chemical potential of A atoms, and ΛA denotes the thermal de Broglie wavelength of each atom of type A, given by the expression: ΛA = √
h 2πmA kT
(3)
where h is Planck’s constant, and mA is the mass of an atom of type A. The probability NA →NA +1 of inserting an atom of type A is equal to the probability of choosing to attempt an insertion, Pins , multiplied by the probability of choosing type A, 1τ , where τ is the number of elements involved in the simulation, multiplied by the acceptance probability min(1, χ). χ is the equivalent of the Boltzmann acceptance factor exp(− E kT ) for standard MC moves, and we derive it for the multicomponent GCMC moves here. The probability NA +1→NA of deleting an atom of type A from the (NA + 1, NB , NC , . . .) configuration is equal to the probability of choosing to attempt a deletion Pdel , multiplied by NA +1 the probability of choosing an atom of type A to delete, N+1 where N = tτ Nt (since we select the atom at random), multiplied by the probability of acceptance, min(1, χ1 ). Therefore, detailed balance (Eq. (1)) requires: 1 V µA − E exp NA + 1 Λ3A kT =
N +1 1 Pins (1/τ) min(1, χ) =χ Pdel (NA + 1/N + 1) min(1, (1/χ)) NA + 1 τ
(4)
Therefore,
τV µA − E exp kT (N + 1)Λ3A τ µ ˜ A − E = exp N +1 kT
χ=
(5)
where we have chosen to set Pins = Pdel . In the final step we have incorporated the factor ΛVA into the chemical potential to make it an effective chemical potential µ ˜ A = µA + kT ln( V3 ). ΛA
This is appropriate as long as our simulations are performed at a fixed volume and temperature. Eq. (5) implies that when considering an insertion, the accepτ ins tance probability is calculated as N+1 exp( µ˜ A −E ). For a delekT tion in a configuration, in which the total number of atoms is initially N + 1, the acceptance probability χ1 can be used. However, we must replace the factor of N + 1 with N if there are N atoms
in the initial configuration from which one is to be deleted. Also note that Edel (N → N − 1) = −Eins (N − 1 → N). So we have an acceptance probability for the deletion step N → N − 1 del of Nτ exp( −µ˜ A −E ). kT A similar argument for the transmutation move type leads to the conclusion that the probability of accepting an atom attempting to change from type A to type B should simply be −EA→B exp ( µ˜ B −µ˜ AkT ). These grand canonical steps can change the system’s composition, so there is a natural question about charge conservation in partially ionic systems such as silicon oxynitrides. A system with a net charge would cause major problems if the periodic images of the simulation cell were allowed to interact correctly. Because the ReaxFF uses a variable charge method to calculate its electrostatics, we are entitled to consider the overall system as always neutral, as long as we only ever insert neutral atoms. The charge will then distribute itself to minimise the total energy. Of course the most favourable energies will be found at stoichiometric compositions, so the system should remain near neutrality. 3.2. Chemical potentials We now turn to the task of determining the effective chemical potentials to use in the IGF simulation. When the IGF is at equilibrium with the glassy pockets at the triple junctions, and with the grains themselves, the effective chemical potentials of silicon, oxygen and nitrogen are constant throughout the entire systems. Because it is easier to extract accurate experimental information about the bulk phases, we may tune our chemical potentials to match experimental data for those phases, and apply these same chemical potentials to the IGF. The target concentrations of ([Si] = 18.48 nm−3 , [N] = 1.44 nm−3 , [O] = 34.78 nm−3 ) in the glassy pockets were based on EELS measurements [25] which also determined that the nitrogen was incorporated into the silicate network rather than being present as molecular N2 . In the  silicon nitride crystal, the corresponding concentrations of ([Si] = 40.6 nm−3 , [N] = 53.6 nm−3 , [O] = 0.56 nm−3 ) were based on the crystal lattice parameters, together with experimental measurements of the oxygen concentration using hot-gas extraction method on liquid-phase sintered -Si3 N4 [27]. Simulations were carried out of both the ␣-SiO2 phase containing nitrogen, and the -Si3 N4 phase containing oxygen. The initial ␣-SiO2 configuration was prepared by Burlakov [28] and contained only silicon and oxygen atoms. An algorithm was adopted whereby if the concentration of any one element was drifting away from the target composition, µ ˜ was adjusted slightly to counteract the drift. The amount by which it was adjusted diminished as the simulation went on, so that eventually convergence was achieved. The values established are shown in Table 1. The value of µ ˜ Si could not be evaluated in the -Si3 N4 simulation, owing to the extremely high energy penalty for creating any Si defect. Thus only our value of µ ˜ Si for the glass is of significance. Assuming the experimental target compositions
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
127
Table 1 Effective chemical potentials at 2000 K ␣-SiO2 phase -Si3 N4 phase Used for IGF
µ ˜ Si (eV)
µ ˜ N (eV)
µ ˜ O (eV)
4.5 N/A 4.5
2.9 1.6 1.9
1.5 2.2 1.9
Fig. 4. (a) Starting configuration of run 2 with a 1.325 nm gap between prism faces of a single crystal of -Si3 N4 . (b) Snapshot of final configuration after 1.44 × 107 GCMC moves.
3.3. Interfacial simulations
are in fact at equilibrium with one another, the differences in the values of µ ˜ N and µ ˜ O for the ␣-SiO2 and -Si3 N4 phases indicate the extent to which our model potential is not thermodynamically consistent. Using the effective chemical potentials obtained from the ␣-SiO2 simulations in a simulation of the Si3 N4 phase produces an excess of O atoms. Similarly using the effective chemical potentials obtained from the -Si3 N4 phase in the ␣-SiO2 phase produces an excess of N atoms. This implies that our potential fails to penalise sufficiently defects in one or both phases, indicating that there is further improvement possible for our potential. For our further simulations we used the intermediate values shown in Table 1, which after equilibration in the -Si3 N4 phase produce an anionic ratio for oxygen of 3%, and after equilibration in the ␣-SiO2 phase produce an anionic ratio for nitrogen of 12%.
In order to simulate a grain boundary, a -Si3 N4 crystal was cleaved at the plane most commonly observed experimentally to terminate Si3 N4 crystals. A gap of 0.3 nm was introduced by separating the two halves of the cleaved single crystal, as shown in Fig. 2a. Because the crystal period in this direction ˚ it is not possible to restore this system to a peris 13.25 A, fect crystal by filling this gap with additional material. Thus the system is forced to accommodate a planar interface. Periodic boundary conditions were imposed both parallel and perpendicular to the interface, and the sides of the periodic supercell (5.59820 nm × 1.52946 nm × 0.87513 nm) were not allowed to change. Our GCMC simulation was then run to determine the equilibrium composition and thickness of this interface. Fig. 2b shows a snapshot of the simulation after 2.64 × 107 MC moves. The system quickly filled the gap, but was unable to generate a well connected bond topology, and so gradually expanded the amorphous oxide layer to relieve this misfit. By tracking the energy and composition of the system (as shown in Fig. 3), we established that this system had not reached an equilibrium by the end of 2100 h of CPU time on a Pentium 4 single-processor machine. A range of interfaces can be generated in this manner, each with a different initial gap between 0 and the crystal period. To
Fig. 3. Cell composition throughout run 1.
Fig. 5. Cell composition throughout run 2.
Fig. 2. (a) Starting configuration of run 1 with a 0.3 nm gap between prism faces of a single crystal of -Si3 N4 . (b) Snapshot of final configuration after 2.64 × 107 GCMC moves.
128
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
Despite the failure of run 2 to close the gap and return to perfect crystal, when compared to run 1 this simulation did show significant differences. We note that the IGF did not penetrate any further into the crystal than the width of the initial gap. The concentration of O defects in the surrounding Si3 N4 is also noticeably lower than in the previous simulation, and does not continue increasing, possibly suggesting a weaker elastic field. Thus in contrast to the run 1, we conclude that the film thick-
Fig. 6. (a) Average composition profile of the IGF during the final period of run 1. (b) Average cumulative frequency of oxygen through the IGF during the final period of run 1.
explore the behaviour of a system which began with a wider gap, in run 2 we chose to study a gap of exactly one crystal period (as shown in Fig. 4a). This has the advantage that the global free energy minimum structure is known — it fills the gap with perfect crystal, leaving almost no defects, and the bulk concentration of oxygen. Thus if our GCMC could explore configuration space very efficiently, it should converge to this structure. However, this did not happen. Just as metastable amorphous materials are unable to crystallise on a physical timescale at low temperatures, our GCMC was unable to find the global minimum within the simulation time of 1200 CPU hours. However, the simulation did converge to a metastable configuration within this time (see Fig. 5). Fig. 4b shows a final snapshot of the simulation. The periodic cell had dimensions of 6.62275 nm × 1.52946 nm × 0.87513 nm.
Fig. 7. (a) Average composition profile of the IGF during the final period of run 2. (b) Average cumulative frequency of oxygen through the IGF during the final period of run 2.
Fig. 8. Average coordination of each type of atom as a function of coordination radius during the final period of run 1. A point f on the line m at radius r means that a fraction f of that type of atom have at least m neighbours within the radius r. (a) silicon; (b) nitrogen; (c) oxygen.
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
129
ness is sufficient to accommodate a well connected oxynitride network, and that any further gains that might be obtained by making it thicker are outweighed by the energetic cost of disrupting crystal. The final period of each run was analysed to extract important physical information about the IGF. Fig. 6a shows the average composition profile through the IGF in run 1. Fig. 6b shows the
Fig. 10. Partial radial distribution functions in the IGF only during the final period of run 1.
average cumulative frequency of oxygen. The film thickness is considered to be the distance between the 10th percentile and the 90th percentile of this distribution (after correcting for the oxygen that is clearly in the crystalline phase), and at the end of ˚ We note that there is substantional this run it is calculated at 16 A. oxygen penetration into the -Si3 N4 phase. Fig. 7 is the equivalent data extracted from run 2. This data gives a metastable film thickness of 1.4 nm. The topology of the system was also studied, specifically with regard to the types of defects present. Fig. 8 shows the fraction of atoms with a particular number of neighbours within a given radius in run 1. Fig. 9 shows the equivalent for run 2. The radius at which first nearest neighbours are typically situated is clearly seen by the first rapid increase in the number of neighbours. At a radius slightly larger than this onset, nearly all silicon atoms have four neighbours, nitrogen atoms have three neighbours, and oxygen atoms have two neighbours, consistent with the valency of each element. A small fraction of the silicon atoms have three neighbours at this radius, and another small fraction have five, indicating the presence of a small quantity of over- and under-coordinated silicon defects. Nitrogen does not show any
Fig. 9. Average coordination of each type of atom as a function of coordination radius during the final period of run 2. A point f on the line m at radius r means that a fraction f of that type of atom have at least m neighbours within the radius r. (a) silicon; (b) nitrogen; (c) oxygen.
Fig. 11. Partial radial distribution functions in the IGF only during the final period of run 2.
130
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
under-coordinated defects, but there are a small fraction of these atoms which have four neighbours and are thus overcoordinated. Oxygen is an interesting case, showing a large fraction of the oxygen atoms to have three neighbours. The distance to the third neighbour varies. The large quantity of over-coordinated oxygen atoms is consistant with the observation that many of them substitute at nitrogen sites in the -Si3 N4 lattice (Figs. 2b and 4b). We have compiled average radial distribution functions (RDFs) for each pair of elements. Because RDFs within the Si3 N4 crystal are not of interest, only atoms within the IGF have been included in the analysis. Normalisation based on average density is performed as usual, however, because we analyse only a thin film, the large-radius asymptote is not to 1 as per normal. Figs. 10 and 11 show the results of these calculations for run 1 and 2, respectively. As expected, these results show a network primarily consisting of bonds between silicon and nitrogen, or silicon and oxygen. Some bonds between pairs of anions are also apparent, which are the primary defects within the IGF layer (bonds between silicon and silicon are not observed).
defects which would show little difference when directly imaged in an electron microscope. Thus unless elementally sensitive techniques such as EELS are used, care should be taken when making conclusions about interfaces that appear to be sharp. Acknowledgements The authors wish to thank Victor Burlakov for providing a simulated amorphous SiO2 cell. Thanks also to Stephen Garofalini and Wai-Yim Ching for their demonstration that our earlier ReaxFF parametrisation [13] needed improvement. Financial support from the European Community Growth Program (Contract No. G5RD-CT-2001-00586) is gratefully acknowledged. There are no competing interests in this work. Appendix A. Reaction force field functional form The terms used in the ReaxFF developed and used in this study are given by Esystem = Ebond + EvdW + ECoulb + Eqdist + Eval + Eover/under + Elp
4. Conclusions We have provided an improved parametrisation for the Si N O ReaxFF. Although not transferrable to systems including molecular N2 , this potential has been successfully applied to condensed Si N O systems. There is room for improvement however, as demonstrated by the thermodynamic inconsistency shown in Section 3.2 when determining chemical potentials in phases which appear experimentally to be in equilibrium. We suggest that further improvements could be obtained by including O ↔ N substitutional defect energies in the fit. We note that the requirement of thermodynamic consistency, in the sense that the chemical potential is constant for each element in phases that are in equilibrium with each other at a given temperature and pressure, is an extremely demanding test of any description of atomic interactions. It is a test that is rarely applied. We have developed and applied a GCMC scheme to simulate thermodynamic equilibrium in IGFs in the Si N O system. However, after many hours of simulation there remains some doubt as to whether ground-state equilibria have been reached, and in the case of run 2, despite reaching a metastable configuration, we are certain that the thermodynamic ground state had not been reached. This suggests that the attainment of true thermodynamic equilibrium in these network ceramics may still be out of reach of modern computation. Nevertheless, we have observed an interface at which an intergranular film formed and thickened into the adjacent crystals. We have also observed an initially wide interface at which an IGF formed and showed no tendency to thicken. The resulting thicknesses of the two simulated IGFs are comparable, and are consistent with experimentally observed IGF widths. Both simulations showed somewhat diffuse interfaces between the IGF and the prism faces of the crystals. The diffuseness was associated with the substitution of O atoms at N crystal sites up to a couple of layers into the crystal. These are substitutional
(A.1)
This is based, with minor alterations, on the form of the ReaxFF potential outlined by van Duin et al. [19], and we have followed the original notation as often as possible. Ebond , EvdW , and ECoulb are two-body1 energy terms representing the bond, van der Waal’s, and Coulomb interaction energies, respectively. Eqdist represents the energy required to ionize neutral atoms to the charges calculated by the electronegativity equalization method (EEM) [29]. Eval is a three-body bond bending term. Eov/und is a many-body term designed to penalise over or under coordination. Finally there are small contributions to the energy, Elp , due to the presence or absense of lone pairs. The differences when compared to Eq. (1) of the overview of ReaxFFs by van Duin et al. are as follows: Epen has been removed from the notation because it was not used in those potentials anyway; Etorsion has been removed to reduce the number of parameters in our potential, chosen because the magnitude of this four-body term is significantly smaller in these systems than other terms; Econj , energy due to conjugation, is not physically relevant to the Si N O system; Eqdist was already included in the practical implementation of the ReaxFF potentials, but was not explicitly mentioned in that article. A.1. Bond energy Two-body (“raw”) bond orders are calculated as follows: pbo,2 pbo,4 rij rij BOij = exp pbo,1 + exp p bo,3 r0 r0 rij pbo,6 + exp pbo,5 (A.2) r0 1
The bond energy is based on bond order. Although it is initially a two-body parameter, it is rescaled to an extent depending on the under or over coordination of both atoms (see Section A.1 for details).
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
where the three terms correspond to the bond, first bond, and second bond (referred to as ), bond orders, respectively, rij the distance between the bonded atoms, r0x the onset radii of each type of bond, and pbo,x are physically related to the onset shape and width of each type of bond. There is a subtlety in this calculation concerning the minimum fraction of a full bond order that is considered to constitute a bond. Many of the other terms in the ReaxFF rely on knowing which pairs of atoms are bonded, and it would be inefficient to consider atoms with a tiny bond order between them to be bonded at all for the purposes of those other terms. A cutoff in the bond order is thus introduced, with a value of 0.01λ30 . However, to ensure continuous potentials for molecular dynamics (MD), bonds with a total bond order of more than the cutoff are stretched such that the bond order becomes 0 as they approach the cutoff. Any bond order less than the cutoff is considered to be 0. These two-body bond orders are then adjusted (Eq. (A.3)) to account for either of the participating atoms being overcoordinated or undercoordinated: BOij = BOij f1 (i , j )f4 (i , BOij )f4 (j , BOij )
(A.3)
where i is the deviation of the sum of the uncorrected bond orders on atom i from its valency Vi (silicon 4, nitrogen 3, oxygen 2), and Vi + f2 (i , j ) 1 f1 (i , j ) = 2 Vi + f2 (i , j ) + f3 (i , j ) Vj + f2 (i , j ) (A.4) + Vj + f2 (i , j ) + f3 (i , j ) f2 (i , j ) = exp(−λ1 i ) + exp(−λ1 j ) f3 (i , j )
1 = − ln λ2
(A.5)
1 [exp(−λ2 i ) + exp(−λ2 j )] 2
where f13 (rij ) =
1 1 + exp(−λ3 (λ4 (BOij )2 − n ) + λ5 )
+
1 γw
λ29 1/λ29 (A.10)
A.3. Electrostatic energy The electrostatic contributions to the total energy of the system consist of two terms, Eqdist , the energy required to place charge on the atoms, and ECoulb , the Coulomb energy of interaction between these charges. The charges themselves are determined by the EEM [29], which requires a large matrix inversion, and is the most computationally expensive part of the energy calculations. Once the charge qi on each atom i is determined, the charge distribution energy can be calculated as follows. Eqdist,i = −χi qi + ηi qi2
(A.11)
where χi is the electronegativity, and ηi is hardness of atom i. The Coulomb interaction is given by ECoulb = T (r) where γcco =
(r 3
C q i qj + γcco )1/3
1 3/2 3/2 γEEM,i γEEM,j
(A.12)
is a screening parameter and may be
different for each pair of interacting elements. A 7th order taper function T (r) is also used to smoothly cut this interaction off in the shell between an inner cutoff radius of λ12 , and an outer cutoff radius of λ13 : T (r) = a0 + a1 r + a2 r 2 + a3 r 3 + a4 r 4
(A.7)
Note some corrections to small typographic errors in Eqs. (3a), (3c), and (3f) of van Duin et al. [20]. Once the bond orders have been corrected in this many-body manner, bond energy is evaluated as: p
+ a5 r 5 + a6 r6 + a7 r 7
(A.8)
a0 = a1 =
a3 =
A.2. Van der Waal’s energy The van der Waal’s energy between each pair of atoms is calculated using: f13 (rij ) EvdW = Dij exp αij 1 − rvdW 1 f13 (rij ) − 2 exp αij 1 − (A.9) 2 rvdW
(A.13)
where the coefficients are given by
a2 =
Ebond = −De BOij exp[pbe,1 (1 − BOijbe,2 )] − De BO ij − De BO ij
λ rij29
and the van der Waal’s interaction radius, rvdW , and γw may be different for each pair of interacting elements. The van der Waal’s energy includes short-range core repulsion as well as long range dispersion. Again note corrections to typographic errors in [20].
(A.6) f4 (n , BOij ) =
131
a4 = a5 = a6 = a7 =
−35λ312 λ413 + 21λ212 λ513 + 7λ12 λ613 + λ713 (λ13 − λ12 )7 140λ312 λ313 (λ13 − λ12 )7 −210(λ312 λ213 + λ212 λ313 ) (λ13 − λ12 )7 140(λ312 λ13 + 3λ212 λ213 + λ12 λ313 ) (λ13 − λ12 )7 3 −35(λ12 + 9λ212 λ13 + λ313 ) (λ13 − λ12 )7 2 84(λ12 + 3λ12 λ13 + λ213 ) (λ13 − λ12 )7 −70(λ12 + λ13 ) (λ13 − λ12 )7 20 (λ13 − λ12 )7
(A.14)
132
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
Eval = f7 (BOij )f7 (BOjk )f8 (j ){ka
A.4. Lone pair energy The energy associated with lone pairs such as those on oxygen and nitrogen atoms is calculated using the expression: Elp =
plp lp,i 1 + exp(−75lp,i )
(A.15)
lp,i is the difference between the actual number of lone pairs on atom i from the optimal number of lone pairs for that element, where the actual number of lone pairs is given by nlp,i = int
ei 2
+ exp −λ16
2 − ei
+ 2 int
ei 2
2
(A.16) ei is the difference between the total number of outer shell electrons, ne,i , and the sum of the bond orders around atom i. Note the corrected signs in comparison to Eq. (4) of van Duin et al. [19].
− ka exp[−kb (Θ0 − Θijk )]}
(A.20)
where p
f7 (BOij ) = 1 − exp(−Aj,26 BOijv,7 )
(A.21)
f8 (j ) = Aj,29 − (Aj,29 − 1) ×
2 + exp(λ15 j ) 1 + exp(λ15 j ) + exp(−pv,2 j )
(A.22)
Note that since previous publications of the ReaxFF, some parameters that were once global have been split into a parameter for each type of bond angle between particular elements (ka , kb , pv,2 , pv,7 ). There is also a subtlety about which pairs of bonds should be counted in this energy. If both bonds are very weak, just above the bond order cutoff, then it would be a waste of computational time to compute bond bending energies. Thus another cutoff and rescaling is included in the bond orders for valence angle energy calculations, similar to that done earlier for bond energies. This only makes a small change to the energy calculations.
A.5. Bond bending energy The bond bending energy calculations must first determine the optimum bond angle, Θ0 , based on bond orders, because atoms with different hybridisations prefer different geometries. This is achieved through the equations: Θ0 = π − Θ0,0 {1 − exp[−λ18 (2 − SBO2)]}
A.6. Over and under-coordination energies The final major term in the ReaxFF energies for this study is due to over or under coordination. The under-coordination energy is expressed as:
(A.17) Eunder,i = −punder
where
SBO2 =
⎧ ⎪ ⎪ ⎪ ⎨
0,
(1 − exp(λ7 i )) (1 + exp(−Ai,25 i ))
f6 (SOVi ) (A.23)
if SBO < 0
SBOλ17 ,
if SBO > 0 and SBO < 1
⎪ 2 − (2 − SBO)λ17 , ⎪ ⎪ ⎩ 2,
where
if SBO > 1 and SBO < 2 if SBO > 2
f6 (SOVi ) =
1 1 + λ9 exp(λ10 SOVi )
(A.24)
(A.18) and
⎡
SBO = ⎣1 −
neighb(j) n=1
+
neighb(j)
and
⎤ exp(−BO8jn )⎦ (j − λ34 lp,j )
BOjn,π
SOVi =
neighb(i) n=1
(A.19)
(A.25)
The over-coordination energy is expressed as: neighb(i)
n=1
Note the corrected bracketing in Eq. (A.19) when compared to Eq. (7c) of van Duin et al. [19]. j is the sum of the bond orders from atom j minus the valency of the atom, minus the optimal number of lone pairs. Once Θ0 has been found, the bond bending or “valence angle” energy can be calculated. For the angle Θijk formed between atoms i, j and k, centred on j, this is expressed as:
(BOn − Vn − lp,n )(BO n + BOn )
Eover,i =
n=1 ov,i
ov,i pbe,3 Dbe,1 BOij + Vi + 10−8 1 + exp(Ai,25 ov,i ) (A.26)
where ov,i =
neighb(i) −lp,i BOni − Vi + 1 + λ33 exp(λ32 SOVi ) n=1
(A.27)
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
133
Appendix B. Reaction force field parametrisation Tables B.1–B.6 Table B.1 General parameters for the reaction force field Parameter
Value
Description
λ1 λ2 λ3 λ4 λ5 λ7 λ9 λ10 λ12 λ13 λ15 λ16 λ17 λ18 λ29 λ30 λ32 λ33 λ34
50.0000 5.5469 21.2839 3.0000 6.5000 1.0159 8.0878 13.0238 0.0000 10.0000 33.8667 25.6125 1.0563 2.0384 1.5591 0.1000 1.6356 5.6937 2.5067
Over-coordination parameter Over-coordination parameter Valency angle conjugation parameter Triple bond stabilisation parameter Triple bond stabilisation parameter Under-coordination parameter Under-coordination parameter Under-coordination parameter Lower Taper-radius Upper Taper-radius Valency undercoordination Valency angle/lone pair parameter Valency angle Valency angle parameter van der Waal’s shielding Cutoff for bond order (×100) Over-coordination parameter Over-coordination parameter Valency/lone pair parameter
These apply to the whole system, and are not element specific. All were previously determined [19,20], and have not been fitted in this study.
Table B.2 Atom parameters Index
Symbol
Si
N
O
Description
Ai,1 Ai,2 Ai,4 Ai,5 Ai,6 Ai,7 Ai,8 Ai,9 Ai,10 Ai,12 Ai,14 Ai,15 Ai,17 Ai,18 Ai,25 Ai,26 Ai,29
r0 V rvdW D γEEM r0 ne ␣ γw −punder χ η r0 plp Ai,25 Ai,26 Ai,29
2.0132 4.0000 2.2203 0.1306 0.8218 1.5629 4.0000 12.1305 1.8838 11.4056 1.8038 7.3852 −1.0000 0.0000 −3.6886 2.9867 2.5791
1.2360 3.0000 1.9333 0.1373 0.8596 1.1602 5.0000 10.0725 8.0467 34.5504 6.8418 6.3404 1.0365 14.7406 −15.0000 2.6491 2.8793
1.1693 2.0000 1.9580 0.0875 1.0804 1.0200 6.0000 10.1761 7.2516 39.2737 8.5000 7.8386 0.9809 18.4195 −7.3000 2.6656 2.9225
cov.r Valency Rvdw Evdw gammaEEM cov.r2 # alfa gammavdW Eunder chiEEM etaEEM cov.r3 Elp ov/un val1 vval4
The index is from the position this variable takes in the structured input to the ReaxFF program by van Duin [21]. Some parameters take this index as a symbol because they were instead denoted as global parameters in previous publications [19,20].
Table B.3 Bond parameters for Si Si, Si N, and Si O Index
Symbol
Si Si
Si N
Si O
Description
Bi,1 Bi,2 Bi,3 Bi,4 Bi,5 Bi,6
De De De pbe,1 pbo,5 ovc13
113.7903 53.9894 30.0000 0.2535 −0.3000 1.0000
144.4831 20.4624 10.0000 − 0.3794 − 0.3000 1.0000
193.1177 41.1424 43.3991 −0.2085 −0.3000 1.0000
Edis1 Dbe1 LPpen Dbe2 n.u. Dbe3 pbe1 pbo5 13corr
134
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
Table B.3 (Continued ) Index
Symbol
Si Si
Si N
Si O
Description
Bi,7 Bi,8 Bi,9 Bi,10 Bi,11 Bi,13 Bi,14 Bi,15
pbo,6 pbe,3 pbe,2 pbo,3 pbo,4 pbo,1 pbo,2 ovc
16.0000 0.0742 0.2586 −0.1963 7.5407 −0.0693 7.9365 0.0000
36.0000 1.5282 3.2297 −0.1012 12.2585 −0.1039 8.1682 1.0000
36.0000 0.7695 0.9220 −0.3683 4.2644 −0.5191 4.4529 1.0000
pbo6 pbe3 pbe2 pbo3 pbo4 pbo1 pbo2 ovcorr
The index is from the position this variable takes in the structured input to the ReaxFF program by van Duin [21].
Table B.4 Bond parameters for N O, N N, and O O Index
Symbol
N O
N N
O O
Description
Bi,1 Bi,2 Bi,3 Bi,4 Bi,5 Bi,6 Bi,7 Bi,8 Bi,9 Bi,10 Bi,11 Bi,13 Bi,14 Bi,15
De De De
147.4287 109.4565 61.7111 0.9507 −0.1616 1.0000 21.7802 0.7500 0.7441 −0.3665 7.2396 −0.2439 4.8277 1.0000
155.9636 0.0000 0.0000 0.2592 −0.1418 1.0000 13.1260 0.3000 0.3989 −0.1336 10.7257 −0.0760 5.2866 1.0000
118.8570 42.8420 167.6132 0.9092 −0.1587 1.0000 10.5117 0.8589 0.8634 −0.2501 6.5155 −0.1613 5.6591 1.0000
Edis1 Dbe1 LPpen Dbe2 n.u. Dbe3 pbe1 pbo5 13corr pbo6 pbe3 pbe2 pbo3 pbo4 pbo1 pbo2 ovcorr
pbe,1 pbo,5 ovc13 pbo,6 pbe,3 pbe,2 pbo,3 pbo,4 pbo,1 pbo,2 ovc
The index is from the position this variable takes in the structured input to the ReaxFF program by van Duin [21]. Table B.5 Off-diagonal parameters Bond
D
rvdW
γw
r0
r0
r0
O Si N Si
0.1233 0.1536
1.7552 1.7871
10.6931 11.0924
1.6096 1.7610
1.2935 1.2495
N/A N/A
If available for a bond type, these override values based on elemental properties.
Table B.6 Bond bending parameters Type O O N O O N Si O O Si O N N O Si N O N
O O O N N N Si Si Si O O Si Si Si N N N O
O N N O N N Si Si O Si Si Si N N Si Si Si Si
θ0,0
ka
kb
pv,2
pv,7
78.4963 79.0601 74.5183 74.8323 77.6183 75.2419 69.3456 70.3016 90.0344 22.1715 73.4663 1.4881 83.1400 77.4641 34.4338 125.3711 69.8728 69.8728
61.0192 38.7432 54.0833 35.2858 29.5041 20.4666 21.7361 15.4081 7.7656 3.6615 25.0761 4.0912 0.0562 4.5724 40.7839 12.9852 32.7155 27.1273
1.0093 2.1471 1.9745 1.4738 1.2928 3.1280 1.4283 1.3267 1.7264 0.3160 0.9143 13.6686 2.2347 1.0849 0.4583 179.7577 1.5875 1.5875
0.4897 0.4897 0.4897 2.7993 2.7993 2.7993 −0.2101 2.1459 0.7689 4.1125 2.2466 12.0668 0.0079 0.7689 6.5339 18.4450 2.2466 2.2466
1.0400 1.0400 1.0400 1.0400 1.0400 1.0400 1.3241 1.0400 1.0400 1.0400 1.0400 1.0400 1.0400 1.0400 1.0400 1.0400 1.0400 1.0400
T.S. Hudson et al. / Materials Science and Engineering A 422 (2006) 123–135
References [1] P.F. Becher, G.S. Painter, E.Y. Sun, C.-H. Hsueh, M.J. Lance, Acta. Mater. 48 (2000) 4493. [2] G.S. Painter, P.F. Becher, H.-J. Kleebe, G. Pezzotti, Phys. Rev. B 65 (2002) 064113. [3] S. Ii, C. Iwamoto, K. Matsunaga, T. Yamamoto, M. Yoshiya, Y. Ikuhara, Phil. Mag. 84 (2004) 2767. [4] C.-M. Wang, X. Pan, M.J. Hoffman, R.M. Cannon, M. R¨uhle, J. Am. Ceram. Soc. 79 (1996) 788. [5] I. Tanaka, H.-J. Kleebe, M.K. Cinibulk, J. Bruley, D.R. Clarke, M. R¨uhle, J. Am. Ceram. Soc. 77 (1994) 911. [6] H. Gu, J. Am. Ceram. Soc. 85 (2002) 33. [7] M. Doeblinger, private communications. [8] H. Gu, in: Interface Science and Materials Interconnection, Silicon Nitride Based Ceramic Materials, Proceedings of JIMIS08. Japan Institute of Metals, Sendai, Japan, (1996), pp. 217–220. [9] A. Ziegler, C. Kisielowski, M.J. Hoffmann, R.O. Ritchie, J. Am. Ceram. Soc. 86 (2003) 1777. [10] G.B. Winkelman, C. Dwyer, T.S. Hudson, D. Nguyen-Manh, M. Doblinger, R.L. Satet, M.J. Hoffmann, D.J.H. Cockayne, Phil. Mag. Lett. 84 (2004) 755. [11] W. McBride, D.J.H. Cockayne, D. Nguyen-Manh, Ultramicroscopy 96 (2003) 191. [12] M. Doblinger, C.D. Marsh, D. Nguyen-Manh, D. Ozkaya, D.J.H. Cockayne, Inst. Phys. Conf. Ser. 179 (2004) 401.
135
[13] D. Nguyen-Manh, D.J.H. Cockayne, M. Doeblinger, C. Marsh, A.P. Sutton, A.C.T. van Duin, Inst. Phys. Conf. Ser. 179 (2004) 139. [14] M. Yoshiya, K. Tatsumi, I. Tanaka, H. Adachi, J. Am. Ceram. Soc. 85 (2002) 109. [15] S.H. Garofalini, W.W. Luo, J. Am. Ceram. Soc. 86 (2003) 1741. [16] X.T. Su, S.H. Garofalini, J. Mater. Res. 19 (2004) 752. [17] D.J. Adams, Mol. Phys. 28 (1974) 1241. [18] D.J. Adams, Mol. Phys. 29 (1975) 307. [19] A.C.T. van Duin, A. Strachan, S. Stewman, Q. Zhang, X. Xu, W.A. Goddard III, J. Phys. Chem. A 107 (2003) 3803. [20] A.C.T. van Duin, S. Dasgupta, F. Lorant, W.A. Goddard III, J. Phys. Chem. A 105 (2001) 9396. [21] A. C. T. van Duin, ReaxFF molecular dynamics program. [22] S. H. Garofalini, private communications. [23] X. Su, S.H. Garofalini, J. Mater. Res. 19 (2004) 752. [24] W.-Y. Ching, private communications. [25] H. Gu, R.M. Cannon, H.J. Seifert, M.J. Hoffmann, I. Tanaka, J. Am. Ceram. Soc. 85 (2002) 25. [26] M.P. Allen, D.J. Tildesley, Computer Simulations of Liquids, Clarendon Press, Oxford, UK, 1987, pp. 39–42. [27] M. Kitayama, K. Hirao, A. Tsuge, M. Toriyama, S. Kanzaki, J. Am. Ceram. Soc. 82 (1999) 3263. [28] V. Burlakov, private communications. [29] W.J. Mortier, S.P. Ghosh, S. Shankar, J. Am. Chem. Soc. 108 (1986) 4315.