. ' . \I Q Ii‘r 'DQ fig; , '53 I‘m" ‘5 filiév; '8‘ : ‘L ”f '1'. . .‘ 5W ' 'I' ' '5‘” 5M?" :33 "I’LL, j'_‘:?.%..’1;I I a: E: «I 94} [ML \ ‘1‘; ' -' ' ' 5 "II , .. | I 'f‘Itfi‘M.‘ 4 I‘M-l". it .I 53. 5?an I "MA-#1.? .l“ ' W'l' || 37‘ .. ‘C". . 9: will“. I ' h. '. . , ‘ 5.,rtiqw‘;;::11hr4b‘ ;Tr:‘b “ ‘ . . I l 2“?" .;::'J.IF “' I :1}- "7:.:: . .. ‘ . .. “Ia-:5»; ."43‘5', . '39 135“ «$1,353 3 . :‘ :JI‘;» 3., " 3‘37 - i \ 7.1.1! ”'1‘ - r: :5 as; ' ‘I hi-V'II'31J 5'11?” ‘ )I‘ I , I I - r1 . WI!“ II‘I :3 ' r4 ‘ W ‘— - :.-_. 'fr,‘ I II '"I‘ I 'I'I‘ bliw MM Alf“: IMIIWI .I'i II” " . .’II(}|':IUI,I ..-M1]","'I’I‘\I£ll 'IN‘I‘ -I.- - ",' l“.:’("‘lll: ': "i. ‘ I. ... III'I' Um :III WI I41“ ' II‘I'Hf'” ".131; I'ng“'l-I :.I\:‘:’II‘L‘IK‘I5: 0": '"WP‘J', }‘ 5‘5"" ~"I,'I[ H¥3$ g” )1," .M ‘3' ' I.|I I [[61] “:19“ A‘F'Q ,. I. ”HM Mi'l '3'. ”mph “i" ‘prfi v I) O D: l I I v I "‘1' ‘9'. U I'J’Il‘ld WM. - I‘M!“ H . 'IIMI Ij5l'I ' II N «I'm - ' 1 I I I HI “ gm M. I‘QE"I.I“‘ )‘:\..:.I MM’I‘I:\;I(|I""}I:?:! 'OI‘:‘ ”h: M 5.35:5; ALI‘f I’ :- .' ;.. ., ‘.:I W. 55.3}. ,. ;‘ W. ,9.)me -. .3. - - . . _ _. M .._.,,'. ' .I II .J;' M ‘Il'l _l ‘I. |"-' | ‘v' I I‘. . ' I .5. I ' P) ', " .73“; ’ I' V 1‘ " ‘J 43:; I III ’— I I9 I“. - M. (M. "ii' ‘.‘;.. t," 'l""'I1I9"-"’;I" 1.4-" jgfiw-Xg’twfizzqrq: 4mg...f}’\t.7:§{§;:cfigz .: ‘I.“ w'f‘v‘ “I.‘l‘ 1"”9'Iy1‘1" I:“" 5".“ I‘r ['II.:h'I|.'I' ' I' :.%:MI I 'IY?)&W§€:.. \‘h, ‘ 1% ' ‘ .5 ”9!..‘H I"; .‘.".| {"315}. U” I“. ' I I I .2 I, (.I._ I . .‘t‘ .I,..'! "Q05. . .irfi‘ . . I I- Ill'l . ”a...“ - ‘- - "J. . . v . ‘ 4' ‘I. I 9 ”'57 ’ I ,1 'Il ‘ :5 M55" ‘ ‘ 'I 'I‘ filmy.” k' " 6‘3! {a ". l I 6"!" "II '.' ' '-. 'II'- 1”"‘5 ‘ '-' "JI' ° ' .k’ 'Pfi"; ’. a» 55:5“ . . I , N ‘ I l "I 1. I,.' ff. ‘ . -. I I} “‘J):€‘V):r . . "} u . ~ :1‘ 3°: :‘, ’- If»; * 1"‘. .q‘; . .‘ 531:! .2...“ I M‘l ,_I _l .. _ “I If?!“ 2:; 'I. .‘f'a III '7.I‘- ‘ I. .. I.».~-‘;,lewn \ - .' .\ ""II‘JC.‘ ,‘j-'~.I5{I:"Mi «5": ‘ ' 'l,‘;'rnl '5 Hit}. : . (I III: I 1'.'.:.!.|'.""‘,|‘3' .5 ..' "I”:"I'I'I . [I l; “h “‘I.-l" . c I. "I‘ In" I :1.""'1f‘ 'HI‘ ' . I;-'.;I "1| 4| | Kim-MILL. . .' . .3 ', M I I ' . .' .M _ I . M. . , It . M .7 , . . 4 .. , x .. M If‘lle‘il I}; . I I v 'I i’.‘ I‘.' . . [ v.9 ‘ ‘ . .M 9, '; , M I' -< I ' , ' .‘I I- 'l_ 1". M I ,|, M .I ‘.I. II " , . ‘II.' $45!; II' Q .;~I;m1;.,~-I..«\ HWII‘I'MA """MIIIIEI—{III' XW' " THE“ ”BRARY Michigan State This is to certify that the dissertation entitled award“ ded “:1ku Since Matt/awe; Mal 514%: Freya/filer 0; 4AM; presented by Hohgxzima HL has been accepted towards fulfillment of the requirements for Mdegteein ER‘ZECJ /{/7be/l£1 Majot [”0me Date 2’2 {' l’g {6785‘ “till...- AM— u... 4 ‘ I" In , I x , 0.12771 MSU LIBRARIES “ RETURNING MATERIALS: Place in book drop to remove this checkout from your record. ‘FINES will be charged if book is returned after the date stamped below. COMPUTER GENERATED VITREOUS SILICA NETWORKS AND ELASTIC PROPERTIES OF GLASSES By HONGXING HE A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Physics and Astronomy 1985 38%805 ABSTRACT COMPUTER GENERATED VITREOUS SILICA NETWORKS AND ELASTIC PROPERTIES OF GLASSES By HONGXING HE A large fully bonded periodic continuous random network model of vitreous silica is generated via a computer simulation. This is done by introducing the defects into a diamond lattice and then decorating with oxygen atoms. The model shows a good agreement in Radial Distribution Function with experiment. Chapter 1 describes the process of generation of the model and compares this work with some other works. The characteristics of the model are also described in this Chapter. The elastic properties of glasses of different mean coordination are quite different. There are two classes of covalent glasses, amorphous solids and polymeric glasses. When the mean coordination number decreases, a transition takes place from the former to the latter. Transition point is close to (r) = 2.”. This is in good agreement with the Constraint Counting Argument. Chapter 2 presents the evidence tOlthiS and describes the behavior of two physical quantities (Fraction of Zero Frequency Modes and the Elastic Moduli) when the mean coordination decreases from four towards the transition point. The implication of the introduction of the dihedral angle forces is also discussed in this Chapter. ACKNOWLEDGEMENTS First of all I would like to thank my advisor Professor M. F. Thorpe for all his patient guidance over last four years. His knowledge and insight of solid state physics helped guide me through this work. I would like to express my appreciation to Professor T. Kaplan for his well taught statistical machannic courses and many useful discussions, to Professor H. McManus for his interesting nuclear physics lectures, to Professor S. Solin for his very intesting and useful third term solid state physics course, and to all three above for serving on my guidance committee. I am indebted to Dr. Simon Deleeuw, Dr. E. Garboczi and Dr. M. Thomsen for their numerous usefull discussions and help over last few years. Dr. D. Heys, Dr. D. Sahu, S. Tang, W. Tang and w. Xia have been my good friend and the source of much help throughout my graduate career. Many thanks to my dear wife Yuqiu and son Zhen, their coming to the United States to Join and help me is a great encouragement to me. Yuqiu's patience and courage to take all the house duty enable me to concentrate to my work and finish this thesis on time. IT List of List of Chapter Section Section Section Section Section Section Chapter Section Section Section Section Section Section TABLE OF CONTENTS Tables Figures 1 Computer Generated Models of Vitreous Silica Networks 1.1 Introduction 1.2 Generation of a-Si Network by the Introduction of Defects into a Diamond Lattice 1.3 Bragg Diffraction Pattern 1.“ Construction of the Vitreous Silica Network 1.5 Radial Distribution Function 1.6 Conclusions 2 Elastic Properties of Classes 2.1 Introduction 2.2 Constaint Counting and Zero Frequency Mode 2.3 Model Description 2.“ Fraction of Zero Frequency Modes 2.5 Elastic Moduli 2.6 Mean Field Behavior VI 12 20 314 "3 us 50 50 62 68 73 78 91 Section 2.? Effective a and 8 section 2.8 Flow Diagram Section 2.9 Introduction of Dihedral Angle Forces Section 2.10 Conclusions Appendix A Relaxation Process Appendix B Extrapolation of Elastic Energy Relaxation List of References 99 105 116 12:4 125 129 131 1.1 1.2 1.3 1.“ 1.5 1.6 2.2 List of Tables The rms deviation of bond length and bond angle and the height of (111) peak in the diffraction pattern vs number of defects introduced. The rms deviation of bond length and bond angle and height of (111) peak in diffraction pattern before and after geometric relaxation. The ring statistics of an a-Si network which has no memory of crystal structure. The ring statistics of an a—Si network after geometric relaxation. Average value and rms deviation of bond length and angle of our a-SiO2 model. Comparison of our model and some other models and experimental value of vitreous silica. Exponent and best power law fit for various value of B/a. Values of elastic modulus before and after extrapolations. 26 28 29 31 38 no 96 130 Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure 1.1 1.2 1.3 107 1.9 1.10 1.11 1.12 List of Figures Comparison of the RDF derived from the Connell and Temkin model and the 519-atom model of Polk and Boudreaux with experimental data for a-Ge. Comparison of RDF calculated from Bell and Dean model and that from neutron scattering experiment. Comparison of the RDF calculated from the Ching's model (averaged over three models) and that from X-ray experiment. 3D pair defect. 2D Pair defect. Four-fold ring is created upon introduction of defect on a side of a five-fold ring. Concave ring gives rise to large angle distortions. A 2D network generated by introducing 250 defects into a 800-site honeycomb lattice. Defect observed by F. H. Sillinger and T. A. Weber in molecular dynamic simulations. Bragg diffraction pattern of the 2D network depicted in Figure 1.8. Bragg diffraction pattern of a 3D network, N-512, 250 defects have been introduced. Bragg diffraction pattern of a generated 3D a-Si random network with no crystal memory left. VI 10 13 1M 15 17 19 2” 25 3O Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure 2.1 2.2 2.3 2.5 2.6 2.7 (a) (b) (a) (b) (a) (b) Bragg diffraction pattern of a 3D random network after geometrical relaxation. A schematic view of the local structure of vitreous silica. Nearest 0-0 distance vs related bond lengths and angles. network. 2 31-0 bond length, O-Si-O angle, and 0-Si-O angle Bragg diffraction pattern of an a-SiO distributions of our a-SiO2 network. Comparison of T(r) calculated from our a-SiO2 model and that from neutron scattering experiment. Partial correlation functions of our a-SiO2 network. Schematic bonding topology of Ge As Se X Y Hybrid bonding orbitals of Si in the solids. 1-x-y glasses. Polymer chain formed by Se atoms. Chains with cross links introduced by As and Ge atoms. Polymeric glass where rigid region does not percolate. Amorphous solid where rigid region percolates. Rigid units connected by a single free hinge. Rigid units connected by a single chain. Fraction of zero frequency modes vs mean coordination from constraint counting. Schematic view of the effect of the closed ring VII 33 35 36 H1 ”2 us In 51 52 53 55 57 66 Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure 2.8 2.9 (a) (b) 2.10 2,11 2.12 2.13 2.1“ 2.15 on the number of independent constraints. 67 Number of three-coordinated sites vs mean coordination for three samples, N-512. 72 Number of three-coordinated sites vs mean coordination of three samples, N=216. Fraction of zero fraction modes vs mean coordination of above three samples. 76 Fraction of zero frequency modes vs mean coordination, averaged over three samples. 77 Chain-like local structure which can be trimmed. 80 Singly connected cluster which can be trimmed. 81 Fraction of bond present after trimming of three samples. 82 C11 vs mean coordination of three samples. 86 C11, Cu”, and B vs mean coordination averaged over three samples. 87 C11, Cu”, and B vs mean coordination for B/a=0.02,0.1,0.5 and 1.0. 88 Comparison of Cu” from our calculations and that from sound velocity measurement of Se Gex 1-x glasses by J. Y. Duquesne and G. Belleesa. 89 Log - Log curve and power law fit of Cij’ B/a-0.1. 92 Log - Log curve and power law fit of Cij’ B/a-O.2. 93 Log - Log curve and power law fit of CiJ’ Elm-0.5. 9M Log - Log curve and power law fit of C filo-1.0. 95 ij’ A vs mean coordination, B/a=0.2,0.1,0.2,0.5,1.0. 100 a * a and B vs mean coordination and CAN fit, VIII Figure Figure Figure Figure Figure Figure Figure Figure Figure 2.2” 2.25 2.26 2.27 2.28 2.29 2.30 2.31 2.32 (a) (b) B/a-0.02, 0.1, 0.2, 0.5, 1.0. It It 8 la vs mean coordination. Diagram of square net with first and second neighbor springs. C / C 11 AA 2D Network P-0.695 before the trimming. and B / CHM flow diagram,3D. 2D Network P-0.695 after the trimming. B and can vs mean Coordination for 2D. C11 / CNN vs mean coordination for 20. C11 vs mean coordination with and without the dihedral angle forces. Cml vs mean coordination with and without the dihedral angle forces. B vs mean coordination with and without the dihedral angle forces. IX 103 10” 106 109 113 1111 115 119 120 121 Chapter 1 Computer Generated Models of Vitreous Silica Networks Section 1.1 Introduction Fundamental understanding of the properties of solids requires that we know the identities of the constituent atoms and their arrangement in space. In crystalline solids, due to the long-range periodicity, we can determine the microscopic arrangement accurately. In.a disordered solid, however, this long-range order does not exist and the structure is usually characterized by short and or intermediate-range order such as nearest neighbor coordinatnon numbers, average bond length and bond angle of bonded atoms and ring statistics etc. Zachariasen (Reference 1.1) proposed the continuous random network model for the glass structure more than fifty years ago. This continuous random network model now is generally believed to be the best description of covalently bonded amorphous solids. Unlike crystalline solids, there is no unique set of atomic coordinates to describe the structure of an amorphous material. There may be many configurations of the material which are degenerate or separated by only a small energy barrier. The experimentally measured Radial Distribution Function (RDF, usually is defined in terms of the atomic density p(r) at distance r from any chosen atom: J(r) = Anr2p(r)) can provide only partial information regarding their structure. Actually many different models can have similar radial distribution function. Model construction is appealing even though it supplies us with only one possible microscopic structure of an amorphous solid. The atomic coordinates so generated may be used to study other properties of amorphous solids via computer simulations, such as electronic, vibrational and optical properties etc. Much work has been done to construct continuous random network models of various amorphous solids, such as tetrahedrally bonded amorphous semiconductors and various amorphous compounds with constituent atoms of different coordination numbers. Building of the model of tetrahedrally bonded amorphous semiconductors has received much attention. Amorphous silicon (a-Si) is a particularly important and interesting material. It is a typical covalently bonded amorphous semiconductor and it is the most promising material for a wide range of application of solar energy conversion. The structure of tetrahedrally bonded amorphous semiconductors have the following main features: (1) Z-M, each atom is four-fold coordinated (Z is coordination number). (2) Narrow distribution of the bond lengths. (3) No dangling bonds. (A) Significant spread in bond angles. (5) No long-range order. It was a model built by D. E. Polk (Reference 1.2) that first demonstrated the possibility of building an extended Zn“ continuous random network with small distortions in bond lengths and bond angles. He built the model by hand, using plastic and metal units. The model contained ““0 atomic centers and was later extended by D. E. Polk and D. S. Boudreaux to contain 519 atomic centers. Paul Steinhardt and others (Reference 1.3) later relaxed the Polk model via a Keating's potential (Reference 1.“): 3 a 2 2 2 v=—-—— {(r- -d) 16 d 1’1 ii +3—L 2(3-3 +92)2 (111) 8 d2 2 1,3 21 ij 3 Here a and B are elastic force constant and d is the desired nearest neighbor distance. Through relaxation, they were able to improve the model by reducing the rms deviation in the bond length to only 0.8 1 and the rms angle deviation to 6.7°. They also built a new 201-atom model starting from a 21-atom seed and serially adding atoms. As each small group of atoms was added, their positions were relaxed by computer to minimize the forces on them. All the above models contain five, six, seven and higher membered rings. G. A. N. Connell and R. J. Temkin later built a 238 atom tetrhedrally bonded amorphous semiconductor model (Reference 1.5) using plastic and metal units with only even membered rings (6,8 . . . ). They did this on the basis of great structure similarity of a-Ge and III - V semiconductor (InAs for example) (Reference 1.6) and that since odd-membered rings would lead to the presence of energetically unfavorable bonds between like atoms. Their model has reasonable distortions with 1.2% and 10.7° in the rms bond length and angle deviations, respectively. There is a main disagreement between Connell- Temkin model and experiment that is the Connell-Temkin model has more features than that the experiment at higher r in the correlation function (Reference 1.7). This indicates that the restriction to even- membered rings imposes an longer range order which does not exist in the amorphous semiconductors. Beeman and Bobbs (Reference 1.8) on the other hand, noticed that there was some experimental evidence for bonds between like atoms in the amorphous III - V semiconductors. They produced a sequence of models containing an increasing number of odd-membered rings, using the Connell-Temkin model as a starting point. The technique was to randomly (r displace atoms in the structure by up to 0.22 r is nearest neighbor 1 distance) and then adjust neighbor table for each atom if appropriate. 1 The resulting model was relaxed to reduce the distortions. The model turned out to have better agreement with experiment in RDF. All the above hand built models had small deviations for bond lengths and bond angles and showed agreement with experimentally obtained radial distribution functions (Figure 1.1 displays some of the comparisons). However, they all had free boundaries. The good agreement with experiment in the radial distribution function is a necessary but not sufficient condition for a good model. Other properties of solids based on the constructed model are more discriminating. In this respect, models with periodic boundary conditions are more valuable because they represent an infinite array of network atoms which are free from surface effect. 50'— s 6 I 1 J ( r ) (Atoms/A) 21g Figure 1.1 RDF, dotted line is from even-membered ring model of Connell and Temkin, dashed line is from 519 atom model of Polk and Boudreaux, heavy line is from experiment (Reference 1.9). A few periodic models of tetrahedrally bonded amorphous semiconductors have been built. The first such models was built by D. Henderson and F. Herman (Reference 1.10). Their first model with 6“ atoms was built via a computer simulation. Their method was as follows. They first placed N atoms in a diamond structure which forms a unit cell. Then they disrupted this crystalline structure by moving the atoms at random. Each atom was considered in turn and moved 10 percent of the distance to the centroid of its four nearest neighbors. Then the four nearest neighbors and the twelve next nearest neighbors were moved radially ten percent of the distance to the nearest neighbor and next nearest neighbor distances of the crystal. Finally each atom was reconnected with the four closest atoms, thus generating a random network. This 6“ model had too large distortions, then Henderson and Herman also hand built a continuous random network with periodic boundary conditions containing 61 atoms, with a rms bond angle deviation of 12° which is slightly higher than the value experimentally suggested (10°). Some other periodic networks of similar size followed (Reference 1.11). L. Guttman generated a periodic random network with a large unit (512 atoms) using a computer algorithm. His method was as follows. He started with a part of crystalline lattice with n)“ (FCC lattice for example, n is number of nearest neighbors). Then he selected four out of n nearest neighbors to bond to each atom at random. Success of such proccess to match all the bonds and periodic boundaries was about one every 1000 attempts. Finally the network was relaxed. The model generated this way turned out to be poor, with a density 35 percent higher than the density of the crystal and the angle deviation was too large, being about 19-20 degrees (Reference 1.12). Model construction for vitreous silica based on the continuous random network model began as early as the mid-19603 (References 1.13 and 1.1“). The structure of vitreous silica are believed by experiment to have the following main features (Reference 1.15): (1) Z(Si)-“ and Z(0)=2, four-fold and two-fold coordinations for silicon and oxygen, respectively, with complete chemical ordering. (2) Almost constant Si-O bond lengths and O-Si-O angles. (3) No dangling bond in the bulk. (“) There is no long range order. (5) The Si-O-Si angles have an average angle of about 150° with a 15-20° in deviation. R. Bell and P. Dean built the first continuous random network model of vitreous silica( Reference 1.1“). The radial distribution function was in good agreement with experiment (Figure 1.2). They built the model by hand with polystyrene spheres and steel wires.‘Their model had free boundaries and contained 188 tetrahedral units. In terms of atoms, this comprised 188 four-coordinated silicon atoms, 326 doubly coordinated or bridging oxygen atoms and 100 singly coordinated non- bridging oxygen atoms at the boundary of the model, giving 61“ atoms in total. This model was characterized by a large number of four-membered rings. We define an n-fold ring as a closed path of bonds in the network which passes once, and only once, through n different silicon atoms. There were “8 four-fold rings, 55 five-fold rings, 9“ six-fold rings and 3 165 seven fold rings. The density of the model is 1.99 g / cm which is 3 lower than the experimental value of 2.20 g / cm . The deviation of 31-0 bond lengths and O-Si-O angles were 3.5% and 6.2°, respectively. I I l T I (d) Computed ' ‘ l I I I r (b) Experimental r I I I I I fl ° 02 04 08 08 rCnm) —.- Figure 1.2 RDF for neutron scattering from vitreous silica: (a) Computed for the Bell and Dean model. (b) Experimental curve. (Reference 1.9) w. T. Ching (Reference 1.16) generated a few periodic a-SiO2 networks by decorating the existing periodic Si networks with oxygen atoms. He breaks each Si-Si bond and inserts an oxygen in between to bond to these two silicon atoms. Finally he relaxed the network with a Keating like potential: 8 Si 3 a -> ~> 2 2 3 1 + + 2 2 v =————— Z(r.-r -d ) +--— 2 (r -r -dcos¢) t 16 d i 21 21 8 (12 1,3 21 23 v 0 - ——-3 -——°‘ 2 (F o? - c12 )2 + -3- :2 2 (F .; - d2COSG)2 n 16 d 1 hi hi 8 (12 i J ni nj Where 521 is the radius vector from atom i to its nearest neighbor atom 1, ¢ is the ideal 0-Si-0 angle, and 0 is the average Si-O-Si angle (was taken to be 1“7°). Ching used Blla - 0.17, 82/a - 1/3. In this manner he produced three periodic random networks with unit cells containing 162, 183, and 186 atoms respectively. These models have correct experimental density,(L7 to 1.7 percent deviation in 81-0 bond lengths, “.8-8.2 degrees deviation in O-Si-O bond angles. The Si-O-Si bond angles have 13-1“ degrees variation. The radial distribution functions calculated from the models agree with the experiment (Figure 1.3). Some defects of these models are that they are small in size and contain only five and higher membered rings since they were generated from existing a-Si networks which do not contain four-fold rings. However, some analysis of the experimenal results do suggest the existence of the small rings in vitreous silica (Reference 1.17). 10 RDF Figure 1.3 RDF of vitreous silica. Solid line is averaged RDF of the three periodic networks generated by W. Y. Ching. Dashed line is experimental X-ray RDF. There is a more recent report (Reference 1.18) of computer generated large a-SiO2 network containing 10,000 atoms, but the network has free boundaries and a large fraction of Si dangling bonds, ring 11 statistics are unrealistic, and three-fold rings dominate providing 80 percent of the total ring population. Much work has been done to generate more or less realistic models in order to understand the structure of amorphous solids from the microscopic level and to produce objects for studying various properties of amorphous solids. What we tried to do is to generate a good model of vitreous silica via a computer simulation. Since previous models generated to simulate the structure of vitreous silica were not quite satisfactory, being either too small or having free boundaries. Our ainl was to generate a model which was fully bonded, with periodic boundary conditions and a large unit cell (at least a few hundred atoms). The model must have the main structural features of vitreous silica so far as we know them from the experiment, the correct radial distribution function for example. The next sections describe how this was done. 12 Section 1.2 Generation of a-Si Network by the Introduction of Defects into a Diamond Lattice There are generally two ways of generating continuous random networks. (hue is generating by hand using plastic and metal units, for example. Second is generating random network via computer simulations. 'There are still two methods in computer generating random networks. One method is to let computer substitute the functions of our hands and brain to simulate the hand-building process. Another method is to disrupt the crystalline system to obtain random networks. For this method, many previous works were imposing random move of the atoms and rearranging the connectivity table and finally relaxing the system. There is another way to create the random network via a computer simulation by disordering a crystal. 0. Weaire (Reference 1.19) has suggested a defect structure which can be introduced into a crystalline stucture to disorder the crystal and eventually generate a model of an amorphous solids. Figure 1.“ depicts such a defect which will be introduced into the diamond lattice to generate model of a-Si. Figure 1.“ (a) is a part of a diamond lattice. The defect is introduced by moving atoms 1 and 2 to new positions and breaking the bonds between 1 and 3 and 2 and 8. We then connect 1 and 8, and 2 and 3 with new bonds. The new structure is shown in the Figure 1.“ (b). As a consequence of this defect, four five membered rings are created, twelve sixfold rings are removed (converted into five-and seven-fold rings), sixteen seven-fold rings are created 13 and the bond lengths and bond angles involved are no longer equal to the values of the diamond lattice. (a) . (b) Figure 1.“ 30 pair defect. (a) is a part of diamond lattice. (b) is the structure after the introduction of the defect. After a large number of such defects have been introduced into the diamond lattice, the system will hopefully lose its crystalline memory and become totally disordered. There will be rings of all different sizes and variations in the bond lengths and the bond angles. This random network may then be a good-candidate to model the structure of tetrahedrally bonded amorphous solids, like a-Si or a-Ge. 111 Il.is.easier to visualize this whole process in two dimension. Similar defects can be introduced into the 2D honeycomb lattice to generate a 2D continuous random network. (a) (b) Figure 1.5 20 pair defect.(a) is a part of honeycomb lattice. (b) is the structure after the introduction of the defect. Figure 1.5 shows such a defect construction in the 2D honeycomb lattice. By changing the positions of atoms 1 and 2, breaking the bond between 1 and “ and 2 and 5, making the new bonds to connect 1 and S and 2 and “, the local structure in Figure 1.5(a) is then transformed into 15 Figure 1.5(b). Consequently, two five-fold rings and two seven-fold rings are introduced into the honeycomb lattice and four six~fold rings are removed. Bond length and bond angle deviation from their crystalline values are also introduced into the network. Figure 1.6 Four-fold ring is created upon introduction of a defect on a side of a five-fold ring. A large number of such defects need to be introduced to disorder the crystal to create a random network. When a small number of defects are introduced at random, only five and seven membered rings can be created. But when fairly large number of defects are introduced, they 16 are run longer separated, rings of different size can be introduced. For example, the four-fold ring can be created when a defect is introduced on a side of five-fold ring (Figure 1.6). Figure 1.7 Concave ring is not allowed, since it gives rise to a large angle distortion. We have generated the 20 random networks by introducing a large number of such defects into the honeycomb lattice. The process is as follows: choose one atom at random, then choose one of its three nearest neighbors at random and introduce the defect using this pair (1. e. let these two atoms play the role of 1 and 2 in the Figure 1.5) under the following restrictions: (1) No four-membered or smaller rings can be introduced, because this would cause a large bond angle distortions. (2) All rings have to be convex. Concave ring would give large angle distortions in the local structure (Figure 1.7). 17 #1 I /\l \ v ‘ a I91 waw Figure 1.8 A 2D periodic network generated by introducing 250 defects into an BOO-site honeycomb lattice. 18 Periodic boundary conditions are always maintained and local relaxation follows each successful introduction of a defect. The entire system is globally relaxed every time a certain number of defects have been introduced. The relaxation is via a Keating-like potential. The details of the relaxation are described in the Appendix A. Figure 1.8 is a 2D network generated in this way. There are 800 atoms in this network, and 250 defects have been introduced into the honeycomb lattice to get this network. The network has periodic boundary conditions and contains 33.5% five-fold rings, 38% six-fold rings, 2“% seven-fold rings and “.5$ eight-fold rings. The rms bond length distortion is 5%, rms angle distortion is 16.5°. The network looks very disordered, and all memory of the honeycomb lattice seems to have been lost. We have also generated 3D networks by introducing defect shown in Figure 1.“ to the diamond lattice. The process is exactly same as in 2D but the restrictions on accepting the defect are slightly different. The minimum ring size is usually chosen to be four to make an ideal underlying a-Si lattice for an a--SiO2 network. One additional condition is that the next nearest neighbor distances must always be greater than the nearest neighbor distances. F. Wooten and D. Weaire have done similar model building, using the same kind of defects to disrupt the diamond lattice to generate the a-Si networks (References 1.20 and 1.21). We will compare their results to ours in the next section. This kind of defect was proposed merely as a means to disrupting crystalline structures in order to get continuous random networks. There was no physical basis to the particular atomic rearangement chosen. 19 However, it is quite interesting that a recent report on molecular dynamic simulations reveals that this may be a real defect structure when a liquid is cooled down to the crystallization point (Reference 1.22). F. H. Stillinger and T. A. Weber used-molecular dynamics to simulate the melting and freezing process of a diamond lattice containing 216 atoms with periodic boundary conditions. They found a local topological defect in the final crystal structure. The defect they found (Figure 1.9) is exactly the one used to disrupt the diamond lattice! This suggests a physical Justification for our model building process. (0) PERFECT CRYSTAL 1b) eowomc DEFECT Figure 1.9 Defects observed by F. H. Sillinger and T. A. Weber in molecular dynamic simulations. 20 Section 1.3 Bragg Diffraction Pattern Introducing atomic rearrangement into a crystal to disorder the system is a novel procedure to generate random networks. But one shortcoming of this method that might be crucial is that the networks are generated from a crystal lattice. There is a question whether the introduction of defects can indeed eliminate all memory of the crystal. A good criteria is needed to check whether the memory of the crystal lattice still exists in the network or not. The radial distribution function is not a good candidate, as many different models can produce similar radial distribution function. We must calculate the diffraction intensity function to judge whether the network is really random. We define the diffraction intensity functdtnill(q) which is proportional to the square of the structure factor. NE) "13'1' X e (1.3.1) For a crystal network, 1(a) should be a series of Bragg peaks on some particular a values. For a random network 1(a) should not have any outstanding Bragg peaks, rather it should have only some features of the given amorphous material and fluctuations on the uniform background. This criteria is very sensitive to any remaining crystalline features. The introduction of each defect gives substantial displacement to only two atoms (1 and 2 in the Figure 1.“) from their crystal 21 positions. A few neighboring atoms move slightly from their crystal positions after the relaxation of the network. If we can introduce defects on each pair that move each atom far enough from its original positions we should be able to eliminate the memory of the crystal entirely. But in the real process this may not be possible. We have imposed several conditions on accepting the defect in order to avoid introducing too large a strain energy, such as the minimum ring size, and the nearest neighbors being always inside the next nearest neighbor shell. When a certain number of defects have already been introduced, some local areas may become very resistant against accepting any further topological change, such as large number of small rings or large angle and length deviations etc, so that any new defect in this area may violate some of conditions. Consequently even if a large number of defects have already been introduced, a small number of atoms may still be situated close to their positions in the crystal. The phase difference between radiation scaterred by any two lattice points is easily shown to be E-fi where R is the lattice vector separating the two sites. The phase coherence of ER gives rise to the enhancement or cancellation in the diffraction pattern, this gives rise to Bragg peaks for crystalline solids. There is no such phase coherence in a random system, therefore there are no Bragg peaks in diffraction pattern of a random system. Thus we have a semi-quantitative measure of the randomness of a solid. Since our system was constructed starting from a crystal structure, there is always a danger that the deviation of the sites from their crystalline positions is not great enough to sufficiently randomize the phases of the scattered radiation and hence eliminate the Bragg peaks. This is especially true for small a which 22 gives rise to small phase differences of the order of 36R where GR is the deviation from the crystalline positions. We first take a 2D network as an example. Figure 1.10 is the Bragg diffraction pattern of the 2D network shown in Figure 1.8. The network itself appears random, but the peaks marked 1 and 2 all have a crystal counterpart. These peaks are much lower in height than that would be in the honeycomb crystal, but they are still higher than the background noise. The memory of the crystal clearly has not been eliminated entirely. The situation is even more severe in three dimension than that in two dimension. Figure 1.11 shows the diffraction pattern for a 3D network we have generated. Two hundred and fifty defects have been introduced into the diamond lattice to build this structure with a unit cell of 512 atoms. Figure 1.11 shows the diffraction pattern along the (100),(110), and (111) directions, where each is averaged over all the possible equivalent directions. Each component of ‘5 is taken to be an integer times 2n/L (L is the box length of the supercell), the x axis of the Figure 1.11 is the magnitude of the q. We can see that there are still pronounced Bragg peaks which have crystalline counterparts. The (111) peak is most pronounced here, and remains so even when more defects are introduced. There are two competing factors here, between the rms deviations in bond lengths and bond angles (a measure of the strain energy), and the height of the Bragg peaks (a measure of crystal memory left in the network). Generally speaking, introducing more defects will bring down the Bragg peaks or disorder the system more, but at the same time will make the bond length and angle distortions larger. 23 Table 1.1 gives the relationship between the rms distortions of the bond lengths and the bond angles, and the height of the most pronunced (111) Bragg peak vs the the number of defects per atom introduced, for a particular series of networks, where each time we introduce more defects on the 512-atom network previously obtained. Local relaxation followed each successful defect and a global relaxation is done after every 10 defects introduced. Only five or more fold rings are allowed, and the ratio of the two elastic force constants B/a=0.2. (dezquH 2“ a 1-~_mrmc*non PATTERN FOR 20 RANDOM NETWORK Nasaao‘1 7 d H d. ruse b 6 ~ ‘3'"? NZ-B '3' h S l " a 15$? n . T a - .°li- (Tjil“rllii. ¢redhon d; 72%;.“ D --a mo r-fl -—e——> B 18 20 3O 40 0(2'E/(N. N. DISTPNCED '1 g A A 1=%('&x"§3) - A A L.‘%{'(‘§X+‘a3) Figure 1.10 Diffraction pattern of the 2D network depicted in Figure 1.8. $8 68 58 4B 18 25 I ‘ ——-————————.-o l ’D 9‘ 3 q J I l 13 Q '2 U .31: 3 3 o '1. 3 fl 1 1 8(0) 54-51 N-512 C=B.S gmuéiea DIRECTION °_”£iiB DIRECTION. Iii DIRECTION Figure 1.11 Diffraction pattern of a 3D network, N-512, 250 defects have been introduced into the diamond lattice. 26 Table 1.1 The rms deviations of the bond length and angle and the height of the (111) peak in th Bragg diffraction pattern vs the number of defects introduced. Defects/atom rms A0 (degree) rms AL (1) Height of (111) Peak 0.098 A 11.15 2.60 1“1.0 0.195 13.86 3.39 79.0 0.293 1“.55 3.63 68.5 0.391 15.05 3.77 51.2 0.“88 15.61 3.95 35.6 0.586 15.85 “.03 3“.3 0.68“ 16.0“ “.06 27.9 0.781 16.53 ”.18 23.0 0.379 16.59 “.23 21.1 0.977 17.27 “.35 16.“ 27 We can see that the difficulty lies in how to eliminate the crystal memory without raising the bond length and bond angle distortions to an unacceptable level. Even when five hundred defects are introduced, the (111) peak is still much higher than the background noise (about 1 ~ 2), while the distortions are already significantly higher than the experiment. F. Wooten and D. Weaire did not realize the stubbonness of the crystal memory when they were generating their a-Si network using the defect method. When only a small number of defects are introduced, the bond length and angle distortions are already high. In order to reduce the strain energy, they adopted two-step annealing process. The first step is: try a defect attempt, if it decreases the energy accept it certainly, if it increases the energy, accept it with probability e-AE/KT (Reference 1.23). The second step is: introduce the defect only when it decreases the energy. This second step is a way to relax the system via the defect introduction. We shall call this step Geometrical Relaxation. This is a highly risky procedure, since the defect construction is a reversible process. From Figure 1.“ we can see that if we introduce a defect twice on the same pair of atoms consecutively it will bring them back to original structure. Wooten and Weaire really found that when only a small number of defects are introduced the system can go back to the diamond lattice through the geometrical relaxation. And even if the whole system does not go back to the crystal a small part of the network may still return to crystalline positions. We have tried the same process to reduce the strain energy and found that if the memory Of the crystal structure is not eliminated, any 28 reduction of the strain energy will sacrifice the randomness of the network. For example, on one sample we introduced five hundred defects first, then introduced additional defects only when it decreased the energy, one hundred more defects were introduced in this way. The deviations.of bond lengths and angles and the height of the Bragg peaks before and after the geometrical relaxation process are in Table 1.2: Table 1.2 rms A0 rms AL Height of the (111) peak before 16.7° 6.1% 2“.0 after 1“.2° “.7% 26.0 ‘We have also tried the step 1 process; it turns out that it always helps the distortions but hurts the randomness, or vice versa. Wooten and Weaire claimed that they obtained an a-Si network which produced a radial distribution function in good agreement with the experiment, and had a rms bond deviation 2.7% and a rms bond angle deviation of 10.9°. However they did not check the Bragg diffraction pattern. And from a private communication we know that they have only introduced one hundred defects in 216 atom system before taking steps 1 and 2 to reduce the elastic energy. From our experience, the Bragg diffraction pattern of their network should be like the Figure 1.11 or worse, implying that their model is not a random network but a quasicrystal! 29 Chn~'ultimate goal is to generate a random network. We must get rid of memory of the crystal entirely before trying to reduce the strain energy. Otherwise, any step which reduces the energy willtndng the system back to a structure which has more crystal memory than before. We finally succeeded in generating one random network with no remaining crystal memory. We did this by introducing a very large number of defects ( > 3 per atom), and allowing four membered rings. Figure 1.12 is the Bragg diffraction pattern of this network. There are no Bragg peaks. But we had to sacrifice the deviations to get this. The angle deviation is 19.8° and the bond length deviation is “.71. These about double the proper values. The ring statistics is listed in Table 1.3: Table 1.3 Size of Ring Number Per Atom “ 0.09 5 0.36 6 0.68 7 1.0“ A similar attempt to generate a random network with no four-fold rings was unsuccessful as we were unable to completely eliminate the memory of the crystal. Since we have already eliminated the crystal memory in above system, we were able to proceed with the geometrical relaxation (needed to reduce the distortions of bond length and angles) with little risk of 3O 5(0) 63-51 N-SlL _ 6-6 188 DIRECTION . 9“”9 11o DIRECTION °“'° 111 DIRECT ION A v q I ‘fi...-..-o.-C Figure 1.12 Diffraction pattern of a generated a-Si network with no crystal memory. ‘s __,,,,.....-.g —-’:.:} $23 31 returning back to the crystal structure. We successfully introduced 180 additional defects into our 512 atom system under the condition that only those which decreases the energy be accepted. This reduces the angle deviation to 17.1°, and the bond length deviation to 2.9%. Using a Boltzman factor, like Wooten an Weaire did, did not improve the deviations. Therefore, this is about the best we could do. Finally our a-Si network which has the following characteristics: (1) Large size, there are 512 atoms in the unit cell. (2) Fully bonded, with no dangling bond. (3) Periodic boundary conditions. (“) There is no long range order. There is no memory of the crystal, Figure 1.13 shows the Bragg diffraction pattern of the network after geometrical relaxation. (5) Ring statistics is listed in the table 1.“. Table 1.“ Size of the Ring Number Per Atom “ 0.1“ 5 0.“2 6 0.71 7 1.06 (6) Deviations of bond length and bond anglestare 2.9% and 17.1°, respectively. 32 This is not a good a-Si network, because the bond lengths and especially the bond angles have large deviations. However, it can supply us with a underlying a~Si network to decorate with oxygen atoms to build a good vitreous silica network. Since the two-fold coordinated oxygen atoms will bring more flexibility to the system, we can perhaps build a less distorted a-SiO2 network. The difference between our model and many other models in terms of ring statistics is that our model has a certain number of four-fold rings but the others do not have. We do need four-fold rings to make our model a good underlying a-Si network for vitreous silica network. -(-lHU)2l"'I-'12H 33 q———BRQGG DIWIW PQTTER’N SI N-512 m W‘L———— -1 H Ema we DIRECTION _ 118 DIRECTION °' '° 111 DIRECTION C--- ---- .a 5 A "v Figure 1.13 Bragg diffraction pattern of a 3D a-Si random network after geometric relaxation. 3“ Section 1.“ Construction of the Vitreous Silica Network We start with our best a-Si continuous randon network mentioned in section 1.3, which has periodic boundary, large size and no Bragg peaks. Now all we need to do is to break each Si-Si bond, insert an oxygen atom in between and bond the oxygen atom to the two silicon atoms. After this procedure is completed, the network will become fully bonded periodic a- 3102 network. To produce the experimental density we need to choose our length scale properly. We have a fixed number of atoms (512 810 units, or 2 chemical units in our unit cell), so the length of our cubic box should be chosen so that the network will have the experimental density of vitreous silica. The experimental density is 2.20 g/cm3 2 3 , or 2.206x1o’ c.u./A for vitreous silica, so our box length is determined by: -!§ - 2.206 x 10-2 L N 2.206.10- 1/3 L-( ) 2 N=512 for our particular system, so L should be chosen: 512 2.206x10- L - ( )1/3 A = 28.52 A 2 Now we decorate the network with oxygen atoms. We break all the 35 Si-Si bonds, insert an oxygen atom in between each pair and bond the oxygen atoms to the Si pairs. As is well known that the main feature of the structure of vitreous silica is that each silicon atom is situated in the center of an almost perfect tetrahedron formed by its four nearest neighbor oxygen atoms. From one oxygen atom's point of view, it should have six nearest neighbors which are divided into two groups, with three oxygen atoms plus. the‘given one making an almost perfect tetrahedron with a silicon atom in its center. The angle between two such tetrahedran can vary quite a bit.due to the wide distribution of Si-O-Si angles (Figure 1.1“ is a schematic view of the this local structure). Figure 1.1“ A schematic view of the local structure of vitreous silica. Motivated by the above feature, we adopt the following process to get a good a--SiO2 networks. We store the six nearest neighboring oxygen atoms of each oxygen atom in a coordination table, then take away all 36 the silicon atoms for the time being and relax the all oxygen system vuth.central forces only. If all oxygen-oxygen pairs had equal bond length, the four oxygen atoms would make a perfect tetrahedron, so that when we bring back the silicon atoms, all Si-O bOnds would be of equal length and all O-Si-O angles would be equal to the tetrahedral angle 109.“7°. This procedure is novel because we are now dealing1with a much simpler system, and the potential will have no adjustable ratio of force constants. The potential is written as: - N 6 + vac) )(r1.-r) ' (1.11.1) i=1 j=1 3 ' 11 is the number of oxygen atoms, and the force constant a is arbitrary. We know the Si-O distance from experiment ought to be 1.61 A (Reference 1.17) and the O-Si-O angle is 109.“7°, so that the oxygen-oxygen distance r0 should be( Figure 1.15): Figure 1.15 Nearest 0-0 distance vs related bond lengths and angles. 0 O-Si-O g 109.“7° p0 = 2 rSI“0 SIN 2 2 X 1.61 x SIN ——§—- ‘ 2.63 A 37 We relax this oxygen system until the energy saturates, finally getting rms deviation of the oxygen-oxygen bond length to be 1.5%. Then we put the silicon atoms back in the center of four nearest oxygen neighbors to get an a-SiO network. We next relax the 810 system 2 2 with a Keating-like potential: “ B 2 3 a 2 2 2 3 1 ~ + + d 2 V - -—— Z Z ( r - d ) + — -—— Z l (r -r + ) 16 d2 2 1:1 21 8 d 2 1’3 21 23 3 + 3 —E3 2 (F .; - d2 cos 8 )2 (1 u 2) 8 d2 2 i1 22 o ' ° Where a, B1, 82 are the three elastic force constants. 9. in the first and second terms is summed over all the silicon atoms, 1,3 are nearest neighbor oxygen atoms of silicon atom 11., 2. in the the third term is summed over the all oxygen atoms, 1 and 2 are two of an oxygen's nearest neighbor silicon atoms, d-1.61 A, and 00 is average Si-O-Si angle. There are three elastic force constants to be chosen. The relative ratios of d, 8,, 8 will determine the rms deviations of 81-0 bond 2 lengths, the O-Si-O angles and the Si-O-Si angles. We take the value of B,/u.- 0.2, and 82 . 0, based on the wide distribution of Si-O-Si angle and the narrow distribution of O-Si-O angles in vitreous silica. Finally we get our a-SiO2 network following relaxation of the network with the potential (1.“.2). Our network has following characteristics: 38 (1) Large size - there are 512 810 units in the systemn giving 512 2 silicon atoms and 102“ oxygen atoms, the total number of atoms is 1536. (2) The network is fully bonded with periodic boundary conditions. (3) The network is random, there are no Bragg peaks in the diffraction .y pattern ( Figure 1.16 ). Here I(q) is defined as: + 1 N - I(q) = 'fi 2 b e (1.11.3) Where 61 is taken to be the isotropic average neutron scattering lenghth of the i'th atom, b / 60 - 0.715. Si (“) This network has the experimental density p . 2.206 x 10_2c.u./A3. (5) The distributions of the bond lengths and angles are listed in the table 1.5: Table 1.5 Average Value rms Deviation Si-O Bond Length 1.62 A 1.2 1 O-Si-O Angles 109.“° 5.8° Si-0~Si angles 1“7.“° 18.8° (6) The ring statistics are same as that listed in the table 1.“, since we did not change the topology of our a-Si network. We made our a-SiO network based on our a-Si network, which was 2 made from diamond lattice by introducing a large number of defects while limiting the smallest ring size to be four. The above ring statistics 39 comes naturally out of this procedure. Although four is likely to be the smallest ring size from energy considerations, some people do suggest the existence of a small number of three-fold planar rings (Reference 1.17). We have actually generated an a-SiO2 network based on this conjecture to see the effect of the three-fold rings (Reference 1.23). Figure 1.17 shows the distribution of the 81-0 bond lengths and O- Si-O angles and Si-O-Si angles. We now compare our a-SiO model with some previously existing 2 models and experiment data in the Table 1.6. HT stands for our model, 80 is the hand-built Bell and Dean model, GS“, H61 and Y62 are models built by W. Y. Ching by decorating already existing periodic a-Si networks. no Table 1.6 Comparison of our HT model, some other models and the experiment. Model HT G5“ H61 Y62 BD No.Of Atoms 1536 162 183 186 61“ Type Periodic Periodic Periodic Periodic Cluster Density(g/cm3) 2.20 2.20 2.20 2.20 1.99 RSi_O 1.62 1.62 1.62 1.6“ 1.62 %5 x1001 1.2 0.7 1.0 1.7 3.5 6(0-31-0) 109.“° 1o9.3° 1o9.3° 109.2° 1o9.3° A0 5.8° “.8° 6.8° 8.2° 6.2° $(s1-o-SI) 1“7.7° 1“7.2° 152.3° 151.5° 153.3° A0 18.8° 13.8° 1“.2° 13.1° 9.6° Exp 2.20 1.61 109.“7° 1uu°-152° 15°~2o° <4Hm2n42” 1—-—-—-—— m DIFTRRC'I’ICN PRTTERN N-1536 5192 Bragg diffraction pattern of our a-SiO “1 AF—fi was 180 DIRECTION . no DIRECTION °" ‘° 111 DIRECTION T 1 ‘1 1l 9 I 1% “ I?l’ l I T E T 1’ k 1 L 11111 1 11.111111111111111 8519152825 V ”—0 W as 35 40 45 so 01%: Figure 1.16 2 network. 429103”? «znnanb 421-10:an . 8? .83 .81 “2 q‘SI-O m LDCTH DISTRIBJTIG'C $102 "-155“ I l I I 1 l l 1.4 1.45 1.5 1.55 1.6 1.65. 1.7 1.73 1.9 '8 11 111 . 81$ .01 (A) MST-O m: DISTRIach smz muslin“fl A. I I I I I I ‘1 ' I B 23 43 60 ea 1% m 140 360 IN 3% .HI-O-SI “(LE DISTRIIJTIM SIQ "-1 q 0 29 48 60 a In 123 16 m Figure 1.17 16. 18 81-0 bond length, O-Si-O angle, and Si~0~Si angle distributions of our a-SiO network. 2 “3 Section 1.5 Radial Distribution Function We now compare Radial Distribution Function of our model with the experimental results. The experimental data we are comparing to are from A. C. Wright's neutron scattering experiment. Function T(r) is used in the experiment, so we also calculate the function T(r) from our model to compare with the experiment data. T(r) is defined b! A. C. Wright as (Reference 1.2“): T(r) - “nrJZbJ {bi gji (r) (1.5.1) Where DJ is isotropical average of neutron scattering length of the j'th atom. BJ - 0.“1“9 x 10"12 cm if the j'th atom is a silicon atom DJ - 0.580“ x 10‘.12 cm if the j'th atom is a oxygen atom gJ1(r) is the atomic distribution of i type atoms about the origin (3) atom. j is summed over a chemical unit, i is summed over all atoms. For large r, if p° is the number density, we get the asymptoth: form T() r""'°+11 (266 +266 +66 +66 ) r n r o sipsi o 090 si 0 po 31 31931 .0 I 0 1131 p . 11° 211 1111””111112'5610126221012661016218 o si o si o si 1111 . “ 11 r p° ( 3: 63) Figure 1.18 shows the comparison of the T(r) calculated from our model and that from neutron scattering data of vitreous silica. They both have the same asymptotic behavior for large r. The first and second peaks of the model are sharper than those of experimental, because we do not include the effects of the thermal broadening in our T(r). The first peak of T(r) is the 81-0 peak, this peak is very sharp, because the length of 81-0 bond has a very narrow distribution. The position of the first peak is slightly shifted to the right from the experimental peak, because the average Si-O distance in our model is 1.62A, which is slightly higher than the experimental value of 1.61A. The second peak is the 0-0 peak. Our model has a narrow distribution of 81-0 bond length r and the narrow distribution of 0- 81-0 81-0 angles 90-Si-O° Thus the 0-0 distance r0_O is given by eO-Si-O rO‘O = 2 rSi’O ° SIN( —-—2——-— ) (1.5.2) ought also to have a narrow ditribution, giving a sharp second peak in T(r). The third peak is mainly decided by the nearest Si-Si pair. This distance rSi-Si is decided by: e Si-O-SI rSI_Si . 2 rSi-O SIN 2 (1.5.3) “5 the value of Si-O-Si angle 0 has a wide distribution, so that we Si-O-Si will not get a sharp Si-Si peak. Our model has a very small third peak. The experimental T(r) does not have a peak but only a shoulder. We can see all the above features more clearly from the partial pair correlation functions. Figure 1.19 are 81-81, 31-0, and 0-0 pair correlation functions. These correlation functions are defined so that the area under the curve gives the number of atoms in the corresponding value of r. We can see the sharp first peak of Si-O andand the broad first peak for 81-81. Further 81-0 and 0-0 neighbor may also fall in the region of the 81-31 distance which will broaden the 81-31 peak in the T(r) further. From Figure 1.18, we can see that T(r) of our model reproduces all the experimental features, with overall good agreement. “6 T(R) SIOZ N=512 BET1-0.Z BETZ-O.B as as - ""n«»1mm£L """ FFDNIDGERDEJW 15'- .- ImrnSJ 2' Io ‘ 5 .. a “n’7 a 1 2 3 14 5 6 7' a 9 Rdh Figure 1.18 Comparison of T(r) calculated from our network and that from neutron scattering experiment. 15 13 4G 10 Partial correlation function of our 810 -—§SI-SI MIR WTIW WIN SHE N'1 I ‘ I I I I I I I O 1 2 3 4 5 6 7 B 9 R18) l—AI-O PAIR WTIU‘I WK!“ SIM I'I'-IE'I$-—___1 d 1 I I I I I I I o i z 3 4 s 6 7 a 9 R15) .1—__o-o PAIR CORRIIATION NCTION 5102 was 1 I I I j I I I a I 2 3 A s 6 7 a 9 R15) Figure 1.19 2 model. “8 Section 1.6 Conclusions From our work we would like to make the following conclusions: (1) The introduction of defects into a crystal structure to generate the structure of an amorphous solid is a novel procedure, which can easily be done on the computer while maintaining periodic boundry conditions. (2) The main difficulty of this method is that it is very hard to eliminate the memory of the crystalline structure. Calculation of the diffraction pattern of the network is a reliable way to judge if the memory of the crystal structure is still retained or not. The RDF can not easily distinguish between a quasicrystal model and a continuous random network model. (3) To eliminate the memory of the crystal and to keep the distortions of bond lengths and angles under control are two competing requirements. So far we are unable to reach both goals to get a really good a-Si network. We have generated an a-Si network which has no crystal memory, but the angle distortion is much larger than that observed in experiment (17° for our model, while 10° for experimental value and for relaxed hand built models with free boundaries). (“) Based on our a-Si model, we have generated a good model of vitreous silica, which has a large size, periodic boundaries, is fully “9 bonded, has no Bragg peaks in the diffraction pattern, has reasonbaly small Si-O bond length and O-Si-O angle distortions, and whose RDF agrees with experiment. Our model contains four and higher fold rings. (5) Our fully bonded large periodic random network is a good model to be used to study the various properties of vitreous silica. S. W. Deleeuw and M. F. Thorpe have used the model to study vibrational properties of vitreous silica (Reference 1.25). Chapter 2 Elastic Properties of Glasses Section 2.1 Introduction The introduction of a small amount of impurity (Phosphorus, Boron for example) in the crystalline tetrahedrally bonded semiconductors can change the conductivity a great deal. In contrast, the electrical properties of chalcogenide glasses do not exhibit great sensitivity to impurities, and are in fact often insensitive in compositional change of even a few percent. The insensitivity of the electrical properties of the chalcogenide glasses led N. F. Mott to propose the "eight-minus-n rule" (Reference 2.1). The valence of each constituent atom is satisfied by every atom in the covalently bonded amorphous solids. The number of electrons needed to fill the outermost shell of each atom is 8 - n, where n is the number of s and p electrons of the outermost shell. Therefore each atom in the covalently bonded amorphous solids form Z 8 - n covalent bonds with its 2 nearest neighbors. In other words, Z - 8 - n is the coordination number, which is satisfied by each atom in the glasses. Consider a ternary system like GexAsySe . Since respectively 1-x-y the valence of Si, As, Se are “, 5 and 6, the coordination numbers are “, 3, and 2 respectively, when we form a compositionally and topologically disordered system out of these three species, they forbIaa 50 51 continuous random network with a definite coordination number for each species. Figure 2.1 schematically illustrates this features of the networks (Reference 2.2). There each Ge atom is connected to four neigboring atoms, each As atom is connected to three, and each Se atom is connected to two neighboring atoms. \ \ SC Se Se] 52 _ \ SE 1 \S X X \ 59 58- —Se e Se- \ A S Se Se / 1],; S‘“——s 1 Se Se -- Se \ 1 e \s ‘ \/ /Se 1 Se SE \ S’e ‘Se Se se/ 68 58/1 ~Se \ Se 5’ /5e /Se 5 9 Se .e’\/S€\Se/ 59 Se— -—-Se\se e Sf Se Se\ Se Se Se / /\Se/\ \ Se 58"" Se \Se/AS /Se \ 58/59 Se; 58 I Se / / . / .Se \ Se ‘ / Se 9 S Se Se '5 Se \ Se ‘\ e \\ \, e Se Se 59/ s \ /s \ " 8 Se ‘ Se e 5 Ge \Se— Se\/ \ Se T / A SC 2 ‘5 Se \ e S \S ’1 / Se/ ‘Se\5e Se K s \ / 8 Se \ / $8 52 Se -Se Se. Se S I 58 \Se \' 58s / e L Se\|‘~se Se-Z—AS . I se/ \ 98 Se be Se Se \\ I S; Se\ Sex \ / Se Figure 2.1 Schematic bonding topology of GexAs Se Glasses. 1-x-y . As is well known, when isolated Si atoms come together to form a crystal, the four electrons in the outermost shell (two “3 and two “p electrons), occupy the four hybrid bonding orbitals to bond each atom to its four nearest neighbors, forming a tetrahedral structure (Figure 52 2.2). The bonds are very strong and directional. This feature is also present in amorphous silicon with only some deviations in bond lengths and bond angles when continuous random networks are formed. This structure can be generalized to all covalently bonded glasses. The dominant forces are covalent forces, so that any Change in bond lengths or bond angles causes a significant increase in strain energy. \ Figure 2.2 Hybrid bonding orbitals of Si in the solid. More distant forces, van der Waals forces (by van der Waals forces we mean all other forces other than strong covalent forces), are rather weak in comparison with the covalent forces. Therefore, it is reasonable to neglect those weak forces in a first approximation. 53 (b) Chains with cross links introduced by As and Ge atoms. Figure 2.3 If we take this approxihmtion, that is we include only bond stretching and bond bending forces and neglect the dihedral angle forces and other forces, we may find two different kinds of regions in a 5“ GexAsySe1__x__y system. One kind is Se-rich regions, where Se atoms form long polymer chains (Figure 2.3 (3)). There are two constraints corresponding to each Se atom: one is from bond stretching (there are two bonds connecting each atom, but it is shared by two atoms), another one is from bond angle bending. However, since there are three degrees of freedom per atom in three dimension, we can twist these chains freely with no cost in energy because the Se atoms can always adjust their positions to relieve any external stress. We call this kind of region Floppy Region. Obviously, a few cross links, which interconnect the chains, stabilizes the structure or makes it rigid ( Figure 2.3 (b)). The higher coordinated atoms ( Ge and As ) introduce these cross links to otherwise isolated chains to stablize the structure. Therefore the Ge or As rich regions are stable and resist any external stress; we call this kind of region Rigid Region. The floppy regions are characterized by low mean coordination. The rigid regions, on the other hand, are characterized by high mean coordination. There are essentially two different situations when rigid regions and floppy regions coexist in the same system. One possibility is as shown in Figure 2.“ (b), where the rigid regions percolate through the system and form an infinitely large rigid backbone (if system is infinite). Second possibility is depicted in Figure 2.“ (a), where the rigid regions just form a few isolated islands in the matrix of the floppy regions. This picture corresponds to the basic notion that there are two classes in the covalent glasses: Amorphous Solids with high Mean Coordination and Polymeric Glasses with low mean coordination. 56 The questions of whether or not the rigid regions percolate for a given mean connectivity defines the Rigidity Percolation problem. These two different systems are very different in terms of elastic properties. The systems where the rigidity percolates are hard, resist any external forces, and their elastic moduli and longitudinal and transverse sound velocities are finite. Those systems where rigidity does not percolate are soft, does not resist any external forces, i.e. the system can be continuously deformed with no cost in energy, the longitudinal and transverse sound velocities are zero (or very small in a real situation). Rigidity percolation is different from regular geometric percolation in terms of Percolation Threshold. One example is the central forces only model. The potential is expressed: 1 + “ 2 V - 5’0 Z.( uij r13) (2.1.1) 1.3 Here rij is.unit vector along the bond (i,j) in the network, fiij-Ji-GJ is relative displacement of the atom on the either side of the bond from their equilibrium positions. Figure 2.5 (a) shows the situation when the two rigid pieces are connected by a single bond, forming a free hinge. Although the two rigid parts are connected geometrically, they are elastically disconnected, since a free hinge can not resistIan external stress. If there is a voltage imposed on the two sides and the system is a resistor network the electrical current is certainly nonzero. But if a stress is exerted on the two sides the system does not resist,implying 57 (a) Rigid units connected by a single hinge. / (b) Rigid units connected by a single chain. Figure 2.5 58 that the elastic moduli are zero. The rigidity percolation threshold is therefore different from the conductivity or geometrical percolation threshold. Another example is a chain-like structure in the bond stretching plus bond bending forces model in three dimension (Figure 2.5(b)). The potential energy can be written as 2 2 1 -> 1 v-—ai(u - )+—82(8) (2.1.2) 2 1’3. 13 13 2 ijk ijk ”3+ When an external stress is imposed on the structure in question, the system would try to relax to reach the minimum energy state. The question is whether the five atoms connecting the two rigid pieces can relax to the proper positions to make the strain energy zero. Mathematically we have 5 x 3 = 15 unknowns and 6 (to fix six bond lengths) plus 5 (to fix five angles in between) plus “ (to fix four angles on the sides) - 15 equations, so that there is always a solution in general. This structure therefore cannot resist an external stress. Again we have the situation where the geometrically connected regions are not necessarily elastically connected. However, if we include more forces than we did in above two cases, life can be quite different. For instance, if there are angle bending forces in addition to the central forces in the two dimensional problem or dihedral angle forces in addition to the bond stretching and bond bending forces in the three dimensional problem above, these two systems would become able to resist external forces. Their elasticity percolation threshold would become identical to their conductivity percolation threshold. In other words, whenever the parts of the 59 networks are geometrically connected they are always elastically connected as well. In summary: when sufficient forces models are adopted, the elasticity percolation problem becomes identical with the corresponding; conductivity problem in regards to percolation threshold, and we call this kind of elasticity percolation problem a class 1 problem. Although for the class 1 problem, the elasticity percolation systems and their conductivity percolation counterparts have identical threshold, these two problems are different in physical nature (A partial explanation of this difference is that conductivity is a tensor of rank 2 in contrast to the elastic tensor which is a tensor of rank ‘1). Various studies by a number of workers (References 2.3-2.6) suggest that the elasticity percolation problem belongs to a different universality class from the conductivity problem although they have the same percolation threshold. We define t as the critical exponent of conductivity t o c ( P - Pc) (2.1.3) Here p is the fraction of bonds present. We can also define the critical exponent for the elastic moduli as follows: c .1 1 p - PC)T (2.1.11) 60 The coductivity exponent t has been found in various numerical and real experiments to be approximately (References 2.3-2.6) t . 1.2 - 1.3 for 20 t . 1.9 - 2.0 for 3D The exponent T for the elastic moduli has been suggested to be substantially higer than the above values for t (References 2.3-2.6): T - 3.6 for 2D T . 3.5 Ifor 3D Insufficient force models make the elasticity percolation threshold different from regular bond percolations. We call this class 2¢problem, or Rigidity Percolation problem. J. C. Philips (Reference 2.7) and M. F. Thorpe (Reference 2.8) have developed the constraint counting argument to locate the rigidity percolation threshold. S. Feng, M. Thorpe and E. Garboczi (Reference 2.9) have studied the central forces only model. They developed an Effective Medium Theory and have done computer simulations on various 20 and 3D networks. Effective medium theory, constraint counting and computer simulations all agree that the rigidity percolation threshold is very close to: Pc-2d/Z (2.1.5) 61 Where d is dimensionality, Z is coordination number. This is substantially higher than the bond percolation threshold Pc~%:% (Reference 2.2). The aim of this research is to study the elastic properties of glasses via computer simulations through a model which is much simpler than the real glasses but we believe it captures the essential physics (Reference 2.10). We use the Keating potential in three dimensions. The Keating potential essentially includes the bond stretching and bond angle bending forces. We neglect the dihedral angle forces at first. This defines a class 2 problem, so that we expect a percolation threshold that differs from geometrical percolation threshold. As we will later on consider the effect of the introduction of the dihedral angle forces. We know that the covalent glasses can be divided into two classes, the amorphous solids with high mean coordination and polymeric glasses with low mean coordination, we are led to study the transition between these two classes of glasses by probing the behavior of two physical quantities such as fraction of zero fraction modes and elastic moduli which couple most directly to the phase transition. The rigidity percolation threshold (class 2 problem), which is different from the geometric percolation threshold, is not known exactly so far. The numerical computation which involves a relaxation Of the system with large floppy regions is also very difficult. We will not try to probe the critical region in any detail. 62 Section 2.2 Constraint Counting and Zero Frequency Modes The notion of Overconstrained and Underconstrained glasses was introduced by J. C. Phillips in 1979 (Reference 2.7). This idea was later put on a firm foundation by M. F. Thorpe (Reference 2.8), who envisaged a glass consisting of Rigid and Floppy regions and a phase transition taking place as the mean coordination is increased and rigidity percolates through the network. The number of zero frequency modes Mo is equal to the number of ways in which the network can be continuously deformed with no cost in energy. In general, M0 is given by the number of degrees Of freedom minus the number of linearly independent constraints NC: M =dN-N (2.2.1) o c ,, .. Here d is the dimensionality, N is the total number of atoms in the system. One way of seeing the validity of equation (2.2.1) is that NO is equal to the rank of the dynamical matrix D (R), so that the nullspace of D(R) has just d N - Nc elements. Mo should be non-negative from its physical meaning. When Mo = 0, there is no way that the system can be continuously deformed with no cost in energy, this is equivalent to the case when the rigid region percolates through the system, so that the system is rigid, the longitudinal and transverse sound velocities are finite, elastic moduli are finite. 0n the other hand, however, if M() > 0, there are some ways 63 in which the system can be deformed with no cost in energy. This case corresponds to the situation when the rigid regions are isolated, system can not resist any external forces, the elastic moduli are zero, longitudinal and transverse sound velocities are zero. Nc is the number of independent constraints. Nc depends upon the choice of the elastic potential expression and the structure of the network. One simple example is the central forces only model, where the potential can be written as ( 2.2.2): 2 v -1/2 a Z (1 -1 ) (2.2.2) 1 i O ( all bonds ) Each bond gives rise to one constraint in the system. If we have a lattice where each site has 2 nearest neighbors, there will be 2 x N/2 bonds altogether so Nc - Z x N / 2, from (2.2.1) the number of zero frequency modes should be approximately (see below): Mo = d N _ Z N / 2 (2.2.3) We can define fo as the total number of zero frequency modes divided by the total number of normal modes: f .. .____9__ (2.2.“) 6“ The computations of various lattices of 2d and 3d systems show good agreement with the above mean field estimates (E. J. Garboczi, Thesis and Reference 2.9). We now consider the model we are interested in. We consider a 3D random network with N atoms in it, composed of a few different species of different coordinations. We have up atoms having r bonds where r . 2, 3, or “ for covalent networks with the sum rule: 2 n - N (2.2.5) We include both the bond stretching forces and angle bending forces. So each bond gives rise to one constraint related to the bond stretching. There are r/2 such constraints to each r coordinated atom, with the factor 1/2 due to the sharing of the bonds with its nearest neighbors. There is one angle for each two coordinated atom, three angles for each three coordinated atom and five angles for each four coordinated atom. There are five angles instead of six for a four coordinated atom because the sixth one is no longer independent. In general there are 2 r - 3 angular constraints for each r coordinated atom. The total number of constraints is then given by: Nc - E nr [ r/2 + ( 2r —3) 1 (2.2.6) The number of zero frequency modes is: Mo . 3 N - Nos 3 N - 2 hr [ r/2 + (2r-3) ] (2.2.7) 65 Define mean coordination as: 2 r nr I" < I" > - fi— (2.2.8) r r The fraction of zero frequency mode will then be: fo . Mo / 3 N . 2 - 5/6 < r > (2.2.9) The transition should take place at Mo = 0, i. e. f}>- 0 which implies that < r > = 2.“ . (2.2.10) The relationship (2.2.9) can be depicted as a straight line (Figure 2.6) where fo - 0 at - 2.“ and fo - 1/3 at -2. At .. 2 every atom is two coordinated so that the number of constraints for each atom is 2 (one from bond stretching and one from the angle formed by the two bonds). But as there are three degrees of freedom for each atom, that implies that 1/3 of the total modes will have zero frequency. 66 “1'; \k «(1') Figure 2.6 Fraction of zero fraction modes fo vs mean coordination from constraint counting. Here we should notice that the expression (2.2.1) is exact as long as NO is the number of independent constraints. However the expressions (2.2.3) and (2.2.7) are not exact, because we did not consider the effect of the formation of closed rings in the network on the number of constraints. Actually, due to the existence of closed rings the constraints we counted above are not all independent (Reference 2.8). This kind of constraint counting is mean field like in nature. One simple example to view the effect of the formation of the closed ring to the number of independent constraints is the 2D network with bond stretching forces plus angle bending forces (Reference 2.8). Mean field like constraint counting would lead to count the number of 67 bond lengths and the angles they formed (notice that the each r coordinated site gives rise to only r-1 independent angle constraints since the last angle is not independent of the others) in the network. However, when the rings are formed, Figure 2.7 for example, the dashed bond and angles indicate that these three constraints are unnecessary as they are determined by specifying the other bond lengths and angles in the ring (Reference 2.8). Figure 2.7 Schematic view of the effect of the closed ring on the number of independent constraints. From above constraint counting argument, we can see that if we neglect all weak forces, there will be a transition between the rigid system where the rigidity percolates, and the floppy system where the rigidity does not percolate. Transition takes place at < r > - 2.“. The implication of this in the real glasses where there are always other weak forces in addition to the covalent forces is: Systems with >2.“ are amorphous solids, the elastic properties of these materials are characrerized by strong covalent forces. Systems with 2 < < 2.“ are polymeric glasses, the elastic properties of the materials are characterized by weak van der Waals forces. 68 Section 2.3 Model Description Our aim is to develop a model to calculate some physical quantities, such as the fraction of zero frequency modes and the elastic moduli, to see if the mean coordination does play the decisive role in determining those quatities and if a transition really takes place at about - 2.“, as predicted by the constraint counting argument. We choose Keating potential (Reference 2.10) for the ealstic potential energy. The form of Keating potential is: 2 2 2 v- °‘ 2 (r -d) T6d2 2,1 9.1 38 -> + 122 + ------—- (r I" .+‘Td) (203-1) 8 (12 2.1.3 ii 13 3 DJ Where :1.) is the vector from the ith atom to jth atom in equilibrium, the d is the nearest neighbor distance in equilibrium. This potential expression includes bond stretching forces and angle bending forces but not the dihedral angle forces. The free choices of dihedral angles make possible the existence of regions which may be geometrically connected but elastically disconnected, i. e. this is a class 2 problem. 69 We choose our system to have periodic boundary conditions so that we do not have surfaces and so that a relatively small unit cell can be used to simulate the infinite system. To choose the actual system on which to do computer simulations we still need to take into account the following problems. Since elastic energies are small deviations from the local minimum of the vibrational potential, in order to calculate the elastic moduli we have to calculate the small change in elastic energies between those before and after external strains are imposed. In order to do this, we need to know the energy of the local minimum to high accuracy. Another factor we have to consider is that we must be able to change the mean coordination continuously over a wide range, so that we can see if is an important parameter, and if changing does lead to a phase transition at about - 2.“. Also each time we change the mean coordination we still need to know the energy Of local minimum of vibrational potential to high accuracy. With the above considerations in mind we chose to work on the diamond lattice, where the the initial coordination is four. We then cut bond at random, so that the mean coordination number decreases from four and can be eventually be decreased to close to two when the network is mainly polymer chain-like. We maintain the periodic boundary conditions so that if a bond is removed in one supercell it is removed in all supercells. At first only three-coordinated sites are created, but with enough bonds removed, two coordinated sites are also created. Although each atom may have different coordination number after the cutting, we still suppose that they all have same position and same interatomic force constants for the remaining bonds. Therefore, the vibrational ground state energies are always zero for the bond-diluted network, 7O implying that all atoms remain in their positions on the underlying diamond lattice. Although this is not very realistic, we believe it can capture the essential physics. The potential energy can be written as follows after the above simplifications. V ' 76"E‘ 2 ( Iii ’ °2 )2 P21 ._ d 1,1 38 2++122 + ._.——— ( r . r , + — d ) P P (2.3-2) 8 d2 2,1,j ii 23 3 ti ij _ . i>j Where P a { 1 if i , j are bonded ij 0 otherwise FTwnn potential (2.3.2), we can see that when a bond is removed, ii for example, all the a and 8 terms associated with it are removed from the potential (2.3.1). We keep the minimum coordination to be two (we do not allow any dangling bond or isolated atoms). so that we can mimic the Gex Asy $61-x-y system. The size of the unit cell we use to calculate the elastic moduli is 8 x 8 x 8 - 512 atoms. After a certain number of bond are removed we have n two- 2 coordinated sites, n three-coordinated sites, and nu four-coordinated 3 sites, and n + n + nu = N (2.3.3) 71 +3n+“n)/N (2.3.5) Mean Coordination < r > - ( 2 n 3 “ 2 Since the equation 2.3.3 has to be satisfied, we have only two independent variables or we have one free parameter. For the same we can have different n2, n3, and nu values (two of them are independent). In order to know if the physical properties concerned are decided by and n only the mean coordination or also the detail of the n n 2’ 3’ “ combinations, we developed three different samples with the same mean coordination but different values of n2, n , and nu. 3 Figure 2.8 shows the differences between the three samples. The x axis is mean coordination, the y axis is the fraction of 3-coordinated sites n3/N. The first two samples are arrived at by totally random cutting. We first pick a site at random, then choose one of its neighbors at random. If we would create a dangling bond by cutting this bond we do not cut it and just move to the next site; otherwise, we always cut it. To get a sample which is substantially different from the first two samples in terms of II (but with the same mean coordination) , 3 we adopt a weighted cutting procedure. We deliberately picked one configuration with an enhanced number Of two-coordinated sites n2, thereby depressing n3, as shown by the circles in Figure 2.9. To do this we first choose a bond at random, if both ends of the bond are 3- coordinated already, we cut the bond always, if one end is 3-coordinated we cut with probability 0.“, if both ends are “-coordinated we cut it with probability only 0.1. This process enables us to create the third sample which has a very different value of n3/N than the other two. 72 076 I l I 1 0.4 — _ “3 ‘. . N . ‘ ‘ U A I“ 0.2 P- ,‘.I' _. 3'. I 6'. an . «500' - O l I l l d 2 2.4 2.8 3.2 3.6 4 < r> Figure 2.8 Number of three coordinated sites vs mean coordination for three samples, N-512. 73 Section 2.“ Fraction of Zero Frequency Modes We calculate the fraction of zero frequency modes by directly diagonalizing the dynamical matrix. The energy of vibrational motion of atoms in the solids can be written as: 1 d 2 1 a a8 8 H - Z -—— P + — X Z u D ,u (2.“.1) 2 a 2 mg i 2 l a 2'8 t it 2' Here, a, B = x, y, z and i = 1,2,3 . . . N, all atoms. 2 D°°,. .__§_I___ (2.“.2) 2“ an“ an“ i i' The equation of motion is: 0d 3 H P2 ' a ’ ml “2 (2'"*3) 3 u i d - i w t For the normal modes uz ~ e Substitute u: in the equation (2.“.3) we have: 2 a dB 8 It w ”I £18 D221 “2' (lfuf3) Take mass of all atoms to be equal and equal to one for simplicity. ‘Now to get the frequencies of the normal modes, all we need to do is to 7“ (:8 2.2." expression of dynamical matrix element as follows: diagonize the matrix D From potential (2.3.2), we can write the Dm'“- 32" 2",!" Bumaun l 2' m n m n ' a I 6221 2 P21 r11 P21 rll'ril'Pil'] m m n n + 1/“ 8 62l'[ 2 ( rg11+ r11)(rli'+ r11) P11 P11. 1,i'¢l m n + Z r. P P P ] i,£*l 15 15 ii 15 + 1/“ 8 [ Z ( rm + rm ) rr1 P P 1‘1, ii 211 It 21 22' n n m + Z ( r I + r. I) r. I P ' P I ] tel I E It I E It i E 1 1’" 8 Z ’21 rTi'Pti P111 (quf‘) 1 Here i, 2' go over the atoms, and m,n go over the x,yn 2 directions (we use m,n other than a, B, because we have already used at and B for the elastic force constant). We have done this on systems smaller than those used for the elastic moduli calculations. We used cells with 6 x 6 x 6 - 216 sites. We developed three similar samples for the 216 atom system. 75 Figure 2.9 (a) show the n3/N vs mean coordination relation for these samples. The process of developing these samples was same as that for 512 atom system. We directly diagonize the above dynamical matrices by using existing subroutine in the library on the CYBER 205 supercomputer and count the degeneracy of the eigenvalue m2- 0. We always subtract out the three so-called Goldstone Modes which correspond to the translational motions of the whole system . Figure 2.9 (b) shows the results of three different samples. It is clear that they are very close despite the large differences in n2, n3, and nu values. From this we can say that the mean coordination does play the decisive role in determining the fraction of zero frequency modes. The constraint counting argument predicted such a result. Figure 2.10 shows fo averaged over the three different samples vs mean coordination. The straight line is drawn from 2.2.9 which is the constraint counting prediction. The calculations show very good agreement with the constraint counting. The small deviation from the straight line when is very close to 2.“ is due to the existence of small isolated floppy regions which persists when the rigid regions percolate and also to finite size effects. As was mentioned before we do not expect the way we count the total number of constraints to be completely accurate. In the network there are always some closed rings, so that some constraints we count are not independent. This is why floppy inclusions can still exist when > 2.“. 76 006 ’I‘ZIE 8.5 ‘ 'QU on. W I: II 94 _ W 8.3 "' I‘ ' 4’ . ‘60.,“ 6.2 .. .VAt/N H . Rmsct- 9 :1: remscT- 87 SHEET-172 00.1 " B I I I I I I I I 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 mm COWDXIRTION (a) Fraction of three-coordinated sites vs mean coordination for three samples, Nn216. aha ;:::F¥WBET- 3 o.___<>RRI~11ssn'1- s? Rmnxx-rnz 7% - 8J2 - 2.1 . g .M A- I 7 Tau-T96 I __ I 2 2.2 2.4. 2.6 2.8 3 3.2 3.4 IERN Figure 2.10 Fraction of zero frequency modes vs mean coordination, averaged over three samples. 78 Section 2.5 Elastic Moduli As we argued before, the system becomes soft when the mean coordination decreases and it will undergo a phase transition when the rigidity ceases to percolate through the system. We estimated this threshold to be = 2.14 for our systems using constraint counting argument. The system we start with is the diamond lattice. This is a cubic system. There are three independent elastic moduli for a cubic system, C11’ Ci moduli from his potential. He writes down the form of the potential 2, and Cu“. Keating (Reference 2.11) has calculated these elastic energy when a general external strain is imposed, assigns an internal displacement u', v', w' between the two FCC sublattices along the x, y, and 2 directions, respectively. He solves for the u', v', w' by minimizng the strain energy with respect to these internal strains: Then he substitutes the u', v', w' into the potential energy and reads off the value of C C 2, and CM from the general expression 11' 1 for the strain energy of a cubic system (Reference 2.12): 1 2 2 2 V - 2 C11 ( exx + eyy + ezz ) + C12 ( exxeyy + exxezz + eyyezz ) 1 2 2 2 + 2 CH” ( exy + exz + eyz ) (2°5'1) Keating's result for C11, C12, and Cu“ for the diamond lattice is: 79 1 CH-Ta(a+38) (2.5.2) 1 C12 = r5 ( a "’ 8) (2.5.3) (2 = “8 (2.5.1:) 414 a(c+8) Here a / J3 is the nearest neighbor distance. The bulk modulus can be obtained from C11 and C12: C +2C B- ”3 ""=(3a+s)/12a (2.5.5) When we remove the bonds so that the mean coordination is close tc> the percolation threshold, some very floppy regions may occur, so that it may take a long time to relax the system completely. We use trimming process to save the computer costs. The trimming process is the following We first identify the part of network, which we know will now contribute to the strain energy, when the network is fully relaxed, in other words we know we can always relax that part to a zero energy state. Then we remove all the bonds in that part of the network when we compute the elastic energy. The following structures can always be trimmed: 80 Figure 2.11 Chain-like local structure which can be trimmed. 1. Chain like structures (Figure 2.11). If there are five sites in the middle part which are two coordinated, and both ends are three coordinated, from the argument in section 2.1, all bond lengths and bond angles in this part of system'can always be satisfied, so that the elastic energy can be relaxed to zero exactly. The dashed curved bonds can be all trimmed away. 2. All isolated clusters can certainly be trimmed, as they always relax to zero energy. 3. After trimming of the chain-like structures some parts of the network can be connected to the backbone through a single bond (Figure 2.12). This dangling part can be trimmed away as well. 81 Figure 2.12 Singly connected cluster which can be trimmed. When the mean coordination is sufficiently low, all bonds may be trimmed, thus giving us a lower bound for the rigidity percolation threshold. The Figure 2.13 shows us‘how many bonds are still present after the trimming process. The y axis is the fraction of the bonds present after trimming. The x axis is the mean coordination and the curves stand for three different samples which we mentioned in the section 2.3. We can immediately see that at about - 2.2 all bonds are gone through trimming, so that the elastic moduli are certainly zero even without any calculations. 82 1. W________.nmx310Ncalculate the elastic moduli we impose an external strain 6, relax the system with the Keating's potential (for detail of the relaxation see the Appendix (A)), and calculate the elastic energy after the energy converges. The final energy can be written as: 2 C 6 (2.5.6) va- 13 NI—b C is the elastic modulus. The appropriate external stnahlis iJ imposed to extract the desired elastic modulus. Since C can not be 12 measured directly, we calculate C Cu“, and the bulk modulus and 11’ extract C12 via equation (2.5.5). 2 Unlike the crystalline cubic system, for the random system, C11and can do rely on the choices of the external strains. To get C we take 11’ 6 to be e , e , e in turn and do the average over these three xx yy zz directions to get better statistics. To get C1m we also do average over three different directions. The actual process is: For C take x direction (exx - 6) for example, change the box 11’ length along x direction from L to L x (1 - 6) and change the x coordinate proportionally: x (i) - x (i) x ( 1 - 6 ) i - 1, 2, . . . N Then relax the system and calculate the final energy to extract the value of C11. For CNH' also take the x direction for example, change the shape of the box and change the coordinates proportionally: SH y (i) - y(i) + x(i) x 6 i - 1, 2, . . . N Then relax the energy and get Cu”. For the bulk modulus we change the box lenghth along all three directions from L to L x (1 - 6) and change the coordinates accordingly: x (i) - x (i) x (1 - 6) y(i)-y(i)x(1—6) z (i) - z (i) x (1 - 6) i - 1, 2, . . . N After relaxation the bulk modulus can be obtained from: V - 1/2 x B x 9 x 6 2 Small value of 6 was used to make the anharmonicity negligible, 6-10—5 was used throughout this work. For the network whose mean coordination is close to the percolation threshold the system becomes very floppy. To save the computer time we used an extrapolation technique, which is described in Appendix B. The size of our system was 8 x 8 x 8 - 512 atoms. The three samples described in section 2.3 were used and periodic boundary conditions were always carefully maintained. Figure 2.1“ shows the result for C for three different samples. 11 B/a - 0.2 is adopted here. The three curves are very close to each other. Similar results were obtained for the other elastic moduli, 85 showing that the elastic moduli depend mainly on the mean coordination instead of the details of the n2, n3, and nu distributions. Figure 2.15 shows CH, C“, and bulk modulus which are averaged over three samples for an... 0.2. All these curves go to almost zero at about - 2.11, implying that the rigidity percolation threshold must be very close to the value 2.1-I, predicted by the constraint counting argument. We have calculated the elastic moduli for various B/c ratios. Figures 2.16 shows the results for B/a equal to 0.02, 0.1, 0.5, and 1.0. They all show that rigidity percolation threshold is near to - 2.“. 86 I I ‘I 1 a G 1.0 — _ g ._ C11 8 Q 0.5 — 8 8 _ . 88 Q E! f? 0 wags” ' l I 2 2.4. ' 2.8 3 2 3 6 ‘4 (T) Figure 2.1” C11 vs mean coordination for three samples. ELASTIC MODULI 1.5 1.0 0.5 87 I I I I m H C11 A—A B A .. D ._ u D (D C] A _. a g 0 _ D 8 DD 0 6 00309 whz'iflfiol I l '—————u..x.L-—: , . 2 2.4 2.8 32 36 4 Figure 2.15 C11, Cu", and B vs mean coordination averaged over three samples, B/c=0.2. 88 .J.§_b.8 n 53 g «3.: “.07 . {\m 95.29» you soaumafiouooo .- v .6 .90 in a.» a ad. a.» cd ~.u u some m> m vow .++u . :0 win ouswwm 28228.8! . .v .6 m,n c.» «a n ad o.» v.» ad a J log...“ a Manta, 0 532.886 [v .6 0.0 v.3 ufi n ad ed v6. Nu u .1 I1 w . a N... v.. 0.. 0.. u ‘u v..- a.v.z ~55.th 'l 89 There is a recent report of the sound velocity measurement on Se1_xGex glasses by J. Y. Duquesne and G. Bellessa (Reference 2.13). They measured the sound velocity at h°K for Se1_x0ex glasses with x - 0, 0.1, 0.25, 0.4, and 1. They found (respectively): VT - 1.05; 1.2“; 1.HM 906910 (longitudinal waves); VR-3.2 in Ge (Rayleigh waves). The elastic modulus in Se, Se7SGe25, SeGOGeno (shear waves); VL - 2.0; 2.10 in Se, Se C1‘“ is pvi (p is the volumic mass). :3 IE“) 1(:) qfi‘ H8 x , :5 U) 33:: Oiolg/OL 22 <: r~ :> ‘4 Figure 2.17 Cu” vs mean coordination from experiment Squares are data from the experiment, the full line is from our calculation (Reference 2.13). The Figure 2.17, which is transplanted from their paper, shows the experimental variation of can as function of the mean coordination number - 2 x + 2 of Se1_xGex glasses (in the case of Se Ge they 90 10 9O assume: V - 1/2 V T L’ in the case of Ge, V - 1.11 V T is assumed, R Reference 2.13). The overall agreement between the experiment and our calculations is good. Unfortunately there were not any data available for the glasses with mean coordination between 3 and ‘1. The elastic properties of the glasses with mean coordination in this region should be best characterized by the strong covalent forces. Below . 2.”, the elastic modulus is not zero, this is because the dihedral angle forces and interchain interactions always exist in the real glasses but are not considered in our calculations. 91 Section 2.6 Mean Field Behavior Because of the complexity of the system, we have not been able to work out a proper theory, like Effective Medium Theory for example adopted in a simpler cases (Reference 2.9 and 2.13). However, we have analised our data and tried to extract the main features of the behavior of the elastic moduli. We did not try to probe the region which is very near the critical point, because the rigidity percolation threshold (class two problem) is not known exactly and the data in the region which is very close'to the critical point is not highly reliable due to the existence of very floppy regions which give rise to relaxation difficulties. Finite size effects are also a problem in this region. We have analised the intermediate region, for between 2.11 and 3.2, and found the data fits a power law fairly well. We can see this by plotting the data on a log - log graph. Figure 2.18-2.21 (a) are log-log curves for various B/a values; they all seem to be quite linear. We have tried to extract the exponent from log - log curve, and get the prefactor by using a least-square fit. Figures 2.18 - 2.21 (b) Show the fit of the power law to our data. Table 2.1 is the list of the exponent and the best power law fit for various B/a ratios and various elastic moduli (Has1 is taken here, a/J3 is nearest neighbor distance). 92 a . LN-LN (my: Bum-0.1 -1 - .2 .1 —3 q 4 - -5 a -6 I I I I I —3 -2.5 -2 -1.5 -1 -o.s a LNtR-2.4) (a) Log - Log curve of elastic moduli vs mean coordination, B/a-0.1. annnlpwmianxmanafe 1.4 b - a—a 1.2 Bug E2 e-uo B 1 - :: 0.398*(r-2.4)**1,5 j» ___,n B.288*(r-2.4)M1,5 l' B.22$*(r-2.4)m1,5 , 2 2J! aha 215 2J3 3 :12 3H4 2L6 SJ! 4 (Y) (b) Power law fit of elastic moduli, B/c a 0.1. Figure 2.18 93 a u+u1cuwzrrnvqjhoa -1 _ H men) " A’ .2 _, ,fl _3 -. _4 -4 -5 - ”air out" a”" 0"": ..... a" -6 ‘ 13"" -7 l l I I ‘5 - ‘4 -3 -2 -1 a UflRihd) (a) Log - Log curve of elastic moduli vs mean coordination, B/u a 0.2. 2 .___aETA/nLP-a.z EXPONENT-1.5__l L75 ‘ H C11 I q G---D 1.5 9-"0 (8:44 :: B.69*(r-2.4)uc2 / 1-25 _ B.35*(r-2.4)**2 - ’,, I .. )> 1 ’/ ,1 835 _ .4 /.‘§" as - 8&5 - a ‘ I l l I I 2 2.2 2.4 2.6 2.8 3 3.2 3.4 PEP" WIMTIW -I l 3.6 3.8 4 (b) Power law fit of elastic moduli, B/a = 0.2. Figure 2.19 911 a -I ‘1 -4 .2 -I _3 - .4 -l "5 I I I I I -3 -2.5 -2 -1.5 -1 —a.s LN(R-2.4) (a) Log - Log curve of elastic moduli vs mean coordination, B/c - 0.5. 3 wear-0.5 man-1.4 ’1 2.5 _ - I ‘ ‘ C11 I “We C44 ’ O-"O B ’1 2 _ :1 1.345t($2.4)fl1.4 ,/ _"., a.715:. The other elastic moduli go to zero as soon as we cut one single bond. We can understand this in terms of the trimming process. This is equivalent to the central forces only model. When one bond is removed, the two sites on the both ends become three coordinated. We can always adjust their positions to make all the bonds have their equilibrium length. We therefore can trim them all. Then the sites connected to those bonds which we Just trimmed away become three coordinated and can be trimmed again. This chain reaction makes the whole structure fall apart, so that the removal of a single bond makes the whole network floppy. The value of B/a - 0.02 is very close to this limiting case, (we from Figure 2.16), Cij can not fit can see this sudden drop of the C 13 98 the power law at all in the region we mentioned above. For very small value of 8/0: (nonzero), the removal of first few bonds causes a large drop of the elastic moduli C and bulk modulus. Thereafter they 11 decrease very slowly to zero at - 2.”. Therefore in the whole region from u - e to 2.“, curves of elastic moduli vs mean coordination are very flat, this correspons to very large exponent (0 in the limit). For larger value of B/a, (B/a - 0.1 or 0.2), the first drop is less dramatic, that will leave more dropping to occur when the mean coordination is further decreased. This gives rise to smaller exponent. This is consistent with the general trend that the exponent is decreasing via an increase in the value of B/a. 99 Section 2.7 Effective c and B S. Feng, M. Thorpe and E. Garboczi developed an Effective Medium Theory for the central forces only model and found a good agreement with their computer simulation results (Reference 2.9). The Keating's potential we have used is far more complicated and we have not been able to work out a similar theory to understand the simulation results. But the data we get clearly show that there are effective force constants like 01* and 8* which govern the behavior of elastic moduli in analogy with the 01* of the elastic moduli of the central forces only model. In the nearest neighbor central forces only model, the elastic moduli can be expressed in terms of a single parameter a (force constant) for the pure crystal case. When a certain amount of bonds are removed, the elastic moduli can be expressed with the same formalism, but the a should be replaced by 01*. 01* can be expressed as a function of P (the fraction of bonds present). In our model, for the diamond lattice, we have three independent elastic moduli C11, Cu”, and C but only two force constants c and 8. 12 ’ This implies that-there is an identity that the elastic moduli must satisfy. This identity is found in Reference 2.10 to be the following: A. 2C""(C”+C‘2) =1 (271) (C11-C12)(CH+3C12) ,. 100 Q-2C44(C11+C12)/(C11-C12)(C11+3C123 1.4 _ EH" 1: Bum-2.02 ,___° saver-0.1 *__* same-2.2 1.2 . . owe BET/ALP-os BET/RLP-1.0 0.6 I l 1 I I 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 HERN COORDINRTION Figure 2.22 A vs mean coordination for various values of B/a. 101 ‘We can replace a and 8 by a* and 6* to make (2.5.2 - 2.5.“) still valid when decreases if A-1 is still right when we remove the bonds. This is indeed the case. Figure 2.22 shows the A value for different values of for various ratios of B/c. A is very close to 1 except in the critical region near. - 2.”. Now we can write formulae 2.5.2 - 2.5.4 as follows: < * + 3 3*) a C11 u a (2.7.2) 5 8* a C11 " a ( a . a) (2.7.3) ( * 3*) a + C12 " u a , (27”) and the bulk modulus can be expressed as: .5 * 4. Ba 3 °‘ 3 (2.7.5) * x The values of a and B can be formed from the computed quantities C11 and B (h a . 1 is taken): a - 8 1 (2.7.6) B = 8‘ (2.7.7) * x We can get Cu" from a and B and compare this directly with the simulation results. The two quantities agree very well for all different 102 ratios of B/a. Figures 2.23 shows such agreement and the behavior of the effective a and B as a function of . We found that when is near the critical region, a fixed point for the value of B*/ 0* is reached no matter what initial value of B/a was chosen to be. This universal ratio is close to 0.5. Figure 2.211 shows this universal critical behavior quite clearly. 103 .uHm 2.0 com coaumcfivuooo Emma m> cm mam yo m~.N 053.5 0 3353585.! ed v.~ ~.~ ~ v ad a...” In ~.n n O.N 0.N '.N N.N N 2853.5 5.... nu; \\ and a std ”$480.5 :9: «M HIM «nu growth . . . .1. 0.-.. m u \ HERE”. EB \\ std 9450.5 . a \. “ESE ad "Han . u a \ . \ u m... a _\ \ x r v.“ ll. 33. ~fl-r{ ~ Em... an 99.5259 B~§§..§ .. v ad Rn v.» «d n ¢.~ o.~ in. ~.~ u v Rn ad in u.» n 9.... u.~ tu ~.~ ~ Bums 30 .11... Sum .- «Ed 928.5 5E 30 9.... 98-28983 ”stout... 9...... 2.20.3755. ”Scout“. I . a... 8.... E at ad ”ESE Weigh! ll . u 6&3 ~.m...1 . 8.8 8.6 8.2 1011 —8‘FECTIVE BUR xEn-‘ECTIVE ALP e a I / 3“ £1" BETA/FILP-B.Bz //° +-_o Egg/Mga.1 / *__* mama-2.2 //° o_*,BEW%¢LP4LS lzx' BETA/nLP-1.a _ r. \' 2.2 2.4 2.5 ,2.8 3 3.2 3.4 3.5 FERN COORDINRTION Figure 2.2“ 1: x B / a vs mean coordination. 105 Section 2.8 Flow Diagram It was proposed that k/u (k is bulk modulus, 11 is shear modulus) would approach a universal value at the depletion transition. This ratio flows to a value that does not depend on the choice of the initial parameters specifying a given model. But this fixed point depends on the microscopic structure of the model (Reference 2.6, 2.1“ and 2.15). For the nearest neighbor central forces only model, there is only one elastic force constant, effective medium theory works very well all the way towards the critical point, the elastic moduli can be expressed in terms of single parameter (spring force constant a for pure system, 0* effective force constant for. random system). It is obvious that the ratios of elastic moduli would be a constant value all the way to the critical region (E. J. Garboczi, Thesis or Reference 2.9). E. J. Garboczi and M. F. Thorpe have also studied this fixed point problem for first and second neighbor central force model (Reference 2.111). They studied the 20 square lattice with the bonds cut at random using the potential: 2 19+ V - 1 (2.8.1) NI—‘ a 2 ( li - l .l 2 1 1 2 1 E ( lj 12) PJ Here 11 and 12 are are desired first and second neighbor distances, respectively. a, Y are the force constants for first and second neighbor, respectively (Figure 2.25). _ { 1 if bond 1 present - 0 otherwise 106 Figure 2.25 Diagram of square net with first and second neighbor springs (from E. Garboczi, Thesis). They used P1 and P2 as independent probabilities of the first and second neighbor bond present, respectively. They applied both effective medium theory and numerical simulations to calculate the elastic moduli. They agree very well and give the fixed point of C11/Can near the critical region. This ratio does not depend on the choice of value Y/u , but only depends on the relative weightbf cutting first and second neighbor bonds. 107 D. J. Bergmann (Reference 2.6) has studied the fixed point probleul for 2D honeycomb lattice with potential (2.8.2): V . NI-II 2 1 2 kE(6bi) +-2-mE(6¢i) (2.8.2) He used the trasfer matrix approach to arrive at the fixed point at the critical point: = 3.5 + 0.2 M. F. Thorpe and P. N. Sen (Reference 2.15) have studied a continuum model of elliptical holes cut randomly in an isotropic elastic continuum. They found the fixed point as: __; . ______.. (2.8.4) Here a, b are semi-axes of the ecllipses. The fixed point behavior seems to be very general in the elasticity percolation problems. We have studied the fixed point problem on our model. Although effective medium theory has not been able to apply to this particular model. From the result shown on section 2.7 that the elastic moduli can be expressed in terms of 01* and 8* using the expression for the pure system. It was also shown that the value of ratio c*/B* goes to a universal value at the critical point. It is obvious from above that the ratios of elastic moduli must also flow to a fixed point in the critical 108 region. Figure 2.26 shows these flow diagrams. These fixed points are approximately: We have also studied the bond-diluted 2D honeycomb lattice, with the Keating-like potential: 1 a 2 2 2 V - ._ Z [-——— Z ( x - a ) P . 2 2 Za2 1 ii 11 B +.+ _1_22 + 2 2 2.( r11 P23 + 2 a ) Piipzj] (2'8‘3) a 1.3 Here a is desired nearest neighbor distance. P21 = 1 when the bond ii present and 0 otherwise. This is same percolation problem studied by B. J. Bergman (Reference 2.6) except that he used the potential (2.8.2), we use potential (2.8.3). 109 41 11C“ mag 14 18.. -_ 9'..- 8 .4..-" " .."=" “gag-cs1; 8 I F I I 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 PER! MIMTIM (a) C11 / Cu“ vs mean coordination, 3D. 14 3444 "-512 ' I 12 - nHuo awn/Moe 02 .mo sown-0.1 19 « .1»-.. enema-a 2 ,mo saw-e s BETWFLP'LB a .1 6 .1 4 - 2 - a I I 1 i I I T I T 2 2.2 2.4 2.6 2.9 3 3.2 3.4 3.6 3.8 4 I‘m MMTIOI (b) B / Cm4 vs mean coordination 3D. Figure 2.26 110 We first calculte the elastic moduli for P-1 (with all bonds present). We use the method by Keating in calculating the elastic moduli for diamond lattice (Reference 2.11). Results are as follows: 1 M 02 + 13 c 8 + 82 1 9 c 8 CM ’ 1 J3 2 a + 3 (2.83.6) 2 2 _ 1 H a - 5 a B + 8 C12 ’ 11 I3— 2 a + B (2’8'7) We can get bulk modulus B either from direct calculation of the energy change when a hydrostatic pressure is imposed or from relation (2.8.8): B - ( C + C ) (2.8.8) 11 12 NI—b Both give the expression: 1 Balm-(Zai’fi) (2.8.9) Since this is an elastically isotropic lattice, the relationship (2.8.10) must be satisfied: - C + C - 0 (2.8.10) This can be easily checked using equations (2.8.5-2.8.7). 111 Then we cut the bond at random and impose the proper external strain for desired elastic moduli, relax the system and calculate the energy difference and extract the elastic moduli. Since this is sufficient forces model for 2D network, it belongs to class 1 problem, i. e. the elasticity percolation threshold is equal tc> the geometric percolation threshold. For honeycomb lattice this is exact: Pc = 0.6473 (2.8.11) We chose an 11148 atom honeycomb lattice and maintain the periodic boundary when we cut the bonds. We still use trimming process to remove the part of the network which does not resist any external forces. Unlike the class two problem, here all parts of structure which are geometrically connected are elastically connected also. All we can trim are completely isolated clusters and dangling clusters which are connected to the rigid backbond by a single bond. The single dangling bond can not be trimmed due to the existence of the angle bending forces. We have developed an algorithm to identify and trim the isolated or dangling clusters. Figure 2.27 shows one example of this trimming process, here P - 0.69u9, (a) is the network before the trimming, (b) is the network of the bonds present after the trimming. We calculate two independent elastic moduli, B and Cu“. The C can 11 be obtained from: + B (2.8.11) 112 (2.8.11) can be obtained from (2.8.8) and (2.8.10) We calculate these two elastic moduli for various value of ratio of 8/0. We averaged over 12 configurations to get final results. Figure 2.28 shows the results for two C“ and B. The percolation threshold is (r) . 3 x 0.6”73 ~ 1.9“. Figure 2.29 shows C11 / C1111 vs mean coordination, this clearly approaching a fixed point near the critical region inspite of the very different values to start with due to the choices of the value of B/a. This suggested the the fixed value to be approximately: C 11 can ~ 3 This is similar to the value ( - 3.5 ) obtained by D. J. Bergman with bond streching and bond bending forces model. 113 \\I I I I I I I I I'I (b) P=0.695, bond present after trimming. Figure 2.27 1111 “-5 4 mm "-443 an: wank I “-4 ' 0H“ sum-1.0 .m‘.‘ cam-9.2 scrap-0.: I 0.3 . 0.2 ' 0.1 'l a l 1.6 1.3 (a) B vs mean coordination of 20 system. 8.5 4 Wm I‘d-448 CL I 8.4 " are [Tm-1.9 .___. Kim-8.2 Elm-9.1 8.3 .. 8.2 ‘ 0.: "I 8 l 1.6 1.8 (b) can vs mean coordination for 2D system. Figure 2.28 115 2D HONEYCONB N-448 C11/C44 § ~\.4 12 18 ‘ aHufl stray-1.2 ' °_-_° BET/M'B.Z ' . h". BET/PLP-BJ ’1 8 - BET/RLP'B.GS ( - 1’ 9 I I I * 1.6 . 1.8 Figure 2.29 1/ Cu” vs Flow diagram of 20 system, C1 mean coordination. 116 Section 2.9 Introduction of Dihedral Angle Forces As was mentioned before, the problem we have studied belongs to class two, that is the rigidity percolation threshold is different from the regular geometric percolation threshold. This is because of the free choice of the dihedral angles, 1. e. we can twist a chain with no cost in energy. What happens if we include the dihedral angle forces in addition to the bond stretching and angle bending forces? All chain-like structures will become rigid due to the extra constraints. There will not be any continuous deformation that costs no energy. We have studied the above model by computer simulation. We introduce the dihedral angle forces in a way similar to how Keating did for the angle bending forces, in the form of a dot product. The potential is as follows: 3 a 2 2 2 v=————-—-2 (r -d > p 16 d2 1,1 ii ii 3 5 2 + + 1 2 2 +_____ (r! or.+—d) P. 8 d2 1,1,1 13 £1 3 ii 13 (DJ) + §.Yn Z ( F . 3 - d2 C080 )2 P P P (2 9 1) 8 d2 l,i,j,k kJ ii n £1 13 jk . n - 1, 2; n - 1 for angle 91: 60°, n . 2 for angle 0 = 180, ‘Y, and 2 Y2 should be different in general, but we will take them same for simplicity. 117 ‘The third term corresponds to dihedral angle forces, where £,i,j,k are connected by a path of bonds. This term is non-zero only if the all bonds 1J1 between are present. Here P - 1 when the 1,3 are bonded, and 13 zero otherwise. Introducing a general external strain, we can solve for the internal displacement u', v', and w' between the two FCC sublattices by minimizing potential energy (2.9.1). Then by substituting u', v', and w' into expression (2.9.1), we can extract the independent elastic moduli C11, C12, and C1‘” for diamond lattice. They are: a + 3 8 + 9 Y 011 “a (2.9.2) o - B + Y C .. 0(1“5)2+B(1+€)2+2Y(1+E)2+3Y(1-€)2 N” M a (2.9.4) a - 8 + 0.5 Y where E c + B + 2.5 Y The bulk modulus becomes: 1 1 B -‘§ ( C11 + 2 C12) ' T23 ( 3 G + B + 11 Y ) (2.9.5) We calculated the elastic moduli in a same way we did before for the network without the dihedral angle forces. Three same samples were used to see the change in elastic moduli solely due to the introduction 118 of the dihedral angle forces. Averaging was taken over three samples. We did calculation for Y - 0.0” and 0.1 respectively. Figure 2.30-2.32 show our data in comparison with previous data for Y-O, B/c - 0.2, B/as0.2 is also used here. Figures 2.30-2.32 are just the blow up of the part(a) of the Figures 2.30-2.32. We can see that the introduction of the dihedral angle forces really washes out the transition at -2.“. Elastic moduli are not zero at -2.11. This is because of the additional contraints introduced by the dihedral angle forces. But from the Figures 2.30-2.32 we can see that the elastic moduli are very small when is near two, this seems to suggest that the elastic moduli may vanish at . 2. In the limiting case (r) - 2, the whole system will consistIof disconnected chains in 3D space- Yacov Kantor and Itzhak Webman have studied the elastic behavior of a single 2D chain formed by a set of N vectors (Reference 2.5). They fix the one end of the chain, exert a force on the other end, and minimize a function W: W-H-F-(R -§) (2.9.6) 1 Where RN - RN is displacement of the end point which is not fixed. H is the potential: (a) C1 8.25 8.2 8. 15 119 1“— C11 "-512 am even 3 m , I 3.2 3.4 3.6 3.8 4 1 vs mean coordination for Y a 0, 0.0”, 0.1. —— C11 "-512 macs: OVER 3 SR‘PI_1"_°~ (b) Blow up of (a). Figure 2.30. 120 1‘ c“ "'5‘? “VF-RAGE OVER 3 snags 2 2.2 20‘ 2.6 208 3 ' 302 .30‘ 3.6 308 l ‘ (a) Cu” vs mean coordination for Y - 0, 0.0", 0.1. (b) Blow up of (a). Figure 2.31. 121 (a) B vs mean coordination for Y a 0, 5H12movmasans 22 24 :26 as 3 :32 34 nnqaumnnnar. .82.. s 105:2 m och 3 my; 1 '3.6 I . 3.8 ‘ 0.0”, (b) Blow up of (a). Figure 2.32. 0.1. 122 N N H - “2 2 6b? +2 2 Mi (2.9.7) 2 a i=1 1 1 Where 31.nsiflth vector, oi is the i'th angle, and a is the equilibrium length of each vector. They solve for the 6b and “1 that minimize w, substitute them 1 into the H, and find the effective force constant of the system to be: 2 2 KRw B / a N (2.9.8) aw here stands for the random walk. This suggests the 1/N2 dependence of the elastic moduli for a single 2D chain, clearly it will go to zero for long chain. A.simtuu'argument can be used to study the chain made up of N vectors in 3D space with dihedral angle forces in addition to the two terms in the potential (2.9.7), which gives the behavior: K ~ Y / a N Y is the dihedral angle force constant, here we suppose that Y is smaller than 8. Therefore the elastic moduli should be very small when - 22. We may understand this in terms of constraint counting. There is one constraint for bond stretching plus one constraint for angle bending plus one constraint for dihedral angle changing per atom in the chain, and this is equal to the number‘of degrees of freedom. The only 123 constraint which prohibits the system from relaxing to the zero energy configuration is the boundary conditions. We have actually done some computations on 30 random chains with dihedral angle forces and found the elastic modulus to be very small indeed. For 3D chains as short as N-20 will make the elastic modulus less than 1 percent of the Y value. This result confirms the above conjecture. From above constraint counting we can also see that the N bond lengths, N angles and N dihedral angles, these 3N variables can always decide the whole structure of a 3D single chain made of N vectors. Any form of inner chain interactions between any further neighbors will not give rise to additional independent constraints. To make the chains rigid in 30 networks the inter chain interactions are needed. In the real polymeric glasses for - 2 (Se for example), the ielastic moduli are not zero. Figure 2.17 displays the result of elastic modulus Cu" from sound velocity measurement. C1M of Se glasses is small but not zero. This indicates that the inter chain interactions certainly exist in the real glasses. These interactions may be due to the van der Waals interactions. Here we should also notice that this does not correspond to the 3D bond percolation where the conductivity vanishes at about - 1.5. The difference is that we restrict the minimum coordination to be two and do not allow dangling bonds or isolated sites. But for the bond percolation problem we cut bonds at random with no restrictbon, so even at =- 1.5 on the average a long thin chain may still percolate to make a nonzero contribution to the conductivity. 1214 Section 2.10 Conclusions It is believed that covalent glasses can be divided into two classes: those with high average coordination (amorphous solids) and those with low average coordination (polymeric glasses). We give the first conclusive evidence through this research that this is correct. We have calculated the elastic properties of random networks with different average coordinations . Our results show that the elastic moduli depend predominately on (r) and go to zero at - 2.” when the dihedral angle forces are not included. We have also calculated the fraction of zero frequency modes for random networks with different average coordinations. Our results show that the fraction of zero frequency modes also depends predominately cu: the average coordination and the results show a very good agreement with constraint counting prediction. The results of calculated elastic moduli and the fraction of zero frequency modes suggest that there is a phase transition between two classes of'covalent glasses. amorphous solids and the polymeric glasses at a mean coordination close to 2.14, which is predicted by the constraint counting argument. The introduction of the dihedral angle forces which give additional constraint to the random network washes out the transion at - 2.14. Interchain interactions are needed to stablize the random network formed by isolated chains. Apppendix A Relaxation Process We have used a local relaxation process to bring a network to a minimum elastic energy state. We move only one site to a local minimum of strain energy at a time with every other sites fixed. We do many cycles to relax the system, here by one cycle we mean that every site has been moved once. For a fairly rigid system as few as 50 cycles are enough to bring the system to a minimum energy state. But for systems with floppy regions, a few thousand cycles may not be enough to find the minimum energy state. At this point we use an extrapolation technique to save some relaxation steps (Appendix B). The individual movement used to bring the given site to local minimum is decided by following steps. First expand the potential, keeping only up to quadratic terms. The energy can be written as follows: 2 “(9.)” v0 + 2 ( 3" )oérm + 2 (—9-——Y-—-) arm 5r“ (3-1) m m n o m 3r 2 m,n argarl Here a, B are the x, y, 2 directions. 2. is a given site, and we take V -0. o The displacement needed to bring site 9. to a local minimum is determined by the equation: ]25 126 3 V m 3 V m n _ E(arm)°6rl+mzn( armarn)°6r9' 6r£=0 (82) i ’ 2 2 3 V 32V n ( 3 rm )0 + E ( arm 3rn )0 6 r1 : 0 (8-3) 2 2 2 ( m ’ 19293 ) From these three equations ( m-1,2,3 ) we can solve for Grn ( n - x,y;z ), then we take step x + x + 6x 6x, 6y, 6z are solutions from equations ( B-3 ). For potential (2.3.2), first and second derivatives are: 3 V _ 1 g 2 _ 2 m m m ' 2 2 Z ( r21 d ) ’21 P21 3 r d i 2 __3_£ *." 1.2 m m u d2 123 (rm1 r9,J + 3 d ) ( r£J+ r21) Plipij i>J -ifi *." 12‘“ u (12 iEJ ( r21 r31 + 3 d ) r31 Plipij 127 2 3 V 2..E_ m n 2 _ 2 m n ' 2 2 Z [ 2 r21 r21 + 5 ( ”21 d ) P21] Br Br d i 2 2 + ;-E- 2 [(rn + rn )( rm rm ) + 26 (F - r d2)] P u (121‘1 9.3 9.3 9.1 mm 9.1 9.3 3 9.1 9.3 i>J _3__§_.. + u d2 H2 rm 13 P11 PiJ j#2 For the potential (2.9.1), first and second derivatives are: 3 V _ g ‘_g_ 2 2 m '""E ' 2 2 Z (r21 d ) r21 P11 3r d i 2 - §._§_ * + .1 2 m m u 2 2 ( ”13 ”11+ 3 d ) (r23+ r21) P2ip2j d 1.J 1>3 — §._1_ 2 ( F .3 d2 0039 ) rm 4 d2 i,J,k kJ 21 kj H21 13 ij -—-—- r r +—d )r P 3 d2 1,J ij i2 ji 2J Jfi2 3 B Z " + 2 + — ——— ( r - r - d COSG ) rm P P P A d2 i,J,k 12 kJ kJ 2i 23 Jk 32 V - 3—3—9— 2 [ 2 rm r“ + 5 ( r2 — d2)] P m n 2 2 21 21 mn 21 21 Br Br d . 2 2 ;‘_g_ m m n n . l 2 + u d2 iZJ[(r23+ r21) (rlj+ r11) + 2 5mn(r2i ”23+ 3 d )] Pzipij i>j ‘ 128 m n g rkj rkj P21 P13 ij :13: r r .P P d i,j,k kJ kj 9.1 9.3 Jk Appendix B Extrapolation of Elastic Energy Relaxation Elastic energy relaxation is a very fast converging process for a perfectly rigid system. But this converging is very slow for a system with floppy inclusions when the mean coordination is very close to the rigidity percolation threshold. In order to save relaxation steps, we have studied the behavior of elastic energy for a long term relaxation process. We found the asymptotic form of elastic moduli near the final convergent value can be written as: - n C - C + a b (B-H) Where constant co is final value of elatic modulus C, n is the number of large relaxation steps. By large relaxation step we usually mean 50 relaxations per site. a > 0 . b > 1 . a, b are constant. If N of these large relaxation steps are used at in relaxing a system to the minimum energy state. Then by fitting the exponential form (B-h) to last three steps, we can solve Co from Cn-Z , Cn—1’ and Cn from: I C C - C C _ n n 2 n 1 o c - 2 cn_ + cn_ ]29 130 Co will only be an upper bound for the value of elastic modulus Ch In effect. extrapolation is just a way to save some relaxation steps. Table 2.2 gives a example where the extrapolated modulus can be compared (x) fo with actual further relaxation steps. This data is taken from C11 r Elm-0.2 and - 2.u3. Table 2.2 Extrapolation vs. Relaxation Relaxation steps Before Relaxation Value After 2000 0.00558 0.00528 ”000 0.00528 0.00513 1.1 1.2 1.3 1.“ 1.8 1.9 1.10 References (Chapter 1) W. H. Zachariasen, J. Am. Chem. Soc. 25, 38M1 (1932). D. E. Polk, J. Non-Cryst. Solids g , 365 (1971). P. Steinhardt, R. Alben and D. Weaire, J. Non-Cryst. Solids 12, 199 (197“). P. W. Keating, Phys. Rev. 1&2 , 637 (1966). G. A. N. Connel and R. J. Temkin, Phys. Rev. B‘g, 5323 (197A). N. J. Shevchik Phys. Rev. Lett. 31, 12H5 (1973). G. Etherington, A. C. Wright, J. T. Wenzel, J. C. Dore, J. H. Clarke, and R. N. Sinclair, J. Non-Cryst. Solids fig, 265 (1982). N D. Beeman and B. L. Bobbs, Phys. Rev. B 12, 1399 (1975). W. A. Paul and G. A. N. Connell in "Physics of Structurally Disordered Solids" edited by S. S. Mitra (Plenum, New York, 1976) P. #5. D. Henderson and F. Herman, J. Non-Cryst. Solids §;19, 3S9 (1972)f L. Guttman, Bul. Am. Phys. Soc. 22, 6M (1977). L. Guttman in "Tetrahedrally Bonded Amorphous Semiconductor" editted by M. H. Brodsky, S. Kirkpatrick and D. Weaire (American Institute of Physics, New York, 197“) pzzu. F. 0rdway , Science 15;, 800 (1964). R. J. Bell and P. Dean , Philos. Mag. 22, 1381, (1972). R. Zallen, The Physics of Amorphous Solids, J. Wiley and Sons, ]3] 1.20 1.21 1.25 132 New York (1983). W. Y. Ching, Phys. Rev. B 26 6610 (1982). F. Galeener, J. Non-Cryst. Solids '52, 53 (1982). A. Tadros, M. A. Klein, and G. Lucovsky, J. Non-Cryst. Solids 65, 215, (198”). D. Weaire in " Excitations in Disorded systems", edited by M. F. Thorpe (Plenum, New York, 1981) P. 579. F. Wooten and D. Weaire, J. Non-Cryst. Solids ‘65, 325 (198”). F. Wooten, K. Winter, and D. Weaire, Phs. Rev. Lett. 53 1392 (1985). F. H. Stillinger and T. A. Weber, Phy. Rev. B 31, 5262 (1985). S. W. Deleeuw, H. He and M. F. Thorpe, to be published. J. A. Erwin, A. C. Wright, Joe Wong, and R. N. Sinclair, J. Non-Cryst. Solids 51, 57, (1982). S. W. Deleeuw, and M. F. Thorpe, to be published. 2.1 2.3 2.“ 2.5 2.6 2.7 2.8 2.9 133 References (Chapter 2) N. F. Mott, Phil. Mag. 13 , 835 (1969). R. Zallen, Physics 0f Amorphous Solids, J. Wiley and sons, New York (1983). S. Feng, P. N. Sen, B. I. Halprin, and C. J. Lobb, Phys. Rev. B‘39 .5386 (198A). S. Feng and M. Sahimi, Phys. Rev. 8.31 , 1671 (1985). Y. Kantor and I. Webman , Phys. Rev. Lett. 23 , 1891 (1984). D. J. Bergman , Phys. Rev. B 31 , 1696 (1985). J. C. Phillips, J. Non-Cryst. Solids 35 , 153 (1979) and J. Non-Cryst. Solids, 23, 37. (1981). M. F. Thorpe, J. Non-Crystalline Solids 31 , 355-370 (1983). S. Feng, M. F. Thorpe, and E. J. Garboczi, Phys. Rev. B 31, 276, (1985). H. He and M. F. Thorpe, Phys. Rev. Lett. 33, 2107 (1985). P. N. Keating, Phys. Rev. 133 , 637 (1966). C. Kittel, Introduction to Solid State Physics , Nth edition. J. Y. Duquesne and G. Bellesa, to be published. E. J. Garboczi and M. F. Thorpe, Phys. Rev.B 31, 7276 (1985). L. M. Schwartz, S. Feng M. F. Thorpe and P. N. Sen, To Be Published.