Computational Study of Random Boron-Carbon-Nitrogen Systems By SIFISO MUSA NKAMBULE A dissertation submitted to the Faculty of Science, University of the Witwatersrand, In fulfillment of the requirements for the degree of Master of Science. School of Physics, and DST-NRF Centre of Excellence in Strong Materials, University of the Witwatersrand, Johannesburg, South Africa. 2010 Declaration I declare that this dissertation is my own work. It is being submitted for the degree of Master of Science in the University of the Witwatersrand, Johan- nesburg. It has not been submitted before for any degree or examination in any other university. Sifiso Musa Nkambule:.......................... Date:......................2010 i Acknowledgements Firstly I would like to thank my supervisor, Professor J.E. Lowther for giving me an opportunity to work under his supervision. I greatly appreciate his mentorship and courageous support which made me to believe in myself even at times when it felt impossible. Sincere thanks to my research focus area colleagues Onyekwelu Okeke, Thabo Letsoalo, Vuyokazi Bangani, and Mahlaga Molepo for the handful contribution to my research and for the light moments of laughter we shared, which really kept me going. I would also like to acknowledge the Dean of the Faculty of Science in the University of Swaziland, Professor M.D. Dlamini for giving me a good recom- mendation. I would also like to thank the people who contributed much to my high school education: Samuel Nkambule (my father), Sabelo Nkambule ii (my uncle), Sabelo Dlamini, Siphosini Mazibuko and Junior Mkhabela (my spiritual brothers). I greatly appreciate the valuable contribution they have made in my life. Special acknowledgement goes to the special people in my life, my wife Diana and my son Junior. I specially thank my wife for believing in me and her ?I know you can do it ?statement, which made me to try harder. I thank my son for keeping me constantly in reminder that he misses me which made me to work very hard so that I can be with him. Lastly, the financial support received from DST-NRF Centre of Excellence in Strong materials and the University of the Witwatersrand is greatly ac- knowledged. !!!!!NgeSiSwati batsi: Kwandza kwaliwa batsakatsi, akwandze lapho nit- setse khona. !!!!! iii Dedication To my wife, Diana Bekisiwe Magagula. With love and gratitude !!!!! Ulubhambo lwami Mtfombeni, Gutjwa, Komtwako, Msutfu!!! iv Citation to previously published work Portions of this work have appeared in a published article: S.M. Nkam- bule and J.E. Lowther, Crystalline and random ?diamond-like? boron-carbon structures, Solid State Communications, 150,133-136,(2010). Part of this work have been presented in the South African Institute of Physics (SAIP) conference, in July 2009, Durban, South Africa, by a poster presentation, entitled; Computational study of random Boron-Carbon systems. v Contents 1 Introduction 1 1.1 Strong Materials . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Hardness Properties . . . . . . . . . . . . . . . . . . . 2 1.1.2 Vickers Hardness . . . . . . . . . . . . . . . . . . . . . 3 1.1.3 Bulk Properties . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Diamond . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.3 B-C complexies . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.4 B-C-N complexies . . . . . . . . . . . . . . . . . . . . . . . . . 18 2 Computational Details 22 2.1 Potential Functions . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2 Interatomic potentials . . . . . . . . . . . . . . . . . . . . . . 23 2.2.1 Tersoff potential . . . . . . . . . . . . . . . . . . . . . . 23 vi 2.3 Previous studies using Tersoff potentials . . . . . . . . . . . . 29 2.4 Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . 32 2.4.1 Molecular Interactions . . . . . . . . . . . . . . . . . . 35 2.5 The Time Integration . . . . . . . . . . . . . . . . . . . . . . 36 2.6 Long-range Interaction Algorithms . . . . . . . . . . . . . . . 37 2.6.1 Ewald technique . . . . . . . . . . . . . . . . . . . . . 38 2.7 The GULP code . . . . . . . . . . . . . . . . . . . . . . . . . 38 3 Application of the Tersoff Potential to B-C Complexies 42 3.1 B-C Crystal Structures . . . . . . . . . . . . . . . . . . . . . . 42 3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2.1 Graphical representation of the results for the B-C crystal structures . . . . . . . . . . . . . . . . . . . . . 52 4 Tersoff Potential Application to B-C-N Complexes 64 4.1 Comparison of the predictions of the Two sets of Tersoff Po- tential Parameters for Nitrogen . . . . . . . . . . . . . . . . . 68 4.2 Further Results for B-C-N with Parameter Set 2 . . . . . . . . 77 4.2.1 Crystalline Phases results . . . . . . . . . . . . . . . . 77 vii 4.2.2 Random Phases Results . . . . . . . . . . . . . . . . . 82 5 Discussion and Conclusion 89 Bibliography 94 viii List of Tables 2.1 Tersoff interatomic potential parameters for Carbon, Boron and Nitrogen . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.2 Tersoff interatomic potential parameters for Silicon and Ger- manium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1 Properties of crystalline phase structures predicted using OPT. Values in brackets are experimental values. . . . . . . . . . . . 46 3.2 Properties of crystalline phase structures predicted using MD. Values in brackets are experimental values . . . . . . . . . . . 47 3.3 Properties of the random phase structures predicted using MD. Values in brackets are experimental values . . . . . . . . 48 3.4 Properties of the random phase structures predicted using OPT. Values in brackets are experimental values. . . . . . . . 49 ix 3.5 Cell parameters predicted for the B-C structures. Values in brackets are experimental values . . . . . . . . . . . . . . . . . 51 4.1 Predicted properties of crystal phase using Parameter Set 1 of B-C-N. Values in brackets are experimental values. Energy here refers to the ground state total formation energy . . . . . 69 4.2 Predicted properties of random phase using Parameter Set 1 of B-C-N. Energy here refers to the ground state total formation energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.3 Predicted properties of crystal phase using Parameter Set 2 of B-C-N. Energy here refers to the ground state total formation energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.4 Predicted properties of random phase using Parameter Set 2 of B-C-N. Energy here refers to the ground state total formation energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.5 Table showing the B-C-N structures?s number of bonds as a percentage of the total number of bonds of a structure . . . . 74 4.6 Table of properties of crystal phase of B-C-N structures, Val- ues in brackets are experimental values . . . . . . . . . . . . . 78 x 4.7 Cell properties of crystal phase of B-C-N structures. The Elas- tic properties are expresse in units of GPa. Energy refers to the total formation energy. Values in brackets are experimen- tal values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.8 Properties of random, 2744 atoms, phase of B-C-N structures. The Elastic properties are expressed in units of GPa. Energy refers to the total formation energy. . . . . . . . . . . . . . . . 83 4.9 Properties of, 2744 atoms, random, phase of B-C-N structures, using Parameter set 2. Values in brackets are experimental values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 xi List of Figures 1.1 Types of cubic lattice symmetries, (a) sc, (b) bcc, (c) fcc . . . 9 1.2 Tetragonal lattice symmetries, (a) simple tetragonal, (b) body- centred tetragonal . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3 Trigonal lattice symmetries [1] . . . . . . . . . . . . . . . . . . 11 2.1 The interatomic bond vector[2]. . . . . . . . . . . . . . . . . . 23 2.2 A schematic diagram showing the bond angle ?[3]. . . . . . . . 26 2.3 A graph showing the relationship ? and g(?) taken from[3]. . . 27 2.4 A graph showing the dependace of ? on rij ? rik taken from[3]. 27 3.1 B-C structures considered in this work. Dark yellow and white represent Boron and Carbon atoms, respectively. . . . . . . . 44 3.2 Percentage of the number of bonds versus the composition for a random phase . . . . . . . . . . . . . . . . . . . . . . . . . . 53 xii 3.3 Energy for crystal structures. (a) is for OPT, while (b) is for MD. The horizontal axes are showing the composition . . . . . 55 3.4 Energy for random structures. (a) is for OPT, while (b) is for MD. The horizontal axes are showing the composition . . . . . 56 3.5 B, G, and  for crystal structures. (a) is for OPT, while (b) is for MD. The horizontal axes are showing the composition . 57 3.6 B, G, and  for random structures. (a) is for OPT, while (b) is for MD. The horizontal axes are showing the composition . 58 3.7 Elastic constants for crystal structures. (a) is for OPT, while (b) is for MD. The horizontal axes are showing the composition 59 3.8 Elastic constants for the random structures. (a) is for OPT, while (b) is for MD. The horizontal axes are showing the com- position . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.9 Cell constant against Composition . . . . . . . . . . . . . . . . 61 3.10 Bulk(?), shear( ) and Young(5) moduli for various sizes of the supercells with random B and C distribution. (a), (b), (c) and (d) represents BC, B3C5, BC3 and BC7, respectively. The horizontal axes are showing the composition . . . . . . . 62 xiii 3.11 Energy difference of the crystal and random phase against B-C trend . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.12 Ratio of the random and crystal lattice cell constants . . . . . 63 4.1 B-C-N structures considered for the composition B4N4?xCx, (x = 0 . . . 4.) , dark yellow, white and yellow represent Boron, Carbon and Nitrogen atoms, respectively. . . . . . . . . . . . 66 4.2 B-C-N structures considered for the composition B4?xN4?xC2x, (x = 0 . . . 4.) , dark yellow, white and yellow represent Boron, Carbon and Nitrogen atoms, respectively. . . . . . . . . . . . 67 4.3 Bulk Moduli and elastic constants in units of GPa as shown on the vertical axes, with Set 1 parameters. (a) and (b) represent the crystal and random phase, respectively. The horizontal axes represent the composition of the structures . . . . . . . . 75 4.4 Bulk Moduli and elastic constants in units of GPa as shown on the vertical axes, with Set 2. (a) and (b) represent the crystal and random phase, respectively. The horizontal axes represent the composition of the structures . . . . . . . . . . . 76 xiv 4.5 (a) Elastic constants, and (b) Moduli for the crystal phase structures. The horizontal axes represent the trend in the B- C-N composition . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.6 The total formation energies( in part (a)) and the cell con- stants (in part (b)), for the crystal structures with Parameter Set 2. The horizontal axes in both graphs represent the com- position in the B-C-N trend . . . . . . . . . . . . . . . . . . . 81 4.7 The elastic constants (shown in part (a), and the Moduli (shown in part (b)) for the 2744 atoms, random structure, using Set 2. The horizontal axes represent the B-C-N compo- sition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.8 Total energy (in part (a)) and cell constants for the random, structures, with 2744 atoms, using Set 2, plotted against the B-C-N composition. . . . . . . . . . . . . . . . . . . . . . . . . 86 4.9 Poisson ratio ? against the B-C-N composition. . . . . . . . . 87 4.10 Percentage number of bonds against the B-C-N composition . 88 xv Abstract Elastic properties, total formation energies and cell parameters of B-C struc- tures are studied by molecular dynamics (MD) simulation. A well tested empirical Tersoff interatomic potentialfor simulation using a General Utility Lattice program (GULP). The results are compared with the available ex- perimental results for Diamond and BC5. The work is extended to the study of B-C-N materials. We compare with available experimental results for BN and BC2N . The Bulk Modulus, Shear Modulus, Young?s Modulus and Elastic con- stants are obtained strongly implying hardness comparable to that of dia- mond. The Novel B-C and B-C-N material are proposed to be ultra-hard ma- terials with different chemical and electrical properties from those of diamond and with potentially very useful application such as advanced abrasives, high temperature electronics and semiconductors. Chapter 1 Introduction 1.1 Strong Materials On daily basis we rely on strong materials for many applications. Materials used in manufacturing cars, aeroplanes, powerstations, spacecraft, and vari- ous power tools need to be strong enough to protect to prevent an accident. Tools of all sorts, from fork and knives, to a wide range of precision instru- ments used in industry and surgery, have to stand up to different degrees of heat, abrasion and pressure. There is a great gain in different sectors if there is an improvement in the quality of strong materials because efficiency increases and costs go down. 1 2 Semiconducting devices also can be greatly improved if , for instance, the temperature at which they can function is increased. 1.1.1 Hardness Properties Hardness is defined as resistance of a material to deformation. The term can also refer to stiffness to temper, resistance to scratching, abrasion, or cutting. It is the property of a substance, which gives it the ability to resist being permanently deformed (bent, broken, or have its shape changed), when a load is applied. The greater the hardness of the substance, the greater re- sistance it has to deformation. Hardness is a very complex property in the sense that it is dependant on many physical properties. It depends on vari- ous parameters like pressure, temperature, impurities, porosity, dislocations and other defects. However, it is correlated to various physical properties including ionicity, melting point, band gap, cohesive energy, etc, and can be then studied indirectly. The hardness of a material is measured in several ways. The simplest test for non-metals is the scratch (Mohs) test[4]. Substance A is harder than substance B if A will scratch B, but B will not scratch A. There are various 3 other hardness test available. These include the Rockwell hardness[5], Brinell hardness[5], Vickers[6], Knoop hardness[5] and the Shore[7] test. Only the Vicker?s hardness will be described in more detail. 1.1.2 Vickers Hardness Modern scales of hardness, such as the Vickers Hardness numbers (VHN), are based on indenter tests in which an indenter is pressed into the surface of the material and the size of the impression is measured. The Vickers hardness test was developed in 1924 by Smith and Sandland at Vickers Ltd as an alternative method to measure the hardness of materials[6]. Most often the Vickers test is much easier to use than other hardness tests since the required calculations are independent of the size of the indenter. Hence the indenter can be used for all materials irrespective of hardness. Like all common measures of hardness, the basic principle is to observe the material?s ability to resist plastic deformation from a standard source. The test can be used for all metals and has one of the widest scales among hardness tests with the unit of hardness known as the Vickers Pyramid Number (HV). 4 1.1.3 Bulk Properties Elastic constants: The harmonic elastic constants represent the second derivatives of the energy density with respect to strain: Cij = 1 V ( ?2U ?i?j ) (1.1) where V, U and  denotes the volume, energy and strain tensor, respectively. Hence the harmonic elastic constants describe the mechanical strength of the material with respect to deformation, Since there are 6 possible strains within the notation scheme, the elastic constant tensor is a 6 x 6 symmetric matrix. The 21 potentially independent matrix elements are usually reduced considerably by symmetry[8]. For a cubic material there are only three unique elements, namely C11, C12, and C44. Bulk and shear moduli: Like the elastic constant tensor, the bulk (B) and shear (G) moduli contain information regarding the hardness of a material with respect to various types of deformation. The bulk modulus is much more easily determined experimentally, than the elastic constant tensor. At zero temperature, B is 5 calculated using the following equation[9]: B = ?V dpdV = V d2U dV 2 (1.2) where V , p, and U are the volume, pressure and energy, respectively. Cohen[9] developed a relation of the bulk modulus to the to the distance of nearest neighbour separation, for compounds which are near the centre of the peri- odic table. His relation was; B = 1761d?3.5 (1.3) Where B is in units of GPa and d, the nearest neighbours distance, is in units of ?A. If the structure of a material is studied as a function of applied isotropic pressure, then a plot of pressure versus volume can be fitted to an equation of state, where the bulk modulus is one of the curve parameters[10]. Typi- cally a third or fourth order Birch-Murgnahan[11] equation of state can be utilised. Also, the bulk and shear moduli are related to the elastic constant elements, but no unique definition for this transformation has been found. Hebbache[12] theoretically asserted that G is correlated to the hardness of a material, which he tested for diamond. 6 The bulk and shear modulus can be obtained using computational meth- ods. Young?s moduli: When a uniaxial tension is applied to a material then the lengthening of the material is measured according to the strain. The stress to strain ratio defines the value of the Youngs modulus for that axis[13]: ?? = ??? ?? (1.4) Where ? denotes the stress tensor. Since a material will always increase in length under tension, the value of this quantity should always be positive[10]. The Youngs moduli in each Cartesian directions can be calculated from the elastic compliances: ?x = 1 S11 (1.5) ?y = 1 S22 (1.6) ?z = 1 S33 (1.7) Poissons ratio: The Poisson ratio measures the change in a material at right angles to the uniaxial stress. Hence, the Poisson ratio and Young?s modulus are com- 7 plementary properties. It is defined as the ratio of lateral to longitudinal strain under a uniform, uniaxial stress. Assuming an isotropic medium, this property can be calculated as given below [8]: ??(?) = ?S?????? (1.8) Because most materials naturally shrink in the direction orthogonal to an applied tension this leads to mainly positive values for Poisson ratio, with a theoretical maximum of 0.5[14]. values for many materials lie in the range 0.2-0.3[10], Athough negative values are also known[14]. This quantity is also be related to the bulk modulus in an isotropic material[8]: B = 13 ? (1? 2?) (1.9) where B, ? and  are as described above. Calculations of bulk ground-state properties, such as lattice constants, atomic positions, bond lengths, and bulk modulus, play an important role in the physics of materials. Such calculations allow us to understand, char- acterize, and predict mechanical properties of materials in our surroundings, under extreme conditions, such as in geological formations and settings, and for industrial applications. 8 Crystalline materials occur in many different structures and, in contrast to isotropic materials, the description of the ground state of crystalline materials may in general need multiple lattice parameters and an atomic basis[15]. Owing to the cubic symmetry of the diamond structure, there are only three independent elastic constants, C11, C12, and C44. Cubic lattice The cubic (or isometric) crystal system is a system where the unit cell is in the shape of a simple cube. It is the most common and simplest shapes found in crystals and minerals. There are three main varieties of these crystals called simple cubic (sc), body-centred cubic (bcc), and face-centred cubic (fcc). The lattices are shown in figure 1.1. A diamond has a fcc cubic symmetry and its primitive basis has two identical atoms[4]. Tetragonal lattice Tetragonal crystal lattices result from stretching a cubic lattice along one of its lattice vectors, so that the cube becomes a rectangular prism with a square base (a ? a) and height (c 6= a). There are two tetragonal Bravais lattices: the simple tetragonal (from stretching the simple-cubic lattice) and 9 (a) (b) (c) Figure 1.1: Types of cubic lattice symmetries, (a) sc, (b) bcc, (c) fcc 10 the centred tetragonal (from stretching either the face-centred or the body- centred cubic lattice). This is shown in figure 1.2[1]. (a) (b) Figure 1.2: Tetragonal lattice symmetries, (a) simple tetragonal, (b) body- centred tetragonal Trigonal lattice The rhombohedral system can be thought of as the cubic system stretched diagonally along a plane. a = b = c; ?, ?, ? 6= 90?, as shown in figure 1.3. 11 Figure 1.3: Trigonal lattice symmetries [1] . The tetragonal lattice can be very close to being cubic if a ? c. A trigonal lattice can be almost cubic if ?, ?, ? ? 90?. Knowledge of crystal structures of different materials at atomic level may lead to improved or new materials. These are refined with the help of com- puter modelling, which simulates the behaviour of the material under differ- ent conditions. This kind of fundamental work helps to indicate, for instance, what processes could be applied to create the crystal structures that, in turn, will most likely give rise to the desired properties and behaviour of the ma- terial. 12 1.2 Diamond Diamond is the hardest known material with many fascinating properties. It has found a wide range of applications in modern science and technol- ogy owing to its unique properties, such as extreme hardness, high thermal conductivity, wide band gap and high electron and hole mobility[16]. It has many unique properties owing to the tetrahedrally coordinated carbon cova- lent bonds. The primitive basis has two identical carbon atoms. Each carbon atom has 4 nearest neighbours and 12 next nearest neighbours. Other prop- erties of diamond include very large elastic modulus, high sound velocities, and high Debye temperature (2340 K)[17]. The carbon atoms form very short covalent bonds. The high atomic density of diamond (1.76? 1023 atom/cm3) accounts for its record elastic and mechanical properties, but at the same time limits the possibilities for doping it[18]. Although there are few prospects of diamond-based microelectronics oust- ing silicon totally, diamond devices could function in situations when silicon electronics fail. For example, diamond chips could potentially still work at temperatures of several hundreds degrees, whereas silicon devices generally 13 fail above 450 K. On the other hand diamond has some pitfalls. It is non-resistant to oxi- dation. While diamond stands out as the best material in cutting and pol- ishing metals and other hard solids, its superiority fails in cutting certain tough materials such as ferrous metals owing to the detrimental formation of iron carbide during high-speed machining[19]. Another pitfall with dia- mond is that its shape can not be easily changed and it wears easily at high temperature. Synthetic materials with a comparable hardness, are potential candidates to mitigate this limitation. There is an ever increasing demand for diamond- like materials in electronic applications. The growing demand for advanced superhard materials for cutting and shaping hard metals and ceramics, as well as in electronic and electrochemical applications, has stimulated the research for novel diamond like phases that are more thermally and chemically stable than pure diamond[16]. 14 1.3 B-C complexies Boron is an important impurity in diamond as a point defect[20]. It is effec- tive in enhancing the hardness of diamond[21]. It makes the material exhibit p-type behaviour with very nearly shallow defect characteristics[22, 21]. It readily make a substitution of carbon atom in the diamond lattice about which there is very little lattice relaxation[20, 23]. The surface is the easi- est place to make substitution for the carbon atoms without disturbing the strong sp3- bonding network of the diamond lattice. Each boron atom re- place a carbon atom on the diamond surface to form covalent bonds with three carbon atoms below it giving the boron 3-fold coordination. Surface boron atoms form planar sp2 bonds with carbon atoms and the compressive stress from the underlying bulk diamond lattice reduces the bonds signifi- cantly. Bulk and shear moduli in the surface region increase substantially, surpassing those of bulk diamond[24]. However, only a very limited amount of boron ( 2 at%) can be introduced into diamond structure by common meth- ods of thermal chemical vapour deposition and high pressure synthesis[20]. The limiting concentration of boron that can be incorporated into diamond structure has not been established yet, although under extreme pressure- temperature conditions, this content could be significantly increased. Since 15 all B-C materials show higher resistance to oxygen and ferrous metals than similar carbon materials[16], the diamond like B-C phases with high boron content are expected to combine the best properties of the elements includ- ing advanced electrical and optical properties, very high hardness, and high thermal and chemical stability. Boron-doped diamond is a very interesting material for the investigation of impurity conductivity[25, 26, 27].This is a result of a wide temperature region of hopping conductivity and to unusual behaviour of the Hall coefficient[28]. Point defect calculations for large concentration of boron in diamond have highlighted the fact that the near-neighbour pairing of boron atoms is un- likely, suggesting that the complex is metastable[29]. However it has been suggested[30, 28, 31] that the presence of such a near-neighbour defects is deduced from the conductivity of heavily doped boron in diamond. The dis- covery of superconductivity[32] following extreme conditions of high pressure treatment further suggest this near-neighbour pairing. Recent studies[33] suggest that the near-neighbour boron (B) pairs may be formed in diamond and even act as traps for hydrogen, implying that large concentration of boron can be induced into diamond, however in a metastable form. Metasta- 16 bility is further implied by the fact that extreme conditions are apparently a prerequisite for the synthesis of diamond in which superconductivity has been found. A typical amount of boron to be doped in diamond that has been success- fully attained is 1017?1021atoms/cm3[34, 35, 36, 37]. Boron doping causes a diamond to exhibit semiconducting properties, until, at a critical concentra- tion of nB ? 1020atom/cm3 the doping results in a metallic-type conductivity at a moderate temperatures [38, 39]. Ekimov et al have pointed out that a further increase in boron concentration to nB ? 1021atom/cm3 produces su- perconductivity in diamond[40] ,in which the boron concentration gives the upper limit of the carrier (hole) concentration. The latest research on the dynamics of the boron-doped diamond lattice by Raman spectroscopy[36, 41] and inelastic x-ray scattering[42] as well as the study of the critical temperature dependence of the doping level[43] indicate some advances in the experimental study of superconducting diamond. Recently BC5, the diamondlike B-C phase with highest boron content ever achieved has been synthesised[16]. It was found to have low compressibility 17 (bulk modulus of 335 GPa), to be conductive, to exhibit extreme Vickers hardness of 71 GPa and high thermal stability (up to 1900K). This makes it an exceptional superbrasive and a promising material for high temperature electronics. A study by Brazhkin et al [18] shows an abnormally high thermal expan- sion coefficient for heavily boron-doped diamonds, which could be related to the softening of the phonon spectrum due to the formation of the soft regions located near boron atoms. It further shows that the concentration dependence of the thermal expansion coefficient suggests two different modes of B incorporation, substitutional at low doping levels and aggregation-type at higher concentrations of the B dopant. In our study we look into 6 different structures of the B-C system, in re- lation to the boron concentration. We study cubic structures, which contain eight atom unit cells of B4?xC4+x,(x = 0 . . . 4). BC3 structures are examined and two structures are considered, trigonal (BC3? a) and tetragonal(BC3? b)), whose departure from cubic symmetry is very small. Our focus is to examine the elastic properties of each material. 18 1.4 B-C-N complexies Extensive theoretical and experimental effort has been employed in find- ing the possibility of new low-compressibility materials with bulk modulus and hardness that is comparable to that of diamond. Ternary boron-carbon- nitride (B-C-N)compounds have attracted much attention due to their poten- tial physical and chemical properties. B-C-N compounds have been of great interest to material scientists since the prediction of superhard ??C3N4[44] compound. Diamond like B-C-N compounds can be used as super-hard materials and high temperature semiconductors[45]. Carbon, nitrogen and boron in their strong covalent bonding configuration have an unusually high hardness and wear resistance and they have been widely used as coating materials[24]. Nitrogen and Boron atoms can be incorporated into the diamond lattice fairly easily[46]. Nitrogen and boron are very important impurities in dia- mond. Boron atoms become acceptor impurities while nitrogen atoms be- come donor impurities[28]. Furthermore, if nitrogen or boron atoms are used to replace surface carbon atoms on the diamond (111) surface, each N or B atom forms covalent bonds with three carbon atoms below it and such a 19 bonding environment is compatible with the tendency of N or B to have 3-fold coordination[47]. Substitutional nitrogen, in the diamond structure, yields a deep donor level resulting from the localization of the unpaired electron on a carbon dangling bond[48]. Considerable efforts have been devoted to the synthesis of the B-C-N com- pounds with different stoichiometric compositions such as BC2N,BC4N, and BC6N , and with different structures such as the hexagonal and cubic structures[49, 50, 51, 52, 53]. Since Diamond and cubic Boron Nitride (c?BN) are well known materials exhibiting outstanding properties such as high hardness, melting point, and bulk modulus, a ternary system of cubic B-C-N, is expected to be another superhard material and to show high thermal and chemical stability. Zhaoet al have highlighted that c? BC2N has a higher Vickers hardness (70 GPa) than c?BN (45-50 GPa)[54]. Extensive experimental work has been done in the study of B-C-N com- ponds. Solozhenko et al [55] obtained a bulk modulus of 282?15 GPa for c?BC2N using traditional axial x-ray diffraction (AXRD) under quasihydro- 20 static compression to 30 GPa. Employing a pressure of 18 GPa, temperature of 2200 K, and a graphitic BC2N starting material, they[55] have obtained a cubic structure that has a much lower density. The lattice constant and bulk modulus of this low-density c?BC2N are found to be about 3.642 A?and 282 GPa, respectively. Brillouin scattering measurements on the nanocrystalline cubic phase of BC2N have been successfully performed using the emulated platelet scatter- ing geometry[56] and they found a bulk and shear modulus of 259?22 GPa and 238?8 GPa, respectively. Metastable complexes of diamond-structured B-C-N complexes have been established with the synthesis of the metastable BC2N material[57]. Different computational approaches have also been employed in the study of B-C-N materials. Sun et al [58] used a first principle study of heterodia- mond BC2N structures and noted that all the structures are metastable. A first principle computational approach has also been applied by Luo et al [53] in the study of three possible configuration of wurzite BC2N . They noted the phase stability of one of the wurzite structures as an indication that it should be experimentally more feasible to be synthesize than zinc-blende-structured 21 BC2N . Azevado et al [59], using first principle calculations, have confirmed that the stable structure of boron ternary monolayers (BCN) is formed by increasing the number of both C-C and B-N bonds, and is independent of the unit cell size. Claims that a B-C-N complex structured material is the second hardest material, after diamond[55], have inspired further studies of hard materials that are potentially diamond-like in structure. In this work we explore prop- erties of cubic B-C-N structures, aiming to ascertain the claims of hardness, by studying the bulk and elastic properties of the system. We further look at the relative energies. In our work we consider trends in the B-C-N system. Each trend is sim- ulated with B and N atoms being positioned at the various lattice sites of the diamond lattice. We look at the trends of B4N4?xCx, B4?xN4?xC2x, and B4?xC4+x, (x = 0 . . . 4.). The approach considers both crystalline and random structures and in order to do this we use a bond order potential formalism employing the Tersoff approach[60, 61]. Chapter 2 Computational Details 2.1 Potential Functions A molecular dynamics simulation requires the definition of a potential func- tion, or a description of the terms by which the particles in the simulation will interact[62]. Molecular interactions are classified as intra-atomic or in- teratomic[63]. Examples of intramolecular potential functions include bond potentials, distance restraints, valence angle potentials, and dihedral angle potentials, where they are defined by a vector rij between two atoms i and j as shown in the figure 2.1[2]. 22 23 Figure 2.1: The interatomic bond vector[2]. 2.2 Interatomic potentials Interatomic potentials are potentials that involve more than one body. An important distinction between these and intramolecular (bonded) forces is that they are mainly specified by atom types rather than atom indices[2, 64]. examples of intermolecular interactions include Buckingham, Lennard-Jones, Stillinger-Weber, Axilrod-Teller[65] and the Tersoff potential[66] which is of particular interest to our study. 2.2.1 Tersoff potential The Tersoff potential is a density dependent potential, which has been de- signed to reproduce the properties of covalent bonding in systems containing carbon, silicon, germanium and their compounds of these elements[67, 2, 61, 68]. Tersoff[68] developed a pair potential whose strength depends on 24 the environment. This idea is similar to that of the ?glue model? in met- als, which uses the coordination of the atom as the variable controlling the energy[69]. This potential allows bond breaking and associated changes in bond hybridisation and it has 11 atomic and 2 bi-atomic parameters[2]. The Tersoff potential is a bond-order potential and it is composed of a two- body expansion which depends on the local environment[70]. E, the total energy is computed by summing the energy of site i (Ei) over all n atoms; where E = n? i=1 Ei = 1 2 n? i=1,i 6=j Uij (2.1) Where the term Uij represents the interaction energy between atoms i and j and is a combination of repulsive and attractive terms[68]. Uij = fc(rij)[fR(rij)? ?ijfA(rij)], (2.2) In equation 2.2, fR and fA are the repulsive and attractive pair potentials, respectively[67], and they have the following form: fR(rij) = Aij exp(?aijrij), fA(rij) = Bij exp(?bijrij) (2.3) 25 The term fc, on the other hand is the smooth cutoff function[68], which can be described as: fc(rij) = ? ??? ?????? ?????? ??? 1 : if rij < Rij 1 2 + 1 2 cos[pi (rij ?Rij) (rij ? Rij) ] : if Rij < rij < Sij 0 : if rij > Sij (2.4) The term ?ij expresses the dependence of each bond upon the local environ- ment and is lowered when the number of neighbours is relatively high[67]. ?ij = ?ij(1 + ??ii ??iij ) ?1 2?i (2.5) The effective coordination number of atom i is defined by the term ?ij[2],where; ?ij = ? k 6=i,j fC(rik)?ijkg(?ijk) (2.6) ?ijk here defines the bond angle between the vector rij and rik, as shown in figure 2.2.1. , and ?ijk can be further computed as [3], ?ijk = {~rij ? ~rik rijrik } (2.7) g(?ijk) = 1 + C 2 i d2i ? c2i d2i + (h? cos ?ijk)2 (2.8) The function g(?) has a minimum for h = cos ?, and the parameter d de- termines how sharp the dependence on the angle is, whilst c expresses the 26 Figure 2.2: A schematic diagram showing the bond angle ?[3]. strength of the angular effect[67]. The radial force which measures the re- sistance to stretch[3], is represented by the term ?, which can be expressed as; ?ijk = exp[? 3(rij?rjk)3] (2.9) The dependence on ?ijk and rij ? rik of the functions g and ? can be repre- sented by the graphs of figure 2.3 and 2.4. The further mixed parameters for a binary system can be defined as[2]; aij = ai + aj 2 , bij = bi + bj 2 (2.10) Aij = (AiAj)1/2, Bij = (BiBj)1/2 (2.11) Rij = (RiRj)1/2, Sij = (SiSj)1/2 (2.12) Here i, j, and k label the atoms in the system, rij is the length of the ij bond, 27 Figure 2.3: A graph showing the relationship ? and g(?) taken from[3]. Figure 2.4: A graph showing the dependace of ? on rij ? rik taken from[3]. 28 which can be calculated as; rij = |ri ? rj | (2.13) The 11 single subscripted parameters (e.g ai and ?i) depend on the atom type[2]. The two sets of bi-atomic parameters, ?ij and ?ij, define the chemistry between different types of atoms[2]: ?ii = 1, ?ij = ?ji, ?ii = 1, ?ij = ?ji (2.14) The Tersoff potentials are distinctively short ranged, typically of order 3 ?A[71]. These potentials have been found to give a convenient and relatively accurate description of the structural properties and energies of carbon, in- cluding elastic properties, phonons, polytypes, and defects and migration barriers in diamond and graphite [66, 72, 61, 68]. The Tersoff potentials work very well for zincblende group IV and group III-V semiconductors and their compounds [3]. 29 2.3 Previous studies using Tersoff potentials Studies have developed Tersoff potential parameters for C[60, 61] and B[72]. There are also two sets of N potentials developed by Matsunaga et al [72] and Fazzio et al [73]. The parameters are listed as shown in Table 2.1: The biggest challenge of the Tersoff potentials is that the fit is difficult. With 6 functions to fit and angular terms, it is difficult to find a good parametrization[69]. Carbon,Silicon and Germanium are some of the ma- terials that have been extensively studied[60, 70, 61] by fitting the tersoff potentials. In Table 2.2 we list the Tersoff parameters of Silicon[60, 70, 61], and germanium[61]. Using the Silicon parameters listed in Table 2.2 and ni- trogen set 2 parameters in Table 2.1, de Brito Mota et al [70] have shown the reliability of the potentials by testing them for a? SiNx for a wide range of nitrogen content(0 < x < 5), and the results are in agreement with available experimental data. Tersoff[68] also developed and improved[60] the Silicon parameters and they yield excellent results for the elastic constants of silicon, when compared with experiments. He also developed parameters for germanium[61], which he used to propose an empirical potential for multicomponent systems, in 30 Elements B[72] N(set1)[72, 74] N(set2)[73] C[60, 61] A(eV ) 2.7702? 102 1.10000? 104 6.36814? 103 1.3936? 103 11 B(eV ) 1.8349? 102 2.1945? 102 5.11760? 102 3.467? 102 ?(?A?1) 1.9922 5.7708 5.4367 3.4879 ?(?A?1) 1.5856 2.5115 2.7000 2.2119 ? 1.500? 10?6 1.0562? 10?1 5.2938? 10?3 1.5724? 10?7 n 3.9929 12.4498 1.33041 0.72751 c 5.2629? 10?1 7.9934? 104 2.03120? 104 3.8049? 104 d 1.5870? 10?3 1.3432? 102 2.55103? 101 4.384? 100 h 0.5000 ?0.9973 ?0.99734 ?5.62390? 10?1 R(A?) 1.8 2.0 1.8 1.8 S(A?) 2.1 2.3 2.1 2.1 Interactions(i? j) B-N B-C C-N ?i?j 1.1593 1.0025 0.9685 ?i?j 1.000 1.0000 0.6381 Table 2.1: Tersoff interatomic potential parameters for Carbon, Boron and Nitrogen 31 Elements Si[60, 70, 61] Ge[61] A(eV ) 1.8308? 103 1.769? 103 B(eV ) 4.7114? 102 4.1923? 102 ?(A??1) 2.4799 2.4451 ?(?A?1) 1.7322 1.7047 ? 1.1000? 10?6 9.0166? 10?7 n 7.8734? 10?1 7.5627? 10?1 c 1.0039? 105 1.0643? 105 d 1.6217? 101 1.5652? 101 h ?5.9825?1 ?4.3884? 10?1 R(?A) 2.7 2.8 S(?A) 3.0 3.1 Interactions(i ? j) Si-Ge C-Si ?i?j 1.00061 0.9776 Table 2.2: Tersoff interatomic potential parameters for Silicon and Germa- nium 32 particular C-Si and Si-Ge system. 2.4 Molecular Dynamics Computer simulations are carried out in the hope of understanding the prop- erties of assemblies of molecules in terms of their structure and the micro- scopic interactions between them[75]. This approach complements conven- tional experiments and enables us to study materials properties that are quite difficult (or even impossible) to find in other ways[8, 64]. They allow detailed investigations and enables us to gain a deeper understanding of the relevant materials processes with the materials[76]. It represents a very convenient and reliable interface between theory and laboratory experiments and can be thought of as a ?virtual experiment?[77]. Because molecular systems gener- ally consist of a vast number of particles, it is impossible to find the properties of such complex systems analytically[8], hence the need for approximation. These accuracy of the predictions is subject to the accuracy of the model used Computer simulations are mainly used as an intermediate between micro- scopic length and time scales and the macroscopic world of the laboratory[8] 33 and provides a guess of the interactions between molecules. There numerous simulation techniques, such as Molecular Dynamics (MD) and Monte Carlo (MC), in addition, there are hybrid techniques which com- bine MD and MC[8]. In our study the MD technique is used. MD is one of the principal tools in the theoretical study of molecules[77]. This is a computational method which simulates the time dependent be- haviour of molecular systems[8]. It is a computer technique which solves numerically a multi-body problem of mechanics[67]. This is a form of com- puter simulation in which atoms and molecules are allowed to interact for a period of time by approximations of known physics, giving a view of the mo- tion for the coordinates components of the N atoms in the assembly[77, 78]. Researchers[79] surmise that the molecular dynamics algorithm in most com- mon use today may even have been known to Newton. MD simulation circumvents the difficulty of finding an analytical solution in the complex systems by using numerical methods[80]. It solves numerically the 3N simultaneous equations of motion[80, 81], by using time integration algorithms (TIA). MD is a multidisciplinary method which probes the rela- 34 tionship between molecular structure, movement and function[77]. The basic laws and theories of MD are from mathematics, physics, and chemistry, and it employs algorithms from computer science and information technology. It was originally conceived within theoretical physics in the late 1950?s[80]. MD can be applied in effectively describing the velocities and positions of atoms disrupted in a structure, as a function of time[67]. The microscopic behaviour of the system can be computed by solving differential equations of motion, from the initial position and velocity of the atoms, and the forces of interaction between them[80]. A MD simulation requires the definition of a potential function[77, 82], or a description of the terms by which the particles in the simulation will interact and is usually referred to as a force field. Those most commonly used are based classical mechanics and embody a classical treatment of particle-particle interactions that can reproduce structural and conforma- tional changes but usually cannot reproduce chemical reactions[80]. The Born-Oppenheimer approximation[83], which states that the dynamics of electrons is so fast that they can be considered to react instantaneously to the motion of their nuclei,and hence the electron and the nuclei may be 35 treated separately. The nuclei may then be treated as point particles which follow classical Newtonian dynamics[77]. these are the two key factors that leads to the reduction from a fully quantum description to a classical po- tential. As a result the effect of the electrons is approximated as a single potential energy surface, usually representing the ground state[80]. 2.4.1 Molecular Interactions Molecular dynamics simulation consists of the numerical, step-by-step, solu- tion of the classical equations of motion, which for a simple atomic system may be written as[62]: mir?i = fi (2.15) fi = ???riU (2.16) For this purpose we need to be able to calculate the forces fi acting on the atoms, and these are generally derived from a potential energy U(rN), where rN = (r1, r2, . . . , rN .) represents the complete set of 3N atomic coordinates. Well known MD softwares includes: ? VASP[84] ? CASTEP[85] 36 ? DL POLY[2] ? NEWTON-X[86] ? AMBER[87] ? GULP[65] ? NAMD[88] ? CHARMM[89] The first two software packages do Quantum MD. In our work we use the General Utility Lattice Program(GULP)[65]. 2.5 The Time Integration Time integration algorithm (TIA) is the engine of MD. It is basically required to integrate the equations of motion and follow the trajectory of interacting particles[69, 78]. TIA is based on a finite difference method[90], which em- ploys the discretization of time on a finite grid, and the distance between two consecutive points on the grid being the time step, ?t. 37 There are a range of TIA available. There is the Verlet algorithm[91], whose basic idea is to have to Taylor expansion of third order for the posi- tions r(t). One is forward, whilst the other is backward in time. The are two versions of the Verlet algorithm namely the Verlet leapfrog(LF) or Ve- locity Verlet(VV) versions[2]. Another algorithm is the Predictor-Corrector algorithm[92], which consist of 3 steps, namely the predictor, force evaluation and the corrector. 2.6 Long-range Interaction Algorithms These methods compute the interaction energies of periodic systems, like crystals, particularly the electrostatic energies. These are[78]: ? Ewald summation ? Particle Mesh Ewald ? Particle-Particle Particle Mesh ? Reaction Field Method ? Central Difference method ? Euler-Cauchy Method 38 ? Nordsieck Method for Newton?s Equations 2.6.1 Ewald technique The Ewald technique represents the most efficient solution for crystalline materials[93] and it is a well-established technique for the evaluation of electrostatic interactions which are necessary for many computer simula- tions. Ewald summation is a special case of the Poisson summation for- mula which replaces the summation of interaction energies which are in real space with an equivalent summation in Fourier space. Its advantage is rapid convergence[93]. 2.7 The GULP code Gulp is a molecular dynamics code capable of handling up to multi-thousand atoms[10]. It is designed to perform a variety of tasks based on force field methods. It is a flexible tool with great variety and diversity for simulation of solids based on interatomic potential models[65], including the Tersoff potentials. In GULP many algorithms have been embodied for the simulation of three-dimensional periodic systems[65]. It is suitable for the treatment of 39 both inorganic and organic systems with fully flexible molecules[10]. One of the principal applications of the program has been to the derivation of empirical potential parameters through least-squares fitting. The ability of GULP to treat multiple structures within the same run[65] is of particular use in our study. Algorithms for symmetry-adapted energy minimisation of solids using an- alytical first and second derivatives have been included and implemented in GULP[65]. The new methods lead to improved computational efficiency of up to an order of magnitude over the standard algorithm which takes no account of symmetry. The largest improvement is obtained from the use of symmetry in the generation of the hessian[10]. Accelerated convergence techniques for the dispersion energy are found to be of great benefit in im- proving the precision at little extra computational cost, particularly when a one centre decomposition is possible or the Ewald sum weighting towards real-space is increased[65]. 40 Energy Minimisation As a prerequisite for any subsequent evaluation of physical properties, ef- ficient minimisation of the energy is an essential part of the simulation of solids[65]. This normally represents the computationally most demanding stage. Gulp is is able to optimise at constant pressure or at constant volume, where the unit cell remains frozen[65]. Calculating bulk properties, GULP uses equations due to Voight, Reuss and Hill[94]. Below are the Reuss and Voight definitions, while the Hill values are defined as the average of the Reuss and Voight values: BV oight = 1 9 (C11 + C22 + C33 + 2(C12 + C13 + C23)) (2.17) BReuss = (S11 + S22 + S33 + 2(S12 + S13 + S23))?1 (2.18) GV oight = 1 15 (C11 + C22 + C33 + 3(C44 + C55 + C66)? C12 ? C13 ? C23) (2.19) GReuss = 15 4(S11 + S22 + S33 ? S12 ? S13 ? S23) + 3(S44 + S55 + S66) (2.20) In this work we employ the Voight convention. The bulk modulus B and first derivative B? can also obtained by monitoring the ratio VV0 of volume, V, and equilibrium volume, V0, over a range of pressures. A numerical procedure 41 can then be used to fit the values of B and B? to the Birch[11] equation of state which is given by; P = 3 2 B0 [( V V0 )?7 3 ? ( V V0 )?5 3 ][ 1 + 3 4 (B? ? 4) ( V V0 )?2 3 ? 1 ] (2.21) Chapter 3 Application of the Tersoff Potential to B-C Complexies An approach using Tersoff potentials is employed in the GULP code, to determine Bulk and Elastic properties of the B ? C complex materials. 3.1 B-C Crystal Structures A possible diamond-related phase may be formed from a simple eight-atom unit cell of the diamond lattice with B or C atoms placed at various lattice points according to the nominal lattice stoichiometry[20]. Figure 3.1 shows 42 43 the B-C structures considered in the study. Supercells of various sizes can be generated from these structures as discussed later in the thesis. 44 BC B3C5 BC3 ? a BC3 ? b BC7 Diamond Figure 3.1: B-C structures considered in this work. Dark yellow and white represent Boron and Carbon atoms, respectively. 45 3.2 Results The study has been conducted for the crystalline and random phases of these structures at an absolute zero temperature (0k). The random structures are created by generating a random distribution of atoms over the crystal lattice sites then letting the geometry relax. The properties simulated are elastic constants, bulk modulus, shear modulus, Young?s modulus, and total formation energies for the structures. A number of supercells were used as GULP input, namely; 1?1?1(8 atoms), 4?4?4(512 atoms), 5?5?5(1000 atoms), and 6?6?6(1728). We run MD for both the crystalline and random phases. A structural optimisation (OPT) procedure[2] was also performed to improve the quality of the starting structure prior to a dynamical simulation. The properties predicted by MD and OPT are in agreement for a smaller supercell( 8-atoms) and differ as the size of the supercell increases. The OPT results is bigger for elastic constants and smaller for the total energy. Tables 3.1 to 3.4 show the predictions of MD and OPT of the total ground state energy (referred to as energy in the table), elastic constants( C11, C12, and C44), bulk modulus (B), shear modulus (G), Young?s modulus () and the Poisson?s ratio (?). 46 Supercell Compound Energy (eV/atom) C11(GPa) C12(GPa) C44(GPa) B(GPa) G(GPa) ?(GPa) ? 111 BC -7.151 598 108 344 271 304 565 0.1524 B3C5 -7.202 690 112 399 305 355 658 0.1400 BC3(b) -7.175 786 115 432 339 393 754 0.1285 BC3(a) -7.252 825 101 456 342 417 799 0.1031 BC7 -7.313 932 107 545 382 492 910 0.1028 Diamond -7.370(?10.125)a 1074 102 642(576)c 426(433 ? 442)b 579 1056 0.0865 444 BC -7.145 598 108 344 271 304 565 0.1524 B3C5 -7.195 690 112 399 305 355 658 0.1400 BC3(b) -7.174 786 115 432 339 393 754 0.1285 BC3(a) -7.252 825 101 456 342 417 799 0.1031 BC7 -7.307 932 107 545 382 492 910 0.1028 Diamond -7.371(?10.125)a 1074 102 642(576)c 426(433 ? 442)b 579 1056 0.0865 555 BC -7.145 598 108 344 271 304 565 0.1524 B3C5 -7.196 690 112 399 305 355 658 0.1400 BC3(b) -7.196 786 115 432 339 393 754 0.1285 BC3(a) -7.252 825 101 456 342 417 799 0.1031 BC7 -7.308 932 107 545 382 492 910 0.1028 Diamond -7.370(?10.125)a 1074 102 642(576)c 426(433 ? 442)b 579 1056 0.0865 666 BC -7.145 598 108 344 271 304 565 0.1524 B3C5 -7.173 690 112 399 305 355 658 0.1400 BC3(b) -7.196 786 115 432 339 393 754 0.1285 BC3(a) -7.252 825 101 456 342 417 799 0.1031 BC7 -7.308 932 107 545 382 492 910 0.1028 Diamond -7.370(?10.125)a 1074 102 642(576)c 426(433 ? 442)b 579 1056 0.0865 BC5 (335 ? 8GPa)d (a) Referecce[95] (b) Referene[96, 97] (c) Reference[98, 99] (d) Reference[23] Table 3.1: Properties of crystalline phase structures predicted using OPT. Values in brackets are experimental values. 47 Supercell Compound Energy(eV/atom) C11(GPa) C12(GPa) C44(GPa) B(GPa) G(GPa) ?(GPa) ? 111 BC -6.9125 538 122 255 260 236 490 0.1866 B3C5 -7.202 690 113 399 305 355 658 0.1403 BC3(b) -7.175 786 115 432 339 393 754 0.1284 BC3(a) -7.252 825 101 456 342 417 799 0.1033 BC7 -7.313 932 107 546 382 492 910 0.1031 Diamond -7.375(?10.125)a 1074 102 642(576)c 426(433 ? 442)b 579 1056 0.0865 444 BC -7.07 619 129 352 292 309 574 0.1731 B3C5 -7.150 692 121 403 312 352 651 0.1469 BC3(b) -7.151 777 115 425 336 386 745 0.1299 BC3(a) -7.227 820 102 455 343 414 792 0.103 BC7 -7.297 927 108 543 381 489 905 0.1042 Diamond -7.367(?10.125)a 1072 102 640(576)c 425(433 ? 442)b 575 1054 0.0868 555 BC -7.026 611 132 350 292 306 564 0.1778 B3C5 -7.127 696 127 406 319 354 654 0.1526 BC3(b) -7.137 776 116 422 337 384 742 0.1314 BC3(a) -7.212 818 104 455 345 413 788 0.1038 BC7 -7.29 924 109 541 380 488 901 0.1054 Diamond -7.362(?10.125)a 1069 102 639(576)c 424(433 ? 442)b 577 1052 0.0871 666 BC -7.742 418 148 126 234 145 226 0.2094 B3C5 -6.950 577 142 298 289 265 527 0.1953 BC3(b) -7.172 760 121 431 335 386 734 0.1379 BC3(a) -7.172 766 119 430 335 387 734 0.1338 BC7 -7.289 919 110 537 380 484 894 0.1071 Diamond -7.370(?10.125)a 1074 102 642(576)c 426(433 ? 442)b 579 1056 0.08651 (a) Referecce[95] (b) Referene[96, 97] (c) Reference[98, 99] Table 3.2: Properties of crystalline phase structures predicted using MD. Values in brackets are experimental values 48 Supercell Compound Energy (eV/atom) C11(GPa) C12(GPa) C44(GPa) B(GPa) G(GPa) ?(GPa) ? 111 BC -6.9125 541 123 256.4 262 238 493 0.186 B3C5 -7.0375 595 127 370 298 308 642 0.1904 BC3(b) -7.175 785 115.5 432 339 393 753 0.1294 BC3(a) -7.252 759 116 483 342 417 799 0.1371 BC7 -7.313 932 107 546 382 492 910 0.1031 Diamond -7.362(?10.125)a 1074 102 642(576)c 426(433 ? 442)b 579 1056 0.08653 444 BC -6.779 458 151 158 247 163 344 0.2436 B3C5 -6.994 606 136 309 290 280 541 0.1843 BC3(b) -7.149 758 123 418 333 377 712 0.1391 BC3(a) -7.164 761 122 430 334 385 721 0.1388 BC7 -7.277 916 111 533 378 481 886 0.108 Diamond -7.367(?10.125)a 1072 102 640(576)c 425(433 ? 442)b 578 1054 0.0868 555 BC -6.756 430 156 179 247 162 351 0.2698 B3C5 -6.961 589 137 304 287 272 541 0.1878 BC3(b) -7.14 723 124 422 333 378 712 0.1404 BC3(a) -7.128 750 123 416 331 374 711 0.1402 BC7 -7.266 906 112 528 377 476 884 0.1099 Diamond -7.362(?10.125)a 1069 102 639(576)c 424(433 ? 442)b 577 1052 0.0871 666 BC -6.7419 418 148 126 234 145 226 0.2093 B3C5 -6.946 577 142 298 289 265 527 0.1953 BC3(b) -7.120 744 128 419 335 374 715 0.147 BC3(a) -7.122 746 125 417 333 375 713 0.1426 BC7 -7.250 905 113 529 377 476 880 0.1108 Diamond -7.358(?10.125)a 1067 102 638(576)c 424(433 ? 442)b 575 1049 0.08741 (a) Referecce[95] (b) Referene[96, 97] (c) Reference[98, 99] Table 3.3: Properties of the random phase structures predicted using MD. Values in brackets are experimental values 49 Supercell Compound Energy (eV/atom) C11(GPa) C12(GPa) C44(GPa) B(GPa) G(GPa) ?(GPa) ? 111 BC -6.916 542 122 258 262 238 494 0.185 B3C5 -7.043 596 126 370 298 309 643 0.189 BC3(b) -7.175 785 115 432 338 393 754 0.1287 BC3(a) -7.252 759 116 483 342 417 799 0.1372 BC7 -7.309 932 107 545 382 492 910 0.1028 Diamond -7.370(?10.125)a 1074 102 642(576)c 426(433 ? 442)b 579 1056 0.08651 444 BC -6.836 479 143 202 250 192 380 0.2331 B3C5 -7.027 618 132 320 292 289 559 0.1775 BC3(b) -7.170 767 120 426 335 384 725 0.1353 BC3(a) -7.186 769 120 435 336 391 733 0.1357 BC7 -7.287 921 110 536 380 483 893 0.1072 Diamond -7.371(?10.125)a 1074 102 642 (576)c 426 (433 ? 442)b 579 1056 0.08651 555 BC -6.830 469 144 211 254 191 410 0.2339 B3C5 -7.009 608 132 316 291 284 567 0.1791 BC3(b) -7.175 766 120 431 335 387 729 0.1346 BC3(a) -7.159 763 120 425 334 383 728 0.1355 BC7 -7.284 915 111 534 379 481 894 0.1082 Diamond -7.371(?10.125)a 1074 102 642(576)c 426(433 ? 442)b 579 1056 0.08651 666 BC -6.836 459 143 214 249 192 400 0.2407 B3C5 -7.021 609 132 320 291 287 567 0.1785 BC3(b) -7.172 760 121 431 335 386 734 0.1379 BC3(a) -7.172 766 119 430 335 387 734 0.1338 BC7 -7.289 919 110 537 380 484 894 0.1071 Diamond -7.370(?10.125)a 1074 102 642(576)c 426(433 ? 442)b 579 1056 0.08651 (a) Referecce[95] (b) Referene[96, 97] (c) Reference[98, 99] Table 3.4: Properties of the random phase structures predicted using OPT. Values in brackets are experimental values. 50 We note the difference in the calculated values from experimental results in diamond for the total ground state energy and C44, but there is agreement for B. Cell parameters Cell parameters (a, b, and c) and cell angles (?, ?, and ?) for all structures are examined. A conventional cell parameter average (a?) is also calculated, using equation 3.1; a? = a + b+ c 3 (3.1) Taking the average here is justified since the structures are almost cubic, as can be seen from table 3.5. 51 Crystal MD a(?A) b(?A) c(?A) a?(?A) ?? ?? ?? BC 3.648 3.648 3.6485 3.648 90 90 90 B3C5 3.635 3.6346 3.635 3.635 90 90 90 BC3(b) 3.601 3.601 3.601 3.601 89.69 89.69 90.32 BC3(a) 3.603 3.603 3.620 3.609 90 90 90 BC7 3.586 3.586 3.586 3.586 90 90 90 Diamond 3.566 3.566 3.566 3.566(3.567)a 90 90 90 Crystal OPT BC 3.662 3.662 3.662 3.662 90 90 90 B3C5 3.634 3.634 3.634 3.634 90 90 90 BC3(b) 3.602 3.602 3.602 3.602 89.68 89.68 90 BC3(a) 3.603 3.603 3.620 3.609 90 90 90 BC7 3.586 3.586 3.586 3.586 90 90 90 Diamond 3.5656 3.5656 3.5656 3.5656(3.567)a 90 90 90 Random MD BC 3.648 3.648 3.649 3.648333333 89 89 89 B3C5 3.641 3.615 3.615 3.624 89.2 90 90 BC3(b) 3.601 3.603 3.602 3.602 89.7 90.3 89.7 BC3(a) 3.620 3.603 3.603 3.609 90 90 90 BC7 3.586 3.586 3.586 3.586 90 90 90 hline Diamond 3.566 3.566 3.566 3.566(3.567)a 90 90 90 Random OPT BC 3.648 3.648 3.648 3.648 89 89 89 B3C5 3.643 3.614 3.614 3.624 89.2 90 90 BC3(b) 3.602 3.602 3.602 3.602 89.7 90.3 89.7 BC3(a) 3.620 3.603 3.603 3.608666667 90 90 90 BC7 3.586 3.586 3.586 3.586 90 90 90 Diamond 3.566 3.566 3.566 3.566(3.567)a 90 90 90 BC5 (3.635? 0.006)d (a) Referecce[95] (d) Reference[100] Table 3.5: Cell parameters predicted for the B-C structures. Values in brack- ets are experimental values 52 3.2.1 Graphical representation of the results for the B-C crystal structures A graphical analysis of some properties predicted in this work is performed to examine the various trends described below. Plots of the number of each type of bonds (B-C, C-C, and B-B bonds) against the B-C trend are done as shown in figure 3.2. Figures 3.3 to 3.4 show the total energy dependence on the trend while the graphs in figure 3.5 to 3.6 show the dependence on the B-C trend of the elastic parameters, B, G, and . The variation of the elastic constants with the B-C trend is graphically represented in figure 3.7 to 3.8, while the calculated cell constants are also plotted against the trend as shown in figure 3.9. For the Random Phases it is important that the size of the supercell chosen for the randomization process does not affect the final result. To confirm this is the case, various sizes of the supercell were tested to ensure convergence and the value of the moduli are as shown in figure 3.10, where the numbers 0, 1, 2 and 3 on the horizontal axis represent systems containing 8, 512, 1000, and 1728 atoms, respectively. It should be noted that in figure 3.10 BC3 53 refers to BC3 ? b because it is the lower energy of the two BC3 structures. Figure 3.2: Percentage of the number of bonds versus the composition for a random phase 54 Energy difference and ratio of lattice constants The difference in energy (4E) between the random structures and their crystalline countepart is also explored, Where; 4E = Ecrystal ? Erandom (3.2) Ecrystal = energy per atom for the crystalline phase (3.3) Erandom = energy per atom for the random phase (3.4) Figure 3.11 shows the trend in 4E as the number carbon-carbon bonds increases. The effect of the random distribution of boron and carbon on the lattice constants is studied. This is done by taking the ratio of the cell parameters for each phase ( i.e a?random a?crystal ) and is shown in figure 3.12. The ground state energy is seen to be decreasing with the concentration and so are the elastic constants and the Poisson ratio. Values of the cell constants are decreasing with the concentration. 55 (a) (a) Figure 3.3: Energy for crystal structures. (a) is for OPT, while (b) is for MD. The horizontal axes are showing the composition 56 (a) (b) Figure 3.4: Energy for random structures. (a) is for OPT, while (b) is for MD. The horizontal axes are showing the composition 57 (a) (b) Figure 3.5: B, G, and  for crystal structures. (a) is for OPT, while (b) is for MD. The horizontal axes are showing the composition 58 (a) (b) Figure 3.6: B, G, and  for random structures. (a) is for OPT, while (b) is for MD. The horizontal axes are showing the composition 59 (a) (b) Figure 3.7: Elastic constants for crystal structures. (a) is for OPT, while (b) is for MD. The horizontal axes are showing the composition 60 (a) (b) Figure 3.8: Elastic constants for the random structures. (a) is for OPT, while (b) is for MD. The horizontal axes are showing the composition 61 Figure 3.9: Cell constant against Composition 62 (a) (b) (c) (d) Figure 3.10: Bulk(?), shear( ) and Young(5) moduli for various sizes of the supercells with random B and C distribution. (a), (b), (c) and (d) represents BC, B3C5, BC3 and BC7, respectively. The horizontal axes are showing the composition 63 Figure 3.11: Energy difference of the crystal and random phase against B-C trend Figure 3.12: Ratio of the random and crystal lattice cell constants Chapter 4 Tersoff Potential Application to B-C-N Complexes Using the Gulp code, the Tersoff potentials were applied to study the trends of the B-C-N complexies as highlighted in Chapter 1. Each trend consist of 5 structures. The trends are are in relation to the composition of B, C and N in the structures and they are; BC to BN , BN to C, and C to BC. The trend B4?xC4+x, (x = 0 . . . 4.) is as outlined in the previous chapter. Only BC3?b is considered for the BC3 structure since it is of lower energy. Figure 4.1 shows the structures of the trendB4N4?xCx, (x = 0 . . . 4.), while figure 4.2 shows the structures of the trend B4?xN4?xC2x, (x = 0 . . . 4.). BN and 64 65 Diamond have been shown previously and are thus omitted here. 66 BC B4C3N B2CN B4CN3 BN Figure 4.1: B-C-N structures considered for the composition B4N4?xCx, (x = 0 . . . 4.) , dark yellow, white and yellow represent Boron, Carbon and Nitrogen atoms, respectively. 67 B3C2N3 BC2N BC6N Figure 4.2: B-C-N structures considered for the composition B4?xN4?xC2x, (x = 0 . . . 4.) , dark yellow, white and yellow represent Boron, Carbon and Nitrogen atoms, respectively. 68 4.1 Comparison of the predictions of the Two sets of Tersoff Potential Parameters for Nitrogen Since there are two sets of Tersoff parameters for Nitrogen as tabulated in Table2.1,The quality of the predictions of the two sets is first looked into. Considering a 512 atoms supercell, we employ the two diffent Tersoff poten- tials, for the simulation of the structures, using only OPT. Tables 4.1 and 4.2 show the results for the total ground state formation energy, C11, C12, C44, Bulk modulus (B), Shear modulus (G), Young?s modulus (), and the Poisson ratio (?) using parameters Set 1. This was also done for Set 2 and the results are tabulated in Tables 4.3 and 4.4. 69 Structure Energy (Ev/atom) C11 (GPa) C12 (GPa) C44 (GPa) B (GPa) G (GPa) ? (GPa) ? BC -7.14 598 108 344 271 304 565 0.152 B4C3N -7.02 560 164 288 296 252 485 0.23 B2CN -6.89 541 230 260 322 207 403 0.33 B4CN3 -6.77 506 273 205 351 169 314 0.35 BN -6.65 494 331 176 386 138 228 0.40 (820)a (190)a (480)a (400)a B3C2N3 -7.04 1662 755 959 1056 752 846 0.28 BC2N -7.38 955 61.4 1152 817.5 656 946.4 (355)a ,(282)b BC6N -7.10 979 169 531.2 439 480 924 0.15 Diamond(C8) -7.37 1074 102 641 426 579 1056 0.087 BC7 -7.31 932 107 545 382 492 910 0.103 BC3 -7.25 825 101 456 342 417 799 0.103 B3C5 -7.20 690 112 399 305 355 658 0.140 BC -7.14 598 108 344 271 304 565 0.152 (a) Referecce[101] (b) Referecce[55] Table 4.1: Predicted properties of crystal phase using Parameter Set 1 of B-C-N. Values in brackets are experimental values. Energy here refers to the ground state total formation energy 70 Structure Energy (Ev/atom) C11 (GPa) C12 (GPa) C44 (GPa) B (GPa) G (GPa) ? (GPa) ? BC -6.84 466 139 220 249 197 408 0.23 B4C3N -6.96 342 152 235 261 184 332 0.281 B2CN -7.20 508 235 248 346 211 374 0.31 B4CN3 -7.26 657 227 325 395 261 562 0.26 BN -6.95 380 217 140 399 121 201 0.364 B3C2N3 -7.20 695 329 362 485 304 549 0.29 BC2N -7.41 818 242 470 445 393 670 0.24 BC6N -7.25 970 210 552 449 492 900 0.16 Diamond(C8) -7.37 1074 102 642 426 579 1057 0.086 BC7 -7.29 918 110 537 380 483 900 0.107 BC3 -7.21 781 117 440 338 397 746 0.13 B3C5 -7.02 627 121 308 291 289 573 0.17 BC -6.84 466 139 220 250 197 408 0.23 Table 4.2: Predicted properties of random phase using Parameter Set 1 of B-C-N. Energy here refers to the ground state total formation energy 71 Structure Energy (Ev/atom) C11 (GPa) C12 (GPa) C44 (GPa) B (GPa) G (GPa) ? (GPa) ? BC -7.15 598 108 344 271 304 565 0.15 B4C3N -6.84 663 101 380 288 341 636.88 0.132 B2CN -6.53 721 94 431 307 382 715 0.11 B4CN3 -6.23 813 83 471.1 327 429 798 0.092 BN -5.92 896 73.2 526 347 480 885 0.076 B3C2N3 -6.20 911 93.6 534 366 484 894 0.093 BC2N -6.53 939 82.9 555 379 497 934 0.097 BC6N -6.93 989 101 580 397 526 970 0.092 Diamond(C8) -7.37 1074 102 642 426 579 1057 0.087 BC7 -7.31 932 107 545 382 492 910 0.103 BC3 -7.25 825 101 456 342 417 799 0.103 B3C5 -7.20 690 112 399 305 355 659 0.140 BC -7.14 598 108 344 271 304 565 0.152 Table 4.3: Predicted properties of crystal phase using Parameter Set 2 of B-C-N. Energy here refers to the ground state total formation energy 72 Structure Energy (Ev/atom) C11 (GPa) C12 (GPa) C44 (GPa) B (GPa) G (GPa) ? (GPa) ? BC -6.76 448 148 219 250 189 384 0.241 B4C3N -6.43 525 137 222 266 215 466 0.21 B2CN -6.14 535 135 360 290 268 622 0.169 B4CN3 -5.89 745 107 348 319 336 719 0.125 BN -5.85 896 73 526 347 480 885 0.076 B3C2N3 -6.012 859 106 452 356 423 834 0.106 BC2N -6.22 846 114 517 366 449 875 0.109 BC6N -6.88 970 102 569 393 515 953 0.095 Diamond(C8) -7.35 1074 102 642 426 579 1056 0.0865 BC7 -7.25 915 111 537 380 484 893 0.108 BC3 -7.12 759 125 432 335 386 723 0.137 B3C5 -6.95 602 137 319 291 284 543 0.181 BC -6.76 448 148 219 250 189 384 0.241 Table 4.4: Predicted properties of random phase using Parameter Set 2 of B-C-N. Energy here refers to the ground state total formation energy 73 The Moduli values from the four different structural types simulated are plotted against the B-C-N composition, as shown in Table 4.5. The per- centage of number bonds are also examined into. The plotted graphs are as shown in figures 4.3 and 4.4. 74 Trend structure C-C B-N C-N N-N B-B B-C BC to BN BC 30 0 0 0 30 40 B4C3N 16 10 8 0 33 33 B2CN 5 23 12 5 34 22 B4CN3 0 33 11 10 33 13 BN to C BN 0 50 0 25 25 0 B3C2N3 5 27 17 17 17 17 BC2N 36 11 22 6 6 22 BC6N 70 2 14 0 0 14 C to BC Diamond 100 0 0 0 0 0 BC7 86 0 0 0 1 12 BC3 69 0 0 0 4 27 B3C5 49 0 0 0 14 37 BC 25 0 0 0 25 50 Table 4.5: Table showing the B-C-N structures?s number of bonds as a per- centage of the total number of bonds of a structure 75 (a) (b) Figure 4.3: Bulk Moduli and elastic constants in units of GPa as shown on the vertical axes, with Set 1 parameters. (a) and (b) represent the crystal and random phase, respectively. The horizontal axes represent the composition of the structures 76 (a) (b) Figure 4.4: Bulk Moduli and elastic constants in units of GPa as shown on the vertical axes, with Set 2. (a) and (b) represent the crystal and random phase, respectively. The horizontal axes represent the composition of the structures 77 4.2 Further Results for B-C-N with Parame- ter Set 2 Parameter Set 2 for nitrogen appears to be yield better results than Parame- ter Set 1. This work is extended by only considering the Set 2 parameters.Set 1 parameters is unphysical as seen in comparing figure 4.3 with figure 4.4 com- pared with Set . the work was then furthered by only considering Parameter Set 2. The supercell size (number of atoms) is varied in the GULP input. 4.2.1 Crystalline Phases results The 1 ? 1 ? 1 (8 atoms) is considered for the crystal phase. Tables 4.6 and 4.7 shows properties obtained. Two approaches are used for calculating the cell constants. The first is to find the average of the 3 cell constants, as done in the previous chapter, shown in equation 4.1.In the other approach, the square root of the total volume, V, per atom, shown in equation 4.2. A dependance of these properties on the B-C-N trend is clearly depicted in the figures 4.5 to 4.6. We note that the Total formation energy per atom is maximum at BN and minimum at diamond. the plot of the cell constants against the composition is showing an almost parabolic curve, with a turning 78 Structure a1(?A) a2(?A) ?0 ?0 ?0 BC 3.662 3.662 90 90 90 B4C3N 3.651 3.651 90 90 90 B2CN 3.64 3.64 90 90 90 B4CN3 3.63 3.63 90 90 90 BN 3.62(3.617)a 3.62(3.617)a 90 90 90 B3C2N3 3.59 3.59 89.76 90.24 90.24 BC2N 3.58(3.602)a , (3.64)b 3.58(3.602)a , (3.64)b 90 90 90 BC6N 3.568 3.567 90.3 90.3 89.7 Diamond(C8) 3.566 3.565 90 90 90 BC7 3.586 3.585 90 90 90 BC3 3.608 3.609 90 90 90 B3C5 3.634 3.634 90 90 90 BC 3.662 3.662 90 90 90 (a) Reference[101] (b) Reference[55] Table 4.6: Table of properties of crystal phase of B-C-N structures, Values in brackets are experimental values point minimum at Diamond. a1 = a + b+ c 3 (4.1) a2 = ? V (4.2) 79 Structure Energy(eV/atom) C11 C12 C44 B G ? ? BC -7.14 598 108 344 271 304 565 0.152 B4C3N -6.84 663 101 380 288 341 637 0.13 B2CN -6.54 721 94 431 307 382 722 0.112 B4CN3 -6.23 813 83.3 471 327 429 798 0.929 BN -5.9 896 73 526 347 480 885 0.076 (820)a (190)a (480)a (400)a (409)c (973)c B3C2N3 -6.2 911 93.7 535 366 484 894 0.0932 BC2N -6.538 940 83 555 379 497 934 0.0965 (355)a, (282)b BC6N -6.925 989 101 580 397 526 970 0.0924 Diamond(C8) -7.37 1074 102 642 426 579 1056 0.08651 BC7 -7.309 932 107 545 382 492 910 0.103 BC3 -7.25 825 101 456 342 417 799 0.122 B3C5 -7.196 690 112 399 305 355 658 0.14 BC -7.144 598 108 344 271 304 565 0.152 (a) Reference[101] (b) Reference[55] (c) Reference[18] Table 4.7: Cell properties of crystal phase of B-C-N structures. The Elastic properties are expresse in units of GPa. Energy refers to the total formation energy. Values in brackets are experimental values 80 (a) (b) Figure 4.5: (a) Elastic constants, and (b) Moduli for the crystal phase struc- tures. The horizontal axes represent the trend in the B-C-N composition 81 (a) (b) Figure 4.6: The total formation energies( in part (a)) and the cell constants (in part (b)), for the crystal structures with Parameter Set 2. The horizontal axes in both graphs represent the composition in the B-C-N trend 82 4.2.2 Random Phases Results A random 4? 4? 4 (512 atoms) and 7? 7? 7 ( 2744 atoms) are considered. Tables 4.4, 4.8,and 4.9 show the properties obtained. Similarly, a dependance of these properties on the B-C-N composition is clearly shown in the figures 4.4, 4.7- 4.9. 83 Structure Energy( Ev/atom) C11 C12 C44 B G ? ? BC -6.8396 466.12 142 212 249.78 192 394 0.238 B4C3N -6.4972 523 137 225 265 215 463 0.21 B2CN -6.195 527 135 360 290 264 577 0.171 B4CN3 -5.9639 747 106 355 320 341 721 0.123 BN -5.8473 895 73 526 347 480 885 0.0755 B3C2N3 -6.0470 862 105 459 357 427 840 0.108 BC2N -6.4275 849 113 522 368 451 884 0.112 BC6N -6.8965 971 102 564 392 512 951 0.0951 Diamond(C8) -7.3706 1074 102 642 426 579 1056 0.08651 BC7 -7.2521 918 110 537 380 484 895 0.107 BC3 -7.1224 767 120 429 334 387.3 728 0.135 B3C5 -6.9603 620 132 322 293 292 563 0.173 BC -6.8396 466.12 142 212 249.78 192 394 0.238 Table 4.8: Properties of random, 2744 atoms, phase of B-C-N structures. The Elastic properties are expressed in units of GPa. Energy refers to the total formation energy. 84 Structure a1(?A) a2(?A) ?0 ?0 ?0 BC 3.644 3.644 90 89.95 90.07 B4C3N 3.627 3.627 89.99 90.07 89.97 B2CN 3.615 3.615 90.06 90 90.04 B4CN3 3.61 3.61 89.98 90 90.05 BN 3.62(3.617)a 3.62(3.617)a 90 90 90 B3C2N3 3.58 3.583 89.96 90.02 89.99 BC2N 3.57(3.602)a, (3.64)b 3.57(3.602)a, (3.64)b 90.07 90 90.01 BC6N 3.567 3.566 89.98 89.98 89.98 Diamond(C8) 3.57 3.57 90 90 90 BC7 3.585 3.585 90 90 89.99 BC3 3.6038 3.604 89.99 90.01 90 B3C5 3.623 3.623 89.98 89.99 90.01 BC 3.644 3.644 90 89.95 90.07 (a) Reference[101] (b) Reference[55] Table 4.9: Properties of, 2744 atoms, random, phase of B-C-N structures, using Parameter set 2. Values in brackets are experimental values 85 (a) (b) Figure 4.7: The elastic constants (shown in part (a), and the Moduli (shown in part (b)) for the 2744 atoms, random structure, using Set 2. The horizontal axes represent the B-C-N composition 86 (a) (a) Figure 4.8: Total energy (in part (a)) and cell constants for the random, structures, with 2744 atoms, using Set 2, plotted against the B-C-N compo- sition. 87 Figure 4.9: Poisson ratio ? against the B-C-N composition. 88 A graph illustrating the percentage number of bonds against the B-C- N trend is plotted as shown in figure 4.10. Here the horizontal axes are showing the percentage of each type of bond to the total number of bonds in the structure. Figure 4.10: Percentage number of bonds against the B-C-N composition Chapter 5 Discussion and Conclusion The B-C and B-C-N materials have been studied using Tersoff potentials applied in the GULP code. The calculated formation energies, elastic con- stants and moduli were compared with the available results for diamond, BN and BC2N . We found a good agreement of the calculated results with experiments[96] for the Bulk modulus and elastic constants but a slight dis- crepancy in the total energy . We found the cell constant for diamond to be 3.565 ?A which is very close to the experimental value of 3.567 ?A by Tiedje et al [95]. There is a 1.6 % difference between the experimental value of the bulk modulus of diamond by Mayer et al [96] and our calculated value of 425 GPa. Amongst the BC3 structures, the tetragonal structures show better stability, 89 90 as evident from lower formation energy than found for trigonal phase. To our knowledge there are no experimental results for the other B-C and B-C-N stuctures and hence our studied materials are novel. To gain a further insight into the results graphs are plotted as shown in figures 3.2-3.12, 4.3- 4.10. This is mainly done to show the dependence of various properties of these materials on the concentration of B, C or N. All our structures show large values for elastic constants, especially the C11 and C44 components. Large Moduli values, especially Young?s modulus and the bulk modulus which is a reciprocal of compressibility. The bulk modulus can be used as an indication of the materials? strength, chemical bonding and electronic structure[102]. This is somewhat an indication of an ultra- incompressible and superhard material. This work has also indicated that similar trends in the behaviour of the bulk, shear and elastic moduli of diamond-like B-C and B-C-N structures are deduced from the Tersoff po- tential calculations. It is of importance that the size of the chosen supercell for the random- ization process does not affect the final result or the choice of numerical 91 randomization. This, is ascertained by varying the supercell size to assure convergence. Employing structural optimisation ensures better quality of the starting structure[10]. We note the convergence of the moduli is quite rapid as illustrated graphically by figure 3.10. Amongst the B-C-N structures we only have considered an 8 atoms supercell. This is done because there is not much change in the structure even if it was bigger, since the atoms are fixed in their lattice positions. As expected , randomization does not have any effect on the results for diamond, since all lattice points contain identical carbon atoms. This is clearly seen in figure 3.11, with energy differences very close to zero and in figure 3.12, where the ratio of cell constants is very close to 1. It is of interest that the results, especially BC2N results are closer to ab initio results by Sun et al [58]. The calculated value for bulk modulus for BC7 of 377 ? 382 GPa is close to the value of 387 GPa obtained by first principles calculation recently done by Liang et al [103]. Two sets of nitrogen Tersoff parameters have reported in literature by Matsunaga et al [72](Set 1) and Justo et al [73](Set 2), and we have used 92 these parameter sets for simulation. Very good results were obtained using nitrogen parameter Set 2, whereas Set 1 parameters for nitrogen yielded poorer results, as they seem to be unphysical. Thereafter, B-C-N properties were only studied using Set 2 parameters. There is a good agreement for our calculated results for BN and BC2N with experimental results found in the literature [55, 101, 18]. Amongst the BC3 structures, the tetragonal structures are observed to have a lower formation energy than the trigonal structures. In conclusion we have successfully used the empirical Tersoff interatomic potentials, applied in the Gulp code, to simulate 13 different diamond-like structures and have determined a number of structural properties. The Gulp molecular dynamics code yielded very good results for B-C and B-C-N ma- terials. The fitted Tersoff interatomic potentials give very good results for diamond in comparison with experiment. We have studied the mechanical properties and formation energies of B-C and B-C-N structures. The values obtained for the B-C and B-C-N structures confirm that they are superhard materials with hardness comparable, though slightly lower, to that of diamond. They are promising hard materials with 93 very interesting properties which may even work when diamond fails, since they have a different chemical and electrical properties from diamond. The Tersoff parameters for nitrogen reported by Matsunaga et al [72] do not give desirable results in this work, in the sense that the properties obtained are unphysical. Hence we can conclude that they are not suitable for B-C-N materials. Since use of the parameter Set 2 nitrogen reported by Justo et al [73] yielded very good results, they were acceptable for B-C-N modelling materials. Experimental values for the lattice constant(3.635? 0.006?A) and bulk modulus(335?8GPa) of BC5[23] are closer to our calculated values for trigonal BC3. The bulk modulus of the B-C materials increases as then number of boron- carbon bonds decreases and the number of carbon-carbon bonds increases. This is expected since there is a weaker bonding between B and C bonds in comparison to the C-C bonding. The experimental values for lattice pa- rameter and bulk modulus of BC5[16] are closer to our calculated values for trigonal BC3. Relative to diamond and boron, it is expected that these structures should be metastable, as seen in previous results[18]. In addition the B-C structures are lighter than diamond and may shape more easily than 94 diamond. With regard to Poisson ratio, elastic constants and moduli, BC6N is second in hardness after diamond. With the exceptional hardness, implied by the high elastic constants, es- pecially C44, and the high moduli values, the B-C and B-C-N materials are attractive alternatives for advanced abrasives and high-temperature electron- ics. Since boron is a metalloid which implies intermediate conductivity, we expect the B-C materials to make good semiconductors. Further since boron has a lower atomic mass than Carbon, hence we anticipate the B-C materials to be lighter than diamond. In addition the atomic radius of N is smaller than that of C, hence substitutional of N atoms for C atoms is expected to result in significant strain in diamond crystal. The pair of nitrogen and boron should be the perfect co-doping candidates for making a n-type diamond-like material. This work could be furthered by looking into the electronic structure of such materials, using ab initio molecular dynamics codes. The electronic structure is of interest since it may reveal some further important properties of these novel materials, such as their use for jewellery which is of great importance in life. Bibliography [1] K.Cornelis H.Cornelius. Manual of Mineralogy. Wiley, 20th edition, 1985. [2] T. Forester W. Smith and I. Todorov. The DL POLY 2 Uer Manual. STFC Daresbury Laboratory, 2008. [3] R. Garg I. Drouzas J. Ulloa P. Koenraad M. Steer H. Liu M. Hopkin- son D. Mowbray V. Haxha, M. Migliorato. The use of abel-tersoff potentials in atomistic simulations of InGaAsSb/GaAs. Online at http://www.nusod.org/conf08/ThA4.pdf, September 2008. [4] C. Kittel. Introduction to Solid State Physics. Wiley, 6th edition, 1996. [5] Online at http://www.gordonengland.co.uk. 95 96 [6] R.L. Smith G.E. Sandland, editor. An Accurate Method of Determining Hardness of Metals, With Particular Reference to Those of High Degree of Hardness,. Institute of Mechanical Engineers, 1922. [7] Hardness testing -the test methodology guide, part-2. online at www.hardnesstesters.com/hardness-method-2.htm, 2008. [8] M. Allen. Comp.Soft Matter, 23, 2004. [9] M. L.Cohen. Phys.Rev.B, 32:12, 1985. [10] J. Gale. General Utility Lattice Program 3.1. Nanochemistry Research Institute, Department of Applied Chemistry, Curtin University of Tech- nology,. [11] F. Birch. Phys.Rev., 71:908, 1947. [12] M. Hebbache. Solid state Comm., 113:427?432, 2000. [13] Online at http://www.engineeringtoolbox.com. [14] A. Schmidt O. Boresi and R. Sidebottom. Advanced Mechanics of Materials,. Wiley, 1993. [15] E.Ziambaras and E. Schro?der. Phys.Rev. B, 68:064118, 2003. 97 [16] V. Solozhenko, O.Kurakevych, D. Andrault, Y. Le Godec and M. Mezouar. Phys.Rev.Lett., 102:015506, 2009. [17] G. T. Gao, K. Van Workum, J. D. Schall and J. A. Harrison. J.Phys.Cond.Matter, 18:S1737?S1750, 2006. [18] R.Hemley V.Brazhkin, A.Lyapin. Philosophical Magazine A, 82(2):231?253, 2002. [19] H-Y. Chung. Science, 316:436, 2007. [20] J. E. Lowther. J.Phys.Cond.Matter, 17:3221?3229, 2005. [21] R. Crenko. Phys.Rev.B., 7(10), 1973. [22] J. Custers. Nature, 176:173?174, 1955. [23] V.Solozhenko, N.Dubrovinskaia, L.Dubrovinsky. Appl.Phys.Letters, 85:1508, 2004. [24] S. Louie M. Cohen S. Han, J. Ihm. Phys. Rev. Lett., 80:995, 1998. [25] A. Collins A. Williams, E. Lightowlers. J. Phys. Condens. matter, 3:1727, 1970. [26] A. Collins A. Williams. J. Phys. Condens. Matter, 4:1789, 1971. 98 [27] R. Chrenko B. Massarani, J. Bourgoin. Phys. Rev. B, 17:1758, 1978. [28] R. Mamin and T. Inushima. Phys. Rev. B, 63:033201, 2001. [29] I. Snook A. Barnard, S. Russo. J, Chem. Phys., 118:10725, 2003. [30] W. Wilson. Phys.Rev., 127:1549, 1962. [31] H.Sagai K.Haruna H.Maeta N.Matsumoto T.Sato, K.Ohashi and H.Otsuka. Phys.Rev.B, 61:12970, 2001. [32] E. Bauer N. Mel?nik N. Curro J. Thompson S. Stishov E. Ekimov, V. Sidorov. Nature, 428:542, 2001. [33] S. Sque R. Jones J.Goss, P. Briddon. Phys. Rev. B, 69:165215, 2004. [34] A. Voronov and A. V. Rakhmanina. Inorg. Mater, 29:707, 1993. [35] V. Torishni A. Vishnevski, A. Gontar and A. A. Shulzhenko. Sov. Phys. Semicond., 15:659, 1981. [36] A. Deneuville F. Fontaine L. Abello P. Gonon, E. Gheeraert and G. Lucazeau. J.Appl.Phys., 78:7059, 1995. [37] P. Germi M. Pernet F. Brunet, A. Deneuville and E. Gheeraert. J. appl.Phys., 81:1120, 1997. 99 [38] H. U. Baerwind E. Obermeier L. Haase W. Seifert A. Ringhandt C. Johnston S. Romani H. Bishop M. Werner, O. Dorsch and P. Chalker. Appl. Phys. Lett., 64:595, 1994. [39] E. Gheeraert E. Bustarret and K. Watanabe. Phys. Stat. Solidi A, 199:9, 2003. [40] E. Bauer N. Melnik N. Curro J. Thompson E. Ekimov, V. Sidorov and S. Stishov. Nature (London), 428:542, 2004. [41] Yu. A. Timofeev A. Utyuzh and A. Rakhmanina. Inorg. Mater., 40:926, 2004. [42] M. McCluskey M. Plano J. Ager III, W. Walukiewicz and M. Land- strass. Appl. Phys. Lett, 66:616, 1995. [43] C. Marcenat E. Gheeraert C. Cytermann J. Marcus E. Bustarret, J. Kae`mare`ik and T. Klein. Phys. Rev. Lett., 93:237005, 2004. [44] L. Cohen A. Liu, M. Wentzcovitch. Phys.Rev. B, 39:3, 1989. [45] D. Roundy M. Cohen S.Louie H. Sun, S. Jhi. Phys. Rev. B, 64, 2001. 100 [46] A.Lyapin S. Popova A.Rakhmanina S. Stishov V. Lebedev Y. Katayama V.Brazhkin, E.Ekimov and K. Kato. Phys.Rev.B, 74(R):140502, 2006. [47] S.Louie M.Cohen S.Han, J.Ihm. Phys.Rev.Lett., 80:5, 1998. [48] R. Farrer. Solid State Comm., 7:685, 1969. [49] A. Badzian. Mater. Res .Bull, 16:1385, 1981. [50] M. Hubacek and T. Sato. J. Solid State Chem., 114:258, 1995. [51] V.Turkevich V. Solozhenko and T. Sato. J. Am. Ceram. Soc., 80:3229, 1997. [52] Y. Zhao. J. Mater. Res, 17:3139, 2002. [53] B. Xu Q. Wu Q. Hu Z. Liu J. He D. Yu Y. Tian H. Wang X. Luo, X. Guo. Phys.Rev. B, 76(094103), 2007. [54] L. Daemen T.Shen R. Schwarz Y. Zhu D.Bish J. Huang J. Zhang G. Shen J. Qian Y. Zhao, D. He and T. Zerda. J. Mater. Res, 17:3139, 2002. [55] G. Fiquet M. Mezouar V. Solozhenko, D. Andrault and D. Rubie. Appl. Phys. Lett., 78:1385, 2001. 101 [56] P. Zinin M. Manghnani S. Tkachev, V. Solozhenko and L. Ming. Phys.Rev.B, 68:052104, 2003. [57] M. Gubachek N. Borimchuk V. Zelyavskii N. Ostrovskaya V. Yarosh .V. Kurdyumov, V. Solozhenko. Powder Metall. Met. Ceram., 93:629? 637, 2000. [58] Y-X. Fan J. Chen H-T. Wang X. Guo J. He J. Sun, X-F. Zhou and Y. Tian. Phys.Rev.B., 73:045108, 2006. [59] S. Azevado and R.De Paira. Euro.Phys.Lett., 75:126, 2006. [60] J. Tersoff. Phys.Rev.lett., 61(25), 1988. [61] J.Tersoff. Phys.Rev. B, 39(8), 1989. [62] A. Hinchliffe. Molecular Modelling for Beginners. Wiley, 2003. [63] J. Ogivile. Dynamic Aspects of Intermolecular Interactions, chapter 5, pages 49?56. Plenum Press, New York, 1998. [64] M.Z. Bazant. Interatomic Forces in Covalent Solids. PhD thesis, De- partment of physics, Havard university , Cambridge Massachusetts, July 1997. 102 [65] J. Gale. J. Chem. Soc., Faraday Trans., 1:629, 1997. [66] M. Gillan W. Smith, G.N. Greaves. J. Chem. Phys., 3091, 1995. [67] D. Saada. Comparison between Tersoff Potential and ab initio results for Surface Graphitazation of Diamond. PhD thesis, Technion, Israel Institute of Technology, Department of Physics, 2000. [68] J.Tersoff. Phys.Rev.B, 37(12), 1988. [69] Furio Ercolessi. A molecular dynamics primer. International School for Advanced Studies, Trieste, Italy, 1997. [70] J. Justo F. de Brito Mota and A. Fazzio. Phys.Rev. B, 58(13), 1998. [71] D. Brenner. Phys.Rev.B, 42(15), 1990. [72] K. Matsunaga, C. Fisher and H. Matsara. Japan J. Appl. Phys., 39:L41, 2000. [73] A. Fazzio F. de Brito Mota, J.F. Justo. J.Int. Quantum Chem., 70:973, 1998. [74] Ho J. Hwang W. Ha Moon. Phys.Lett. A, 315:319?324, 2003. 103 [75] D. Rachedb K. Amaraa, B. Soudinib and A. Boudalia. Comp.Mater.Sci., 44(2):635?640, 2008. [76] P. Erhart J. Nord, K. Albe and K. Nordlund. J. Phys.Cond. Matter, 15:5649?5662, 2003. [77] J. Haile. Molecular Dynamics Simulation, Elementary Methods. Wiley Interscience, 1992. [78] J.R. Beeler Jr. and G.L. Kulcinski. Agenda discussion:computer tech- niques. In J.R. Beeler Jr. P.C. Gehlen and R.I. Jaffe, editors, Inter- atomic Potentials and Simulation of Lattice Defects, pages 735?751. Plenum Press, 1972. [79] C. Lubich E. Hairer and Gerhard Wanner. Acta Numerica., 12:399?450, 2003. [80] B. J. Alder and T. E. Wainwright. J. Chem. Phys, 31(459), 1959. [81] W. Vetterling and B. P. Flannery W. Press, S. Teukolsky. Numerical Recipes. Cambridge University Press,, 2nd edition, 1992. 104 [82] Mark E. Tuckerman and Bruce J. Berne. Molecular dynamics algorithm for multiple time scales: Systems with long range forces. J.Chem Phys., 94(10), 1991. [83] M.Born and R.Oppenheimer. Ann. Physik, 84(457), 1927. [84] J. Furthmuller G.Kresse, M. Marsman. VASP the GUIDE, 2009. [85] Online at http://www. castep.org. [86] Online at http://www.univie.ac.at/newtonx. [87] C. Bayly I. Gould K. Merz D. Ferguson D. Spellmeyer T. Fox J. Cald- well W. Cornell, P. Cieplak and P. Kollman. J. Am. Chem. Soc., 117:5179?5197, 1995. [88] A. Gursoy A. Dalke L. Kale R. Skeel K. Schulten M. Nelson, W. Humphrey. Namd- a parrallel, object -oriented molecular dynamics program. 1996. [89] B.Olafson D. States S. Swaminathan B. Brooks, R. Bruccoleri and M. Karplus. J. Comp. Chem, 4:187?217, 1983. [90] D. Burnett. Finite Element Analysis: From Concepts To Applications. Addison-Wesley, 1987. 105 [91] L. Verlet. Phys. Rev, 159(98), 1967. [92] A. Rahman. Phys. Rev, 136(405), 1964. [93] B. Montgomery Pettitt Paul E. Smith. Efficient ewald electrostatic calculations for large systems. Comp. Phys. Comm., 91:339?344, 1995. [94] D.Chung and W.Buessem. J. Appl.Phys., 39(2777), 1968. [95] T. Tiedje K. Myrtle M. Kasrai J.R. Dahn, B.M. Way. Phys.Rev.B, 46:1697, 1992. [96] L. Mayer Y. Baskin. Phys.Rev., 100:44, 1955. [97] D. Roundy M. Cohen S. Louie S. Jhi, H. Sun. Phys.Rev.B, 64:755, 2008. [98] A. Pollan E. Zouboulis, M. Grimsdditch. J. Appl.Phys., 76(882), 1994. [99] S.Kim L.Kleinman D.Bylander, S.Lee. Phys.Rev.B, 45:3248, 1992. [100] D. Andrault Y. Le Godec V. Solozhenko, O. Kurakevych and M. Mezouar. Phys.Rev.Lett., 102:015506, 2009. [101] R. Jeanloz E. Knittle, R. Kanner and M.L. Cohen. Phys.Rev. B, 51(18), 1995. 106 [102] C. Yoo H. Cynn, J. Klepeis and D. Young. Phys.Rev. Lett., 88(13), 2002. [103] J. Zhao Y. Liang, W. Zhang and L. Chen. Phys.Rev.B, 80(113401), 2009.