a z 3 x: , u 1%“. as.” 1 3i 1“ .m ‘7‘?! 1 q r,. 3 U I 1 " t‘.§.l.! , Ei\Nv\.t\ i, .l :1} 3.51\ sifi‘tsii . :tt}. . . ‘ .. P)», ._ .5. g .rémmq... x p {\x ‘\ 4 .5 Jr; yfiflmfifigfi, . _ . . . :Kuflfi. memwu . . ‘ $2. {a ) 7 , V v.\. 1 1.. bug 131; .11. 1 3.. .4. , «SvISIiwbvu‘ n~pNuluwi ’1 tifiVrlgs‘yth‘éiufli 7323“,.1): ‘i (’31.. 3412‘: »>\) I)". 334‘ .0}: kn EK‘ ui.¥s}yk§2(i.hilz§xl§{§ \ v 53 ax ‘ a 2134.)}! a )5! 1‘! All“ ‘1‘ inz‘ivfl. «Hung. Nun”; .. *hai_}1)§s 111%,.“5 2 n 33‘ 3. ‘ S . yltlt ) i .l: a .. . 2!; 1.151511! Iul-l’ultn‘lli {51% $31.17!)! it“s-1‘32?! . 1‘4. . iii. A . Ly??? fl m Vi; mum;TIMITITQMTTIIMTTm « rHEG'B Maw... .. 1:: state -‘-.J Umvermty This is to certify that the thesis entitled BAND STRUCTURE AND REFRACTIVE INDEX CALCULATIONS FOR SEMICONDUCTORS presented by Mark Douglas Ewbank has been accepted towards fulfillment of the requirements for Master degreein Physics % IVE“ r Major professor Date 4/é/76/ / / 0-7639 OVERDUE FINES ARE 25¢ RER DAY PER ITEM Return to book drop to remove this checkout from your record. BAND STRUCTURE AND REFRACTIVE INDEX CALCULATIONS FOR SEMICONDUCTORS By Mark Douglas Ewbank A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Physics 1979 ABSTRACT BAND STRUCTURE AND REFRACTIVE INDEX CALCULATIONS FOR SEMICONDUCTORS By Mark Douglas Ewbank A "chemical bond" approach was employed to provide a means of theoretically predicting material properties of semiconductors which can be used to delimit the performance characteristics in specific device applications. First, a set of empirical parameters used to represent the interaction between atomic orbitals was investigated by performing numerical band structure computations following the LCAO or tight binding method and utilizing the Slater-Koster 2-center integral approximation. Band structures were acquired for zinc- blendes (Si, Ge, GaAs, ZnSe), chalcopyrites (ZnSiAsZ, CdGeAsZ and CdGePZ) and a sulfosalt (Tl3AsSe3). Then Harrison's bond-orbital model (BOM) was engaged to predict the principal refractive indices of the chalc0pyrites through a similar parameterization. Both numerical calculations were implemented with "generalized" Fortran computer programs (listings included). To Laurian, Shawn, and Cathy ii ACKNOWLEDGMENTS The author would like to express his appreciation to E.A. Kraut and w.A. Harrison for the vast amounts of consultation they offered. Also, P.R. Newman deserves recognition for many fruitful suggestions concerning the subject matter of the work. Finally, the author wishes to thank J. Tracy and P.A. Schroeder for making this effort possible. TABLE OF CONTENTS Page LIST OF FIGURESO0.00...0..OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOVii LIST OF TABLESOOOOOOOOOOOOOOOOOO0.00.00.00.00.0.0.0000... X1. INTRODUCTION.......... ................................... l ENERGY BAND CALCULATIONS... .............. .... ..... . ...... 3 FORMULISM. O O O O O O O O 0 O O O O O O O O O O O O O O 0 O O O O O O O O O O O O O O O O O O O 3 IMPLEMENTATION. O O O O 0 O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O 9 RESULTS. 0 O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O 0 O O O O O O 11 Zincblendes...................................... 11 Chalc0pyrites.................................... 25 T]3Asse300000000000OOOOOOOOOOOOOOOOOO ........... O 28 TlgAsSe3 with Group Theory ....................... 28 Symmetry Considerations... ................ ....... 37 CONCLUSIONOOOOOO00.0.0000...O...OOOOOOOOOOOOOOOOOOOOO 40 DENSITY OF STATES CALCULATIONS........ ..... .............. 41 RESULTSOOOOOOOOOOOOOOOOO00.000.00.000.000000000000000 44 SUMMARY.OOOOOOOOOOOOOOOOOOOOOO0.000000000000000000... 44 iv TABLE OF CONTENTS (Cont'd) Page INDICES OF REFRACTION AND BIREFRINGENCE CALCULATIONSOOO0.00000000000.000000000000000000000.... 49 FORMULISM........................................... 49 IMPLEMENTATION....... .......... ...... ............... 56 RESULTS FOR CHALCOPYRITES........................... 58 SYNOPSIS............................................ 62 CONCLUSION .............. . .......... . ...... ..... ......... 65 APPENDICES A. ATOMIC ORBITALS.. ...... ........ ....... . ........ ..... 67 B. SLATER AND KOSTER 2-CENTER INTEGRALS OF ATOMIC ORBITALS...OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 73 MEANING OF Z-CENTER INTEGRALS... ................ 73 COMPONENT NATURE OF ATOMIC ORBITALS............. 74 INTERACTION BETWEEN TWO ATOMIC ORBITALS......... 78 C. ENERGY BAND COMPUTER PROGRAM........................ 87 D. DENSITY OF STATES COMPUTER PROGRAM................. 107 TABLE OF CONTENTS (Cont‘d) E. PRINCIPAL INDICES OF REFRACTION COMPUTER PROGRAM...................... .................... F. CHALCOPYRITE CRYSTAL STRUCTURE....... ..... ......... Go TI3ASS€3 CRYSTAL STRUCTUREOOOOOO00.00.00.00000.0000 H. SUBROUTINE FOR CHADI AND COHEN'S FIRST NEAREST NEIGHBOR PARAMETERS...’OOOOOOOOOOOOOOO0.0.0....O. I. SIMPLIFICATION OF THE FIRST ORDER BOND POLARIZATION FORMULA............................. REFERENCESOOO0.0000000000000000 OOOOOOOOOOOOOOOOOOOOOOOO vi Page 118 132 140 145 148 150 LIST OF FIGURES Figure Page 1 Identical interaction between pairs of atomic orbitals from 2-D lattice periOdiCitYOOO0.0.0000000000000000000000000.0.0....O... 8 2 Energy bands of GaAs using Chadi and Cohen's first nearest neighbor parameters; (Top) from r to L, (Bottom) from r to X.......................... 15 3 Energy bands of ZnSe using Chadi and Cohen's first nearest neighbor parameters; (T0p) from F to L, (BOttom) fromr to XOOOOOOOOOOOOOOOOOOOOO00000 16 4 Energy bands of Si using Chadi and Cohen's first nearest neighbor parameters; (Top) from r to L, (Bottom) from r to X.......................... 19 5 Energy bands of Ge using Chadi and Cohen's first nearest neighbor parameters; (T0p) from r to L, (Bottom) from r to X..................... ..... 20 6 Energy bands of GaAs using Harrison's parameters for first nearest neighbor parameters; (Top) from r to L, (Bottom) from r to X............... ...... 21 7 Energy bands of ZnSe using Harrison's parameters for first nearest neighbor parameters; (Top) from 11 to L, (Bottom) from r to XOOOOOOOOOOOOOOOOOOOOO 22 8 Energy bands of Si using Harrison's parameters for first nearest neighbor parameters; (Top) from r to L, (Bottom) from r to X..................... 23 9 Energy bands of Ge using Harrison's parameters for first nearest neighbor parameters; (Top) from r to L, (Bottom) from r to X............ ...... ... 24 10 Energy bands of GaAs using Harrison's parameters for second nearest neighbor parameters; (Top) from r to L, (Bottom) from r to X..................... 26 vii LIST OF FIGURES (Cont'd) Figure 11 12 13 14 15 16 17 18 19 20 21 Energy bands of GaAs using Harrison's parameters for third nearest neighbor parameters; (Tap) from r to L, (Bottom) from r to X..................... Energy bands of ZnSiAsz using Harrison's para- meters for first nearest neighbors.................... Energy bands of CdGeAsz using Harrison's para- meters for first nearest neighbors.................... Energy bands of CdGePz using Harrison's para- meters for first nearest neighbors.................... Energy bands of Tl3AsSe3 using Harrison's para- meters for first nearest neighbors.................... Atomic orbitals with triangular symmetry................ Energy bands of two irreducible represent- tations of Tl3AsSe3 using Harrison's para- meters for first nearest neighbors; (Top) A1, (BOttOfll) AZOOOOOOOOOOOOOO0.00000000000000000000000 Energy bands of Si along c-axis; (Tap) one degree of symmetry omitted, (Bottom) all symmetry included..................................... Wave—vector, k, which falls outside a cubic Brillduin zone is related to wave-vector, k, within the Brillouin zone by reciprocal lattice vector, G..................................... Density of states of Si for first nearest neighbors; (Top) Chadi and Cohen parameters, (Bottom) Harrison's parameters........................ Density of states of Ge for first nearest neighbors; (T0p) Chadi and Cohen parameters, (Bottom) Harrison's parameters........................ viii Page 27 29 3O 31 32 34 38 39 43 45 46 LIST OF FIGURES (Cont'd) Figure 22 23 24 A1 Bl B2 B3 B4 B5 B6 C1 C2 D1 D2 Density of states of GaAs for first nearest neighbors; (Top) Chadi and Cohen parameters, (Bottom) Harrison's parameters........................ Density of states of ZnSe for first nearest neighbors; (Top) Chadi and Cohen parameters, (Bottom) Harrison's parameters........................ Theoretical birefringence of CdGeAsZ as a function of lattice constant, x................................ Directional nature of atomic orbitals................... Two dimensional transformation of a p orbita]0.0.00..0.000000000000000000000000 ..... O. 00000 0 Three dimensional transformation of a p orbita]0.0.0.0....COOOOOOOOOOOOOOOOOOOOO0.0.0....0.... Slater and Koster's matrix element symbols.............. Transfer integrals between atomic orbitals which are zero.............. ..... .. ....... . ........... Interaction between an s orbital and a randomly oriented p orbital........................... Interaction between two randomly oriented p orbitaISooooooooooooooooooocoo-0000.00.00.00 00000 0.. Flowchart for energy band program....................... Computer listings for energy band program............... Flowchart for density of states program................ Computer listings for density of states program.............................................. ix Page 47 48 63 72 75 77 79 8O 82 83 88 89 109 LIST OF FIGURES (Cont'd) Figure Page E1 Flowchart for principal indices of refraction program...OOOOOOOOOOOOOOOOOOOO0.00.000000000000000... 119 E2 Computer listings for principal indices of refraction computer program........... ............... 120 F1 Computer listings for chalcopyrite crystal structure program.................................... 135 G]. CT'YStaI StY'UCtUY‘e 0f T13Asse3ooooooooooo0000000000ooooo 144 H1 Computer listing of subroutine for Chadi and Cohen's first nearest neighbor parameters............ 146 Table A1 A2 Bl F1 G1 G2 LIST OF TABLES Selected Energy Eigenvalues for GaAs Using First NeareSt NeighborSOOOOOOOOOOOO00..000......000000.00... Selected Energy Eigenvalues for ZnSe Using First Nearest NeighborSOOOO...00.00.000.000.0000000000000000 Selected Energy Eigenvalues for Si. ................... Selected Energy Eigenvalues for Ge .............. ....... Indices of Refraction and Birefringence for Various Chalcopyrites--A Comparison of Theory and Experiment.0.00.0000...O...0.0.0.000...OOOOOOOOOOOOOO. Theoretical Ordinary Index of Refraction and Bire- fringence for Various Chalcopyrites................... Spherical Harmonics... ..... . ............................ Atomic orbita] FunctionSOOOOO0......OOOOOOOOOOOOOOOOOOO. Correlation of Harrison and Slater-Koster Notation...... List of Atom Positions for the Chalc0pyrite, CdGeASZOOO...0.0.0000....OOOOOOOOOOOOOOO0.0.0.0 000000 Lattice Parameters for Tl3AsSe3. ............ . .......... List of Atom Positions for the Sulfosalt, T13Asse3oooooooo00000000000000.0000oooooooooo....... xi Page 13 14 17 18 59 61 7O 71 86 139 141 143 INTRODUCTION Sophisticated space surveillance requires a "solid state" elec- tronically tunable filter which functions in the infrared region of the electromagnetic spectrum. In other words, the optical filter should exhibit, for example, a spectral passband of approximately ten nanometers about a central wavelength that is remotely tunable from three to five or eight to twelve microns. Two possible devices which can fulfill these types of requirements are Acousto-Optical Tunable Filters (AOTF) and Electra-Optic Tunable Filters (EOTF). Also, future applications may arise in the visible and even ultraviolet portions of the Spectrum which could be satisfied by similar devices. There is a variety of physical problems to resolve in the de- velOpment of such devices--crystal growth, material pr0perties char- acterizations, acoustic wave pr0pagation, infrared imaging through an optically anisotropic medium, etc. The primary purpose of this work is to initiate the evolution of working theoretical models that are capable of predicting some pertinent material properties which are necessary for tunable filter design. In particular, by qualitatively ranking ternary crystalline semiconductors with respect to a specific physical property (e.g., birefringence, elasto-Optical coefficients or electro-Optical coefficients), the numerous candidates can be screened in an attempt at optimizing crystal selection for tunable filter development. More generally, the underlying intention of this 1 endeavor is to furnish theories that can benefit the needs of eventual, yet unforeseen, device applications. The emphasis of the present work involves predicting the princ- iple refractive indices using a "chemical bond" approach. The empirical parameters used to characterize a "bond interaction" are initially investigated, in detail, with energy band structure cal- culations. And then, subsequently, the principal indices of refrac- tion are computed by employing an equivalent set of parameters. Both numerical calculations are implemented on a computer, incorporating the concept of general applicability to any semiconductor which helps to fulfill the criterion of providing a theoretical foundation for future requirements. ENERGY BAND CALCULATIONS To test the numerical accuracies of the parameters used to des- cribe the interaction of electronic orbitals between two atoms (transfer integral), energy band calculations have been performed and compared with similar band structures that can be found in the liter- ature ad infinitum.1 The approach taken here utilizes the LCAO (Linear Combination of Atomic Orbitals) or tight binding method to determine the quantized energy states of electrons at various points in the Brillouin zone.2’3 FORMULISM First, define and specify the notation, following that of Zunger.4 Atomic orbitals (see Appendix A) are symbolized by ¢u(F-Rg), where n denotes the type of valence orbital (s, px, py, pz, etc.) and R: is the vector distance to the ath atom in the nth unit cell. Bloch functions are formed by taking linear combinations of ' "equivalent" atomic orbitals in each unit cell where the expansion coefficients are determined by the periodicity of the lattice, a #- —)- = +-+a °u (k,r) ¢ (r Rn) (1) The summation is over unit cells and approaches an infinite sum for real crystals. "Equivalent" atomic orbitals imply one type of atomic orbital on the same atom in each unit cell. Crystal wave functions or eigenfunctions are linear combinations of the above Bloch functions: wit?) = z z c . (i?) a (TC?) (2) 3 u=1 a=l where j specifies which band (j = 1,2...., M and M = on = number of bands). The summation over u is a sum over all valence atomic orb- itals and the summation over a is a sum over all atoms in the unit cell. Each wavefunction |¢j>, is a solution to the time independent Schrodinger equation, To determine the coefficients, C35(R), minimize the energy eigen- values, Ej, with respect to each coefficient. In other words, J = o (4) for each u and a, where <¢.|H|¢.> E. = 3 3 (5) J (wjle> This results in a set of M = on "Secular Equations" (see Equations 47 through 64 for an example) n 1 1 a=l "MG 3 -+ B + + a + + B -+ _ _ .k< , H k,>- . E.- , u qu( ) ¢v(k r)! |.u( r) CvJ(k) J 0 (b) one for each v = 1,....,o and B = 1,....,n. These M equations con- tain M unknowns -- the coefficients, C3j(R). If a nontrivial solu- tion exists for a system of M linear, homogeneous equations in M unknowns, then the determinant of the coefficient of those M unknowns must be zero.5 For the sake of illustration, let the double indices u and a be represented by a single index. Then the Bloch function ¢:(R,F) goes to 61(E,F) and the "Secular Determinant" can be written as L<¢1|H|¢l>-Ej] <¢1|H|¢2> <¢1|H|¢3> .... <¢1|HI¢M> (7) <¢2|Hlo1> I<¢2|HI¢2>-Ej] <¢2lHlo3> .... <¢2|HI¢M> <¢3|Hlol> [-Ej] .... = 0 <¢M|Hl¢1> <¢MIH|¢2> <¢M|H|¢3> oooo [<¢M|H|©M>-Ej] The solution of this secular determinant, which results in M energy eigenvalues, E ¢j(J = 1,2,.... Hamiltonian, H. js that correspond to each of the eigenfunctions, That is, perform a unitary transformation, H' U'lHU, which puts the Hamiltonian, H', into diagonal form. words, ,M), is obtained by "diagonalizing" the matrix of the In other [-Ej] 0 0 ... 0 (8) O [Ej] 0 ... 0 ? 9 [<@.3|H|¢l3>-Ej] 9 = 0 0 0 0 ... [<¢.M|H|¢'M>-Ej] o The roots of this new secular determinant are Ej = <4 jIH lo 3) . (9) Also, the new eigenfunctions associated with these eigenvalues, Ej, I j. transformation, the secular determinant of H and the diagonalized are wj = 6 Since the eigenvalues are invariant under unitary secular determinant of H' yield the same eigenvalues. (This diago- nalization is done numerically on a computer by a Fortran subroutine, provided by the International Mathematical and Statistical Libraries (IMSL), which requires each of the matrix elements, , as input and returns the energy eigenvalues, Ej, and Optionally the eigenfunctions, wj.) Consequently, to determine the eigenvalue energies at a given wave-vector, k, in the Brillouin zone, the matrix elements, <¢3(R,F)IHIQ€(R,F)> must be evaluated -- that is, the matrix element of the Hamiltonian between two Bloch functions of atomic orbitals. Writing out each Bloch function as an expansion of atomic orbitals from Equation 1, ¢v(?—Eg. )]> . (10) Since the Hamiltonian only Operates on functions of F, N ik°(R T-R ‘* ”*a ‘* “*8 A n n <¢u(r-Rn)lHl¢v(r-Rn.)>. (11) This double sum can be reduced to N times a single sum, as illus- trated in the following example. Consider the interaction between atomic orbitals PAS and ¢Bl in the two-dimensional lattice with nine unit cells, as shown in Figure 1. This interaction corresponds to only one term in Equation 11. But an identical interaction occurs between the pairs of atoms (A6, 82), (A8, B4), and (A9, B5). In general, for N unit cells, there are N such identical interactions, neglecting edge effects. Therefore, Equation 11 becomes 8 a B __ ) +30; +38 <¢UIH|¢v> — 1 <¢u(r Rn)|H|¢v(r R1)> (12) "N2 (D I —lo x . A m :9 I m where Rf is arbitrarily taken to imply the unit cell "at" the origin. Now recall that the sum over N unit cells is essentially an in- finite sum. To reduce this to a finite sum, consider only "near neighbors". In other words, assume that the matrix element, <¢u(F—R:)IHI¢V(F-R€)>, is negligible when the overlap is small. A Figure 1 o o o O p /O A9 B 40’ 350’ o 0 A5 A6 B o’ d” o 1 B2 0 A1 0 O Identical interaction between pairs of atomic orbitals from Z-D lattice periodicity 9 small overlap between the two atomic orbital occurs when (F—Rfi)- (FLA?) = Rfi-Rf is large. To insure that all interactions are weighted equally by distance, choose the cutoff such that only atoms (specified by n and a) inside a spherical shell of radius, Dmax’ about the central atom (specified by R?) are included. That is, the matrix element, <¢u(?-Rfi)lHl¢v(F-Rf)>, contributes to the sum in Equation 12 only when the two atoms involved are "near neighbors" or +0: +8 an - R1! < Dmax . (13) Then, <¢Q|H|¢§> = ze'1R (RnR1)<¢ (r R“ Ml'Hli (F Rf)> . (14) “ near neighbors n IMPLEMENTATION In view of Equation 14, the energy band structure can be obtained by specifying the crystal structure (i.e., knowing the position, R?) of each atom in the unit cell along with the positions, Rfi, of all "near neighbor" atoms in real coordinate space) and numerically eval- uating the matrix element of the Hamiltonian between pairs of atomic orbitals. It is these same matrix elements that enter into calcula- tions of index of refraction, birefringence, elastic constants, etc., and whose validity we wish to investigate via the band structure calculations. The implementation of Equations 7 and 14 to determine the energy eigenvalues as a function of wave-vector has been accomplished by 10 writing a Fortran program, which is flowcharted and listed in Appen— dix C. (Note: the listings are commented in detail for readability without an indepth Fortran background.) The problem of obtaining the crystal structure was handled separately -- i.e., communicate the crystal geometry via a "list" of atoms with vector coordinates gener— ated by a separate computer program (see, for example, Appendix F). The evaluation of matrix elements of the Hamiltonian between two atomic orbitals was done by using the Slater and Koster "two-center integral" approximation2 (see Appendix B). One generalized method of numerically specifying the parameters in this approximation was sug- gested by Harrison.6’7 For two 5 orbitals, an s and p orbital, and two p orbitals, re- spectively, Harrison's "Solid State Interatomic Matrix Elements" are given below. <52|H|sl> = VSS (15) /\ /\ <52|H|p1> = vSp (p1 - d) (16) A o A AoA AOA = Vppn (p1 p2) - [Vppo + Vppfl](p1 d)(p2 d) (17) 2 _ h VSS - -1.4O REF? (18) 2 _ h VSp - -1.84 2 (19) 11 h2 v = -3.24 h2 V = "0081 _ 2 ppn mdz ( l) h2 = 7.62 eV - A2 (22) m '6 is the vector distance from atom #1 to atom #2 and 61 and 6% are unit vectors describing the orientations of the positive lobes of the p orbitals. Note that the geometrical factors of the matrix elements involving 5 and p orbitals are compatible with those of Slater and Koster as shown in Appendix B. For matrix elements involving two atomic orbitals on the same atom, it is convenient to use the Herman and Skillman atomic term values for non-orthogonal orbitals.6’8 For the sake of simplicity, d orbitals were neglected in this work since, for most semiconductors of interest, the d shell is completely filled. RESULTS Zincblendes The first step in performing numerical calculations with a com— puter is to show that the computer program being used generates valid results. This was accomplished by reproducing the energy bands for GaAs and ZnSe from Chadi and Cohen,1 using their first nearest neigh- bor parameters (Table 5) for the matrix elements of the Hamiltonian between two atomic orbitals. A special subroutine, incorporating these parameters, was written to evaluate the matrix elements and is 12 listed in Appendix H. Energy eigenvalues of GaAs and ZnSe, at sym- metry points r, L and X in the Brillouin zone, are compared in Tables 1 and 2, respectively. Note that r implies R'= (0,0,0), L implies T? = n/& (1,1,1), and x implies T? = 2n/a (1,0,0) where "a" is the lat- tice constant. The calculated band structures for GaAs from r to L and from r to X are shown in Figure 2. These can be compared with 1 which is mislabeled as Ge instead of Figure 7 of Chadi and Cohen, GaAs. Similarly, the energy bands for ZnSe from r to L and r to X, plotted in Figure 3, are equivalent to portions of Figure 9 in Chadi and Cohen.1 For the sake of completeness, the energy bands of Si and Ge, using Chadi and Cohen's first nearest neighbor parameters, are given in Tables 3 and 4, and in Figures 4 and 5, respectively. Since Harrison's formulas in Equations 15 through 21, in conjunc- tion with the Herman-Skillman atomic term values, provide a more gen— eral way of numerically evaluating matrix elements of the Hamiltonian between two atomic orbitals, a comparison of the energy band struc- tures for the Zincblendes using these parameters was made against the calculations with the Chadi and Cohen parameters. The results for GaAs, ZnSe, Si and Ge, assuming only first neighbor interactions, are given in Tables 1, 2, 3 and 4, and in Figures 6, 7, 8 and 9, respec- tively. The qualitative shape of the energy bands from Harrison's parameters agree very well with the bands of Chadi and Cohen, but the quantitative pr0perties such as the band gap for GaAs are not as accurately predicted. To investigate Harrison's formulas further, the inclusion of higher order near neighbors was attempted. One would expect the band 13 Table 1 Selected Energy Eigenvalues for GaAs Using First Nearest Neighbors Symmetry Chadi & Cghen . Calculated Degeneracy P01nt Results Chadi & Cohen/Harrison Parameters r -12.4 eV -12.43 eV -22.07 eV 1 r 0.0 0.00 -9.54 3 r 1.6 1.63 -6.63 1 r 4.8 4.78 -3.27 3 L -10.7 -10.69 -20.20 1 L -6.2 -6.24 -15.58 1 L -1.2 -1.19 -11.44 2 L 1.7 1.70 -6.19 1 L 6.0 5.97 -1.37 2 L 9.21 0.47 1 X -9.7 —9.71 -19.35 1 X -6.8 -6.76 -15.31 1 X -2.8 —2.82 -13.43 2 X 2.2 2.16 —3.97 1 X 7.60 -2.88 2/1 X 8.29 0.62 1/2 14 Table 2 Selected Energy Eigenvalues for ZnSe Using First Nearest Neighbors Calculated Symmetry Chadi & Cghen Degeneracy Point Results Chadi & Cohen/Harrison Parameters r -12.1 eV -12.11 eV -23.62 eV 1 r 0.0 0.00 -10.57 3 r 2.9 2.91 -5.10 1 P 7.5 7.54 -2.34 3 L -11.0 -11.02 -22.36 1 L -4.7 -4.71 -15.22 1 L -0.75 -O.75 -12.14 2 L 3.9 3.88 —4.96 1 L 8.3 8.29 -O.77 2 L 10.19 0.91 1 X -10.6 -10.58 -21.88 1 X -4.8 -4.81 -14.37 1 X -1.9 -1.93 -13.94 2 X 4.7 4.65 -3.56 1 X 9.08 -1.82 1 X 9.47 1.03 2 15 2/ 5/1979 ENERGY BANDS --- ENERGY vs UAVEVECTOR 18 I I l I l I l I h ‘l 5 _. E - ‘1 N - . S a r G - A Y - _ é .. .1 V -S r- - - -( )— .— ~10 - ._ _15 ' 1 1 1. L 1 I 1 1 1 ‘ 0.0 0.2 0.4 0.6 0.8 1.0 BSEKJ (ixeN STR N5) K-VECTOR was 11 c .00. .55. .950) T0 c .55. .55. .55) 2/ 5/1575 ENERGY BANDS --~ ENERGY vs UAVEVECTOR 1° ' T ' r ' T r r '* l 1 - 1 5 =a-________I__.g5=_.=_—=:::::::::::::::::::::.':'_'::':Z:::::::::: .1 L _ S - 1 E ._ R o _ _ G Y b q 5 -s _- - -13 _. .. _15 ’ 1 I 1 1 1 1 1 L 1 ) 1 ‘ G.G 9.2 0.4 a.G 0.3 1.6 1.2 ABSEK] (l/ANGSTRO 3 K-VECTOR VARIES FROM C .20. .08. . 0) TO (1.11. .88. .80) Figure 2. Energy bands of GaAs using Chadi and Cohen's first nearest neighbor parameters; (Top) from F to L, (Bottom) from F to X 16 2/ 5/1979 ENERGY BANDS --- ENERGY VS UQVEVECTOR 15 _ T I l r T I l r l C I E s —- l N b q E - R b _ G L - Y 0 ._ E - _ V - J - _ 4 -SF— ——( )- -( -1a —- i r 1 _15 ' 1 l 1 1 1 1 1 1 1 ‘ 0.0 0.2 0.4 0.5 0.8 1.0 BSCK] Cl/AN STRO 3 K-VECTOR VfiRlES FRdE C .00. .05: .03% T0 C .55: .55. .55) 2/ 5/1979 ENERGY BRNDS --- ENERGY VS UAVEVECTOR 10 1 1 1 I T—l i I r 1 L )- -1 p- —( 5 _- ‘7 E ‘ 2 E a- R _ . G Y _ _ 5-5— ~ -10 —- _-——_—-—'_—_———______________ “ L I _15 1 1 1 1 1 1 1 1 1 1 1 ‘ 0.0 0.2 0.4 0.5 0.8 1.0 1.2 ABSEKJ Cl/ANGSTROHS) K-VECTOR VARIES FROM C .00. .00. .00) TO (1.11. .00. .00) Figure 3. Energy bands of ZnSe using Chadi andCohen's first nearest neighbor parameters; (Top) from F to L, (Bottom) from F to X 17 Table 3 Selected Energy Eigenvalues for Si Symmetry Chadi & Cohen* Calculated** Degeneracy Point Results Chadi & Cohen/Harrison Parameters P -12.16 eV -12.16 eV -21.27 eV 1 r 0.00 0.00 -9.50 3 r 3.42 4.10 -5.83 1 r 4.10 6.34 -3.54 3 L -9.44 -9.49 -18.82 1 L ~7.11 -6.64 -16.34 1 L -1.44 -2.17 -11.73 2 L 3.92 —5.82 1 L 8.51 -1.31 2 L 10.49 0.84 1 X -7.70 -7.32 -16.86 2 X ”2085 '4034 ”13096 2 X 6.46 -3.21 2 X 10.68 0.92 2 Notes: * for first and second nearest neighbor parameters ** for first nearest neighbor parameters only. 18 Table 4 Selected Energy Eigenvalues for Ge Calculated** Symmetry Chadi & Cohen* Degeneracy Point Results Chadi & Cohen/Harrison Parameters P -12.57 eV -12.57 eV -21.49 eV 1 r 0.00 0.00 -9.10 3 r 0.99 .99 -7.27 l r 3.24 -5.24 -3.62 3 L -10.30 -10.33 -19.13 1 L —7.52 —7.25 -16.03 1 L -1.60 -2.10 -11.16 2 L 1.96 -6.64 1 L 7.34 -1.56 2 L 9.28 0.32 1 X -8.60 -8.36 -17.09 2 X -3.30 -4.20 -13.22 2 X 5.19 -3.65 2 X 9.44 0.50 2 Notes: * for first and second nearest neighbor parameters ** for first nearest neighbor parameters only. Figure 4. «mmzm 1 EOPQM>IX mm Dmhm82¢\fiu Hyummc 8.8 m.8 v.8 M.8 mw.8 .H.8 .NH 8.3m.E F_FL______________ _—_____ _ _ _ mOH8w>w>¢3 m> >8mmzm III mwma\h \N mozcm >8mwzm .8 mwml .8NI man. 8H1. ILLI>I UJZLLJOCL’D- 3O .mgonsmwm: pmmgmm: umgwm Low mgmmemcma m.:0mwggmz mcwm: Nm mo+um>1¥ nmzomkmoz¢\fiu Axummc mw.® mw.® EV.® mn.® mw.® 8.8 mw.® I _ _ _ _ _ _ _ a _ _ _ _ _ _ fl _ _ _ _ _ _ _ _ q _ _ _ q ... mmml l L 1 lAHHlIIIII/ 11 911.. I mHI II nmmn111nnHuHnHuunnHHunnHHHHHHHHHHHHHHHHHHHHHII. 11 881. H 1 I IIflIIIIIIllIIIIIIIIIIIIIIIII 1 TI 11 m1 1 111.1111 M 111111 1111111 a I 1 l_.___________E_______________lm moeum>m>¢3 m> >ommzm 111 mozcm >ommzm mnma\w \N IL|J>I LUZUJDCL'D- 31 mgongmwm: pmmgmwc mewm Low mcwpwEmLma mcomwcgm: mcwm: Namwuu mo mvcmn xmgmcm .efi m.=mEE nmm. .88. .88. 8 8H n88. .88. .88. 8 [88¢ mm_m¢> m8+88>1¥ nmzomkm82¢\~u nXummq m.m. mw.® E..® m.mw m.mw “.3 “w.® — — _ T _ d _ _ . _ _ _ _ _ fl 8 _ _ 8 fl _ _ . _ _ _ _ _ mm... H1111nnHHHHH11. H I IIIIIIIIIIIIIIII I I1 8N1 I IIIIIIIIIIIIIIIIIIIIIIH I IIIIIIIVII‘II‘II‘I‘IIIII I I lIflIIIIIIIIIIIII 1 1. 11. .1 m_1 INK 81 ll Ill lit 11 W... H 11.... a 11.111 Iluufldllufluhnfluu III I _ _ _ _ _ . _ _ EL _ _ _ _ _ _ _ _ p b _ _ _ E _ _ _ _ J m mapum>m>¢3 m> >omwzm 111 mozam >ommzm mwmfi\w \N ILLI>| uzmmo> 32 mgonsmwmc pmmgmmc pmgvm Lek mgwmemLma m.:0meLmI mcwm: mmmm mOka>nx nmzomhmwz¢\~u n¥umm ¢.H N.“ ®.~ m.& m.& ¢.® N.® 8.8 _ _ _ _ _ — _ _ _ _ m _ m _ m r‘ m _ _ _ m _ «Ohom>m>¢3 w> >ommzm 1|: mozqm >ommzm mwmfi\m \N mmml .wml mfi.. 8H1. |L1J>| UJZUJ01(D>- 33 ticular, the threefold screw axis symmetry was used to reduce the 84 states to 28, or to seven atoms, as above. But in addition, these seven atoms (one formula unit, Tl3AsSe3) exhibit triangular sym- 7 metry; i.e., three reflections and two rotations. The irreducible representation for the group 03v are given below.10 E 30 3o 30 2C 2C A1 1 1 1 1 1 1 112 1 -1 -1 -1 1 1 _1_ -3 1J3 -1 3 -1f3 (23) A3 (10>(-10 2 —2_ 72-7 77 77 01 01 —\/§ ;1_ fl ;1_ 431 _._/_3_'_;1_ 2 2 2 2 2 2 2 2 The character table for the irreducible representations of C3v is written as follows E 30 2C A1 1 1 1 A2 1 -1 1 (24) A3 2 o -1 By considering a linear combination of atomic orbitals that is com- patible with the triangular symmetry (see Figure 16) and examining how these states transform under symmetry operations, the following character table for TAS was obtained. 34 s orbitals pr orbitals pt orbitals Figure 16. Atomic orbitals with triangular symmetry 35 Therefore, A = 8A1 + 2A2 + 9A3, (25) (26) which implies that in the three irreducible representations A1, A2 and A3, there are eight, two and nine basis states, respectively. These basis states are7 A1: |1> = s(As) |2> = pZ (As) |3>= 2f; [s(Se)a + s(Se)b + |4> = 7%_[pP(SE)a + pr(Se)b |5> =711pZ (Se) + p213)»b |6> =«-\—/r:1§-[s(Tl)a + s(T|)b + 17> =7l1p (TI) + prunb 18> =71—[p m) + plumb A2: |1> = Ué—[pt(5e)a + pt(Se)b |2> =-3%-[pt (Tl)a + pt(T|)b s(Se)C] + pr(Se)C] + pz(Se)C] s(T|)C] + + + + pr(T|)C pt(Se) pt(T|) ] pZ(T|)C] C] J (27) (28) (29) (30) (31) (32) (33) (34) (35) (36) 36 A3: l1> :'U%'[px(AS) + ipy(As)] (37) 12> = f1s15e1a + 1.513121b + ¢*S(Se)c] (38) |3> = 7%—[pz(Se)a + <1>pz(Se)b + o*pZ(Se)c] (39) 14> = % [pr(Se)a + opr(Se)b + <1>*pr(Se)C] (40) 15> =7§1pt15e1a + opt(Se)b + *pt(Se)C] 1411 |6> = 7;- [5111161 + «151111b + ¢*s1111c1 (42) 17> =7§— [pZ(T|)a + <1>pZ(T|)b + *pZ(T|)C] (43) |8> = 7:1, [pr(T|)a + 1111,1111b + 1*pr1111c1 1441 19> =7; [pt(T|)a + opt(T|)b + o*pt(T|)C] (45) where o = exp(2ni/3). Nine additional basis states are obtained by complex conjugation of the A3 states, thus giving the desired 28 states. Now, by following the formularization set forth in Equations 1 through 7, one only needs to replace the atomic orbitals in the Bloch function (Equation 1) by the above basis states. Since these basis states are in the irreducible representations, Equation 7 turns out to be in "block-diagonal" form. Therefore, the diagonalization of the 28x28 matrix is reduced to four independent diagonalizations: an 8x8 for A1, a 2x2 for A2, and two 9x9's for A3. This analytical 37 "simplification" was considerably more cumbersome to code as a Fortran program when compared with the program listed in Appendix C. Nonetheless, the resultant energy band structures for A1 and A2, using Harrison's formulas in Equations 15 through 21 and the Herman- Skillman atomic term values to numerically evaluate the matrix ele- ments for first nearest neighbors, are shown in Figure 17. Upon close comparison of this figure with Figure 15, it is obvious that the diagonalization of the 28x28 matrix yields identical eigenvalues to the method utilizing the basis states in the irreducible represen- tations, provided that the interaction parameters and number of nearest neighbors are the same. Symmetry Considerations The above examples of energy band structure calculated for TAS indicate that the symmetry or group theory portion of the energy bond problem can be set aside by the programmer and handled solely by having the computer diagonalize a larger matrix (at correspondingly longer computer execution times). In particular, a graph similar to Figure 15 for TAS would be generated, in a reduced Brillouin zone, by a diagonalization of the 84x84 along the c-axis. To make this point more explicit, silicon was treated as a chalcopyrite -- that is, eight atoms per "unit cell" with a twofold screw axis symmetry. The resultant band structure along the c-axis is shown in Figure 18 and, for direct comparison, the bands for the two atoms per unit cell are also plotted. The existence of an additional symmetry can be Figure 17. 38 2/ 5/1979 T95 ENERGY BANDS a I I f I I I I I I I r I I _ J - J E I R _ _ G 5 -15 —- _ 1— -1 F- f d r : -25 P 1_ i L l 1 l 1 I L l 1 1 1L 0.0 0.2 0.4 0.5 0.8 .0 1.2 1 K2 Cl/HNGSTRONJ 2/ 5/1979 T93 ENERGY BANDS e r I 1 I I r I I I I l I l ‘5 j\ _ E L ' E 10 F _ R L _: g _ 1 - ‘i (I. -1 5 -15 _ __ - )" q '- "l )- .4 - 4 -29 L- _ 1' -1 _25 1 . 1 1 L 1 1 1 1 1 1 . 0 G 6.2 0.4 . 6.5 6.8 1.9 1.2 ' 4 K2 (l/ANGSTRCH) Energy bands of two irreducible representations of T13AsSe3 Harrison's parameters for first nearest neighbors; usin (Top 3 A 19 (Bottom) A2 :«fifih 39 ' 2/ 5/1979 ENERGY BANDS --- ENERGY vs UAVEVECTOR S I I I I I r I T I I I I I 9 L— .2 F q E -s - - N 1 E _ R d $ a - -1 E ..J V .1 _25 r 1 1 1 1 1 1 1 1 1 1 1 1 9.9 9.2 9.4 9.5 9.3 1.9 1.2 BSEK] (l/A GSTRO sn K-VECTOR VARIES Enda c .99. T99. .99) T0 c .99. .99.1.15) ‘" ' ‘”“fl" 2/ 5/1979 ENERGY BANDS --- ENERGY vs UAVEVECTOR S I I I I I I I I I T I 9 —- ._ I 1 E -s ;;______...__—-——E*”"’"~——'——» _ N —1 E .. .1 R - CI! 9 _ Y -18 — _..1 E ~ - v - - - 1 - ~15 1- ._ )— .1 -2a —- -1 ” 1 _25 C 1 1 1 1 1 1 1 1 1 1 1 i 9.9 9.2 9.4 9.5 9.8 1.9 1.2 ABSEKJ c1/nNGsrRoNS) K-VECTOR VARIES ERGN c .99. .99. .99) T0 c .99. .99.1.1s) Figure 18. Energy bands of Si along c-axis; (T0p) one degree of symmetry omitted, (Bottom) all symmetry included 40 recognized from the mirror reflection symmetry about the reduced Brillouin zone edge at k2 = n/(c/Z) in the energy bands. More importantly, by using the larger matrices, the diagonal- ization is not restricted to wave-vectors that lie on symmetry axes, but can be done for any wave-vector in the Brillouin zone. CONCLUSION To summarize the section on energy bands, a computer program has been written and tested to perform energy band structure calculations on_any_crystal. The foundation of the approach lies in the LCAO model and the Slater and Koster two-center integral approximation. The only input requirements are the crystal geometry and a set of parameters which characterize the interactions between pairs of atomic orbitals. DENSITY 0F STATES CALCULATIONS Since it is now possible to calculate the energy eigenvalues for any crystal at an arbitrary wave-vector in the Brillouin zone as shown above, the corresponding energy density of states may be ob- tained by performing an integration over a constant energy surface in the Brillouin zone:11 dS 9(E) ~ I VIE . (46) Considerable effort has been put forth to increase the computational efficiency of algorithms evaluating this type of integral.12'14 But if the calculation is done on a minicomputer, which is idle more than 50% of the time, simplification in programming is traded for maximum efficiency of execution. Therefore, the statistical histogram or root sampling method15 can be used to greatly simplify the density of states calculation. This method involves generating the energy eigenvalues at a large number of random wave-vectors in the Brillouin zone and then sorting these eigenvalues into "buckets" that subdivide the energy range. This procedure suggests that the shape of the Brillouin zone must be known. Because different crystals possess Brillouin zones of various shapes, the sought after generality of applicability to all crystals, as exhibited by the energy band pro- gram, would be lost. In order to retain this generality feature, the following fact was recognized and employed: choosing a random wave- 41 42 vector in a single Brillouin zone is equivalent to choosing a random wave-vector in any random Brillouin zone. Consequently, in order to simplify the geometry, a wave-vector that falls within an "arbitrarily large sphere" was selected. This "sphere" must be large enough that the edge effects do not interfere with randomizing. The edge effect problem can be readily understood with an example. Consider a cubic Brillouin zone and let the "sphere" circumscribe the cube. A wave-vector, ?, that falls outside the cube is equivalent to the wave-vector, E“, shown in Figure 19, where K" =‘? + 6 and 3'15 a reciprocal lattice vector. Therefore, for a random wave-vector inside the sphere, the probability that this wave-vector (or its equivalent) falls between the dotted circular arcs of Figure 19 and the Brillouin zone edge is enhanced by a factor of two, due to the portion of the sphere outside the cube which is folded back into the Brillouin zone. By choosing a sphere of suffi- ciently large radius, the perturbation on the randomizing process, from fractional Brillouin zones at the edge of the sphere is small because: i) the probability that a random wave—vector falls within a partial Brillouin zone is given by the ratio of the volume of all fractional Brillouin zones to the volume of the sphere; and ii) the fractional Brillouin zones, when combined, will form a number of "complete" Brillouin zones which do not contribute to the edge effect. 43 '— x 5") Figure 19. Nave-vector, E, which fallsfioutside a cubic Brillouin zone is related to wave-vector, kg within the Brillouin zone by reciprocal lattice vector, 44 RESULTS The flowchart and listing for a Fortran program to compute the density of states, using the above method, are given in Appendix D. The density of states of various zincblende compounds -- Si, Ge, GaAs and ZnSe -- are shown in Figures 20, 21, 22 and 23, respectively, for 1 and Harrison's para- first nearest neighbors and both Chadi-Cohen meters (Equations 15 through 21). The plots for Ge, GaAs, and ZnSe can be compared with Figures 4, 7 (mislabeled), and 9 of Chadi- Cohen,1 respectively. To carry out the density of states calcula- tions on the chalcopyrites on Tl3AsSe3, the symmetries along the c- axes are no longer useful. That is, the diagonalization of the full 64x64 or 84x84 matrix for each random wave-vector is required. These calculations will be performed when the program is transferred onto a CDC computer. SUMMARY In summary of the density of states calculation, a computer pro- gram is available to determine the density of states for any crystal, requiring only the specification of atomic numbers and atomic coordi- nates of the atoms in the crystal, along with a set of parameters which characterize the interactions between pairs of atomic orbitals. Figure 20. 45 1/31/1979 ENERGY BAND DENSITY 0F STATES 1.5 I IfI I I I I I I I I I I I IT I I I I I I I I I r L _ E - 1 N ? 108 _ .— T - Y 1 0 _ F - d S T - ‘ ‘1‘ 9 s g . I- -1 r - 0.0 L I 1 LI I l [J I I I I l I l L l I l l I l L l 1 l -15 ~10 -5 D S 10 IS ENERGY (EV) / 1/1979 ENERGY BAND DENSITY OF STATES 1'5 I—rT I I I I I I I I I I I I’I I I I FI 7 T I I I I L— ..1 D _ E F T b ‘ Y 0 I K _ F 1— _1 S T __ _ A g 305 '— —1 S P . '- -1 .- I 1 6.9 I I l l I l I I I I 1 L1 1 I l I I I I I l l I ' l L l -25 -20 -15 -1@ -S D 5 Density of states ENERGY (EV) of Si for first nearest neighbors; (T0p) Chadi and Cohen parameters, (Bottom) Harrison's parameters Figure 21. 1.5 1.0 (om—11> -IU) 'no (fi—(DZITIU 1.5 1.0 (om-1:1) -I(.0 TIO K-I—(DZI'TIU 46 1/30/1979 ENERGY BAND DENSITY 0F STATES III I I IlllI’IIil I I l l l FTI I IIII I ll — -—-1 LIIIIIIIIIILIJ IIJIIIILLIIII -15 -1@ —S 0 S 10 15 ENERGY (EV) 2/ 1/1879 ENERGY BAND DENSITY 0F STATES IITTIIIIFIVIIIIIIIIIIIIIITIll 1 _ "“ 1 Ll'lIllllIIlJiJ IIIJILLIILI -2S -2@ -15 -1@ -S G S ENERGY (EV) Density of states of Ge for first nearest neighbors; (Top) Chadi and Cohen parameters, (Bottom) Harrison's parameters 47 1/25/1379 ENERGY BAND DENSITY 0F STATES 2'5 T I I I I I I I I I [—F I I I I I I I T I I I I 2.0 L— —' D b _ E - - N _. -1 s _ . 0 ' q r: .. _ 13' 1 9 1— — A — - T E u— -1 S - _ ans _— '— e'e _ l 1 L44 I I I I Ll I I I I I I I I I ‘ ‘15 '19 '5 G 5 10 ENERGY (EV) 2/ 1/1979 ENERGY BAND DENSITY 0F STATES 2‘0 I I T I II I I FI T I II I III I I I I I I I I IT D 1.5 —- '— E N F .1 S L _ I T — d Y _ - 0 F 1 9 — __ S ._ — T _ _ A _ 1 T E _ _ S 9-5 "'— —4 h- J 8.8 L L I I l I I I l L l l l I I I I I l l I I I 1‘1 ‘25 -2@ ‘15 '19 ‘5 0 5 ENERGY (EV) Figure 22. Density of states of GaAs for first nearest neighbors; (Top) Chadi and Cohen parameters, (Bottom) Harrison's parameters 48 1/29/1979 ENERGY BAND DENSITY 0F STATES 2'5 I I I I I I I I I I I I I I I I I T I I l I I I I I I I I '- A 2.0 h- _. o t - E P A : 1 T . —- —~ Y 1 5 _ - D b d F - II J S P ... ‘T 1.0 —- _. A .. - T E ” ‘ S - _ 0.5 - .4 ~ J 0.0 “.ng 1 I 1 1 1 I 1 1 LL 14 1 I 1 1 1 1 14 1 1— -15 -10 -S G S 10 15 ENERGY (EV) 2/ 1/1979 ENERGY BAND DENSITY 0F STATES 105 Ij I I I I I I I I I T I I I I I I r1 I fTI I I l I v- -1 D _ - E N ? 1'9 — _—‘ T _ _ Y O _ _ F _ - S T _ _ A E 9.5 E. .— S _ 0.8 I I I LI I I I I l I I I I I I I I I I I I I I I L4 -25 -2@ -15 -10 -S D S ENERGY (EV) Figure 23. Density of states of ZnSe for first nearest neighbors; (Top) Chadi and Cohen parameters, (Bottom) Harrison's parameters INDEX OF REFRACTION AND BIREFRINGENCE CALCULATIONS A systematic method for numerically describing the interaction between two atomic orbitals has been presented in the preceding sections and appendices. Now this "chemical bond“ approach will be applied to the prediction of the electronic susceptibility in semiconductors. An equivalent set of interaction parameters, as was utilized in the band structure calculations, will enter into the formularization. Harrison's Bond-Orbital Model (BOM),16’17 which provides the theoretical foundation, is reviewed below as an introduction to the concepts and notation. FORMULISM Consider a linear combination of an anion sp3 hybrid orbital and a cation sp3 hybrid orbital, lw> = ualha> + uclhc> (47) where 2+u2=1 (48) for normalization. The coefficients, ua and uc, are determined by minimizing the energy, E, of the hybrid orbital with respect to these coefficients. That is, 49 -—- = 0 and ——— = 0 (49) where _ < IHI > E — 'Tgififig‘ . (50) By defining the following parameters: v2 = — = - (51) ea = (52) cc = , (53) the energy can be written as E = uazea +2uC25C2- Zuaucv2 . (54) ua + uC Then, from Equations 49, .§§__ = 2uaea - 2uCV2 - E(2ua) = 0 aua uaZ + UCZ uaZ uC2 or Ua(ea-E) - ”CV2 = 0 . (55) Similarly, from (BE/auc) = O, 51 ”C(EC-E) " UaVZ = 0 o (56) If a nontrivial solution for ua and uc, in Equations 55 and 56, is to exist, then the determinant of the coefficients must be zero. That is, (ea-E) -V2 = 0 (57) "V2 (EC-E) OY‘ ca€C+E2 -E(ea+ec) -V22 =0. (58) The energy solution to this quadratic equation in E is 8 +8 a c 1 E =-———————- i 2 2 (59) 2 72"/(ga+cc) -4€acc+4V2 , Now, defining v = -1—(e -e) (60) 3 2 c a and E=L(e +8) (61) 2 c a yields I'T'I II —- 2 2 s: t ‘fihz +‘v3 . (62) 52 The "minus" sign corresponds to the bonding state and I'plus" sign gives the antibonding state. By substituting the bonding energy back into Equations 48 and (55 or 56), solutions to u and uc are obtained a for the bonding orbital: ua = __7?_JL and uC = -—7TJ1 where 2 2 V u - u 3 a c a = = . (64) p 2 2 2 2 \«yz + V3 ) ua + uc Note that V2, V3 and a are called the covalent energy, the polar P energy, and the polarity, respectively. Now, since the probability of a bonding electron being "on" the anion or cation is ua2 or ucz, the charge per bond (two electrons for each bond) associated with the anion is -26 u 2 ‘ -e(1 + up) (65) and with the cation is -2e u 2 - -e(1 - up) . (66) The electronic dipole moment of the bond is therefore 9 = % qixi = -e(1 + up) - e(1 - up) (67) 53 —). where and <§ a C> represent the average vector positions of the bonding electron distributions about the anion and cation, respec- tively. If the origin is chosen to be at the center of the bond and . +~. . . . . if d 15 the internuclear vector distance from the cat10n to anion, then (Ya) = +y d/2 (68) and <78 = -y 6/2 (69) where y scales the "center of gravity" of the electron charge distri- bution in both anion and cation hybrid orbitals. Therefore, the di- pole moment per bond can be written 3 = -e(1+ap)y’6/2 + e(1-ap)y‘d’/2 which reduces to 3 = we up 3. (70) A change in this dipole moment occurs in response to an external electric field, 3. Upon application of this external electric field, there is an interaction energy of -p°€ = weap d°Z (71) which corresponds to a change in polar energy, V3,16 given by 54 V3 + V3 ‘Yéfi'E/Zo (72) The polarization per bond, when an external electric field is present, is derivable from Equations 64, 70 and 72 as follows: —> + + dV -yed°g/2 #(V2 + (V3-yed°g/2) ) Using a binomial expansion,19 n x2 (1+X) =1+nX + "(n-1) —2"|—+ 00000, (74) on the denominator for small electric field, and collecting terms with similar powers of E yields Y2e2v 2 3': ’Yeafia + g 2 372 3(3.2) + 2(V + V ) 2 3 (75) 3y3e3V 2v 2 3 + +o+ 2 2 25/2 d(d E) + .. 8(v2 + V3 ) As an alternative, the first and second order terms of the polar- ization per bond can also be separately obtained from —) (76) Q) Q) 0113+ H A O) O) 3 3.1: 3.. 3P. __3_ = 3 p 3 2 3 3V3 respectively. In order to compute the overall polarization, the contributions from each bond in a unit cell are added together and normalized to volume. That is, 3 = %. X (36 + Bi + Bé + ....) (78) bonds where 36, Si and 32 are the zeroth, first and second order dipole moments per bond as given in Equation 75 and V is the volume of the unit cell. For non-ferroelectric materials, the zeroth order term averages to zero. (This makes a nice consistency check to verify that the apprOpriate bonds are all included in the sum.) And for zincblende crystals, with a regular tetrahedral structure, the first order geometric factor, 3(3‘3), averages to dzacosze = dZE/B (see Appendix I), as is found in the literature.16’18 By evaluating Equation 78 for each of the three orthogonal elec- tric fields, Ex, 5 and $2, the nine components of the susceptibility y tensor are obtained since Px X11 X12 X13 5ix 0 0 py = x21 x22 x23 0 or 5y or 0 (79) P2 X31 X32 X33 0 0 52 - Diagonalization of this 3x3 susceptibility matrix results in the three principal susceptibilities which are denoted by x1, X2 and 56 x3. Consequently, the principal indices of refraction can be calculated from n. = + nxi for i = 1,2,3. (80) For all three principal indices of refraction equal, the crystal is 20 If only two of the principal indices of re- optically isotropic. fraction are equal (e.g., no = n1 = n2 # n3 = ne), then the crystal is uniaxial with the birefringence given by An (81) ll 3 I 3 where subscripts e and "0" denote extraordinary and ordinary. Finally, when all three principal indices of refraction are unequal, the crystal is biaxial. IMPLEMENTATION Therefore, the principal indices of refraction can be numerically predicted by incorporating Equations 75, 78, 79 and 80 into a com- puter program. The only requirements to perform this calculation are the crystal geometry, a procedure for summing over bonds, and the set of parameters V2, V3 and 7, which describe the interaction between 3 hybrid orbitals. A Fortran program (flowcharted and each pair of sp listed in Appendix E) was written utilizing the above formulari- zation, specifically for the tetrahedrally coordinated chalc0pyrites (see Appendix F). The geometry was again handled separately by the generation of a list of atoms and their coordinate positions. (In 57 fact, both this calculation and the band structure calculations utilized the same geometry file list.) Two different sets of parameters describing each hybrid bond were tried. The first is a set of empirical parameters based on experi- mental dielectric constants of the III-V and II-VI binary semiconduc- tors.17 V2 and 7 were found to correspond to the row of the periodic table to which each atom involved belongs. And the set of V3 was directly related to a set of hybrid energies via Equation 60. A second set of parameters, based on atomic orbitals rather than sp3 3 hybrid orbitals can be decomposed hybrids, were also used. Since sp into a linear combination of atomic orbitals, it is possible to ex- press V2 and V3 in terms of Harrison's "Solid State Interatomic Matrix Elements" (Equations 15 through 21) and the Herman-Skillman atomic term values.6’8 To be more explicit, the hybrid orbitals for the anion and the cation are written - 1 -J3 and Iha> _.7453> +"z-a'mx'pax) + dylpay> + dzlpaz>) (82) _ 1 ‘J3 IhC> -'2'ISC> -‘2'a- (dxlch> + dylpqy) + dZ|pC2>) (83) where d = dxx and dyy + dz? is the vector distance from cation to anion. [Notesz i) the minus sign in Equation 83 makes the cation hybrid orbital point in the Opposite direction of the anion orbital; ii) the factors di/d permit the directions of the hybrid orbitals to deviate from the (i1, t1, :1) directions of a "regular" tetrahedron -- that is, these hybrid orbitals can reflect the distorted tetra- hedral angles of the chalcopyrites; iii) the orthogonality of the 58 hybrid orbitals is lost except, for the "regular" tetrahedron, when dX = dy = dZ = d/-J3]. Then, from Equation 51, 4 4 3 . v = ) ) d.d. (84) 2 4d? i=1 j=l ‘ J C‘ 33 where aci and aaj are the cation and anion atomic orbitals (i and j = 1,2,3,4 denotes s, px,py,pz), di and dj equal -d/\/3 and +d/~J3, dx, dys dz for i and j = 1929394 and (aC‘ilH'aaj) is a "501'id State Iriter- atomic Matrix Element“. Also, from Equations 52, 53, and 60, V : 3 116;). 36g) - (a: + 363)] (85) and e are "atomic term values" for the cation and where 65, 8p, 55 p anion. RESULTS FOR CHALCOPYRITES The calculations were performed on the chalc0pyrites because: i) the lattice parameters and Optical pr0perties have been experi- mentally established and ii) the chalc0pyrites, being closely re- lated to the zincblende semiconductors, can be treated by simply extending some well-develOped empirical theories for Zincblendes. The indices of refraction as a function of wavelength for many chal- copyrites were measured by Boyd et. al.21‘25 The predicted static (long wavelength limit) indices of refraction and birefringence, using the above BOM theory with parameters from the binary semicon- ductors and only the first order polarization, are compared to experiment in Table 5. Also listed in Table 5 are the results of 59 TABLE 5 INDICES OF REFRACTION AND BIREFRINGENCE FOR VARIOUS CHALCOPYRITES--A COMPARISON OF THEORY AND EXPERIMENT Experimental* Theoretical (Boyd et. al.) Lines & lst Order Binary Waszczak BOM Parameters Compound no An n 110 he An CuGaS2 2.4275 -.0154 2.51 2.525 2.484 -.0414 AgGaS2 2.3266 -.0550 2.45 2.524 2.305 -.2186 CulnS2 2.5166 -.0179 2.53 2.519 2.522 +.0032 ZnGeP2 3.0552 +.O397 3.10 3.114 3.084 -.0302 CdGeAs2 3.5004 +.0868 3.36 3.414 3.277 -.1368 CdGeP2 3.1165 +.0167 3.08 3.123 2.990 -.1330 ZnSiAs2 3.1685 +.0268 3.22 3.248 3.178 -.0704 CuAlSe2 2.4659 —.0126 2.64 2.667 2.605 -.0615 AgGaSe2 2.5731 -.0327 2.61 2.684 2.485 -.1985 CuGaSe2 2.6872 +.0026 2.69 2.698 2.655 -.0430 AglnSe2 2.6167 +.OOl7 2.63 2.655 2.573 -.0822 * at longest measured wavelength (generally ~12pm) 60 Lines and Naszczak,18 who used an equivalent BOM theory, but treated only the average index of refraction and did not attempt to include the birefringence. As is evident from Table 5, the computer program of Appendix E is capable of reproducing the number of Lines and Naszczak for the average index of refraction which shows reasonable agreement with the experimental values. However, the predicted birefringence appears to have_ng_correlation if compared to the corresponding measured quantity. More particularly, even the sign of the birefringence differs when contrasting theory and experiment for all of the II-IV-V compounds and some of the I-III-VI compounds. As an attempt to resolve the above conflict with regard to birefringence, the same calculation was repeated using a different set of parameters (Harrison's "Solid State Interatomic Matrix Elements" and Herman-Skillman's atomic term values). The resultant first order ordinary index of refraction and birefringence are given in Table 6. For comparison, the first order predictions with the binary compound BOM parameters and calculations including the second order polarization term are also presented in Table 6. It is evident from Table 6 that the relative qualitative arrangement of birefrin- gence from compound to compound is, in general, unchanged by using a different set of BOM parameters or including the second order polarization. In addition, slight variations in the assumed chalc0pyrite geo- metry were tried in order to account for the discrepency in bire- fringence between experiment and theory. Specifically, the calcu- lated birefringence is plotted as a function of lattice constant, x, 61 TABLE 6 THEORETICAL ORDINARY INDEX OF REFRACTION AND BIREFRINGENCE FOR VARIOUS CHALCOPYRITES lst Order lst Order lst & 2nd Polarlzation Polarization Order Polarization Decomposed sp3 Binary BOM Binary BOM Parameters Parameters Parameters Compound no An no An no An CUGaSZ 10620 “0018] 20525 '00414 30307 ’00587 AgGaS2 1.603 -.0928 2.524 -.2186 3.374 -.3308 CulnS2 1.590 +.0014 2.519 +.0032 3.384 +.OO45 ZnGeP2 1.823 -.0176 3.114 -.0302 3.862 -.0449 CdGeAs2 1.844 -.O666 3.414 -.1368 4.405 -.2271 CdGeP2 1.827 -.0736 3.123 -.1330 3.988 -.2125 ZnSiAs2 1.839 -.0349 3.248 -.0704 4.071 -.1077 CuAlSe2 1.630 -.0245 2.667 -.0615 3.533 -.0902 AgGaSe2 1.621 -.0780 2.684 -.1985 3.682 -.3075 CuGaSe2 1.640 -.0173 2.698 -.O430 3.623 -.0622 AglnSe2 1.591 -.0330 2.655 .0822 3.721 .1242 62 in Figure 24, using the first order polarization with binary semi- conductor BOM parameters for CdGeAsZ. As illustrated in this figure, the dependence on lattice constant is not strong enough to accomodate a change in sign of the birefringence. Synopsis The simplistic "chemical bond" approach to calculating the re- fractive indices and birefringence is very attractive when contrasted with the more common "oscillator strength" method. Consequently, a computer program employing the former scheme has been written in which, again, the crystal geometry and a set of parameters that de- scribe the interactions between orbitals are supplied. Upon appli- cation to the chalc0pyrites, the indices of refraction are predicted quite well (for the appropriate set of BOM parameters) but the bire- fringence is not. In trying to infer a qualitative aspect of the birefringence calculation, it appears that the geometric factors 3(3-3)" in Equation 75 are the dominant contributions to the sign of the birefringence as calculated using the theory presented above. The fact that the calculations do not agree with experiment implies that something in the theory is missing. Speculation as to possible origins is currently in progress and has led to reformulation and generalization of the BOM in terms of coupled states between valence and conduction bands which yields a susceptibility related to, in part, the energy difference between the states.7 In addition to not being able to account for the birefringence, the Bond-Orbital Model theory, as it now stands, assumes a static 63 x .ucmpmcoo mowppm_ to cowpoczw m mm Nm = EI¢> (A1) or 2 -h 2 _ W V (I) ‘1' VII) - E11) (A2) where w is the wavefunction for a single electron in the crystaI. In spherical coordinates, with vector components written as x = rcos¢sine y = rsin¢sine (A3) 2 = rcose , Equation A2 becomes 2 2 -h 1 a 2 a l a a l a -———-E—— -——(r -———) + -—— (Sine-——-) + ——-] w 2m r2 Br Br r251.“e 36 36 r251n29 302 + Vw = Em . (A4) Assume that all interactions between the electron and the rest of the crystal are described by the potential, V(?) and, furthermore, assume that this potential is a central potential; i.e., V(r,6,¢) = V(r) . (A5) 67 68 Then, in order to solve the differential equation, Equation A4,27 separate angular and radial coordinates, (20.6.6) = R(r)Y(e,¢) . (A6) Consider only the resultant angular differential equation, where B is a separation constant, [ .1 ~31 (sine-—§-) + 2 2 ] Y = -8Y (A7) and perform another separation of variables, 6 and 6, Y(e,¢) = O(e)<1>(¢) (A8) which yields 5fi:% = -a2¢ (A9) d4 and [sine-ag-(sine-ag-) + Bsinze] e = 029 (A10) where 62 is the second separation constant. The normalized solutions to this last pair of differential equa- tions are the set of spherical harmonics, YT(6.¢) = (-1)( m+|ml)/2 (il?%l);1;?l)l PT(cose) eim¢ (A11) 69 where B = l(l+l) for l = 0,1,2,..., a = m for m = 0, 11,12,..., and P1(cose) = P1 (cose) are the associated Legendre functions given by28 m _ 2 m/2 dm Pn(x) - (l-x ) (BET-P"(X) (A12) for m>O and n 2 n Pn(X) = 1 d (X -1) . (A13) 2"n1 dx" A few Of these angular solutions are written out in Table A1.3 By taking the apprOpriate linear combinations of these solutions, the "standard" atomic orbitals, which are real, of elementary chemistry are obtained, as shown in Table A2. Note that the subscript notation in Table A2 utilizes the spherical to rectangular coordinate con- version of Equations A3. The directional characteristics of some of the atomic orbitals, as shown in Figure A1, are determined by making three dimensional polar plots of each orbital in Table A2 against the angular direction specified by e and 6. The sign associated with each "lobe" Of an atomic orbital is also indicated. The primary purpose of this appendix was to present a review of atomic orbitals while specifying the notation. In addition, the assumption of a central potential seen by any electron in an atomic orbital is again explicitly emphasized. 70 Table A1 Spherical Harmonics m Symbol Y¥(e,o) O s ‘L%; .. ___ . i6 p+1 ‘(8311 San e )3 0 p0 T C059 . -i¢ 8M San e .p I ...-x 'C I ....) W 51 n26 8214) N O. + N (AI-J N01 —S'I (10 C056 e11) ....I/151T 1 dtl 8n l)5 2 0 do 16" (3COS O ’ l) 15 . —i -1 d_1 -—§;—51ne cose e 15 2 -2i¢ -2 d_2 32 sin e e 71 Table A2 Atomic Orbital Functions _ l S - 34H _ 1 _ )3 . pX — 7[p_1-p+1] - Tfls1necos¢ _ ~1 _ 3 . — 17 [p_1+p+11 - Wsmecosa __ _ 3 pZ _ p0 — 4n C059 _ _ 15 .2 2 .2 d 2 2 — 7[d_2+ d+2] - T6? Sln 8(cos 6-51n 6) d:y-:l i %[d_2 - d+2] =‘/—%—%— sinzecos¢sin¢ dyZ = 1 -%-[d_1 + d+1] =‘%:%§T sinecosesin¢ dxz = %[d_1 - d+1] =‘/—_[1T_§ sinecosesinq; d 2 = do = -T%% (3cosze-1) 72 Figure A1. Directional nature of atomic orbitals APPENDIX B SLATER AND KOSTER 2-CENTER INTEGRALS OF ATOMIC ORBITALS Slater and Koster2 have provided a reasonable approach toward handling the geometrical considerations of a tight-binding calcu- lation; that is, a convenient and concise way of separating the "interaction parameters" from the orientation or geometry of the atomic orbitals involved in a "chemical bond". MEANING OF 2-CENTER INTEGRALS From Equation 14, the matrix element of the Hamiltonian between two atomic orbitals is written —+ —> (F—R:)H¢v(F—Rf)dv (B1) <¢u(?-Eg)|u|¢v(F—E§)> = f¢u where the Hamiltonian, H, is H =<fi—- v2 + v(?—R:) + V(T—fif) + 2 v (F-Eg). (32) (m.v) # (n,a) or (1.8) and the summation over (m,y) includes all atoms in the crystal except the two at R: and Rf. Therefore, the matrix element in Equation 81 can be broken up into separate terms. The "Three—Center Integrals" are 73 (‘F-fifmv (B3) where (m,Y) # (fl,a) or (1.8); that is, the matrix element involves two atomic orbitals on two distinct atoms interacting via a potential on yet a third atom. These are assumed to be negligible. The "Two- Center Integral", on the other hand, is a matrix element of atomic orbitals and potentials on only two atoms. COMPONENT NATURE OF ATOMIC ORBITALS Before proceeding with the discussion of the interaction (trans- fer integral) between two atomic orbitals, the transformation prop- erties Of an atomic orbital will be explored. First consider a two- dimensional transformation of coordinates for a p orbital, as illus- trated in Figure Bl. The p orbital is written in primed coordinates (see Table A2) as, p = c sine'cos¢' =-—$- x' . (B4) x r A transformation from the unprimed tO the primed coordinates, corres— ponding to a rotation of angle or about the z-axis, is represented in matrix notation by x' coser siner O x y' = ~siner coser 0 y (BS) 2' O 0 1 z . By expreSsing px' in terms of the unprimed coordinates, 75 Figure Bl. Two dimensional transformation Of a p orbital. 76 _ c . pX - -FT- (xcoser + ys1ner). (B6) And since r' is equal to r for a pure rotation, px = pX coser + py Sanr . (B7) In other words, a p orbital oriented in some arbitrary direction in the x-y plane can be written as a linear combination Of px and py orbitals. Now generalize the above result to three dimensions; that is, consider the transformation of coordinates for a p orbital in three dimensions, shown in Figure B2. Writing the p orbital in primed coordinates (again see Table A2), I _ I _ C I pZ - c cose - -Fr- 2 . (B8) Also, by transforming coordinates with a rotation or ¢r about the z- axis followed by a rotation of 9r about the y'-axis, x' coser O -siner cos¢r sin¢r O x y' = O 1 O -sin¢r cos¢r 0 y (89) z' siner O coser O O 1 z or x' c056r coser sin¢r coser -siner x Y' = -sin¢r cos¢r 0 y (810) z cos¢rsiner sin¢rsiner coser z . Figure B2. 77 Three dimensional transformation of a p orbital 78 Finally, as was done in the two-dimensional case, expressing pz' in terms of the unprimed coordinates with r'=r yields, __ _£; . . . pZ - r (xcos¢r51ner + y51n¢rs1ner + zcoser) (811) O)“ p z = pX cos¢rsin6r + pysin¢rsiner + pzcoser . (B12) Consequently, any randomly oriented p orbital can be written as a linear combination of three mutually orthogonal p orbitals, as in Equation 812. Furthermore, the coefficients of this linear combi- nation are seen to be identical to the spherical coordinate compon- ents of a unit vector (see Figure 82) oriented in the direction of the positive lobe of the p orbital: A A , /\ , , A p = x cosorsmer + y $1n¢r51ner + z coser . (813) INTERACTION BETWEEN TWO ATOMIC ORBITALS Define the Slater and Koster symbols2 (ssO), (Spa), (p50) and (ppn) as indicated in Figure 83. For certain orientations of s and p orbitals, the matrix element of the Hamiltonian is zero, as shown in Figure 84. Now specify the interaction between two atomic orbitals, oriented in arbitrary directions, in terms Of these symbols. In other words, express the two arbitrary atomic orbitals as a linear combination of atomic orbitals, where the coefficients Of each linear 79 \m U1 (550) = Bonding (pso) = (spa) = Bonding Antibonding X X Z p2 p1 p2 p1 (ppo) = Antibonding (ppm) = Bonding Figure 83. Slater and Koster's matrix element symbols 80 Figure 84. Transfer integrals between atomic orbitals which are zero 81 combination is determined by a rotation into a coordinate system with the orbitals oriented as in Figures 83 and 84. Consider first a p orbital interacting with an s orbital, a shown in Figure 85a. By expressing the p orbital as a linear combination of px, py and p2, as done in Equation 812, = cos¢rsiner + sin¢rsiner + (814) <52|H|plz>coser. This is illustrated in Figures 85b, c, and d. Comparing Figures 3, 4 and 5 suggests that only the last term of Equation 814 is nonzero and consequently Equation 814 becomes <52|H|p1> = (pso)coser . (815) Next, consider the interaction of two randomly oriented p orbi- tals, shown in Figure 86a. As before, both p orbitals can be written as linear combinations of px, py and pZ in the same manner as Equation 812. Therefore, = <(p2xcos¢zsin62 + pzysin¢zsin62 + pZZcosez)| H|(p1xcos¢1sinel + plysin¢lsinel + plzcosel)>. (816) Upon comparison of Figures 3, 4 and 6, it is evident that only three matrix elements are nonzero, as illustrated in Figures 6b, c and d. ,2 s s b) 2 + c) 2 + > - m p1 p1x + y y <521HIP1X> 0 <521HIP1y> = 0 x x z s d) 2 T <521Hlp12> : (p50) Figure 85. Interaction between an s orbital and a randomly oriented p orbital a) Z 2 b) c) p ' > p2y + 2x - +

(ppm)

(ppw) 2x 1x 2y 1y P1y _ + plx + Y x x 2 C + p22 _ d) (ppo) C9 D1 Figure 86. Interaction between two randomly oriented p orbitals 84 In terms of the Slater and Koster symbols, Equation 816 can be written = (ppn)cos¢zsinezcosolsinel + (ppn)sin625inezsin615inel + (817) (ppo)cosezcosel . Now introduce the "direction cosines" to put Equations 815 and 817 in the notation of Slater and Koster. Note that Slater and Koster required that the 21 and 22 axes of Figure 86a be parallel or perpendicular. The parallel case implies 41 = 42 and 81 = 82. Then the direction cosines can be defined as l = C050 = cos¢1sinel m = coss = sin¢lsinel (818) n = cosv = cosel . Hence, Equation 815 becomes = n(p56) (819) and Equation 817 becomes (using 12 + m2 + n2 = 1) = n2 (ppo) + (1-n2) (ppn). (320) Both Equations 819 and 820 appear in Table 1 of Slater and Koster.2 85 It is also possible to express the interaction between two atomic orbitals in vector notation.7 Define 3 as the unit vector from orbi~ tal #1 to orbital #2. In addition, for a p orbital, define 8 as the unit vector in the direction of the positive lobe of the orbital. Thesevectors are indicated in Figures 85a and 86a. From Equations 815 and 817, it is evident that <52|H|p1> = (61.3) (psO) (821) and = L61°6é - (61‘8) (62°8)] (ppn) + (61-8) (6; a) (ppo) (322) since 6108 = cosel, 6§.8 = c0592, and fiiofié = cosolsinelcosozsinez + sin¢lsinelsin¢zsin62 + coselcosez. Equations 821 and 822 are similar to those of Harrison‘s listed in Equations 16 and 17. The corre- lation between Harrison's parameters (Equations 18 through 21) and Slater and Koster's symbols (Figure 83) is given in Table 81. To conclude this appendix, the Slater and Koster method permits the numerical evaluation of the matrix element of the Hamiltonian between two atomic orbitals with the specification of a set of well- defined parameters. 86 Table 81 Correlation of Harrison and Slater-Koster Notation Harrison Slater and Koster VSS +(SSO) VSp +(p56) = -(Spo) Vppo -(ppo) V +(ppn) PP"I APPENDIX C ENERGY BAND COMPUTER PROGRAM An overall flowchart for the energy band program is shown in Figure Cl. The three software functions performed in this energy band computer program -- interactive input, calculation, and graphic output -- are handled by separate, stand-alone routines. The list- ings for these three programs, including non-library subroutines, are given in Figure C2. This division was made in an attempt to minimize the memory requirements of the calculation portion. The program is designed to run both on 16-bit minicomputers (e.g., DEC, DG, etc.) and the large number crunchers (e.g., IBM, CDC, etc.). The sub- routines listed here which evaluate the matrix elements Of the Hamiltonian between two atomic orbitals utilize Harrison's formulas (Equations 15-21) and Herman-Skillman's atomic term values. Replacing these subroutines with the subroutine in Appendix H yields the Chadi and Cohen calculations. 87 88 Specify crystal Start ‘ geometry Iterate wave-vector to R max Save energy eigenvalues wave-vector N0 —+ _ —> 7 Yes k _ kmax' Plot E I I .7 Figure Cl. Flowchart for energy band program u——>-—-—‘ from zero +' Determine each matrix element in Equation 7 from Equation 14 A at each -—"""— Diagonalize matrix to (Obtain energy eigenvalues 4—1 Stop ) 89 :UDD:EUBRNK:EBRND:DIR:EBINPUT.FR 0000000000000000000000000 000 158 168 00000 Figure C2. EBINPUT.FR ROUTINE TD RLLOU INPUTTING OF CONTROL PQRQMETERS USED BY PROGRRM ’EBRND.FR‘ FROM CONSOLE 9ND GENERRTES THE CONTROL PQRRMETER FILE. INPUT PBRRMETERS: FROM CONSOLE OUTPUT PRRRMETERS: I) CONTROL PRRBMETER FILE --- CONTRINING: Q) NRPUC -*- # RTOMS PER ’UNIT CELL’ B) NQTMNUM --- RTOMIC NUMBER OF EQCH RTOM IN ’UNIT CELL’. (ORDERED RCCORDING TO ‘ENUMERRTION’ (SEE BELOU)) C) UVKMRX --- 3D K-VECTOR (KX.KY.KZ) REPRESENTING MRXIMUM K-VECTOR D) KPTS --- # PTS IN K-SPRCE RT UHICH TO CRLCULQTE THE ENERGY BRNDS E) DMQXZ --- MRXIMUM 'NERR NEIGHBOR’ DISTRNCE SOURRED FOR RTOMS TO BE INCLUDED IN BLOCH FCN F) IFILE --- FILENRME OF INPUT DRTQ FILE (BELOU) FORMRT --- BINQRY EXTERNBL REFERENCES: NONE PRRQMETER IN=11.IOUT=IO.NPRGE=2?*4OOK+12 PRRQMETER KKPTS=IOO :MRXIMUM # KPTS RLLOUED PRRRMETER MRXRPUC=IO :MQXIMUM # QTOMS PER UNIT CELL DIMENSION NQTMNUM(MQXQPUC).UVKMRXI3).KXYZ(3).IFILE(31) DRTR KXYZ/’KX’.’KY’.’KZ’/ ENTER # QTOMS PER ’UNIT CELL‘ URITE (IOUT) NPQGE URITEo(IOUT.168) FORMAT (/." ENTER # RTOMS PER UNIT CELL: ".Z) RERD FREE (IN.ERR=ISO) NQPUC IF CNRPUC.LT.I .OR. NRPUC.GT.MQXRPUC) GOTO 158 ENTER HTOMIC NUMBERS RSSOCIQTED UITH ERCH BTOM IN UNIT CELL. (ORDER MUST REFLECT ’ENUMERQTION' OF DRTR FILE TO BE INPUT) URITE.(IOUT.218) FORMQT (/.“ ENTER RTOMIC NUMBER OF ERCH RTOM IN UNIT CELL'. /." (ORDER MUST CORRESPOND TO “.IH’."ENUMERRTION".2H’)) DO 223 I=1.NRPUC URITE-(IOUT,24B) I FORMAT (BX."RTOM #“.I2.4X.Z) RERD FREE (IN.ERR=238) NRTMNUM(I) IF {NQTMNUM(I).LT.1 .OR. NRTMNUMII).GT.BS) GOTO 238 Computer listings for energy band program 9O ‘ :UDD:EUBANK:EBAND.DIR:EBINPUT.FR C. C C. 25 268 288 298 278 0000 188 118 C. 35 368 378 000000 000 Figure C2. ENTER 3D K-VECTOR URITE-(IOUT.268) FORMAT (/." ENTER MAXIMUM K-VECTOR : “) DO 278 I=1.3 URITE (IOUT.298) KXYZ(I) FORMAT (BX.A2." = “.Z) READ FREE (IN.ERR=288) UVKMAX(I) IF (UVKMAX(I).LT.8.8 .OR. UVKMAX(I).GT.18.8) GOTO 288 ENTER # KPTS URITE-(IOUT.68) FORMAT (/." ENTER # PTS IN K-SPACE AT UHICH TO CALC EPBANDS: '.2) READ FREE (IN.ERR=58) KPTS IF (KPTS.LT.2 .OR. KPTS.GT.KKPTS) SOTO 58 ENTER ‘NEAR NEIGHBOR’ MAXIMUM DISTANCE FOR INCLUSION IN BLOCH FUNCTION SUM URITE-(IOUT.118) FORMAT(/." ENTER ".1H’.“NEAR NEIGHBOR'.1H’.' MAXIMUM DISTANCE '. /.” (ANGSTROMS) FOR INCLUSION IN BLOCH FCN: II.Z) READ FREE (IN.ERR=188) DMAX IF (DMAX.LT.8.9 .OR. DMAX.GT.48.8) GOTO 188 DMAX2=DMAX*DMAX :DISTANCE SOUARED ENTER FILENAME OF INPUT DATA FILE URITE-(IOUT.368) FORMAT(/." ENTER FILENAME OF FILE CONTAINING ATOM POSITIONS: ".Z) READ (IN.3?8.ERR=358) IFILE(1) FORMAT (868) OPEN 6.IFILE.ERR=358 CLOSE-6 URITE INPUT CONTROL PARAMETERS TO FILE. (TO SAVE CORE. READ IN CONTROL PARAMETERS FROM CONSOLE IN THIS STAND-ALONE PROGRAM. THEREBY AVOIDING ASSOCIATED FORMAT STATEMENTS IN EBAND.FR) OPEN 6."?EB6INPUT" URITE BINARY (6) KPTS.DMAX2.NAPUC.(NATMNUM(I).I=1.NAPUC). (UVKMAX(I).I=I.3).(IFILE(I).I=I.31) CLOSE.6 NOU CHAIN TO EBAND --- CALCULATE BANDS CALL FCHAN("EBAND.PR") END (Continued) 91 :UDD:EUBANK:EBAND.DIR:EBAND.FR EBAND.FR PROGRAM DOES THE FOLLOUING: I) GENERATES A SLATER-KOSTER MATRIX CONTAINING ELEMENTS BETUEEN BLOCH FUNCTIONS & ATOMIC ORBITALS II) DIAGONALIZES EACH MATRIX TO OBTAIN THE ENERGY EIGENVALUES AT A GIVEN UAVEVECTOR. (KX.KY.KZ) III) PLOTS THE ENERGIES VS UAVEVECTOR TO CONSTRUCT THE THE ENERGY BANDS. INPUT PARAMETERS: I) CONTROL PARAMETER FILE --- CONTAINING: A) NAPUC --- # ATOMS PER ’UNIT CELL“ B) NATMNUM --- ATOMIC NUMBER OF EACH ATOM IN 'UNIT CELL’. (ORDERED ACCORDING TO ’ENUMERATION’ (SEE BELOU)) C) UVKMAX --- 3D K-VECTOR (KX.KY.KZ) REPRESENTING MAXIMUM K-VECTOR D) KPTS --- 4 PTS IN K-SPACE AT UHICH TO CALCULATE THE ENERGY BANDS E) DMAXZ --- MAXIMUM ’NEAR NEIGHBOR’ DISTANCE SOUARED FOR ATOMS TO BE INCLUDED IN BLOCH FCN F) IFILE --- FILENAME OF INPUT DATA FILE (BELOU) FORMAT --- BINARY II) DATA FILE (’IFILE') --- EACH RECORD OF FILE CONTAINS: A) ATOMIC NUMBER 8) ATOM ENUMERATION IN ’UNIT CELL’ ’ENUMERATION’ IMPLIES THAT EACH ATOM IN ’UNIT CELL’ IS ASSIGN A UNIOUE NUMBER. 1 TO NAPUC (UHERE NAPUC IS THE # ATOMS PER ’UNIT CELL’) C) ATOM POSITION (X.Y.Z) IN RECTANGULAR COORDINATES FORMATTED AS FOLLOUS: (12X.I2.3X.12.3(3X.E14.?)) OUTPUT PARAMETERS: . I) DATA FILE --- "?EB+DATA“ --- CONTAINING CALCULATED EIGENVALUES. ETC BEING PASSED TO PLOTTING PROGRAM EXTERNAL REFERENCES: I) INPUT CONTROL PARAMETER FILE --- l"P‘EBéINPUT" GENERATED BY STAND-ALONE PROGRAM "EBINPUT.FR" II) INPUT DATA FILE --- ATOM POSITIONS B TYPES AJSUBROUTINE ATOMIN 1) COMMON/ATMPOS/.... III) GENERATE MATRIX ELEMENTS A) SUBROUTINE BLOCH 1) SUBROUTINE SSIAME A) SUBR SSIME 1) FCN ATTMV 000000000000000000000000000000000000000000000000000 Figure C2. (Continued) 92 :UDD:EUBAHK:EBAND.DIR:EBAND.FR 000000 000 00(700 C Figure C2. 2) SUBR UVECT 3) FCN DOTPD 2) COMMON/ATMPOS/.... IV) DIAGONALIZE A) SUBROUTINE EIGCH (IMSL LIBRARY) 1) EHBKH 2) EHOUHS.. 3) EORTZS. 4) UERTST V) PLOT ENERGY BANDS --- STAND-ALONE PROGRAM --- EBPLOT A) SUBROUTINE GRIDR B) SUBROUTINE CURVE C) SUBROUTINE TPOT D) SUBOUTINE OPLOT E) PLOT18 LIBRARY (TEKTRONIX) F) AGII LIBRARY (TEKTRONIX) PARAMETER MAXAPUC=18 :MAXIMUM # ATOMS PER UNIT CELLS PARAMETER NOPA=4 :# ORBITALS PER ATOM PARAMETER NDIM=NOPA*MAXAPUC :MAXIMUM # BANDS PARAMETER NNDIM=NDIM*(NDIM+1)/2. NNNDIM=3mNDIM DOUBLE PRECISION COMPLEX HA(NNDIM). VALUE DOUBLE PRECISION UORKSP(NNNDIM). E(NDIM) DIMENSION EE(NDIM).UVK(3).UVKMAX(3).IFILE(31).NATMNUM(MAXAPUC) EOUIVALENCE (EE(1).E(1)) READ INPUT CONTROL PARAMETER FILE. (TO SAVE CORE. READ IN CONTROL PARAMETERS FROM CONSOLE IN STAND-ALONE PROGRAM.THEREBY AVOIDING ASSOCIATED FORMAT STATEMENTS HERE) OPEN 6."?EBéINPUT" READ BINARY (6) KPTS.DMAX2.NAPUC.(NATMNUM(I).I=1.NAPUC). (UVKMAX(I),I=1.3),(IFILE(I).I=1.31) CLOSE.6 READ ATOMIC POSITIONS FROM DATA FILE CALL ATOMIN(IFILE.DMAX2.NAPUC.NATNNUM) TO SAVE CORE. DO PLOTTING IN SEPARATE STAND-ALONE PROGRAM. THEREFORE. NEED TO URITE EIGENVALUES TO A DISK FILE. CALL FDELETE ("?EBéDATA") OPEN 6. "?EB€DATA" NBANDS = NAPUC * NORA URITE BINARY (6) NBANDS.KPTS.(UVKMAX(I).I=1,3) (Continued) 93 :UDD:EUBANK:EBAND.DIR:EBAND.FR C C.. C 358 X378 . (310 (D - - (3)0 (000 Figure C2. ITERATE UAVENUMBER DO 388 KINCR=1.KPTS VARY K-VECTOR FROM (8.8.8) TO (KXMAX.KYMAX.KZMAX) DO 358 I=1.3 . UVK(I)=UVKMAX(I)*(KINCR-1)/(KPTS-1) URITE (12.378) (UVK(I).I=1.3) FORMAT(/.18X.’(KX.KY.KZ)=(’.2(F5.3.1H.).F5.3.1H)) FIND EACH MATRIX ELEMENT A(I.J) --- ON & BELOM DIAGONAL DO 488 IALPHA=1.NAPUC :UHICH ATOM (COL) DO 488 IMU=1.NOPA :UHICH ORBITAL (COL) DO 488 JBETA=IALPHA.NAFUC :UHICH ATOM (ROU) JNUMIN=1 IF (JBETA.EO.IALPHA) JNUMIN=IMU DO 488 JNU=JNUMIN.NOPA :UHICH ORBITAL (ROM) IATOM=NATMNUM(IALPHA) :ASSIGN ATOMIC NUMBER JATOM=NATMNUM(JBETA) :ASSIGN ATOMIC NUMBER CALL BLOCH (IALPHA.IMU.IATOM.JBETA.JNU.JATOM.UVK.VALUE) I=(IALPHA-1)*NOPA+IMU :I B J TAKE ON VALUES J=(JBETA-1)*NOPA+JNU :FROM 1-NBANDS (I.LE.J) A(I.J)=VALUE CONVERT 2-D ARRAY INTO l-D ARRAY USED BY EIGCH (REFERENCE PAGE INTRO‘IO OF IMSL MANUAL) L=(J*(J-1)/2) + I HA(L)=A(I.J) URITE (12.398) I.J.L.VALUE FORMAT(3H A(.I2.1H..I2.5H)=HA(.12.2H)=.2(DI4.7.4X)) HA(L)=VALUE NOU DIAGONALIZE CALL EIGCH (HA.NBANDS.8.E.DUMMY.DUMMY.UORKSP.IER) CONVERT DOUBLE PRECISION REAL TO SINGLE PRECISION REAL DO 588 I=1.NBANDS URITE (12.458) I.E(I) FORMAT(5X.2HE(.I2.2H)=.D14.7) EE(I)=E(I) URITE EIGEN VALUES TO DISK FILE URITE BINARY (6) (EECI).I=1.NBANDS) CLOSE.6 NOU PLOT CALL FCHAN("EBPLOT.PR“) :TC SAVE CORE.CHAIN END (Continued) 94 :UDD:EUBANK:EBAND.DIR:EBPLOT.FR 00000000000000000000000 0000 000 Figure C2. EBPLOT.FR ROUTINE TO MAKE THE GRAPHICS CALLS FOR EBAND.FR IE. PLOT ENERGY BANDS VS UAVENUMBER (K-VECTOR) INPUT PARAMETERS: DATA FILE --- CONTAINING THE FOLLOUING: 1) RECORD 1 --- #BANDS.#PTS/BAND.MAX K-VECTOR(3D) 2) RECORDS 2 THRU N --- EIGENVALUES FOR EACH BAND AT PARTICULAR K-VECTOR. (N = # PTS/BAND + 1) FORMAT --- BINARY OUTPUT PARAMETERS: PLOTFILE.“?EBGPLOTI EXTERNAL REFERENCES: FCN DOTPD SUBR GRIDR SUBR CURVE SUBR TPOT SUBR OPLOT- AGII LIBRARY PLOT18 LIBRARY PARAMETER IOUT=18. NPAGE=27*488K+12 PARAMETER NDIM=58.KKPTS=188 DIMENSION UVK(KKPTS).UVKMAX(3).UVKXYZ(3).EE(KKPTS.NDIM) READ IN DATA FOR BANDS FROM FILE OPEN 6. “?EB6DATA“ READ BINARY (6) NBANDS.KPTS.(UVKMAX(I).I-1.3) GENERATE UAVENUMBERS AND INPUT ENERGY EIGENVALUES DO 388 KINCR=1.KPTS READ BINARY (6) (EE(KINCR.J).J=1.NBANDS) VARY K FROM (8.8.8) TO (KXMAX.KYMAX.KZMAX) DO 358 I=1.3 UVKXYZ(I)=UVKMAX(I)*(KINCR-1)/(KPTS-1) NOU FIND MAGNITUDE OF UAVEVECTOR K UVK(KINCR)=SORT(DOTPD(UVKXYZ.UVKXYZ.3)) CONTINUE CLOSE.6 CONSTRUCT GRID CALL AMNMX (KPTS.UVK.XMIN.IDUMMY.XMAX.IDUMMY) YMIN=+1E12 YMAX=-1E12 DO 488 J=I.NBANDS (Continued) 95 :UDD:EUBANK:EBAND.DIR:EBPLOT.FR CALL AMNMX (KPTS.EE(1.J).EMIN.IDUMMY.EMAX.IDUMMY) IF (EMIN.LT.YMIN) YMIN=EMIN 488 IF (EMAX.GT.YMAX) YMAX=EMAX CALL GRIDR ($999.XMIN.XMAX.YMIN.YMAX."7EB6PLOT".37. "ENERGY BANDS --- ENERGY VS UAVEVECTOR".28. "ABSEKJ (l/ANGSTROMS)".11.“ENERGY -EV-". XMN.XMX.YMN.YMX.8) OPEN 7.“?EBéPLOT”.ATT="AO" URITE (7 .458) (UVKMAX(I).I=1.3) 458 FORMAT(33(/).9X.'K-VECTOR VARIES FROM ( .88. .88. .88) TO (’. & 2(F4.2.’.’).F4.2.’)’) CLOSE.7 9°90? NOU PLOT BANDS 000 D8 588 J=1.NBANDS 588 CALL CURVE ($999.KPTS.UVK.EE(1.J)."?EB¢PLOT".XMN.XMX. & YMN.YMX.8.8.8.8) C. URITE BINARY (IOUT) NPAGE :CLEAR CONSOLE SCREEN CALL TPOT (6."?EBGPLOT".IOUT) :COPY TO CONSOLE CALL OPLOT ("?EBéPLOT".IER) :COPY TO VERSATEC IF (IER.NE.1) GOTO 999 C. STOP.. C. C. 99 URITE (18.998) 998 FORMAT (" PLOTTING ERROR“) STOP.. C. END Figure C2. (Continued) 96 :UDD:EUBANK:EBAND.DIR:ATOMIN.FR Figure C2. SUBROUTINE ATOMIN(FILENAME.DMAX2.NAUC.NATMNUM) ROUTINE TO INPUT ATOM POSITIONS. ATOMIC NUMBER. AND ENUMERATION OF ATOMS IN UNIT CELL FOR BLOCH FCN SUM. ALSO. ASSIGN ’NEAR NEIGHBOR’ ATOMS TO EACH ATOM IN ’UNIT CELL’. INPUT PARAMETERS: 1) FILENAME --- DATA FILE CONTAINING ATOM POSITIONS.ETC EACH RECORD OF FILE CONTAINS: A) ATOMIC NUMBE B) ATOM ENUMERATION IN ’UNIT CELL’ C) ATOM POSITION (X.Y.Z) IN RECTANGULAR COORDINATES FORMATTED AS -~- (12X.I2.3X.I2.3(3X.E14.7)) NATOM --- ATOMIC NUMBER OF EACH ATOM NENUM -~- ENUMERATION OF EACH ATOM IN ’UNIT CELL’ FOR BLOCH FCN SUM. IE. ASSIGN EACH ATOM IN ’UNIT CELL’ A UNIQUE NUMBER. 1 THRU N. UHERE THERE ARE N ATOMS IN ’UNIT CELL’ R --- 3-D ARRAY FOR ATOM POSITION (X.Y.Z) ASSUMPTIONS MADE ABOUT INPUT DATA FILE: I) NEED TO INCLUDE AT LEAST AS MANY UNIT CELLS AS BLOCH SUM ’NEAR NEIGHBOR’ DISTANCE NECESSITATES 2) ALL ATOMS IN THE UNIT CELL AT THE ’ORIGIN’ MUST.BE.LISTED.FIRST. 2) DMAXZ --- DISTANCE SOUARED UITHIN UHICH ’NEAR NEIGHBORS’ UILL BE INCLUDED IN BLOCH FCN 3) NAUC -~- NUMBER OF ATOMS PER ’UNIT CELL’ 4) NATMNUM --- ARRAY OF ATOMIC NUMBERS CORRESPONDING TO EACH ATOM IN ’UNIT CELL’ (ORDER GIVEN BY ENUMERATION) OUTPUT PARAMETERS: COMMON/ATMPOS/AR(..).JENUM(..).NMAX(..) AR(IXYZ.J.K) --- GIVES ATOMIC COORDINATES (IXYZ=1.2.3) OF JTH ATOM IN ’UNIT CELL’ (ENUMERATION=J 8 K=1) ALONG UITH POSITIONS OF ALL ITS ’NEAR NEIGHBORS’ (K=2.....NMAX(J)) JENUM(J.K) --- GIVES ENUMERATION ASSOCIATED UITH EACH ATOM ABOVE NMAX(J) --~ NUMBER OF ’NEAR NEIGHBORS’ + 1 AROUND THE JTH ATOM IN THE ’UNIT CELL’ EXTERNAL REFERENCES: FUNCTION DOTPD DIMENSION FILENAME(I).R(3).D(3).NATMNUM(1) PARAMETER IOUT=18 PARAMETER MAXAUC=18 :MAXIMUM 4 ATOMS/UNIT CELL PARAMETER NNMAX=4O :MAXIMUM # NEAR NEIGHBLRS COMMON/ATMPDS/AR(3.HAXAUC.NNMAX).JENUH(MAXAUC.NNMAX).NMAKCMAXAUS) (Continued) 97 :UDDiEUSAXKIEBAND.DIR:ATOMIN.FR C.. - DO 1 I=1.NAUC 1 NMAX(I)=8 C.. OPEN 6. FILENAME.ATT="I".ERR=988 DO 5 IENUM=1.NAUC NFLAG=1 18 READ (6.28.END=38) NATOM.NENUM.(R(I).I=1.3) 28 FORMAT (12X.I2.3X.I2.3(3X.E14.7)) GOTO (48.58) NFLAG 48 IF (IENUM.NE.NENUM) GOTO 18 REUIND 6 NFLAG=2 GOTO 68 C.. 58 DO 55 I=1.3 55 D(I)=R(I)-AR(I.IENUM.1) D2=DOTPD(D.D.3) IF (D2.GT.DMAX2 .OR. D2.LT.8.81) GOTO 18 68 IF (NATOM.NE.NATMNUM(NENUM)) GOTO 928 NMAX(IENUM)=NMAX(IENUM)+1 DO 65 I=l.3 65 AR(I.IENUM.NMAX(IENUM))=R(I) JENUM(IENUM.NMAX(IENUM))=NENUM GOTO 18 38 IF (NMAX(IENUM).EO.8) GOTO 948 IF (NMAX(IENUM).GT.NNMAX) GOTO 968 REUIND 6 5 CONTINUE CLOSE.6 C.. RETURN C.. C. ERRORS- C.. 988 URITE-(IOUT.918) 918 FORMAT (" ERROR ON GEOMETRY FILE OPEN".Z) GOTO 999 928 URITE (IOUT.938) 938 FORMAT (" ERROR ON ATOMIC NUMBER".Z) GOTO 999 ' 948 URITE (IOUT.958) 958 FORMAT (“ ERROR ON ENUMERATION".Z) GOTO 999 968 URITE-(IOUT.978) 978 FORMAT (" ERROR --- TOO MANY NEAR NEIGHBORS".Z) 999 URITE-(IOUT.1888) 1888 FORMAT (" --- SUBR ATOMIN") STOP... END Figure C2. (Continued) 98 :UDD:EUBANK:EBAND.DIR:BLOCH.FR ("J0(")000000000000000000000000000000000000 C.. C C Figure C2. SUBROUTINE BLOCH (IALPHA.IMU.IATOM.JBETA.JNU.JATOM.UVK.VALUE) ROUTINE TO EVALUATE A MATRIX ELEMENT OF THE HAMILTONIAN BETUEEN 2 BLOCH FUNCTIONS OR (UITH NORMALIZATION) BETUEEN A BLOCH FUNCTION ORBITAL (LCAO) AND A SINGLE ATOMIC ORBITAL --- IE. EVALUATE (BLOCH FCN (STATES 2) EH] BLOCH FCN (STATES 1)) N*(BLOCH FCN (STATES 2) [H] ATOMIC ORB (STATE 1)> N* INPUT PARAMETERS: IALPHA --- ENUMERATION FOR ATOM # I --- IE. UHICH ATOM IN ’UNIT CELL’ TO USE FOR STATE 1 IMU --- ENUMERATION FOR ATOMIC ORBITAL # 1 --- IE. UHICH ATOMIC ORBITAL ON ATOM # 1 TO USE AS STATE 1 (S.PX.PY.PZ) --> (1.2.3.4) NOTE: IALPHA & IMU TOGETHER SPECIFY UHICH ATOMIC ORBITAL IN ’UNIT CELL’ IS STATE 1 IATOM --- ATOMIC NUMBER FOR ATOM # 1 JBETA --- ENUMERATION FOR ATOM # 2 --- IE. UHICH ATOM IN ’UNIT CELL’ TO USE FOR STATE 2 JNU --- ENUMERATION FOR ATOMIC ORBITAL # 2 --- IE. UHICH ATOMIC ORBITAL ON ATOM # 2 TO USE AS STATE 2 (S.PX.PY.P2) --> (1.2.3.4) NOTE: JBETA 8 JNU TOGETHER SPECIFY UHICH ATOMIC ORBITAL IN ’UNIT CELL’ IS STATE 2 JATOM --- ATOMIC NUMBER FOR ATOM 4 2 UVK -*- UAVEVECTOR AT UHICH TO EVALUATE THE BLOCH FCN (3-D ARRAY) OUTPUT PARAMETERS: VALUE --- EVALUATED BLOCH FUNCTION MATRIX ELEMENTS EXTERNAL REFERENCES: I) SUBROUTINE SSIAME A) SUBR SSIME 1) SUBR ATTMV 2) SUBR UVECT 3) FCN DOTPD II) FUNCTION DOTPD III) COMMON/ATMPOS/.... PARAMETER IAXAUC=18 :MAXIMUM # ATOMS/UNIT CELL PARAMETER NNMAX=48 :MAXIMUM # NEAR NEIGHBORS COMMON/ATMPOS/AR(3.MAXAUC.NNMAX).JENUM(MAXAUC.NNMAX).NMAX(MAXAUC) DIMENSION DD(3).UVK(3) COMPLEX CEXP.BLHEXP DOUBLE PRECISION COMPLEX VALUE VALUE = (8.8.8.8) SUM OVER ’NEAR NEIGHBOR’ ATOMS (Continued) :uOD: 000M 000 000 000 000 Figure CZ. E 99 (uANK:EBAND.DIR:BLOCH.FR DO 18 J=1.NMAX(IALPHA) INCLUDE ONLY ATOMS UITH CORRECT ENUMERATION IF (JENUM(IALPHA.J).NE.JBETA) GOTO 18 FIND ’SOLID STATE MATRIX ELEMENT’ OR ’ATOMIC TERM VALUE’ CALL SSIAME(AR(1.IALPHA.1).AR(1.IALPHA.J).IMU.JNU.IATOM.JATOM.SSME) . VECTOR DISTANCE BETUEEN ATOMS DO 28 I=1.3 DD(I)=AR(I.IALPHA.J)-AR(I.IALPHA.1) CALCULATE ARGUMENT OF EXPONENTIAL IN BLOCH FCN BLHEXP = (8.8.1.8) * DOTPD(UVK.DD.3) :DOT PRODUCT DETERMINE ONE TERM OF BLOCH FCN 8 ACCUMULATE TERMS IT! VALUE = VALUE + SSM *CEXP(BLHEXP) :COMPLEX EXPONENT CONTINUE RETURN END (Continued) 100 :UDD:EUBANK=UTILITY.DIR:SSIAME.FR SUBROUTINE SSIAME(R1.R2.NOR1.NOR2.NATOM1.NATOM2.SSME) ROUTINE TO EVALUATE A ’SOLID STATE INTERATOMIC MATRIX ELEMENT’ OF THE HAMILTONIAN BETUEEN THO ATOMIC ORBITALS OF ’STANDARD’- ORIENTATION---IE. S.PX.PY. OR PZ. IF ATOMIC ORBITALS ARE ON DIFFERENT ATOMS. USE THE SLATER-KOSTER 2-CENTER INTEGRAL APPROXIMATION ALONG UITH HARRISON’S PARAMETERS. IF ATOMIC ORBITALS ARE ON SAME ATOM. THEN MATRIX ELEMENT IS: 1) ZERO FOR DIFFERENT (ORTHOGONAL) ORBITALS 2) HERMAN-SKILLMAN ATOMIC TERM VALUE FOR SAME ORBITAL INPUT PARAMETERS: R1-~VECTOR POSITION IN RECTANGULAR COORDINATES OF ATOM #1 (ANGSTROMS) (3 ELEMENT ARRAY) R2--VECTOR POSITION IN RECTANGULAR COORDINATES OF ATOM #2 (ANGSTROMS) (3 ELEMENT ARRAY) NOR1---SPECIFIES ATOMIC ORBITAL #1 =1 IMPLIES S ORBITAL =2 IMPLIES PX ORBITAL =3 IMPLIES PY ORBITAL =4 IMPLIES PZ ORBITAL NOR2---SPECIFIES ATOMIC ORBITAL #2 =1 IMPLIES S ORBITAL =2 IMPLIES PX ORBITAL =3 IMPLIES PY ORBITAL =4 IMPLIES PZ ORBITAL NATOM1---ATOMIC NUMBER OF ATOM # 1 NATOM2---ATOMIC NUMBER OF ATOM # 2 OUTPUT PARAMETERS: SSME--EVALUATED SOLID STATE MATRIX ELEMENT EXTERNAL REFERENCES: SUBROUTINE SSIME DIMENSION R1(3).R2(3).OR1(4).OR2(4) 00 000000000000000000000000000000000000 CHECK INPUTTED ORBITALS IF (NOR1.LT.1 .OR. NOR1.GT.4) GOTO 28 IF (NOR2.LT.1 .OR. NOR2.GT.4) GOTO 28 C. C SENERHTE 4-D ORBITAL VECTORS REQUIRED BY SSIME D0 18 131:4 ORI(I)=8.B 18 UR2(I)=8.8 C. OR1(NORI)=I.8 OR2(NUR2)=1.8 C. Figure C2. (Continued) 101 :UDD:EUBANK:UTILITY.DIR:SSIAME.FR EVALUATE SOLID STATE MATRIX ELEMENT CALL SSIME (R1;R2.0R1,0R2,NATOM1,NATOM2;SSME) RETURN ERROR URITE (18,38) NOR1,NOR2 FORMAT (' ERROR---ILLEGAL ORBITAL SPECIFICATION---'. & uSUBROUTINE SSIAME";/.IBX,"NOR1='.13,5X;‘NOR2-“,13) PAUSE STOP.. END NNUDO (T) O DQ- . . . Figure C2. (Continued) 102 I :ODD:EhBANK:UTILITY.DIR:SSIME.FR (T100(")(")(")DUO0000000OOOOODOUUDOOODOOOOO DOD—- Figure C2. SUBROUTINE SSIME(R1.R2.0R1.0R2.NATOM1.NATOM2.SSME) ROUTINE TO EVALUATE A ’SOLID STATE INTERATOMIC MATRIX ELEMENT’ OF THE HAMILTONIAN BETUEEN TUO ATOMIC ORBITALS OF ANY ORIENTATION. IF ATOMIC ORBITALS ARE ON DIFFERENT ATOMS. USE THE SLATER-KOSTER 2-CENTER INTEGRAL APPROXIMATION ALONG UITH HARRISON’S PARAMETERS. IF ATOMIC ORBITALS ARE ON SAME ATOM. THEN MATRIX ELEMENT IS: 1) ZERO FOR DIFFERENT (ORTHOGONAL) ORBITALS 2) HERMAN-SKILLMAN ATOMIC TERM VALUE FOR SAME ORBITAL INPUT PARAMETERS: R1--VECTOR POSITION IN RECTANGULAR COORDINATES OF ATOM #1 (ANGSTROMS) (3 ELEMENT ARRAY) R2--VECTOR POSITION IN RECTANGULAR COORDINATES OF ATOM #2 (ANGSTROMS) (3 ELEMENT ARRAY) OR1--4 ELEMENT ARRAY TO SPECIFY ATOMIC ORBITAL #1 AND ITS DIRECTION (Sl.PX1.PY1.PZl) (DON’T UORRY ABOUT NORMALIZATION) EG. ’S’ ORBITAL ==> (1.6.8.6) EG. ‘P’ ORBITAL DIRECTED IN X-Y PLANE. 36 DEGREES ABOVE X-AXIS ==> (B.COSSO.SIN38.0) OR2--4 ELEMENT ARRAY TO SPECIFY ATOMIC ORBITAL #2 AND ITS DIRECTION (SZ.PX2.PY2.P22) (DON’T UORRY ABOUT NORMALIZATION) NATOMl--ATOMIC NUMBER OF ATOM #1 NATOM2--ATOMIC NUMBER OF ATOM #2 OUTPUT PARAMETERS: SSME--EVALUATED SOLID STATE MATRIX ELEMENT OR ATOMIC TERM VALUE EXTERNAL REFERENCES: SUBROUTINE UVECT FUNCTION DOTPD FUNCTION ATTMV DIMENSION R1(3).R2(3).OR1(4).OR2(4) DIMENSION D(3).P1(3).P2(3) DIMENSION UD(3).UP1(3).UP2(3) PARAMETER H2M=7.62. C1=-1.40. C2=-1.84. C3=-O.BI. C4=+4.0S SEPARATE ’8’ AND ’P’ ORBITALS SI=OR1(1) SZ=OR2(1) DO 1 I=1.3 P1(I)=OR1(I+1) P2(I)=OR2(I+1) CHECK INPUT PARAMETERS ORl & 0R2 TO VERIFY THAT ONLY AN ’5’ OR ’P’ TYPE ORBITAL IS SPECIFIED (NOT BOTH!!!) NSl=O (Continued) 103 :UDD:EUBANK:UTILITY.DIR:SSIME.FR F3001 - G 0000 OODH DUO 01") DOC—)0 C C 11 C. C 12 C. Figure C2. NS2=8 NP1=8 NP2=8 IF (Sl.NE.8.8) NSI=1 IF (SZ.NE.8.8) N82=1 DO 5 I=1.3 IF (P1(I).NE.8.8) NP1=1 IF (P2(I).NE.8.8) NP2=1 IF ( (NSl+NP1).NE.1 .OR. (NSZ+NP2).NE.1 ) GOTO 28 CALCULATE VECTOR D DO 18 I=1.3 D(I)=R2(I)-R1(I) CALCULATE MAGNITUDE OF VECTOR D SOUARED (=VECTOR DOT PRODUCT OF D UITH ITSELF) D2=DOTPD(D.D.3) FIND UNIT VECTORS CORRESPONDING TO D.P1.P2.Sl.& 82 UNIT VECTOR OF A ‘NULL’ VECTOR IS DEFINED TO BE A NULL VECTOR. CALLVUVECT(D.3.UD) CALL UVECT(P1.3.UP1) CALL UVECT(P2.3.UP2) CALL UVECT(Sl.1.USl) CALL UVECT(S2.1.USZ) DO R1 & R2 SPECIFY 2 DIFFERENT ATOMS (IE. D>.1 ANGSTROM)? IF 2 DIFFERENT ATOMS. CALCULATE "SOLID STATE MATRIX ELEM" IF SAME ATOM. ASSIGN "ATOMIC TERM VALUEll IF (D2.GT.(.81)) GOTO 18 FOR D<.1 ANGSTROM. MAKE SURE ATOMIC NUMBERS SAME IF (NATOM1.NE.NATOM2) GOTO 48 ARE ORBITALS ORTHOGONAL? IF SO. RETURN ZERO. IF NOT. CHOOSE TYPE OF ORBITAL (’3’ OR ’P’) SSME=8.8 SCALE=1.8 NORB=8 IF (N81.EO.1 .AND. NSZ.EO.1) NORB=1 :S ORBITAL IF (NP1.EO.1 .AND. NP2.EO.1) NORB=2 :P ORBITAL GOTO (19.12.11) (NORB+1) FOR P-ORBITALS. INCLUDE POSSIBILITY THAT ORBl & ORB2 ARE NOT STRICTLY ORTHONORMAL. SCALE=DOTPD(UP1.UP2.3) ASSIGN ATOMIC TERM VALUE SSME=SCALE$ATTMV(NATOM1.NORB) GOTO 19 (Continued) 104 :UDD=EUBAHK:UTILITY.DIR:SSIME.FR C. C 18 19 28 36 48 58 54 56 68 Figure C2. 909°? *** CALCULATE SOLID STATE MATRIX ELEMENT SSME=(H2M/D2)*( C1*USl*US2 + C2*(U82*DOTPD(UP1.UD.3)-USl*DOTPD(UP2.UD.3)) + C3*DOTPD(UP1.UP2.3) + C4*DOTPD(UP1.UD.3)*DOTPD(UP2.UD.3) ) RETURN ERROR ON INPUT PARAMETERS URITE (18.38) (OR1(I).I=1.4).(OR2(I).I=1.4) FORMAT (/./.” ONE AND ONLY ONE TYPE OF ATOMIC ORBITAL MUST BE". ” NON-ZERO"./.“ Sl=“.F?.3." PX1=“.F?.3.' PY1=”. F7.3." PZl=".F7.3./. “ 82=".F7.3." PX2=“.F?.3." PY2=“.F7.3.“ PZ2=“.F7.3) GOTO 68 URITE (18.58) FORMAT (/./." FOR ATOM POSITIONS IDENTICAL. ATOMIC NUMBERS ”. IIMUST BE SAME“./." ATOM 1 :".Z) URITE (18.54) NATOM1.(R1(I).I=1.3) FORMAT (“ ATOMIC NUMBER=".13." X=“.F?.3.“ Y=".F?.3.” Z='. F7.3) URITE-(18.56) FORMAT (" ATOM 2: ".Z) URITE (18.54) NATOM2.(R2(I).I=1.3) PAUSE " *** ERROR ON INPUT PARAMETERS *** SUBROUTINE SSIME *Xm' STOP.. END (Continued) 105 " :UDD:EUBANK:UTILITY.DIR:ATTMV.FR 00 00000000000000 (TIUIF) 18 28 38 35 48 45 58 Figure C2. FUNCTION ATTMV(NATOM.NORB) ROUTINE TO ASSIGN AN ENERGY FOR THE DIAGONAL MATRIX ELEMENTS OF THE HAMILTONIAN (ATOMIC TERM VALUES) FROM F. HERMAN & S. SKILLMAN. “ATOMIC STRUCTURE CALCULATIONS" (1963) USING THE HARTREE-FOCH-SLATER SELF-CONSISTENT ATOMIC FIELD APPROXIMATION. INPUT PARAMETERS: NATOM --- ATOMIC NUMBER OF ATOM UNDER CONSIDERATION NORB --- 1-->S ORBITAL:2-=>P ORBITAL OUTPUT PARAMETERS: ATTMV --- ATOMIC TERM VALUE EXTERNAL REFERENCES: NONE DIMENSION ATV(2.89) DOUBLE DIMENSIONED ARRAY STORED BY INCREMENTING LEFT MOST INDEX FASTER A DATA ATV/13.68..8..8..8.5.4B..8.B.1?.4.14.12.54.6.64. $17.52.B.9?.23.84.11.4?.29.14.14.13.35.88.16.99..8..8. &5.13..8.6.86.2.99.18.11.4.86.13.55.6.52.17.18.8.33.28.88.18.2?. $24.63.12.31.2*.8.4.19..8.5.41..8.5.B5..8.14*.8.6.92.1.83.8.4.3.3B. &11.37.4.98.14.3B.6.36.17.33.7.91.28.32.9.53.23.35.11.28..8..8. &3.S4..8.5.88..8.5.53..8.14*.8.6.41.2.85.7.78.3.38.18.12.4.69. $12.58.5.94.14.B8.7.24.1?.11.8.59.19.42.9.97..8..8.3.56..8.4.45..8. &4.86..8.42*.8.6.48.2.3B.?.68.3.4B.9.92.4.61.12.8?.5.??.14.15.6.9?. &16.21.B.19.18.24.9.44..8..8.3.48..8.4.24..8.4.63..8/ IF (NORB.LT.1.0R.NORB.GT.2) GOTO 5 IF (NATOM.LT.8.0R.NATOM.GT.B9) GOTO 18 IF (ATV(NORB.NATOM).EO.8.8) GOTO 18 ATTMV=-ATV(NORB.NATOM) RETURN URITE (18.6) NORB FORMAT (/./.'I ONLY S ORBITAL (NORB=1) OR P ORBITAL (NORB=2)". & "ALLOUED -—- NORB=".I1) GOTO 58 URITE (18.28) NATOM FORMAT (/./.“ ATOMIC TERM VALUE NOT AVAILABLE FOR:'. & " ATOMIC NUMBER=".I3.5X.Z) GOTO (38.48) NORB URITE-(18.35) FORMAT ("S ORBITAL“) GOTO 58 URITE (18.45) FORMAT ("P ORBITAL") PAUSE “ *** ERROR ON INPUT PARAMETERS xxx FUNCTION ATTMV xxx" STOP.. END (Continued) 106 :UDD:EUBANK:UTILITY.DIR:UVECT.FR (T) DOOOUUOOODOOU - (9 OF) 00(2).» ("J 28 C. 38 SUBROUTINE UVECT(X.N.XUNIT) ROUTINE TO GENERATE A N-DIMENSIONAL ’UNIT' VECTOR FROM A GIVEN VECTOR. IF INPUT VECTOR IS A ’NULL’ VECTOR. THEN THE ’UNIT’ VECTOR IS ALSO NULLED. INPUT PARAMETERS: X--VECTOR TO BE ’NORMALIZED’: A N-DIMENSIONAL ARRAY N--DIMENSIONALITY OF VECTOR OUTPUT PARAMETERS: XUNIT--UNIT VECTOR: A N-DIMENSIONAL ARRAY EXTERNAL REFERNCES: FUNCTION DOTPD DIMENSION X(1).XUNIT(1) DO 18 I=1.N XUNIT(I)=8.8 MAGNITUDE OF VECTOR SOUARED IS EOUAL TO THE DOT PRODUCT OF VECTOR UITH ITSELF D=SORT(DOTPD(X.X.N)) NULL VECTOR FOR D=8 IF (D.EO.8.8) GOTO 38 DO 28 I=1.N XUNIT(I)=X(I)/D RETURN END :UDDzEUBANK:UTILITY.DIR:DOTPD.FR DOC—)DUDOUODOO F) - (5) Figure C2. FUNCTION DOTPD(X1.X2.N) ROUTINE TO EVALUATE A N-DIMENSIONAL VECTOR DOT PRODUCT BETUEEN 2 VECTORS INPUT PARAMETERS: X1--VECTOR #1: A N-DIMENSIONAL ARRAY X2--VECTOR #2: A N-DIMENSIONAL ARRAY N--DIMENSIONALITY OF VECTORS X1 & X2 OUTPUT PARAMETERS: DOTPD--EVALUATED DOT PRODUCT EXTERNAL REFERENCES: NONE DIMENSION X1(1).X2(1) DOTPD=8.8 DO 18 I=1.N DOTPD=DOTPD+X1(I)*X2(I) RETURN END (Continued) i” APPENDIX D DENSITY 0F STATES COMPUTER PROGRAM An overall flowchart for the density of states program is shown in Figure 01. The basic structure of this density of states computer program is identical to the energy band program in Appendix C. The listings for this program, excluding subroutines which appear in Appendix C, are given in Figure D2. 107 108 Specific crystal Start ‘ geometry Generate Determine each trix element ._>__‘ random __>___ma , wave-vector 1" Equation 7 from Equation 14 Sort Diagonalize energy matrix to eigenvalues L—"‘—" obtain energy eigenvalues No Enough Yes statistics? Normalize and plot L .7 a) Figure D1. Flowchart for density of states program 109 :UDD:EUBANK:DSTATE.DIR:DSINPUT.FR nonnnnonnnnnnnnn00000000000 1‘) non 168 (1000C) 218 238 248 Figure 02. DSINPUT.FR --- DENSITY OF STATES ROUTINE TO ALLOU INPUTTING OF CONTROL PARAMETERS USED BY PROGRAM ’DSTATE.FR’ FROM CONSOLE AND GENERATES THE CONTROL PARAMETER FILE. INPUT PARAMETERS: FROM CONSOLE OUTPUT PARAMETERS: I) CONTROL PARAMETER FILE --- CONTAINING: A) NAPUC --- # ATOMS PER ’UNIT CELL’ B) NATMNUM --- ATOMIC NUMBER OF EACH ATOM IN ’UNIT CELL’. (ORDERED ACCORDING TO ’ENUMERATION’ (SEE BELOU)) C) EMIN & EMAX --- ENERGY MIN & MAX FOR ’UINDOU’ OF HISTOGRAM D) KPTS --- # PTS IN K-SPACE AT UHICH TO CALCULATE THE ENERGY EIGENVALUES E) DMAX2 --- MAXIMUM ’NEAR NEIGHBOR’ DISTANCE SOUARED FOR ATOMS TO BE INCLUDED IN BLOCH FCN F) IFILE --- FILENAME OF INPUT DATA FILE (BELOU) FORMAT --- BINARY EXTERNAL REFERENCES: NONE PARAMETER IN=11.IOUT=18.NPAGE=2?*488K+12 PARAMETER KKPTS=1E5 :MAXIMUM # KPTS ALLOUED PARAMETER MAXAPUC=18 :MAXIMUM # ATOMS PER UNIT CELL DIMENSION NATMNUM(MAXAPUC).IFILE(31) URITE (IOUT) NPAGE ENTER # ATOMS PER ’UNIT CELL’ URITE-(IOUT.168) FORMAT (/." ENTER # ATOMS PER UNIT CELL: n.Z) READ FREE (IN.ERR=158) NAPUC IF (NAPUC.LT.1 .OR. NAPUC.GT.MAXAPUC) GOTO 158 ENTER ATOMIC NUMBERS ASSOCIATED UITH EACH ATOM IN UNIT CELL. (ORDER MUST REFLECT ’ENUMERATION’ OF DATA FILE TO BE INPUT) URITE-(IOUT.218) FORMAT (/.“ ENTER ATOMIC NUMBER OF EACH ATOM IN UNIT CELL“. /." (ORDER MUST CORRESPOND TO “.1H’.“ENUMERATION'.2H’)) DO 228 I=1.NAPUC URITE-(IOUT.248) I FORMAT (8X."ATOM #“.12.4X.Z) Computer listings for density of states program 110 :UDD:EUBANK:DSTATE.DIR:DSINPUT.FR 228 258 268 278 288 mun—:00 (5)0. . 0000 188 118 C. 35 368 378 0000000 Figure D2. READ FREE (IN.ERR=238) NATMNUM(I) IF (NATMNUM(I).LT.1 .OR. NATMNUM(I).GT.BS) GOTO 238 CONTINUE ENTER ENERGY MIN & MAX URITE-(IOUT.268) FORMAT (/.“ ENTER ENERGY MINIMUM: ll.Z) READ FREE (IN.ERR=258) EMIN URITE-(IOUT.288) FORMAT (“ ENTER ENERGY MAXIMUM: I'.Z) READ FREE (IN.ERR=278) EMAX IF (EMIN.GE.EMAX) GOTO 258 ENTER # KPTS URITE-(IOUT.68) FORMAT (/.' ENTER # PTS 1N K-SPACE AT UHICH TO CALCULATE'./. " THE ENERGY EIGENVALUES: ”.Z) READ FREE (IN.ERR=58) KPTS IF (KPTS.LT.2 .OR. KPTS.GT.KKPTS) GOTO 58 ENTER ’NEAR NEIGHBOR’ MAXIMUM DISTANCE FOR INCLUSION IN BLOCH FUNCTION SUM URITE.(IOUT.118) FORMAT(/." ENTER “.1H’.“NEAR NEIGHBOR“.1H’.” MAXIMUM DISTANCE l'. /.“ (ANGSTROMS) FOR INCLUSION IN BLOCH FCN: ”.2) READ FREE (IN.ERR=188) DMAX IF (DMAX.LT.8.9 .OR. DMAX.GT.48.8) GOTO 188 DMAX2=DMAX*DMAX :DISTANCE SOUARED ENTER FILENAME OF INPUT DATA FILE URITE-(IOUT.368) FORMAT(/.“ ENTER FILENAME OF FILE CONTAINING ATOM POSITIONS: “.Z) READ (IN.3?8.ERR=358) IFILE(1) FORMAT (868) OPEN 6.IFILE.ERR=358 CLOSE.6 URITE INPUT CONTROL PARAMETERS TO FILE. (TO SAVE CORE. READ IN CONTROL PARAMETERS FROM CONSOLE IN THIS STAND-ALONE PROGRAM. THEREBY AVOIDING ASSOCIATED FORMAT STATEMENTS IN EBAND.FR) OPEN 6."?DSéINPUT" URITE BINARY (6) KPTS.DMAX2.NAPUC.(NATMNUM(I).I=1.NAPUC). EMIN.EMAX.(IFILE(1).I=1.31) CLOSE.6 (Continued) 111 :UDD:EUBANK:DSTATE.D1R:DSINPUT.FR C. C. C NOU CHAIN TO DSTATE --- CALCULATE DENSITY OF STATES C TO ALLOU FOR THE OPTION OF EXECUTING DSTATE.PR IN C BATCH AT OFF-HOURS. MAKE CHAINING OPTIONAL. C. 488 URITE (IOUT.418) 418 FORMAT (/." CHAIN? (YES/NO): ".Z) READ (IN.428.ERR-488) NANS 428 FORMAT (A2) IF (NANS.EO.’NO’) STOP IF (NANS.NE.’YE’) GOTO 488 CALL FCHAN(“DSTATE.PR”) C. END Figure D2. (Continued) 112 :UDD:EUBANK:DSTATE.DIR:DSTATE.FR DSTATE.FR --- DENSITY OF STATES PROGRAM DOES THE FOLLOUING: I) SELECTS RANDOM K-VECTORS II) GENERATES A SLATER-KOSTER MATRIX CONTAINING ELEMENTS BETUEEN BLOCH FUNCTIONS & ATOMIC ORBITALS FOR EACH K-VECTOR III) DIAGONALIZES EACH MATRIX TO OBTAIN THE ENERGY EIGENVALUES AT A GIVEN UAVEVECTOR. (KX.KY.KZ) IV) SORTS EIGENVALUES INTO ENERGY RANGE BUCKETS V) NORMALIZE DENSITY OF STATES VI) PLOTS THE DENSITY OF STATES VS ENERGY INPUT PARAMETERS: I) CONTROL PARAMETER FILE --- CONTAINING: A) NAPUC --- # ATOMS PER ’UNIT CELL’ B) NATMNUM --- ATOMIC NUMBER OF EACH ATOM IN ’UNIT CELL’. (ORDERED ACCORDING TO ’ENUMERATION’ (SEE BELOU)) C) EMIN & EMAX --- ENERGY MIN & MAX FOR ’UINDOU’ OF HISTOGRAM D) KPTS --- # PTS IN K-SPACE AT UHICH TO CALCULATE THE ENERGY EIGENVALUES E) DMAX2 --- MAXIMUM ’NEAR NEIGHBOR’ DISTANCE SOUARED FOR ATOMS TO BE INCLUDED IN BLOCH FCN F) IFILE --- FILENAME OF INPUT DATA FILE (BELOU) FORMAT --- BINARY II) DATA FILE (’IFILE’) --- EACH RECORD OF FILE CONTAINS: A) ATOMIC NUMBER B) ATOM ENUMERATION IN ’UNIT CELL’ ’ENUMERATION’ IMPLIES THAT EACH ATOM IN ’UNIT CELL’ IS ASSIGN A UNIQUE NUMBER. 1 TO NAPUC (UHERE NAPUC IS THE # ATOMS PER ’UNIT CELL’) C) ATOM POSITION (X.Y.Z) IN RECTANGULAR COORDINATES FORMATTED AS FOLLOUS: (12X.12.3X.I2.3(3X.E14.?)) OUTPUT PARAMETERS: I) DATA FILE --- ”?DS€DATA“ --- CONTAINING CALCULATED DENSITY OF STATES. ETC BEING PASSED TO PLOTTING PROGRAM EXTERNAL REFERENCES: I) INPUT CONTROL PARAMETER FILE --- '?EB+INPUT“ GENERATED BY STAND-ALONE PROGRAM 'EBINPUT.FR" II) INPUT DATA FILE --- ATOM POSITIONS & TYPES A)SUBROUTINE ATOMIN 1) COMMON/ATMPOS/.... 000000000000000000000000000000000000000000000000000 Figure D2. (Continued) 113 :UDD:EUBANK:DSTATE.DIR:DSTATE.FR 000000000000000000000000000000 0000000 Figure DZ. III) GENERATE RANDOM K-VECTOR A) SUBR TIME (SYSTEM ROUTINE) B) SUBR IRANDOM (SYSTEM ROUTINE) C) SUBR RANDOM (SYSTEM ROUTINE) IV) GENERATE MATRIX ELEMENTS A) SUBROUTINE BLOCH 1) SUBROUTINE SSIAME A) SUBR SSIME l) FCN ATTMV 2) SUBR UVECT 3) FCN DOTPD 2) COMMON/ATMPOS/.... V) DIAGONALIZE A) SUBROUTINE EIGCH (IMSL LIBRARY) 1) EHBKH 2) EHOUHS. 3) EORT2S 4) UERTST VI) NORMALIZE --- INTEGRATED DENSITY OF STATES - Q BANDS A) SUBROUTINE TRAPINT --- TRAPEZOIDAL INTEGRATION VII) PLOT DENSITY OF STATES VS ENERGY --- STAND-ALONE PROGRAM --- DSPLOT A) SUBROUTINE GRIDR B) SUBROUTINE CURVE C) SUBROUTINE TPOT D) SUBOUTINE OPLOT E) PLOT18 LIBRARY (TEKTRONIX) F) AGII LIBRARY (TEKTRONIX) PARAMETER MAXAPUC=18 :MAXIMUM # ATOMS PER UNIT CELLS PARAMETER NOPA=4 :# ORBITALS PER ATOM PARAMETER NDIM-NOPA*MAXAPUC :MAXIMUM # BANDS PARAMETER NNDIM=NDIM*(NDIM+1)/2. NNNDIM=3*NDIM ‘ PARAMETER NEPTS=288 : * BUCKETS FOR HISTOGRAM PARAMETER CONV=1888 :CONVERT RANDOM K-VECTOR FROM :8 THRU 1 TO 8 THRU 1888 (1/ANGST) DIMENSION IFILE(31) DIMENSION NATMNUM(MAXAPUC) DOUBLE PRECISION COMPLEX HA(NNDIM) DOUBLE PRECISION COMPLEX VALUE DOUBLE PRECISION UORKSP(NNNDIM). E(NDIM) DIMENSION G(NEPTS).UVK(3).KSEED(3).IHMS(3) READ INPUT CONTROL PARAMETER FILE. (TO SAVE CORE. READ IN CONTROL PARAMETERS FROM CONSOLE IN STAND-ALONE PROGRAM.THEREBY AVOIDING ASSOCIATED FORMAT STATEMENTS HERE) (Continued) 114 ‘ :UDDzEwBANK:DSTATE.DIR:DSTATE.FR OPEN 6."?DSGINPUT" READ BINARY (6) KPTS.DMAX2.NAPUC.(NATMNUM(I).I=1.HAPUC). & EMIN.EMAX.(IFILE(I).I=1.31) CLOSE.6 C.. C.. C READ ATOMIC POSITIONS FROM DATA FILE C.. CALL ATOMIN(IFILE.DMAX2.NAPUC.NATMNUM) C.. C INITIALIZE DENSITY OF STATES ARRAY C.. DO 188 J=1.NEPTS 188 G(J)=8.8 C.. C.. C GENERATE RANDOM ’SEEDS' FOR K-VECTOR C.. CALL TIME (IHMS.IER) DO 288 I=1.3 288 CALL IRANDOM(KSEED(I).IHMS(I)) C.. C GENERATE RANDOM UAVENUMBERS (KPTS) TIMES C. DO 388 KINCR=1.KPTS DO 358 I=1.3 CALL RANDOM(UVK(I).KSEED(I)) KSIGN=KSEED(I)/3+UVK(I)*CONV :IFIX CALL RANDOM(SIGN.KSIGN) IF (SIGN.GT.8.5) SIGN=+1.8 IF (SIGN.LE.8.5) SIGN=-1.8 358 UVK(I)=SIGN*CONV*UVK(I) C.. C.. C FIND EACH MATRIX ELEMENT A(I.J) --= ON a BELOU DIAGONAL C.. DO 488 IALPHA=1.NAPUC :UHICH ATOM (COL) DO 488 IMU=1.NOPA :UHIcH ORBITAL (COL) DO 488 JBETA=IALPHA.NAPUC :UHICH ATOM (ROU) JNUMIN=1 IF (JBETA.EO.IALPHA) JNUMIN=IMU DO 488 JNU=JNUMIN.NOPA :UHICH ORBITAL (ROU) IATOM=NATMNUM(IALPHA) :ASSIGN ATOMIC NUMBER JATOM=NATMNUM(JBETA) :ASSIGN ATOMIC NUMBER CALL BLOCH (IALPHA.IMU.IATOM.JBETA.JNU.JATOM.UVK.VALUE) I=(IALPHA-1)MNOPA+IMU :I a J TAKE ON VALUES J=(JBETA-1)*NOFA+JNU :FROM 1-NBANDS (I.LE.J) C A(I.J)=VALUE C CONVERT 2-D ARRAY INTO 1-D ARRAY USE- BY EIGCH C (REFERENCE PAGE INTRO-18 CF IMSL MANUAL) L=(J>k(J-I)/2) + I Figure DZ. (Continued) 115 :ODD:EUBANK:DSTATE.DIR:DSTATE.PR HA(L)=A(I.J) 8 HA(L)=VALUE NOU DIAGONALIZE NBANDS = NAPUC * NOPA CALL EIGCH (HA.NBANDS.8.E.DUMMY.DUMMY.UORKSP.IER) ACCUMULATE EACH EIGENVALUE IN APPROPRIATE BUCKET 0000 DO 588 I=1.NBANDS IF (E(I).ST.EMAX .OR. E(I).LT.EMIN) GOTO 588 J=IFIX( (E(I)-EMIN)*(NEPTS-1)/(EMAX-EMIN) + 1 ) C(J) = G(J) + 1 588 CONTINUE 388 CONTINUE NORMALIZE 0000 CALL TRAPINT (G.NEPTS.EMIN.EMAX.GI) DO 688 J=1.NEPTS 88 G(J)=G(J)*NBANDS/GI TO SAVE CORE. DO PLOTTING IN SEPARATE STAND-ALONE PROGRAM. THEREFORE. NEED TO URITE DENSITY OF STATES TO A DISK FILE. 00000001 CALL FDELETE ("?DSéDATA“) OPEN 6. “?DS€DATA" . URITE BINARY (6) NEPTS.EMIN.EMAX.(G(J).J=1.NEPTS) CLOSE.6 NOU PLOT 0000 CALL FCHAN('DSPLOT.PR") :TO SAVE CORE.CHAIN END Figure DZ. (Continued) 116 I:UDD:EUBANK:DSTATE.DIR:DSPLOT.FR 0000000000000000000000 0000 000 188 C. X Figure DZ. & DSPLOT.FR --- DENSITY OF STATES ROUTINE TO MAKE THE GRAPHICS CALLS FOR DSTATE.FR IE. PLOT DENSITY OF STATES VS ENERGY INPUT PARAMETERS: DATA FILE --- CONTAINING THE FOLLOUING: # ENERGY PTS. ENERGY MIN & MAX. DENSITY OF STATES ARRAY FORMAT --- BINARY OUTPUT PARAMETERS: RLOTFILE.“?DSGPLOT" EXTERNAL REFERENCES: FCN DOTPD SUBR AMNMX SUBR GRIDR SUBR CURVE SUBR TPOT SUBR OPLOT. AGII LIBRARY PLOT18 LIBRARY PARAMETER IOUT=18. NPAGE=27$488K+12 PARAMETER NNEPTS=288 DIMENSION E(NNEPTS).G(NNEPTS) READ IN DATA FOR DENSITY OF STATES FROM FILE OPEN 6. '?DS¢DATA" READ BINARY (6) NEPTS.EMIN.EMAX.(G(J).J=1.NEPTS) CLOSE-6 GENERATE ENERGY DATA ARRAY DO 183 J=1.NEPTS E(J)=((J-1)/FLOAT(NEPTS))*(EMAX—EMIN)+EMIN CONSTRUCT GRID CALL AMNMX (NEPTS.G.GMIN.IDUMMY.GMAX.IDUMMY) CALL GRIDR ($999.EMIN.EMAX.GMIN.GMAX.“?DSéPLOT“.38. I'ENERGY BAND DENSITY OF STATES".11.“ENERGY (EV)“. 17."DENSITY OF STATES".XMN.XMX.YMN.YMX.8) NOU PLOT BANDS CALL CURVE ($999.NEPTS.E.G.“?DSGPLOT".XMN.XMX. YMN.YMX.8.8.8.8) URITE BINARY (IOUT) NPAGE :CLEAR CONSOLE SCREEN (Continued) 117 :UDD:EUEANK:DSTATE.DIR:DSPLOT.FR X CALL TPOT (6.“?DSéPLOT".IOUT) :COPY TO CONSOLE CALL OPLOT ("?DS€PLOT".IER) :COPY TO VERSATEC IF (IER.NE.1) GOTO 999 C. STOP.. C. C. 99 URITE (18.998) 998 FORMAT (“ PLOTTING ERROR") STOP.. C. END ” :UDD:EUBANK:DSTATE.DIR:TRAPINT.FR SUBROUTINE TRAPINT (FX.NPTS.XMIN.XMAX.FXINT) ROUTINE TO DO TRAPIZOIDAL RULE INTEGRATION OF FUNCTION FX BETUEEN XMIN & XMAX GIVEN FX AT NPTS EOUALLY SPACED (IN X) POINTS INPUT PARAMETERS: FX --- ARRAY OF DEPENDENT VARIABLE VALUES AT EOUALLY SPACED INTERVALS OF INDEPENDENT VARIABLE VALUES NPTS --- NUMBER OF POINTS IN ARRAY FX XMIN & XMAX --- DEFINITE INTEGRAL ==> EVALUATE BETUEEN XMIN AND XMAX OUTPUT PARAMETERS: FXINT --- EVALUATED INTEGRAL EXTERNAL REFERENCES: NONE DIMENSION FX(1) 0 000000000000000 FXINT=(FX(1)+FX(NPTS))/2 DO 18 I=2.(NPTS-1) 18 FXINT=FXINT+FX(I) XDELTA=(XMAX-XMIN)/(NPTS-1) FXINT=XDELTA*FXINT RETURN END Figure DZ. (Continued) APPENDIX E PRINCIPAL INDICES OF REFRACTION COMPUTER PROGRAM The flowchart for the principal indices of refraction computer program is shown in Figure E1, and the listings are given in Figure EZ. Unlike the general applicability to any crystal exhibited by the energy band and density of states programs, portions of this program were specifically written for the chalcopyrite crystals. In partic— ular, the summation of polarization over each bond in the unit cell assumed a chalcopyrite structure. However, only two routines reflect this assumption. The first is the subroutine to input the crystal geometry (CPATMIN) and the second is the subroutine that computes the net polarization (POLARIZ). Also listed are two different routines (HBOMHYB and HBOMATM) for calculating Bond Orbital Model parameters V2, V3 and y. Some lower level subroutines, which were also used by the band structure programs, are not listed here but appear in Appendix C. 118 119 Specify geometry Select 3 Calculate orthogonal BOM ! ’ " electric , ' dipole moment fields for each bond Calculate 3 Sum to A components of obtain net susceptability """"“"‘ polarization tensor Diagonalize No Final tensor giving , electric , 3 . . 1 field? PI‘UCTP? . susceptibilities Calculate principal indices of refraction Figure E1. Flowchart for principal indices of refraction program Output _...‘__ (O‘DCD CHU(D~= C) (7 Figure EZ. C)F)F)03F)F)F)F)F)F)F)F)F)F)C)F)C)F)FJF)F30)03(301(30)?)F)F)F)F)F)F)F) 120 =UDD: EUEANK: DIE ECTRIC. DIR: PRINDEX. FR FRINDEX.FR MAIN PROGRAM TO DETERMINE THE PRINCIPAL INDICES OF REFRACTION FOR A CRYSTAL USING HARRISON’S BOND ORBITAL MODEL (BOM) APPROACH. THIS METHOD INVOLVES CALCULATING THE POLARIZATION FOR EACE. POND AND SUMMING UP THE CONTRIBUTIONS FROM ALL BONDS. BY CALCU!_ATING THE NET POLARIZATION AT 3 ORTHOGONAL ORIENTATIONS OF THE EXTERNAL ELECTRIC FIELD. THE NINE COMPONENTS OF THE SUSCEPTABILITY TENSOR ARE OBTAINED. DIAGONALIZATION CF THIS 3X3 MATRIX YIELDS THE -RINCIPAL SUSCErTASILITIES AND. CONSEOUENTLY. THE PRINCIPAL INDICES OF .EFRACTION. INPUT PARAMETERS (PPCH CONSOLE): IFILE *-- FILENAME OP FILE CONT.)IMI“G ATOM POSITIONS AND ATOMIC NUTICERS DESCRIBING CRYSTAL A.C -** LATTICE PARAMETERS OUTPUT PARAMETERS (TO CONSOLE): PRINCIPAL INDICES OP REPRAC TION EXTERNAL REFERENCES: SUBROUTINE PIIMPUT SUBROUTINE PIOUTPUT SUBROUTINE CPATMIN UHHOH/HTVPOS/ S'OROUTIHE POLPPIZ SUBR HSOHHYB OR HBOMATM SUER-DOTPD- COHMON/ATMPOS/.... SUBROUTINE EIGRP (IMSL LIBRARY) SUSP EORLAP SUBR EHESSP- SUSR-EQRH3P- SUBP EHOCKP- DUUP -ECOCKF- PARAMETEP PI=3.141592?. N3=3 DIMENSION JA TOM(3).I FILE(1§) DIMENSION E(N3) .P-OL(N3 .CVIEHV(“:.13) PIRfN3) DOUBLE PRECISION CHI(N3.N3 .LA (N3) DOUBLE PP ECISION UWPL-. ..-Piufla. -INPUT CONTROL IFFPPM’T'O! FROM CONSOLE NPLAG=~1 CQLL PI: “PUT (“FLT-1'3), .1" 15E: TIC) WET . ) INPUT ATOM POSITION FROM DATA FILE Computer listings for principal indices of refraction program 121 :UDD:EOS:NK2DIELECTRIC.DIR2PRINDEX. R C.. (7 {—3 (T) I\.\ - - ("5) Flt—)0 - G) F)(_)(')(—)H ('71???) (1000-1 (”)0 Figure E2. ki“ DO 10 I=I.do DO 28 J=I.N3 E(J)=0.0 E(I)=1.0 CALCULATE POLARIZATION F CALL POLARIZ(A:C,E,JATOM.POLJ CALCULATE SUSCEPTABILITY DO IO J=1;3 CHISAVCJII)=POL(J) CHI(J.IJ=POL(J) 3R EACH ELECTRIC FIELD VECTOR TENSOR COMPONENTS DIAGONALIZE SUSCEPTABILITY MATRIX TO GET PRINCIPAL SUSCEPTABILITIES CALL EIERF(CHI.N3.N3.0.PRCHI.DUM CALCULATE PRINCIPAL INDI DO 36 I=I.N3 PRC=PRCHI(I) PIR(I)=SORT(1.0+4*PI*PRC) OUTPUT.. MY,DUMMY,UK.IER) .q-fi C;a OF REFRACTION CALL PIOUTPUTICHISAV.“RCHI.PIR.IFILE(4)) GOTO 1 END (Continued) FIFIF)FICIF)F)F)F)F)F)FJC)F)F) ... t-fi U! f.) (... (T) **L9Ku- 636) - Figure E2. rn W‘LEEI‘.1I\-J; :ELECTP ICoD SUBROUTINE PIINRU ROUTINE TO INPUT INPUT PRRRMETERS: NFLQG --- 122 ER: PIE”: UT FR T(NFLRG;IFILEIR;C»JRTOM) PRRRRETERS FROM CONSOLE FOR PRINDEX.FR .EO.-1 ==> FIRST TIME THRU OUTPUT PRRRMETERS: NFLRG --- IFILE -—- 4 R E L “ JRTON --- EXTERNRL REFERENC NONE .E0.0 => DON’T NEED TO RERD RTOM POSITION FILE .EO.1 => NEED TO RERD RTOM POSITION FILE FILENRME OF FILE CONTRINING RTOM POSITIONS LRTTICE CONSTRNTS RRRRY OF RTOMIC NUMBERS FOR ERCH OF THE 3 TYPES OF RTOMS IN UNIT CELL OF CHRLCOPYRITE ES: DIMENSION IFILE(16).JRTOM(3) PRRRMETER IN=11, IOUT=IO, NRRGE=O?*4OOK+12 IF (NFLRG.NE.-1) GOTO 2 MRITE (IOUT) NRREE URITE- C-IOIJTII) FORNRT (/,251,"CHTLCOPYRITE INDEX OF REF RR MIDI") GOTO 99 URITE-(IOUTp3) FORNRT {/," DONE? RERD (IN 4 ERR =2) FORNRT (92) IF (IRNS.EO.‘YE’) IF (IRNS.NE.’NO’) NFLRS=O (YES/NO}: "32) IRNS STOP GOTO 2 URITE (IOUT) NPRGE HRIi E- ’14,], O) F05 WHT (/: /;/ " CHRIGE CDYSTQL LLQH:TTYT (.’ES/IIU)€ ".Z) PEGD ( TN?B!EI\P=16) IFII‘IS FERHRT LRZ) IF (IRNS.EO.’NO’F J- \,."') IE (IFHIS.IIE. ac INPUT LRTTICE CONS NFL? :G=I URITE-‘ICU .18) FORKRT (/ " ENIEr RERD FREE (IN.:3F= IF (R.LT.4.0 .OR. URZTE QIOUTIISO) FC'I’WPI.\T ‘l-‘E\I:’ ”1': ._. II RE.3D H EE {INIERR= Ir 1" L‘.” ‘ r71 (“D 0 Kb. InLu'.‘ ob‘l" (Continued) GOTO 253 GOTO 10 U3 LRTTICE CCNSTQNTS: "./.5x."9 = ',:) 168) 9 Q.GT.?.O) GOTO 163 Z) - n. 1.0: 3" f‘ C .GT.‘3.0F SOTO ISO 123 :UDD:EUBRHK:DIELECTRIC.DIR:PIINPUT.FR 000 ISO 169 198 195 "1 CD DJUIDFJUt-m GNU!- . (5)63. NUUJ 01- ’\l (3). CS) Figure E2. INPUT RTOMIC NUHSERS URITE-(IOUTIISO) FORMQT (/," SPECIFY 3 TYPES OF RTOMS IN UNET CELL IN ORDER") DO 178 I=1,3 URITE-CIOUT,195) i. FORHQT C" ENTER RTOMIC NUMBER FOR QTOM #",Il,“: ",Z) RERD FREE (IN,ERR=195) JRTOM(I} 1F (JQTOM(I).LT.5 .OR. JQTOM(I).GT.89) GOTO ISO CONTINUE ENTER FILENQHE Cl ‘Tl INPUT DRTQ FIL URITE-CIOUTy368) FORMQT(/," ENTER FILENRNE OF FILE CONTHINUNG RTOM POSITIONS: ",Z) RERD (IN,3?O,ERR=3SO) IFILE(1) FORMQT (868) RETURN END (Continued) 124 :UDD:EUBRN£zDIELECTRIC.DIR:CPRTMIN.FR SUBROUTINE CPQTMIN(IFILE,Q,C,JQTOMa$) C.. C ROUTINE TO INPUT RTOM POSITIONS 9ND QTOMIC NUMBERS C.. C INPUT PQRRMETERS: C 1) IFILE --- DRTQ FILE CONTQINING RTOM POSITIONS;ETC C ERCH RECORD OF FILE CONTRINS: C Q) QTOMIC NUMBER C B) QTOM POSITION (XIY,Z) IN RECTQNGULQR C COORDINRTES C FORMRTTED RS --- (12X;12,5X,3(3XpE14.7)) C NRTOM --- RTOMIC NUMBER OF ERCH RTOM C R --- 3-D RRRQY FOR QTOM POSITION (X’Y,Z) C 2) R & C --- LRTTICE CONSTRNTS C 3) JRTOM --- 3D RRRQY OF RTOMIC NUMBERS C OUTPUT PRRRMETERS: C ERROR RETURN C EXTERNRL REFERENCES: C COMMON/RTMPOS/R,NRTOM C.. PRRRMETER GTONE=1.BI COMMON/RTMPOS/R(3,NRUC),NRTOM(NRUC) DIMENSION IFILE(1),RR(3),JRTOM(3) PRRRMETER NRUC=31 :# RTOMS IN UNIT CELL (INCLUDE C :FULL RTOM RT EDGES) J=O OPEN 6; IFILE,RTT=“I",ERR=SOB IO RERD (6,28,END=40) IRTOM.(RR(I);I=1,3) 2O FORMQT (12X,12;5X’3(3X.E14.7)) IF (RBS(RR(1)).GT.(R/2*GTONE)) GOTO 16 IF (RBS(RR(2)).GT.(R/2*GTONE)) GOTO 18 IF (QBS(RR(3)).GT.(C/2*GTONE)) GOTO 18 DO 25 I=1;3 25 IF (IRTOM.EO.JRTOM(I)) GOTO 27 GOTO 986 2? J=J+1 IF (J.GT.NRUC) GOTO 988 NRTOM(J)=IRTOM DO 38 I=1;3 3O R(I,J)=RR(I) GOTO IO 48 CLOSE 6 C.. RETURN C.. C. ERRORS- C.. 98 RETURN 5 END Figure E2. (Continued) 125 :UDD:EUBRNK:DIELECTRIC.DIR:POLRRIZ.FR UOUOOOUUUUUOUUODUU DOC-11') DOC—)0] (Tint—2‘04 Figure E2. SUBROUTINE POLQRIZ(R:C,E,JRTOM,POL) PROGRRM TO CRLCULRTE THE POLRRIZRTION FOR 9 GIVEN ELECTRIC FIELD DIRECTION FOR 9 CHQLCOPYRITE CRYSTRL USING HRRRISON’S BOND ORBITRL MODEL (BOM) INPUT PRRQMETERS: R & C --- LRTTICE CONSTQNTS E --- 3D RRRQY FOR THE ELECTRIC FIELD VECTOR JQTOM --- RRRRY OF RTOMIC NUMBERS FOR ERCH OF THE 3 TYPES OF RTOMS IN UNIT CELL OUTPUT PRRRMETERS: POL --- 3D QRRRY FOR THE POLRRIZRTION VECTOR EXTERNRL REFERENCES: FCN DOTPD SUBR HBOMHYB OR HBOMRTM COMMON/QTMPOS/... PRRRMETER EE=3.?949; N3=3 DIMENSION PO(N3);P1(N3),P2(N3) DIMENSION JRTOM(1),E(1),D(N3).POL(1J PRRRMETER NQUC=31 :# QTOMS IN UNIT CELL (INCLUDE :FULL RTOM RT EDGES) COMMON/QTMPOS/R(3,NRUCJ»NRTOM(NRUCJ QNION DISPLRCED FROM CENTER OF TETRRHEDRRL BY DELX DELX=(.25-SORT(C*C/32/H/Q-I./16.))*9 NEQREST NEIGHBOR DISTRNCE SOURRED DNN2=(R*R/8+C*C/64+DELXmDELx+anELX/2)*1.3 INITIRLIZE POLRRIZRTION VECTOR DO 5 K=1,N3 PO(K)=0.0 P1(K)=0.0 P2(K)=0.0 CHOOSE RNION --- 8 PER UNIT CELL IRNEXT=1 DO 18 IIQ=1;8 DO 26 I=IRNEXT,NQUC IR=I IF (NRTOM(IR).EO.JRTOM(N3)) GOTO 38 PRUSE " CRNNOT FIND RNION" STOP... IRNEXT=IR+1 CHOOSE CRTION --- 4 PER ERCH RNION ICNEXT=1 DO 48 IIC=1,4 (Continued) 126 :UDD:EUSRNK:DIELECTRIC.DIR:POLRRIZ.FR 45 - (9 00000 000000001 78 48 IO C88 OXXXXCO - (SI Figure E2. DO SO I=ICNEXTINRUC IC=I DO 45 K=IIN3 DfK)=R(K,IC)-R(K.IQJ D2=DOTPD(D;D;N3) IF ((D2.LT.DNN2) .QND. (D2.GT.0.0I)) GOTO SO PQUSE " CRNNOT FIND CQTION“ STOP.... ICNEXT=IC+I CRLCULRTE BOM PQRRMETERS V2;V3 & GRMMQ FOR BOND BETUEEN RNION SP3 HYBRID 9ND CRTION SP3 HYBRID USE DECOMPOSED SP3 HYBRIDS RND HRRRISON’S SOLID STRTE MQTRIX ELEMENTS FOR RTOMIC ORBITRLS RLONG UITH HERMQN- SKILLMRN’S RTOMIC TERM VBLUES CRLL HBOMRTM(R(I,IR),R(1,IC).NRTOM(IR),NRTOM(IC),V2,V3,GRMMR) USE HRRRISON’S EMPIRICRL BOM PQRRMETERS FOR BINHRY CMPDS CRLL HBOMHYB(NRTOM(IC),NRTOM(IR),V2;V3,GRMMR) VECTOR SUM OF POLRRIZRTION FROM ERCH BOND IN UNITCELL SV23=SORT(V2*V2+V3*V3) DO 76 K=1,N3 PO(K)=PO(K)*(—GRMMR*EE*V3)*D(K)/SV23 P1(K)=P1(K)+((GnmnnmEEmv2)**2)*DOTPD(E,D,N3)*D(K)/2/(SV23**3) P2(K)=P2(K)+3*(GRMMR*EE)**3*V2*V2*V3*D(K)*DOTPD(E,D;N3)/8/(SV23**S) CONTINUE CONTINUE CQLCULRTE 8 NORMQLIZE POLRRIZRTION DO SO K=1»N3 IF (RBS(PO(K)).GT.IE-6) URITE FREE (12) " FERROELECTRIC? " POL(K)=(PO(KJ+PI(K)+P2(K))/(R*Q*C) :INCLUDE 2ND ORDER TERM POL(K)=(PO(K)+PI(K))/(Q*R*C) :FIRST ORDER ONLY URITE FREE (12) " PO URITE FREE (12) " PI ",(PI(I),I=I,N3) URITE FREE (12) " P2 ",(P2(I),I=I,N3) URITE FREE (12) " POLRRIZQTION VECTOR = ",(POL(I),I=1,N3) "p(PO(I);I=IIN3) RETURN END (Continued) 127 :UDD:EUBHNK:DIELECTRIC.DIR:HBOMRTM.FR 0000000 0 H (S) Figure E2. SUBROUTINE HBOMQTM(RI,R2,NRTOMIJNQTOMZIVZ,V3IGRMMR) ROUTINE TO CRLCULRTE HQRRISON’S BOND ORBITQL MODEL (BOM) PRRRMETERS FOR THE DIELECTRIC FUNCTION FROM RTOMIC ORBITRLS USING HIS ’SOLID STRTE INTERRTOMIC MRTRIX ELEMENTS’ 9ND HERMQN-SKILLMQN’S QTOMIC TERM VRLUES. INPUT PQRRMETERS: RI~-VECTOR POSITION IN RECTRNGULRR COORDINRTES OF RNION (RNGSTROMS) (3 ELEMENT RRRRY) R2--VECTOR POSITION IN RECTRNGULRR COORDINRTES OF CRTION (RNGSTROMS) (3 ELEMENT RRRQY) NRTOMI---RTOMIC NUMBER OF ONION NRTOM2---RTOMIC NUMBER OF CRTION OUTPUT PQRQMETERS: V2,V3,GRMMQ --- BOM PQRRMETERS NEEDED FOR DIELECTRIC FCN EXTERNRL REFERENCES: SUBROUTINE SSIQME SUBROUTINE SSIME SUBROUTINE QTTMV SUBROUTINE UVECT SUBROUTINE DOTPD NOTES: NORBI---SPECIFIES QTOMIC ORBITRL OF RNION =1 IMPLIES S ORBITRL =2 IMPLIES PX ORBITQL =3 IMPLIES PY ORBITRL =4 IMPLIES P2 ORBITQL NORB2---SPECIFIES NTOMIC ORBITRL OF CRTION =1 IMPLIES S ORBITRL =2 IMPLIES PX ORBITRL =3 IMPLIES PY ORBITRL =4 IMPLIES PZ ORBITAL PRRRMETER SQR3=I.732851 DIMENSION R1(3),R2(3),D(3),SR(4),SC(4) GRMMQ=1.0 DETERMINE V2 FROM HRRRISON’S "SOLID STRTE INTERRTOMIC MQTRIX ELEMENTS“ FOR RTOMIC ORBITQLS. DECOMPOSE SP3 HYBRID ORBITQLS OF RNION & CQTION INTO RESPECTIVE LINEQR COMBINQTIONS OF RTOMIC ORBITRLS DISTRNCE BETUEEN RNION & CQTION DO 18 I=1,3 D(I)=R2(I)-RI(I) DD=SORT(DOTPD(D,D,3)) (Continued) 128 :UDD:EUBONK:DIELECTRIC.DIR:HBOMOTM.FR C.. C SET SCOLE FOCTORS & SIGNS OF ONION O COTION OTOMIC C ORBITOL TERMS THOT MOKE UP HYBRID C.. SO(1)=DD/SOR3 SC(I)=-DD/SOR3 DO 26 I=2,4 SO(I)=D(I-1) 26 SC(I)=SO(I) Cu. C ITEROTE OVER ONION & COTION OTOMIC ORBITOLS C OND OCCUMULOTE MOTRIX ELEMENTS C.. V2=6.6 DO 46 NORBI=II4 DO 46 NORB2=II4 COLL SSIOME(RI,R2,NORBI,NOR82,NOTOMI,NOTOMZISSME) V2=V2+SO(NORBI)*SC(NORBZ)*SSME 46 CONTINUE V2=3*V2/4/DD/DD COLCULOTE V3 FROM OTOMIC TERM VOLUES FOR OTOMIC ORBITOLS. DECOMPOSE SP3 HYBRIDS INTO OTOMIC ORBITOLS 00000 OTOMIC TERM VOLUE FOR ONION S-ORBITOL COLL SSIOME(RI,R1,I,I,NOTOMI,NOTOMI;EOS) OTOMIC TERM VOLUE FOR ONION P-ORBITOL COLL SSIOME(RI;R1,2,2;NOTOMI;NOTOMI,EOP) OTOMIC TERM VOLUE FOR COTION S-ORBITOL COLL SSIOME(R2.R2,I.IINOTOM23NOTOMZIECS) OTOMIC TERM VOLUE FOR COTION P-ORBITOL COLL SSIOME(R2,R2,2;2,NOTOM2,NOTOM2,ECP) 0 0 V3=(ECS+3*ECP-EOS-3*EOP)/8 000 RETURN END Figure E2. (Continued) 129 :UDB:EUSONK:DIELECTRIC.DIR:HBOMHYB.FR SUBROUTINE HBOMHYB (NOTOMC:NOTOMO:V2:V3:GOMMO) ROUTINE TO OBTOIN HORRISON’S BOND ORBITOL MODEL (BOM) POROMETERS FOR O GIVEN COTION OND ONION FROM BINORY CMPDS EMPIRICOL NUMBERS INPUT POROMETERS: NOTOMC --- OTOMIC NUMBER OF COTION NOTOMO --- OTOMIC NUMBER OF ONION OUTPUT POROMETERS: V2 --- COVOLENT ENERGY V3 --- POLOR ENERGY GOMMO --- OVERLOP INTEGROL:INTERNOL FIELD: ETC EXTERNOL REFERENCES: NONE 00000000000000 DIMENSION EH(42) DOTO EH/3.BB:4.74:5.76:7.16:8.12:8.76:II*6.6: & 3.76:3.4B.4.44:5.86:6.86,B.66:9.36:11*6.6: & 3.34.3.36:4.28:5.36:6.47:7.46:8.64/ IF (NOTOMC.LT.I2 .OR. NOTOMC.GT.52) GOTO 999 IF (NOTOMO.LT.I2 .OR. NOTOMO.GT.52) GOTO 999 DETERMINE V2 t’) (”I 0 V2C=V20=2.2 :ROU 3 IF (NOTOMC.GE.19) V2C=2.15 :ROU 4 IF (NOTOMO.GE.19) V20=2.15 IF (NOTOMC.GE.37) V2C=I.?6 :ROU 5 IF (NOTOMO.GE.3?) V20=I.76 V2=SORT(V2C*V20) DETERMINE GOMMO 0(‘J0 GC=GO=1.2 :ROU 3 IF (NOTOMC.GE.19) GC=I.42 :ROU 4 IF (NOTOMO.GE.19) GO=1.42 IF (NOTOMC.GE.37J GC=1.69 :ROU 5 IF (NOTOMO.GE.37) GO=I.69 GOMMO=SORT(GC*GO) DETERMINE V3 (TI00 V3=-(EH(NOTOMC-II)-EH(NOTOMO-II))/2 RETURN ERRORS- @000 SS POUSE " ERROR --- SUBROUTINE HBOMP ...." STOP I n . . END Figure E2. (Continued) 130 :UDD:EUSONK:DIELECTRIC.DIR:PIOUTPUT.FR SUBROUTINE PIOUTPUT(CHI.PRCHI:PIR:IXSTOL) ROUTINE TO OUTPUT THE PRINCIPOL INDICES OF REFROCTION. OLSO: SPECIFIES ISOTROPIC:UNIOXIOL: OR BIOXIOL OND OPTIONOLLY CON OUTPUT THE COMPLETE SUSCEPTOBILITY TENSOR OND THE PRINCIPOL SUSCEPTOBILITIES. INPUT POROMETERS: CHI --- 3X3 MOTRIX FOR THE SUSCEPTOBILITY TENSOR PRCHI --- 3D ORROY FOR THE PRINCIPOL SUSCEPTOBILITIES PIR --- 3D ORROY FOR THE PRINCIPOL INDICES OF REFROCTION IXSTOL --- OSCII CHOROCTER STRING DESCRIBING CRYSTOL OUTPUT POROMETERS: ONY OR OLL OF THE INPUT POROMETERS EXTERNOL REFERENCES: NONE 0000000000000000 POROMETER N3=3: SN=IE-3. IOUT=12 DOUBLE PRECISION COMPLEX PRCHI(N3J DIMENSION CHI(N3:N3), PIR(N3): IXSTOL(I) OUTPUT SUSCEPTOBILITY TENSOR URITE (IOUTI4) ((CHI(I,J).J=1:N3),I=I:N3) 4 FORMOT (/:/:IBX."SUSCEPTOBILITY TENSOR": 3(/:5X:IH(:3(E12.5,3X):IH)):/) OUTPUT PRINCIPOL SUSCEPTOBILITIES "(000XXX000 90 d URITE (IOUT:6) (I:PRCHI(I).I=I.N3) X6 FORNRT (/,/,18X,"PRINCIPRL SUSCEPTOBILITIES", X & 3(/,11X,“CHI¢",I1,5X,E12.S," + I*",E12.S)) C.. C OUTPUT PRINCIPOL INDICES OF REFRQCTION C.. X21=OBS(PIR(2)-PIR(1)) X31=OBS£PIR(3)-PIR(1>) X32=QBSIPIR(BB-PIR(2)) NTYPEaS :OIRXIRL IF (X21.LT.SN) NTYPE=2 :UNIRXIRL IF (X31.LT.8N) NTYPE=3 :UNIOXIRL IF (X32.LT.SN) NTYPE=4 :UNIOXIOL IF ((X21.LT.SN) .RND. (X31.LT.SN) .QND. 'X32.LT.8N)) & NTYPE=1 :ISOTROPIC URITE (IOUT,8) IXSTOLCI) O FORMOT (/,/,22X,824,/) GOTO (18:28:22,24,30) NTYPE C.. IO URITE-CIOUT,15) PIRCI) 15 FORNOT (23X,"ISOTROPIC CRYSTOL"./,1SX:"INDEX OF REFRRC", & "TION = ",F6.3:/) Figure E2. (Continued) 131 :UDD:EUBONK:DIELECTRIC.DIR:PIOUTPUT.FR GOTO 46 C.. 26 NO=I NE=3 GOTO 26 22 NO=I NE=2 GOTO 26 24 NO=2 NE=1 26 BIREF=PIR(NE)-PIR(NO) URITE (IOUT:28) PIR(NO):PIR(NE):BIREF 28 FORMOT (26X:“UNIOXIOL CRYSTOL":/:12X:"ORDINORY INDEX OF REF": & llROCTION = “:F6.3,/,12X:"EXTROORDINORY INDEX OF REFROCTION = ': & F6.3,/:17X,"BIREFRINGENCE = ":F6.4,/) GOTO 46 C.. 3 URITE (IOUT,35) (I’PIR(I):I=I:N3) 35 FORMOT (26X:"BIOXIOL CRYSTOL"./:ISX:”THE THREE PRINCIPOL “: & "INDICES OF REFROCTION ORE:“:3(/:IBX,"#",II:5X:F6.3),/) C.. 46 RETURN END Figure E2. (Continued) APPENDIX F CHALCOPYRITE CRYSTAL STRUCTURE The tetragonal structure of the chalc0pyrites closely resembles the cubic zincblende structure. For example, by replacing the Ga atoms in the binary compound, GaAs, with Cd and Ge atoms, the ternary chalcopyrite, CdGeAsz, is obtained (see Figure 1 of Lines and Naszazak18). The space group is 152d and the point group is 12m.29 The lattice parameters a and c determine the distortion of the tetrahedron; that is, c = 2a implies a "regular" tetrahedron and c>2a means that the tetrahedron is stretched along the z axis or com- pressed along the x and y axes. Lattice parameter x describes the displacement in the x-y plane of an anion from the center of a tetra- hedron, formed by four cations, one at each corner. Specifically, the anion is displaced by an amount, alx-.25l, towards the smaller of the cations; for example, in CdGeAsZ, the As anion is closer to the two Ge cations. The lattice constants for many of the chalcopyrites have been measured and/or predicted by Abrahams and Bernstein.29’30’31 Also, they obtained an empirical equation for lattice parameter x which is given by, 2 x=-§--‘/C2-—% (F1) 32a To generate the crystal structure in crystal coordinates (which are also rectangular), the following formulas32 could be used (with CdGeA52 again as the example): 132 133 Cd: [(o,o, o;%, %, é) + (0, o, o ,0,%, i)] (a, a, c) + (la, ma, nc) (F2) Ge: [(0,0, 0;;,%, i) + (0, 0,%;0 ,%,i)] (a, a ,c) + (la, ma, nc) (F3) . .111 IL—iii 1_- AS. [(0,0,0,-2373'?) + (xaTs'gsxa49 :98 49 x398 49X 3%)]. (a, a ,C) + (la, ma, nc) (F4) where 2 = l-x and l,m and n are integers. These formulas result in four Cd atoms, four Ge atoms, and eight As atoms per unit cell; in other words, there are four formula units of CdGeAsz per unit cell. I However, an equivalent structure can be generated by considering only 9 two formula units. The corresponding relations are Cd: (0,0,0, 0,2,4)° (a, a, c) + ([2l- nJa/Z, [2m+n]a/2, nc/Z) (F5) Ge: (;%—%.0;%.0,%)°(a,a,c) + ([21-n]a/2.[2m+n1a/2,nC/2) (F6) , 1 1-1 3— 3 1 -3 1 3 , AS. (" "X94, 8; :—‘49L—2’ +X], g” gag—4"? 43 [‘T-XJyg) (aga,C) + ([2l-n]a/2,[2m+n]a/2,nc/2) . (F7) This translation vector between 'unit cells', aflQ+flU+nF§§+§9+%2) (W) corresponds to a twofold screw axis. 134 A Fortran program was written (listed in Figure F1) to generate a list of atom positions, compatible with the list required by the energy band program, by using Equations F1, F5, F6, and F7. As an example, a portion of the list of atoms for CdGeAsZ is given in Table F1 for a = 5.94A and c = 11.22A. From this list, it is easy to verify the interatomic distances of 2.629A for Cd-As and 2.429A for Ge-AS o 135 :UDDzEUBONK:GEOMETRY.DIR:CPGEOM.FR 000000000(‘)00000000000000000000000000000000 1 166 116 Figure F1. Computer listings for chalc0pyrite crystal structure program CPGEOM.FR ROUTINE TO GENEROTE CHOLCOPYRITE OTOM POSITIONSIOTOMIC NUMBERS. ENUMEROTION OF OTOMS IN UNIT CELL: ETC IN DOTO FILE TO BE INPUT FOR EBOND.FR CHOLCOPYRITES HOVE O TETROHEDROL STRUCTURE: O(II)B(IV)C2(V) [LIKE O III-V COMPOUND] OR O(I)B(III)C2(VI) [LIKE O II-VI COMPOUND] USE THE FOLLOUING TO GENEROTE THE CRYSTOL STRUCTURE: CONSIDER 2 FORMULO UNITS (ON ’O’ OTOM OT THE ORIGIN. I COMPLETE TETROHEDROL OBOUT O ’C’ OTOM. OND 3 OTHER ODJOCENT ’C’ OTOMS) OND THE TRONSLOTION VECTOR: ((L*O+N*O/2):(M*O-N*O/2).(N*C/2)) UHERE L.M.N ORE INTEGERS OND O G C ORE THE LOTTICE CONSTONTS NOTE: THE ’C' OTOM IS DISPLOCED FROM CENTER OF TETROHEDROL f BY DISTONCE DELX TOUORD THE SMOLLER TUO ’B’ OTOMS ’fi X=FCN(O»C) FORMULO FROM OBROHOMS G BERNSTEIN;J CHEM PHYS; VOL 55.NO 2,(?/15/71):P?96 INPUT POROMETERS: (FROM CONSOLE) 1) LOTTICE CONSTONTS (O & C) 2) OTOMIC SYMBOL & OTOMIC NUMBER FOR-3 TYPES OF OTOMS IN UNIT CELL 3) NCELL --- DETERMINES NUMBER OF UNIT CELLS TO CONSTRUCT 4) DMOX --- MOXIMUM DISTONCE FOR OTOMS TO BE STORED ON FILE OUTPUT POROMETERS: DOTO FILE --- "?GEOMéGOOS“ --- CONTOINS THE FOLLOUING INFORMOTION ON EOCH RECORD: 1) RUNNING INDEX OF NUMBER OF OTOMS IN FILE 2) OTOMIC SYMBOL 3) OTOMIC NUMBER 4) ENUMEROTION 5) OTOMIC POSITION IN RECTONGULOR COORDINOTES 6) DISTONCE FROM ORIGIN FORMOT IS (IX.I3:3X.O2.2(3X:12):3(3X.EI4.7),6X,EI4.7) EXTERNOL REFERENCES: SUBROUTINE CPUC FUNCTION DOTPD SUBROUTINE CPSORT COMMON/OTMPOS/NENUM(866).RR(4.866):INDEX DIMENSION IOTOM(3).JOTOM(3) POROMETER IN=II.IOUT=16:NPOGE=2?*466K+12 URITE (IOUT) NPOGE URITE-(IOUTyl) FORMOT (25x." CHOLCOPYRITE GEOMETRY") URITE-(IOUT,IIO) FORMOT (/," ENTER LOTTICE CONSTONTS: ",/.5x,"g = «,2; 136 :UDD:EUBONK:GEOMETRY.DIR:CPGEOM.FR 136 I46 145 C. 15 166 186 185 18? I96 195 176 16 15 26 25 000M 000 . . . Q . 0 Figure F1. (Continued) REOD FREE (IN:ERR=166) O IF (O.LT.4.6 .OR. O.GT.?.6) GOTO 166 URITE-(IOUTpI36) FORMOT (5X."C = ".Z) REOD FREE (IN.ERR=1@6) C IF (C.LT.9.6 .OR. C.GT.13.§) GOTO 166 URITE (IOUT.145) FORMOT (5X.“X (ENTER 6 IF UNKNOUN) = ".Z) REOD FREE (IN.ERR=146) X IF (X.NE.6.6 .OND. (X.LT.6.15 .OR. X.GT.6.35)) GOTO 146 IF (X.EO.6.6) X=.5-SORT(C*C/O/O/32-1./16) DELX=(X-.25)*O URITE.(ICUT.166> FORMOT mOmZmeuxzcm3munmsn :«o mucmym_o .AN.x.xv mmgw:EULoou .cowmemsscm .LmnE:z oweou< .Fonszm ow50p< .xmccH ”OsmEL op pymfi SOLE vmpmw4_ .Nm Yevz -)-—)--> ->-->-> p 2* 2 372 d(d°€) = Ad(d'€) (11) 2(V + V ) 2 3 where E is the electric field vector and d is the vector distance # from cation to anion. By evaluating the geometric factor when summed over all bonds in the unit cell for one specific direction of elec- tric field (choose E'= £2 3) which gives a principal susceptability, 2+ ai WW n a 2 (dxdzfi + ddeQ + dzzQ) (12) 2 bonds Now, when considering tetrahedrally coordinated crystals, for each term of the form (+ldxldz Q) or (+Idyldz 9), there will be corres- ponding terms of (-|dxldz 9) and (-|dyldz 9) which cancel the former terms. However, similar pairs of the 2 component terms do not cancel; that is, 2 A 2 A - 2 A (+Idzl) z + (-ldzl) z - 2dZ 2. (I3) Therefore, the summation in Equation 12 can neglect components perpendicular to the electric field or 148 149 z IRE-E) = 52 z dzzfz‘ (I4) bonds bonds which can be written as ram-E) = z dza c0526 (15) bonds bonds where 6 is the angle between E and d. It is this form that appears in the literature,16’18 but the result is used before the summation over bonds is performed. Further simplification is realized by considering a "regular" /\ i tetrahedral structure. Then d = (i 9 i 9 t 2)/‘J§ and, again letting ‘ 3 = :2 Q, the geometric factor from Equation 14 reduces to z 8(8-2) = z (a 1 >22. (16) bonds bonds :7?- Consequently, there is a geometric factor of one—third associated with each bond in a "regular" tetrahedron. LIST OF REFERENCES 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. LIST OF REFERENCES D.J. Chadi and M.L. Cohen, Phys. Stat. Sol. 68, 405 (1975). J.C. Slater and G.F. Koster, Phys. Rev._94, 1498 (1954). A. Nussbaum, in Solid State Physics, edited by F. Seitz and D. Turnbull (Academic Press, New York, 1965), Vol. 18, p. 165. A. Zunger, Phys. Rev. B 11, 626 (1978). G. Arfken, Mathematical Methods for Physicists, (Academic Press, New York, 1970), p. 159. w.A. Harrison, Festkorperprobleme XVII, 135 (1977). N.A. Harrison, private communication. F. Herman and S. Skillman, Atomic Structure Calculations, (Prentice-Hall, Englewood Cliffs, N.J., 1963). E.A. Kraut, private communication. N.A. Harrison, Solid State Theory (McGraw-Hill, New York, 1970), pp. 25-40. C. Kittel, Introduction to Solid State Physics (John Wiley and Sons, Inc., New York, 1971). p. 211. G. Gilat and L.J. Raubenheimer, Phys. Rev._144, 390 (1966). J. Rath and A.J. Freeman, Phys. Rev. B 113 2109 (1975). An-Ban Chen, Phys. Rev. B_1§, 3291 (1977). M. Buchkeit and P.D. Loly, Amer. J. Phys._49, 289 (1972). w.A. Harrison, Phys. Rev. 8.8, 4487 (1973). N.A. Harrison and S. Ciraci, Phys. Rev. B_19, 1516 (1974). M.E. Lines and J.V. Naszczak, J. Appl. Phys._4§, 1395 (1977). G. Arfken, Op. Cit. p. 262. 150 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 151 G.R. Fowles, Introduction to Modern Optics (Holt, Rinehart, and Winston, New Ybrk, 1968), Pp. 168-190. G.D. Boyd, E. Buehler and F.G. Storz, Appl. Phys. Lett. 18, 301 (1971). T‘— G.D. Boyd, H. Kasper and J.H. McFee, IEEE J. Quantum Electron. QE-7, 563 (1971). G.D. Boyd, E. Buehler, F.G. Storz and J.H. Nernick, IEEE J. Quantum Electron. QE-8, 419 (1972). G.D. Boyd, H. KaSper, J.H. McFee and F.G. Storz, IEEE J. Quantum Electron. QE-8, 900 (1972). S.H. Nemple, J.D. Gabbe and G.D. Boyd, J. Appl. Phys. 46, 3597 (1975). T’— J.F. Nye, Physical Pr0perties of Crystals (Oxford, London, 1975), pp. 235-259. D.S. Saxon, Elementary Quantum Mechanics (Holden-Day, San Francisco, 1968), pp. 268—275. G. Arfken, op. cit., pp. 554-558. S.C. Abrahams and J.L. Bernstein, J. Chem. Phys. 55, 796 (1971). S.C. Abrahams and J.L. Bernstein, J. Chem. Phys._§23 5607 (1970). S.C. Abrahams and J.L. Bernstein, J. Chem. Phys. 59, 5415 (1973). International Tables for X-Ray Crystallography, Vol. 1, edited by N. Henry and K. Loundale (Kynoch Press, Birmingham, England, 1965), p. 212. H. Y-P Hong, J.C. Mikkelsen and G.N. Roland, Mat. Res. Bull. 9, 365 (1974). International Tables for X-Ray Crystallography, op. cit., p. 266. IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII