.Iv. . if. v .6 , t, ;. . ? 9 {Ii-:3...— -& :t . V. v — .'. thus.) . t I....P\ Xh- . a; .7 ; .Il: .Vipv‘v: . . rlv'.’ t v. I In}... II...- 7p- . .2. 5" .str.!# . X Dull .l!(.rv’. . (-1", ¢ 2...? .201 .1: s... {.f .0... LAY? bk! I {lit 1.: . ,. IquVotlrd tut .57. V . I :67! 71. Iv..‘0nll. (1.3!...Ir'...l E11: \l 213.67... .. l1. .tz...rv)\~.“ v i). LEI-.35: . Illlgt.alul‘v'bsfi. l'z‘vltll I‘IDVQ‘E v'ftlfl.l~'t'|f v. u In}. cl I‘sl .‘tfo‘v‘liv’lll‘o'n'... ’01:.btlcvlvo . Ll‘ .1... A4A‘I’I.‘vv$0l.i., - . . .r I} I 4 . ..Ir. . .- .tliovl. y‘.... I {n.llu‘r-os, . . . : .. L»... .. .lcflvv.» Iln.vv.. not? 3 ll' litivl '05:... l0!!!r-Ivyfirri...‘. .5 ~ 7: . til! . . 1“... :19: . 4 :12.) ..v.‘n. :‘I .0)..nu1?3ttt.u‘o III If): 5.. Ivy Iv. I I... I airr‘! .. 18!!!!tvyltn RI. 1.. it; 1!... I I. 1’. l' AV Uch. 1.61,...5 Lr fl .. 1r .3. 2:. .hl.‘ 91;?! A I I‘ll-IV J.!Ovbl nu‘l‘rl uAitn 0.! v\4\l‘ll".!n\b It. .v»v.AtY. .1: on . ~ .1. I .-. Mr... Vllvo01!.|9~.ixt¢.. .II' I‘M-IN”; .lllo.». ‘ ‘ v1.3 ' “F3“; Ml HIGAN STAT NIVERS IIIllllIHlll ll ll lulwlliilllll‘ll 3 1293 00891 9361 This is to certify that the dissertation entitled Part I: The Simulation of Effective Transport Coefficients in Composite Materials II: Electron Localization: Quantum Molecular Dynamics presented by Sheh-Yi Sheu has been accepted towards fulfillment of the requirements for Ph.D. Chemistry degree in 22431514.; Robert I. Cukier Major professor Date 7-19-90 MSU is an Affirmative Action/Equal Opportunity Institution 0-12771 LIBRARY 7 M'Chlgan State UnlversIty ——_. I I L; 1 PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MSU Ie An Affirmative ActlorVEquel'Opponunity Institution cWMS-nt PART I THE SIMULATION OF EFFECTIVE TRANSPORT COEFFICIENTS IN COMPOSITE MATERIALS PART II ELECTRON LOCALIZATION : QUANTUM MOLECULAR DYNAMICS BY SHEH-YI SHEU A DISSERTATION Submitted to Michigan State Univerlity in partial fulfillment of the requiremente for the degree of DOCTOR OF PHILOSOPHY Department of Chemistry 1990 ABSTRACT PARTI THE SIMULATION OF EFFECTIVE TRANSPORT COEFFICIENTS IN COMPOSITE MATERIALS PART II ELECTRON LOCALIZATION 3 QUANTUM MOLECULAR DYNAMICS BY SHEH-YI SHEU PART I Simulation methods for the investigation of effective transport coefficients in composite materials are an extremely important technique in science. We study composite materials consisting of spherical impenetrable inclusions embedded in a homogeneous matrix. The method is referred 'to as an analytic-simulation method. The effective transport properties can be: diffusion, dielectric behavior, elastic and viscous constants, electric and thermal conductivity, and magnetic susceptibility. For a given configuration of inclusions, a set of coupled algebraic linear equations written in terms of the t-operators and the propagators between inclusions'is solved to obtain the multipole moments of the inclusions. The minimum image convention, a spherical interaction cutoff, is used to evaluate the multipole moments for each inclusion. The resulting total sample polarization is related by macroscopic electrostatics to the sample’s effective dielectric constant. We simulate two classes of properties of a composite material: the static and frequency dependent dielectric constant. For the static dielectric constant, we investigate how the effective dielectric constant, ee’ is controlled by multipolar effects. The investigated systems are conducting inclusions in an insulating matrix, the inverse case, and coated inclusions (composite-composites). Due to the difficulty of obtaining converged results for £6, for conducting inclusions in an insulating matrix, at volume fractions above 0.5, we evaluate Ce by a random walk method. The random walk method permits an accurate evalution of 8e up to volume fractions corresponding to near to close packing of the inclusions. We also consider the frequency dependent effective dielectric constant, sew), of composites with metallic inclusions, modelled as Drude oscillators, in an insulating matrix such as a glass. The optical properties such as line-broadening and line—shift of 'the dielectric constant lineshape of these composites have been studied. We find that the lineshape of saw) is greatly broadened by the allowance for the electrostatic interactions among the inclusions in comparison with the Maxwell-Garnett results. We consider different types of disordered configurations for this problem. Comparison is made between results based on the minimum image convention and the lattice-sum approach. The former method is more efficient than the customary lattice sum approach, which employs Ewald sums, and yields results in good agreement with the latter method. PART II Quantum molecular dynamics (an Adiabatic Simulation Method) has been used to discuss an excess quantum electron which interacts through .pseudopotentials with a fluid of classical molecules. A detailed algorithm for the investigation of the equilibrium and dynamical properties of this coupled quantum-classical system is described. This study focuses on the localization, dynamics, and mode of transport of an excess electron in condensed helium. Properties investigated include the correlation functions and electronic energy of the ground and lowest excited states, and the diffusion coefficient of the ground state electron. Copyright by SHEH-YI SHEU 1990 TO MY FAMILY iii ACKNOWLEDGMENTS I want to thank Professor Robert I. Cukier for his guidance, patience, and support in conducting this research. I would like to thank Professor K. L. C. Hunt for her reading this manuscript in great detail and pointing out the errors. I am grateful to Professors J. Dye, D. Nocera, and T. liaplan to be my committee. Finally I want to thank Center for Fundamental Materials Research in Michigan State University for grant support. iv TABLE OF CONTENTS LIST OF TABLES LIST OF FIGURES PART I CHAPTER 1 INTRODUCTION CHAPTER 2 HISTORY 2.1 DEVELOPMENT OF THE STATIC as 2.2 DEVELOPMENT OP THE FREQUENCY DEPENDENT gem) 2.2.1 DRUDE FUNCTION FOR METALLIC INCLUSION 2.2.2 THE QUANTUM SIZE EFFECT OF gem) CHAPTER 3 ANALYTIC DEVELOPMENT 3.1 THE PRINCIPAL EQUATION 3.2 DIPOLE TENSOR 31m“) CHAPTER 4 NUMERICAL METHOD 4.1 METHODS OP GENERATING CONFIGURATIONS 4.1.1 IMMEDIATE METHOD 4.1.2 METROPOLIS METHOD 4.1.3 THE DISORDER OP CONFIGURATIONS a) RANDOMIZED LATTICE b) VACANCY LATTICE 4.1.4 ' MINIMUM IMAGE METHOD 4.2 COMPUTATIONAL ALGORITHM V 13 13 17 20 21 29 38 38 38 4O 42 42 44 45 45 CHAPTER 5 5.1 THE 5.1.4 5.1.5 5.2 THE 5.2.1 5.2.2 5.2.3 5.2.4 5.2.5 CHAPTER 6 APPENDIX A APPENDIX B APPENDIX C REFERENCES PART II CHAPTER 1 CHAPTER 2 RESULTS AND DISCUSSION STATIC DIELECTRIC CONSTANT CONDUCTING INCLUSIONS AN INVERSE CASE A SPECIAL CASE COATED INCLUSIONS CONCLUSION FREQUENCY DEPENDENT sew) RANDOMIZED LATTICE VACANCY LATTICE RANDOM SIZE-DISTRIBUTED EFFECT SILVER IN GLASS COMPARISON OF MINIMUM IMAGE AND LATTICE SUM METHODS CONCLUSION DERIVATION OF THE LATTICE SUM FOR INDUCED DIPOLAR SYSTEM IRREDUCIBLE CARTESIAN TENSORS LIST OF COMPUTER PROGRAMS INTRODUCTION COMPUTATIONAL METHODS 2. 1 MOLECULAR DYNAMICS 2.2 SCHRODINGER EQUATION 2.3 THE ADIABATIC SIMULATION METHOD vi 80 84 88 97 101 133 2.4 THE MOVING GRID TECHNIQUE CHAPTER 3 RESULTS AND DISCUSSIONS CHAPTER 4 CONCLUSIONS APPENDIX A LIST OF REFERENCE PARAMETERS APPENDIX B LIST OF COMPUTER PROGRAMS REFERENCES vii 167 PART I Table 1 Table 2 PART II Table 2. 1 Table A.1 LIST OF TABLES Comparison of tie/£2 for a special case. Comparison of tie/Em for composite-composites case. List of the diffusion coefficient of helium. List of the reference parameters. \riii 61 154 184 PART I Figure 3.1 Figure 4.1 Figure 5.1 Figure 5.2 Figure 5.3 Figure 5.4 Figure 5.5 Figure 5.6 Figure 5.7a Figure 5.7b Figure 5.8 Figure 5.9 Figure 5.10 Figure 5.11 Figure 5.12 Figure 5.13 LIST OF FIGURES Supercells. Comparison of the radial distribution functions. Plots of Ite/e:2 for conductors-in-insulating matrix versus ¢. Plots of 58/82 for insulators in conducting matrix versus (D. Plots of tie/£2 for composite-composites versus 4). Plots of CPU time versus the dimension of matrix. The 58/15 versus ¢ for the different types of 2 configurations. 1m 68(0)) for fluid and randomized lattice. Im sea») for vacancy lattice. Im 58(4)) for vacancy lattice. The log-normal distribution. Comparison of Im £e((d) for equal-sized and size-distributed fluid. Im 88(0) for size-distributed fluid. Im 58(0) for Ag in a glass. Cole-Cole plot of silver spheres in glass. Im sew) for silver spheres in glass. ix 31 43 82 83 PART II Figure 2.1 Figure 2.2 Figure 3.1 Figure 3.2 Figure 3.3 Figure 3.4 Figure 3.5 Figure 3.6 Figure 3.7 The helium-helium radial distribution function. Plot of (R2) versus time. Plot of the electron radial distribution function. Contours of the excess electron density. The radial distribution function of helium measured from the center of mass of the electron. Plots of the electronic energy versus the number of quench steps. Contours of the electronic configurations for ground and lowest excited states. Plots of the fluctuation of the energy in the molecular dynamics adiabatic simulation. The electron mean-square displacement versus time in helium. 173 174 176 181 PART I THE SIMULATION OF EFFECTIVE TRANSPORT COEFFICIENTS IN COMPOSITE MATERIALS 1. INTRODUCTION Many composite materials are of great scientific and technological interest. The concept of the dielectric constant is of considerable importance in the description of the macroscopic electric behavior of composite materials. Based on the electrostatic point of view, when an electric field is imposed on a perfect conductor, there is no difference between microscopic and macroscopic fields. In contrast, when we apply the same field on a dielectric system, the dielectric will respond to the field and a polarization will develop at a molecular level, and dipoles will reorient in response to the field. If the constituents of the composite are sufficiently large that they can be assigned a spatially-dependent dielectric constant €(r)1, then the electric displacement D(r) in the system is related to the electric ~~ field E(r) as D(r) = 6(5) 51(5). (1.1) The electric field E(r) is given in terms of the potential field as E(r) = - Y Mr). (1.2) ~~ Let us combine Eqs. (1.1) and (1.2) to obtain a first order 2 differential equation, which can be applied equally well to similar phenomenologies such as diffusion, electric and thermal conduction, elasticity, viscosity and magnetizationz. In order to calculate the dielectric constant of a composite material, we must solve Maxwell’s equation V - D(r) = Slr). (1.3) ~ ~~ where S(r) is the source term. Many composite materials are heterogeneous on a scale such that at each space point the material properties obey the macroscopic constitutive equations. Thus, for a macroscopic sample of the material the ensemble average obeys (r) = ee(r) (1.4) which gives a definition for the effective dielectric constant 6e of the composite. The concept for this description of phenomena must be on a length scale that is large compared to the typical scale of the inhomogeneities, or the microscopic correlation lengths. The average means a configuration average over many realizations of the material distribution functions, a so-called statistically ensemble average. In order to evaluate the effective dielectric constant so, the composite material’s structure must be prescribed; this is done 3 statistically. Then an averaging procedure must be carried out over the prescribed distribution of configurations. We study materials consisting of spherical non-overlapping inclusions, which are assumed to be statistically distributed with a probability distribution independent of the position of the external applied field, embedded in a homogeneous matrix. Simulations for the investigation of effective transport coefficients in composite materials have become an extremely powerful technique in science. The properties of diffusion, dielectric, elastic, and viscous constants, electric and thermal conductivity, and magnetic susceptibility can be studied by analogous methods. There are two steps in the study of effective transport properties by this technique: First, an expression for the transport property must be given in terms of some microscopic quantities, which can be obtained from the simulation method. For example, the effective dielectric constant 6e of a composite material is obtained via the polarization of sample and the electric field. The polarization is defined as the total dipole moment per unit volume. Here the microscopic quantity is the dipole moment on each inclusion. Second, a simulation method must be developed to evaluate these microscopic quantities. The method we present is a combination of an analytic calculation and a numerical simulation, which we will refer to as an analytic-simulation method. It is exact in principle. For the given dielectric constants of the inclusion and matrix, the sizes, shapes, 4 volume fraction ¢ and distribution in space of the inclusions it will yield as to the accuracy implied by a finite size simulation of an infinite material. In order to evaluate an accurate Ce, we employ a technique similar to the method used by Lebenhaft and Kapral in their study of diffusion—controlled reactionsa. For a given configuration of inclusions, a set of coupled linear equations for the multipole moments of the polarization of the inclusions can be written in terms of the t-operators of each inclusion and the propagator between pairs of the inclusions. The total polarization of the sample is obtained by matrix algebra. The multipole expansion for the polarization is required to convert the coupled integral equations for the inclusion polarizations to linear algebraic equations. The t-operator gives the response of the ith inclusion to an "arbitrary field", which arises from the external field and all the inclusions excluding the ith inclusion. Correlations among the inclusions are exactly accounted for with this method. For high volume fractions of the inclusions it is necessary to be concerned with the higher multipole contributions for the calculation of an accurate effective dielectric constant. A systematic and direct method is used to solve these linear coupled equations, Eq. (3.19). In the following, numerical results for some special models will illustrate the importance of multipole contributions to the effective dielectric constant. The simulations will be carried out by generating the inclusions in a primary cell (typically using a Monte Carlo method) and 5 periodically imaging the cell. We use a minimum image method to evaluate the dipole moment and higher order multipoles of a given inclusion inside the primary cell. As the multipole moments are defined in terms of the propagator 1 tensors, the dipole field involves the dipole tensor W(—r—). This is a long-ranged interaction and leads to conditional convergence of dipole sums. Due to the long-ranged behavior of the dipole tensor, the dielectric constant is not well-defined i.e. it is sample shape-dependent in the thermodynamic limit. The lattice sum approach, which employs a large number of periodic images of the primary cell and involves an Ewald-summation method‘, has been used to represent the dipole propagator in an efficient and rapidly convergent form. Comparison is made between results based on the minimum image method and the lattice-sum approach. In the investigation of the optical properties of composite materials, it is the frequency dependent effective dielectric constant sew) that is experimentally probed. The absorption lineshape-broadening and line-shift of the frequency dependent dielectric constant for composite 'materials can also be studied by the analytic-simulation method. The types of the configurations we investigate include: randomized lattice, vacancy lattice, and mono- and poly-disperse liquid-like structure. The frequency dependent dielectric constant of the inclusions is given by a Drude model, as appropriate to metallic inclusions. We find 6 that the lineshape of the imaginary part of the effective dielectric constant is greatly broadened by the allowance for the electrostatic interactions among the inclusions in comparison with the Maxwell-Garnett results. The rest of our work is outlined as follows. In Chapter 2 we review some methods for calculating the effective dielectric constant of composite materials. The first section discusses the historical development of the methods of calculation of the. static dielectric constant in random media. The second section contains an introduction of the concept of the Drude model for a metallic inclusion and the development of the frequency dependent dielectric constant in composite materials. We also discuss the Quantum size effect that becomes important for small particles. In Chapter 3 the analytic development of analytic-simulation method is presented in terms of the t-operator or the polarizability of the inclusion and the propagator tensor gal. The relation between the microscopic dipole moments and the macroscopic averaged polarization is discussed. The other part in this Chapter discusses the properties of the dipolar tensors including the basic idea of the lattice 'sum approach. In Chapter 4 we describe the numerical simulation methods and introduce the different methods to generate the different types of configurations. We extend the minimum image method to calculate the dipole moments and higher order multipoles for each inclusion. The rest of this Chapter demonstrates the intent of the computer programs and lists the algorithm for computations. In Chapter 5 we present the results for the different physical models. First, in the static dielectric problem we use the 7 analytic-simulation method to discuss the conducting case, the inverse case and a composite-composite case. Then, the results from the random walk method are compared with those of the analytic-simulation method for the conducting case. The second section contains the calculation of the frequency dependent dielectric constant for the different types of configurations for a hypothetical and more realistic model of silver in glass. Finally, we compare the results of the analytic-simulation method using the minimum image method with the results of the lattice-sum approach. Chapter 6 contains the conclusions and outlines future applications of the method to composites with non-spherical inclusions, electromagnetic , versus electrostatic effects and wavevector-dependent effective transport coefficients. 2. HISTORY 2.1 DEVELOPMENT OF THE STATIC e:e Theories about the static dielectric constants of composite materials have a long history. For this field there is a great deal of literature reviewed by Landauer‘. In the following we will briefly summarize some approaches to these problems. A general expression for the. dielectric virial coefficient for the one inclusion problem was derived by Maxwell5 in 1873. He considered a spherical inclusion of dielectric constant t:1 embedded in a matrix of dielectric constant 62. The dielectric constant Ce of a macroscopic sample of the material depends on the ratio 61/62 and the volume fraction of inclusions, ¢. His formula takes into account the induced dipole moments of the inclusions, and the Ce is calculated by taking the result from one inclusion and multiplying by the number density of inclusions. This can be applied to obtain an accurate 8e for small volume fractions. For the dilute inclusion problem 6e is given by 86/6 = 1 + 3 a 4’ (2.1) £1 - :2 4n 3 where a = and 1» = -—- a N. V is the volume of the system, :1 + 252 3V 3 is the radius of inclusion, and N is the number of inclusions. Maxwell-Garnetts presented an effective medium theory for the calculation of £6. The 8e evaluated by this theory is identical to a result from the Clausius-Mossotti formula7’a. In effective medium theory the effective field applied on a given inclusion is defined as the sum of an applied field and a field due to the polarizations induced on the other inclusions. The matrix is thought of as a homogeneous polarized medium surrounding the inclusions. This leads to an expression ( 5 - E ) ( E - 5 ) ° 2 = 1 2 4. (2.2) ( c + 252 ) ( £1 + 252 ) called the Maxwell-Garnett or Clausius-Mossotti formula. A similar expression was also obtained by Bruggeman 9. He iterated the polarization of the inclusions and the matrix until the average net polarization vanishes. It is a self-consistent field calculation. He employed this method to improve the Clausius-Mossotti result. The formal expression is (8-6) (8-6) [ 1 ° ]¢+[ 2 ° ](1-¢)=o (2.3) ( £1 + 26° ) ( £2 + 256 ) 10 This treats the two phases in a symmetric fashion. These mean field theories neglect the correlations among the inclusions, which become important at high volume fractions of inclusions. Certainly, these theories have the virture of simplicity and give physically meaningful results over the entire range of volume fractions. However, it is necessary to ascertain whether 66 can be calculated accurately by these theories. Ce can be expressed as an expansion in the volume fractions of inclusions. It is a virial series : _ 2 56/5 - 1 + A1¢ + A24 4» (2.4) 2 Batchelorm discussed the calculation of the second virial coefficient. The second virial coefficient involves the two-body distribution function and reflects the local structure of the system. In 1968, Levine and McQuarrien obtained the polarizability of metallic spheres embedded in insulating medium by solving Laplace’s equation in bispherical coordinates to get an exact second dielectric virial coefficient. In 1973 Jeffrey” expanded the work of Maxwell to second order by using Batchelor’sm method, which reduced the problem to a consideration of interactions between pairs of spheres, in order to avoid the cOnvergence difficulties. He used a twin spherical harmonics expansion due to Ross” to solve the two-spheres Laplace problem. In Jeffrey’s work the second dielectric virial coefficient AZ, corresponding to thermal conduction, is evaluated explicitly for all 11 values of the ratio of conductivities of the two phases. Subsequently, Felderhof, Ford and Cohen“ developed a multiple scattering expansion in order to obtain a formal expression for the virial coefficients. They formulated an infinite system as an expansion in cluster terms, each of which is absolutely convergent. By rearrangement of the terms, corresponding to removing the depolarization effects, problems of conditional convergence never arise in their problem. Sridharan and Cukieris developed a t-operator multiple scattering (e -s theory for the Clausius-Mossotti function as a virial expansion in terms of a reference medium dielectric constant am. There are several results obtained from their expressions. When they neglect the correlations among the inclusions and cm is set to be equal to the medium dielectric constant 62, they get the Clausius-Mossotti formula. When cm is set to the effective dielectric constant 66, the Bruggeman result is established. As the correlations among inclusions are considered, the second virial coefficient is able to be calculated exactly when we set emzez. Another result also can be obtained with emu‘3 called pair order effective medium theory. It is not possible to get results without an approximate evaluation of the formal expression. This theory generalizes .Bruggeman’s self-consistent lowest order theory. However, the shortcoming of this theory is the difficulty of evaluating the results even at pair order. 12 Because the virial diagrammatic expansion has not been carried out beyond the second virial coefficient, it is restricted to low volume fractions. Meanwhile, the self-consistent expansion is not well controlled. This theory is limited to spherical inclusions and to static material properties. The evaluation of this theory for arbitrarily shaped inclusions and frequency dependent material properties has not been carried out. A variational approach leading to upper and lower bounds for the effective material properties was developed by Hashin and Shtrikmanm. For a two-phase problem the bounds for 8e are in terms of the ratio 511/62 and 45. When the ratio, 81/82, is extremely large or small, the bounds are far apart for moderate volume fractions. The variational approach is summarized in the literature”. In 1892 Rayleigh18 solved the partial differential equation to second order for periodically arranged conductors and provided results for a wide range of conductivity and the entire range of volume fraction. McPhedran and McKenzie19 extended his method to calculate the conductivity of periodic lattices consisting of conducting spheres in an insulating matrix. They solved the problem by a potential and . field expansion On each sphere and obtained the expansion coefficients of the fields. They were able to include multipolar effects to very high order. Their results for the conductivity of a regular array of spheres yield good agreement between theory and experimental measurements”, even when the volume fraction approaches a lattice close packed value. 13 2.2 DEVELOPMENT OF THE FREQUENCY DEPENDENT 88(0)) 2.2.1 DRUDE FUNCTION FOR METALLIC INCLUSION The development of frequency dependent dielectric properties starts from the assumption of the Drude21 model for the dielectric constant of metallic inclusions. If the inclusions are so small that the quantum size effect comes out, then Drude model can no longer be used. Let us consider a system consisting of spherical inclusions of H N radius a embedded in a homogeneous medium with dielectric constant 61(0)). If the absorption of the medium can be neglected 62 is frequency independent and therefore real. For example, the medium 22,25 23,24 could be a. ceramic , glass , K0125, gelatin25 and so on. The definition of the Drude model for the dielectric constant of metallic inclusions is up? 5 (w) = e -' (2.5) 1 1‘” u ( u + 1 r ) where 610 is the (complex) high-frequency dielectric constant, (61" = I 810 :1: i 81;), up is the plasma frequency, which is the natural frequency of the density fluctuation for the free electron in bulk metal, and 1‘ is the damping constant or width, defined as 14 u l 100 2 I‘:7+I‘.b:ya+e . (2.6) a 1 where (a. is the resonant frequency and 7a is the damping constant for a ' 00 ' I e s 26 radius a' of Inclusion given by 7a 3 10° + VF/a. (2.7) where 10 is the bulk damping constant, VF is the Fermi velocity, and 1‘ib is the interband contribution. The Drude model provides that the dielectric constant of metal inclusion is a function of the frequency and the damping constant 1". The dielectric constant can also be written in terms of the frequency dependent conductivity 0(0) i 0(4)) 5 (u) : s - ————- (2.8) 1 100 m For one inclusion effect, the frequency independent 8e is replaced by 56(0)) in Eq. (2.1). By combining Eqs. (2.1) and (2.5), the effective dielectric constant 66(0) has an absorption peak of width 1' centered at the approximate resonant frequency w u = p (2.9) ( c ’ + 2 £2 )1/2 1m When the Clausius-Mossotti (CM) or Maxwell-Garnett (MG) expression is applied for 86(4)), we can write the CM formula as ( E (0) - E ) e 2 = 4 (2.10) ( sew) + 262 ) ( 61(0) 4» 262 ) Eqs. (2.5) and (2.10) show that the lineshape is similar to a Lorentzian distribution with its absorption peak centered around r 1 1/2 ( 1-¢ ) (2.11) 8| 11 E: P 515(1-¢)+(2+¢)ezl The width is independent of the volume fraction of inclusions. These approximate expressions result from the conditions 1" << up and (08 >> "1'2; The CM results predict that there is a shift of the line to lower frequency due to increasing CD and no line broadening since the width does not change. The description of electromagnetic propagation and the investigation of the optical properties of composite materials in terms of the properties of the dielectric constants of the components and the volume fraction of the inclusions have been done both experimentally and theoretically1'27. One of the key issues involves the broadening of the spectrum due to the effect of interactions among the inclusions. Kreibigza calculated 59(0) for silver in gelatin by using the average-t-matrix approximation (ATA)6. The ATA is also a Maxwell-Garnett approximation, corresponding to obtaining an accurate result at low 4!. Another approach is the Coherent Potential 16 Approximation (CPA). The CPA corresponds to Bruggeman’s effective medium theorys. These approximations have the common properties : 1) They are long-wavelength approximations. 2) The scattering process is only dipole, because this is one sphere problem and the response is to a constant applied field. 3) The structural information is only the volume fraction of inclusions. The interesting problems are related to the strong scattering properties and the higher volume fractions of inclusions. Davis and Schwartz”. employed multiple scattering theory to evaluate 88(6)) for a disordered system with spherical inclusions at quite high ct. They used Lax’s quasicrystalline approximation (QCA)29, which is known as the Maxwell-Garnett approximation, and compared with the application of Roth’s effective medium approximation (EMA)30'_31 . The-EMA provides a complete description of the spatial correlations of the inclusions but not the electromagnetic correlations, and an accurate description of the plasma resonance in composite problems. They found that the EMA is better than the QCA. Felderhof and Joneszr3 calculated send) for silver in glass or water by a cluster expansion technique, which is combined with a spectral representation analysis based on the work by Bergman”. The 58(0) in Bergman representation is an integral of Hilbert type with a spectral density determined by the statistical geometry of the configurations, and is accurate to second order in ¢. The results obtained by them assume dipolar .intersphere interactions, and neglect higher multipolar interactions. From section 2.1, note that it is 17 necessary to include the higher order multipolar effects during the calculation of the static dielectric quantities at high 4’. For example, the contributions from all multipoles to the second dielectric virial coefficient are roughly double that of the dipole contributions. Meanwhile, they” found that around the resonance peak the dipolar interaction term dominates the multipolar contributions. They predict that when the distance between the pairs of inclusions is roughly equal to the mean pair separation, d : (1.49)“3 a {U3 , (a is the radius of inclusions) only pairs of inclusions dominate the line-broadening and line-shift of resonance. The dipolar approximation is quite accurate for 66(0)) at large pair separation or at low 4). At high «t due to the contributions of the electrostatic interactions between inclusions the absorption peak differs from that at low 4’. Kantor and Bergman2 also used Bergman’s spectral representation and an expansion in spherical harmonics in order to describe the interactions among the inclusions. They derived an exact theory for calculating 86(0) for composites. 1 Their theory is not restricted to lower order interactions between inclusions. 2.2.2 THE QUANTUM SIZE EFFECT OF 86(0)) 33, 34 (QSE) What is the quantum size effect? Quantum size effects occur when the inclusion size is so small that the material properties ( e.g. 61(0) ) do depend on size because of the discreteness of the energy levels, yet, it is still useful to describe ,the inclusion in 18 terms of a macroscopic material property. From the expression Eq. (2.6), 7a is larger than 10 due to the V finite size of inclusion or the term ~52. The size of inclusion controls the value of 7a' Consequently, the width 1‘ is dominated by 18. Therefore, the Drude model of 61(0)) is size dependent and 58(4)) is too. The interesting problem is how 68(0) is affected as the inclusion size decreases. "As the particle size continues to decrease to the tens of angstrom range, new phenomena involving quantum physics, electromagnetics, and hydrodynamics are still to be explored."35 For very small inclusions of equal size and shape, the plasma resonance absorption shifts and broadens, and shows fine structure corresponding to transitions between discrete conduction band energy levels. Two main effects on optical properties may occur as a function of size: a shift of the peak and a change in the width. There are two expressions for the dielectric constant: the classical Drude model used with a size limited mean free path and the quantum mechanical theory of Kawabata and Kubo:39 accounting for quantum size effects. If the inclusions are so small that the Drude model can no longer be used, the dielectric constant should take quantum size effects into account. Friihlich37 pointed out in 1937 that the continuous electronic conduction band of a metal should break up into observable discrete states when the dimension of the metal become small enough. Kubo 35 19 formulated this problem quantitatively in 1962 and observed these quantum size effects. Both theoretical and experimental results on these problems have been reviewed by. Kreibig and Genzelaa. There are quantum mechanical and classical models used to explain the changes in the optical properties as the size of particle decreases. In the classical model, the damping constant in the Drude free electron theory, which is the inverse of the collision time for conduction electrons, is increased due to increasing collisions with the boundary of the particles. In quantum mechanical model, Kawabata and Kubo 39 have argued that the optical spectrum should be discrete with the spacing increasing as the particle size decreases. The width of the resonance should be described as due to the plasma mode damping which resultes from the excitation of one-particle modes. The 'peak width predicted by the quantum mechanical model is in better agreement with experiment:33 than that predicted from the size-dependent Drude model, where the size dependence arises from the classical limitation of the electrons’ mean free paths. 3. ANALYTIC DEVELOPMENT We study the problem of the effective transport properties of small material inclusions dispersed in a dielectric host. We use a method similar to that used by Lebenhaft and Kapral in a reaction-diffusion problems. Let us consider a system consisting of N non-overlapping spherical inclusions of dielectric constant 61 and radius "a" embedded in a composite medium of the dielectric constant 52 all within a volume V. If one attempts to impose an external electric field on. such a heterogeneous material, it is important to understand the response of the system and that of the inclusions due to an external field. For instance, in rather dilute media, since the inclusion separation is quite large, there is little difference between the macroscopic electric field of the system and that incident on every inclusion. But in dense media due to the closeness of the inclusions, the polarization of neighboring inclusions induces an internal field Ei at any given inclusion in addition to the average macroscopic electric field of the system. The induced dipole moment is proportional to the electric field acting on the inclusions. The polarizability of an inclusion, or, is defined as the ratio of the average polarization to the total applied field on the inclusions. Thus, the inclusion’s polarizability 20 21 characterizes the capacities of the response to an applied field for inclusions. The approach is: First, the microscopic induced dipole moments of the inclusions must be addressed by an analytic method. Then, the effective dielectric constant has to be obtained by relating the ensemble average of the total polarization and the macroscopic electric field. 3.1 THE PRINCIPAL EQUATION According to the electrostatic problem, the material equation for the electric potential, (Mr), is 2 ° ( 6(5) ENE) ) = 8(5) (3.1) where E(r) = 61 ( £2 ) for 5 inside ( outside ) any of inclusions, and 8(2) is a source term, which is independent of the positions of the inclusions. For convenience, let 65(5) = E(r)-£2. If I; is inside an inclusion, 5dr) equals 61-62; otherwise it is zero. Rewriting Eq. (3.1) by using the definition 65m yields e2 v2 1411;) = - g - 1 55(5) 3 4(5) 1 + 5(5) (32) 22 This alternative permits us to introduce the medium propagator, corresponding to the inverse of the operator EZVZ. Let us take the gradient of Eq. (3.2) and solve it formally to get E(r) = f dr'gdr-r’ I) v v - [6am E(r'n + How (3.3) The electric field in Eq. (3.3) has been defined as E(r) = - V Mr). (3.4) the propagator in the medium is -1 g(r) = (4xezr) , (3.5) and the external source field Eo(r) is Eom = 3] dg'ng-g' I) 815') (3.6) Then, integrating Eq. (3.3) by parts, we have E(r) : I dr’§(|r-r'|) - sew) E(r') + E°(r) (3.7) where 23 1 C(lr|)=VVs(r) =H(r)- 3E 16(r) (3.8) - ., ~~ -~ 2- ~ 1 A. g(r)=-_—'—‘3—'(l-3rr). (3.9) ” 4n€2 r A g(r) is a dipole propagator for r at 0. r is a unit vector of r. ° We only need g(lrl) for r at 0, so that the 5-function in Eq. (3.8) is dropped. Therefore, the electric field is given in terms of dipole propagators and 56(r), the response of the inclusion due to Eo(r). In order to get an expression for the effective dielectric constant, 6e or 88(0), we define a polarization related to the electric field by 66(r) N P : 5cm E(r) = 2 “1"" (3.10) .. .. .. .. i=1 .. .. where we anticipate that in isotropic media the polarization and field are parallel. The last equality is the sum of the polarization of the inclusions over the N inclusions. Because 6£(r) is zero outside the inclusion, the ui(r) is non-zero only inside the ith inclusion. We review the theory developed by Sridharan and Cukier 15. It is an application of the t-operator method expressed as a multipole expansion. First, we express the { ui(r) ) in terms of the ~~ 24 one-inclusion scattering operator ti(r,r’), then we write ui(r) to be ui(r) = I 45' li(r.r') - Ei,eff(5') (3.11) ~ ~ - ~~ Here E. eff(r) is the effective 1 field at r that arises from the ~ ’ ~ presence of all the inclusions excluding the ith. From the perspeCtive of the ith inclusion, all the other inclusions serve to produce a field Einff‘f” which acts on the ith inclusion. Hence, the field Ei,eff(£) is given by Elififflr) = E 0(1') + j£° Idr' g(Ir-r’ I) “ff ) (3.12) We use the solution of the electrostatic problem of one inclusion's response to an external field E°(r) to obtain an expression for Eilrm’) as a multipole expansion and we hav915,* °° c ' t .(r,r’) = 4ue2 1: 2 t ( v ) 6(r-Ri) O ( v ) 6(r’-Ri) (3.13) l * The multipole expansion for the t-operator is defined as: 0 I CI 155. 5') =42 (-1)M (4202)" 1&0 0w ‘Mr- R. ”(W ) «Hy-gin. I where L“ is a tensor of rank ((+0 +2). As E“ is constructed as: 1"" a 5 = 5 . ....5 5 “B’HIHZ. ° wt”?! ' “7/ t a M1”1 "Cut 015 this yields Eq. (3.13). 25 where Ri is the center of the ith inclusion relative to an arbitrary . . C ’ C . . or1g1n. The term ( V ) 5(r-Ri)@( V) 5(r -Ri) 1nvolves the {th order gradient with respect to r of a 5-function centered at Ri fully contracted with the corresponding Cth order gradient with respect to r' to give a scalar. The operator 0 means the contraction of two tensors. For a spherical inclusion, t? is given in terms of the polarizability of inclusion, are, to be (X 2+1 tf’c : (3.14) [ ((+1)! (2C+1)!! ] and f. ( s - s ) 2 1 ac = am”) (3.15) el£+ez(£+1) £i(r,r’) is proportional to the unit tensor ,_1_ and only the terms tf’c exist in Eq. (3.13), as t? = 0 for C at C'. This simplification is only obtained for the spherical inclusion problem and since the Eq. (3.1) is a scalar Maxwell equation. By combining Eqs. (3.11) and (3.12), an important equation can be produced as 26 N u.(r) : I dr' t.(r,r’) ' E (r’) + E I dr’ J‘ dr" t.(r,r’) ~14» 4» -l~~ ~O~ j¢i ~ ~ -1~~ g(l§'-g"l) - (3(5"). (3.16) The first term of the right hand side in Eq. (3.16) describes the single inclusion’s polarization induced by the applied field Eo(r). The second terms imply that the modification of the { ui(r) } arises from the interaction among the inclusions. As an important issue, let us consider the effects on the ith inclusion, as arising from an external field due to the other inclusion. In other words, the inclusions j at i set up some field that the ith inclusion feels. The i(r,r') operator gives the response of the ith inclusion to this .1. field. Eq. (3.16) is an integral representation for the { ui(r) }; given in terms of the one-inclusion t operator, the dipole propagator c; and the external applied field. In order to obtain { ui(r) ), if we assume §i(r,r') is known as a multipole expansion such as Eq. (3.13) and if we express ui(r) as a multipole expansion, a) u.(r) = 2 (-1)é (42)" D? 0 VC 5(r-Ri), (3.17) (1:0 ~ ~ ~ ~ ~1~ then we might be able to obtain a set of coupled algebraic equations for the multipole moments of the ui(r). Here the multipole moment .(If 27 is an (6+1)th rank tensor for each Cartesian component wi’t) (u = x, y, z) of the vector ui(r). It is defined by ~~ (i u. = I dr 1.1.(r) r r ------ r (3.18) Where U - ' - - DC each run on the cartesian indices x, y, z and we have 1 set R1: 0. Note that (J? is independent of r and it is a constant. ~ Substituting the multipole expansion of Eqs. (3.13) and (3.18) in Eq. (3.16), and equating terms of equal orders in ( V )Car-Ri) yields the main analytic equation N e _ _ 00 0° a ‘51 - 4m: [ti 5&0 E + t. 2 (’2' 6' 2 g (IRij|)®£l ] (3.19) The first term on the right hand side involves the assumption that Eo(r) is a constant external field, BC. This is to correspond to the standard exper1mental situatIOn. The quantity G (IRijI) depends only on the center to center distance IRUI = I Hi - le between the ith and jth inclusions and is defined as g“ (3.20) ' C as = v v - ’ I O ”51.)“ ~ .. g(lff " IE-I; I=IR I It is a rank ( C+C'+2) propagator tensor. 28 The reason that we use the multipole expansion in this problem is to reduce the coupled integral equations to coupled algebraic equations which can be solved easily. In practice, the problem becomes numerically large, since for N inclusions, we have to solve the linear algebraic equations, and to get convergence C may have to be quite large. In order to relate ( “i0 ) to the Ee’ we need to define a macroscopic parameter, 9, called the average polarizability tensor, whichconnects the mean polarization to the external electric field. The mean polarization is written as 1 N 0 av = 9.510 (3'21) 2 .. 1 Due to the dipole propagator’s 3 range, the polarization depends on r the macroscopic system’s shape. .9 is a shape-dependent polarizability tensor, but 6e cannot be shape dependent”. Macroscopic electrostatics provides the connection between 9 and 864°. F‘ixman“1 has provided a clear discussion of this connection. For a spherical macroscopic system, the relation is given as : —— (3.22) 29 1 where Q = 3 Tr( Q )- Other macroscopic geometries lead to other functional dependences of 5e on 9, but 9 also changes and the result will always be that Ce is independent of the shape of the macroscopic sample. The value of Se is determined as well by this relation, since we can evaluate Q from Eqs. (3.17), (3.19) and (3.21). 3.2 DIPOLE TENSOR W (:4) When we attempt to calculate the effective dielectric constant of a system, we must solve the set of equations Eq. (3.19), which includes the dipole propagator goo acting between the inclusions. A source of difficulty is the long-ranged nature of the dipole field. Since these equations involve the long-ranged dipole-dipole interaction operator, 1 W (T), summation over the inclusion separations leads to a ~~ conditionally convergent, sample shape-dependent sum. Qualitatively, if a divergence can be subtracted from the given conditionally convergent integral to get another conditionally convergent integral whose integrand has the same asymptotic behavior as the integrand of the original integral, then the ensemble average of the given integral is unchanged, but the divergent behavior of the integrand is removed. The result is well-defined. For a finite system the dipole tensor is well defined. However, 30 in order to coincide with the thermodynamic limit, N —) 00 V -) do and N/V = constant, the present consideration must apply to an infinite system. To include the long range interactions, we calculate the interactions of the dipoles and higher-order multipoles under periodic boundary conditions, where a primary cell is periodically imaged to form a large spherical (or ellipsoidal) supercell surrounded by vacuum (see Figure 3.1). Macroscopic electrostatics can be applied to the supercell. We proceed to sum the interactions by the Ewald method 4, and to obtain the effective dipole propagators by the lattice sum approach“. The dipole propagator evaluated at the positions i and j is l (3.23) I with I fijl : | {j - Si' . As we develop the dipole propagator over all the images, we can express the dipole propagator as the following ' 1 I:Z{_vz[—_—]} (3.24) n ~ “311' Here,|Rij|=|rj-ri+n|,andnznx1+nyj+nzk. Herenfiny, nz are arbitary integers, and n = 9 corresponds to the primary cell. 31 Figure 3.1 Spheroidal macroscopic sample consists of supercells and the center heavy line is the primary cell. 32 The prime on the sum indicates that if i = j, the n = 0 should be omitted. 1 By applying the two grads on —-— in Eq. (3.24), we obtain I 13,, I . I 3(rij+n)(rij+n) T : E - «125) - |r..+n(3 Ir..+n|“S n ~IJ ~13 In such circumstances, there are two issues: (1) The existence of the sum and (2) The shape-dependent effects of the thermodynamic limit. 1 As we can see, I has a singularity at RU: 0 and the integrand 3 ‘ * |R..| ~1J is shape-dependent at Rij -2 the boundary. The lattice sum of the dipole propagator in Eq. (3.25) is conditionally convergent. Both sums in Eq. (3.25) diverge independently. So it is important to use a summation method to adjust I the cancellation of the divergences‘ Let us multiply the terms of the conditionally convergent series by a convergence factor exp(-s|n|2). We may write Eq. (3.25) as ' 1 3 ( Eij + “ ’( Eij + n ’ -s|n|2 I = 2 3 - 5 e «126) ' Ir..+n I |r..+n| . n ~13 ~13 33 Then we expand these sums as power series in s and subsequently take the limit 8 —> O to obtain convergent results. We start with a general form Then Eq. (3.27) multiplied by exp(-s|nlz) is 2 ’ . _ ’ -2x -s|nl ¢(rij,x)-§|rij+nl e . Henceforth, lim 4) ( ri., n, x ) = ¢( 151., n, x ). s--20 ” J “' J We introduce the identities”. ' 1 ' m 2 -22 _ Z-l -R t 0 By using Eq. (3.30), we express Eq. (3.28) as representation ' ' 1 0° x-1 -|r +n|2t -sInI2 ¢(rij,x)'£ F(x) dtt 9 1J e n o (3.27) (3.28) (3.29) (3.30) the integral (3.31) Another important identity is the Jacobi transformationnb. This 34 transformation serves to invert the lattice space sum from the coordinate space to the reciprocal space, wherein the lattice vectors arek=2nn. Z 2 Z . Z e-|;~~+n| t : 3/2 X e( -n n /t + 1 2n n-r L (3.32) Applying the Jacobi transformation, we rewrite the integral of Eq. (3.31) (assuming rm: 0) as 11> I 1 n ¢1 r... no X ) 2 X r(x) J dt tx'l ( t )3/2 D 2 2 . 2 e‘ "' n A I 1 2" “'Elj’] (“I") (3.33) x-5/2 is an analytic function for x > In Eq. (3.33) the integral of t 3/2. For x s 3/2, there is a singularity as t —9 0. We split the integral in Eq. (3.33) into two parts. The integral of the first part is from a2 to 0 , with a2 an arbitrary value. The second one is over the range (0 , <12)- Now, we are able to integrate Eq. (3.33) over ((12, 00) and take the limit s -b 0 without difficulty. The divergence is due to the second integral at t = 0 so that in k space we subtract this term at k = 0 (n = 0) from the integrand and add it on again as a separate term. 35 In other words, the divergent behavior is due to the lattice sum over the whole r space at large r, k —4 0 in reciprocal space, which causes the shape-dependent effect of the problem. At n = O (k = O) the ~ transformed sum yields the shape-dependent term. For instance, for a spherical supercell, the shape-dependent term is 4g 1,. , as r —P 00 . For PU: 0 , then n must have a nonzero value. But it is necessary to include n = 0 term during a Jacobi transformation, and to do that we add a contribution from n = 0 in r space in the second part, and subtract again the same term , which can been evaluated to produce 3 4 O: the self-term operator of the inclusion, 1 , at r = 0 . 3(11' " From Appendix A we obtain the two different dipole propagators: (1) 1‘ diff is the dipole propagator for the different dipoles, rij # O. and (2) 1‘ self is the self-dipole propagator for the dipole and the image dipole, til: 0. l 3 diff =2 { 3 B 0.1, periodic boundary condition must be used to minimize boundary effects. We must: (1) set up a cell with length, L, containing N inclusions. (2) duplicate such a cell -— called the primary cell — in all directions to provide the conventional periodic boundary conditions. The initial configuration is established by placing the inclusions randomly inside the primary cell or choosing a lattice structure, i.e. simple cubic or face-centered cubic, and then equilibrating such a configuration by using the Metropolis algorithm“ to correspond to a hard sphere fluid. The Metropolis procedure consists of moving an inclusion to a new position at random. The new position is accepted only if the ith inclusion does not overlap with any other inclusions already placed in the primary cell and its periodic images. The density is conserved because when an inclusion moves out across the surface of the cell, then another will move into the cell from the opposite surface. Thus if the position of an inclusion’s center is outside the primary cell, it is brought inside the primary cell by translating it by :t L, where L is the length of the cell. This procedure of randomly moving, checking overlap and translating into the primary cell is repeated for each inclusion inside the primary cell. One attempt to move each inclusion is called one Monte Carlo step (MCS). 41 In order to provide equilibration, the acceptance rate for the Monte Carlo moves is chosen to be about 50% and the maximum step size for one Monte Carlo step is determined. If after such a move an inclusion happens to overlap with another inclusion, we place it in its original position. For higher volume fractions ¢ 2 0.3 the configurations are started from a lattice structure, and run for several million 6 x N MCS) to equilibrate. For the volume configurations (typically 10 fractions between 0.1 and 0.2 the configurations are started as generated from the immediate method, and equilibrated for 10:3 x N MCS initially. After such an initial equilibration, these configurations have been selected at intervals of 102 x N MCS apart. This is in order to establish a Markov chain of configurations, thereby eliminating correlations among the configurations. To check that these configurations correspond to a fluid, the pair distribution functions for these configurations have been calculated. A similar case, the pair distribution function for a hard sphere fluid has been evaluated by Alder et al.45 in 1954, and by Henderson“5 to higher precision. The following outline describes how to measure the pair distribution function by the Monte Carlo method: (1) In a generated configuration we choose each inclusion as the central inclusion and evaluate, ni(r,dr), the density (3) (4) 42 of the inclusions in the ith shell of radius r and r+dr centered around the central inclusion, normalized by the total number of inclusions in the bulk. The pair distribution function is defined as 1 m g(r) = 75- .Zlni(r.dr) (4.1) I: where m is the number of configurations. The number of configurations must be rather large in order to reduce the statistical error. Repeat step (1) In times, and sum ni(r,dr). Average the sums over m to obtain g(r). The comparison of our configurations with Henderson’s result“5 is shown in Figure 4.1. This shows that our configurations are correct. 4.1.3 THE DISORDER OF CONFIGURATIONS The effects on 5e of other types of disordered configurations have been investigated (see section 5.2). We will consider: a) RANDOMIZED LATTICE At volume fractions below 0c: 0.524 (simple cubic lattice close packed volume fraction), the inclusions are initially placed on a 43 l ' T I ‘ 1 .8 ‘ 1 100 particles L 1.54 . ‘ 4 vol. fraction-0.2 . J 500 configurations I 1.3-4 1 A < b J 1 v 1 o 4 4 l .0: 1 J 1 0.8-4 .1 I A exact 1 i e Henderson IS Ti V r I o o i 2 3 4 r/Za Figure 4.1 Comparison of the hard-sphere radial distribution functions between our liquid-like structure and Henderson’s results. 44 simple cubic lattice and are then displaced randomly within the Wigner-Seitz primitive cell47 centered on their respective lattice cells. Here, a DS parameter is defined to control the degree of randomness of a configuration. In a unit cell of side 1 let max = 0.5 - a, where a is the radius of inclusion. DS is the actual displacement allowed in the randomization around the lattice position divided by max. DS = 0 corresponds to the periodic lattice and DS = 1 the maximum randomness. When the volume fraction is above 0.524, we replace the simple cubic lattice by a face-centered cubic lattice (fcc) and use the same procedures. Note that the Wigner-Seitz primitive cell of the fcc lattice is a rhombohedron. For convenience, we choose the biggest sphere with 00: 0.74 (fcc close-packed volume fraction) instead of the rhombohedral cell. b) VACANCY LATTICE Here we start with a lattice structure as in a). The inclusions are deleted randomly and the size of the remaining inclusions is expanded to regain the original volume fraction. The original volume fractions must be smaller than the close packed volume fractions of the lattice structures. 45 4.1.4 MINIMUM IMAGE METHOD In the previous sections, the system has either a spherical boundary condition or a periodic boundary condition. In our calculations, we have extended the periodic boundary condition by using a minimum image method to evaluate the induced multipole moments for each inclusion in the primary cell. This method is to generate the inclusions in a primary cell and duplicate its identical replicas throughout space. Then, the mutual interaction between ith inclusion and jth is not restricted to be within the primary cell. The ith inclusion is only allowed to interact with those inclusions in the primary cell or the periodic replicas, which are inside a radius R z L/2 of a sphere centered on the ith inclusion -- image sphere. This method corresponds to setting a spherical cutoff, R, for the pair interaction range. We have measured the effective dielectric constant followed by changing the size of this image sphere: 0.5L, L, and 1.5L. The results show that the effective dielectric constants seem very stable with any size of image sphere. This method provides accurate values for the induced multipole moments. 4.2 COMPUTATIONAL ALGORITHM Now, we have to evaluate the induced. dipole moment and higher multipole moments for each inclusion in the primary cell by using the minimum image method. We repeat this method for each inclusion inside in-.. _ 46 the primary cell to yield a set of 3n N linear algebraic equations for the induced multipole moments )4? on each inclusion. Here n : t’. +1 is the nth multipole moment, and N is the number of inclusions inside the primary cell. The solution of the linear equations, Eq. (3.19), for the multipole moments depends on the shape of the system, due to the dipole-dipole interactions. That is, the second term on the right hand side of Eq. (3.19) is conditionally convergent when the summation its extended to an infinite system, for example, by using periodic boundary condition. We have discused this problem in Chapter 3 and Appendix A. We have evaluated the dielectric constants of a composite material at the dipole level approximation using our minimum image method and compared the results with those of Cichocki and Felderhof 48, who used a version of the lattice sum method, and found that the results are in excellent agreement. We will show this comparison of our results with those of Felderhof and Cichocki"a for the frequency dependent dielectric constants of a composite composed of metallic inclusions in an insulating matrix evaluated at the dipole approximation level in Chapter 5. The minimum image method is considerably faster than lattice sum method; results for the system in Chapter 5 are in accord with the minimum image method. We use it for the simulations presented in our work. In Order to get dipole moments and higher multipoles, the 47 propagators, Ga (Irl), must be determined. We discuss the construction of these tensors in Appendix B. Now we must relate the microscopic dipole moments to the macroscopic polarization of the sample. To do so we inscribe a sphere of radius L/2 in the primary cell. The polarization of this sphere is defined by the volume average of the sum of the dipole moments of the inclusions inside the sphere. The macroscopic polarization is then obtained by averaging the polarization over a suitable number of configurations. Then, the effective dielectric constant can be evaluated with the use of Eqs. (3.21) and (3.22). We have tested the convergence of the simulations by changing the number of inclusions in the primary cell. The dielectric constants from three sets of numbers of inclusions: 25, 125 and 256 are nearly independent of the number of inclusions. The spatial distribution of the polarization was also monitored. The fluctuation of the polarization is quite small throughout the space, even up to the edge of the cell, indicating that the boundary effects are under control. Computer programs have been written for the numerical simulations. A listing of the programs is in Appendix C. These programs are written in FORTRAN. The program STDIEL is used to compute the static effective dielectric constants (used in section 5.1). This program is written at 48 the multipole level, expanding n up to 5. The multipole moments on the inclusions are evaluated by the minimum image method. The matrix equations are solved by using the LINPACK facilities. It is essential to have sufficient memory space for storing the large matrices generated by the multipolar expansion method. The program FRQDIEL was written for the calculations of the frequency dependent effective dielectric constants (used in section 5.2). This program is constructed at the dipole level. There are two subroutines called LSMATRIX and MIMATRIX. LSMATRIX is used to compute the dipole moments by the lattice sum approach. MIMATRIX uses the minimum image method to compute the dipole moments. In the program FRQDIEL complex matrices are created. Therefore, we solve these linear complex equations by calling CSOLVQ, which is one of the routines in the mathematical library (IMSL) of the FPS-164 Floating Point Systems Attached Processor. The statistical error analysis is addressed in both of programs. The principle of the error analysis is the application of the statistical standard deviation. For instance, if the measured quantity and the number of measurements are respectively x and m, the statistical error is defined by“ m 1/2 ABS( 2 x2 -m z ) 5x ”'1 (4.2) m(m-1) 49 In Here (x) = 2: xi /m. The larger the ID value the smaller the error. In i=1 our calculations In is chosen between 100 and 200, and our results are quite stable. The algorithm of the simulation procedures is summarized below: (1) The configurations must be generated with the information on the number of inclusions, the number of configurations, the volume fraction, the radius of the system, the radius of the inclusions and the positions of the inclusions in x, y, z stored in "TININP". (2) The parameters must be set: the number of inclusions, the number of configurations, the radius of the inclusions, the radius of the system, the volume fraction, the s dielectric constants : c 2 and cm, the magnitude of 1) the applied electric field, the convergence factor a, the frequency scanning range ( from BF to FF ), -- - ° and so on stored in "TININZ". (3) Data from (1) and (2) is read. (4) The 6 matrix and the constant parts, B, are created, i.e. Q - u = B and the matrix equations are solved. to obtain all the multipole moments on the inclusions. (5) For each run all of the dipole moments of the inclusions are summed. The above steps from (3) to (5) are repeated In times to accumulate the total dipole moments, then the average of the total dipole moments for each run is taken. (6) as is calculated according to Eqs. (3.21) and (3.22), then the error analysis is performed. For convenience, we set the side of the system to be unity, and fix the number of inclusions to vary the radius of inclusion in accord with the desired volume fractions. In these programs the radius of image sphere is set to be half of the side length. The x-component of applied field is a constant, and the others are 0. The parameters have been set in Chapter 5 according to these considerations. 5. RESULTS AND DISCUSSION In this chapter, we will present results on the effective dielectric constant of a composite material that are obtained by use of the analytic-simulation method. We simulate two properties of a composite material, the static and frequency dependent dielectric constant. To address the former problem, a multipole simulation of the static effective dielectric constant Ce is carried out. We study the effects of multipole moments on the effective dielectric constant for a uniformly-conducting inclusion case, the inverse case of insulating inclusions, and for coated inclusions. Then we compare the results obtained from the analytic-simulation of the uniformly conducting case with a random walk method‘g’ 50. Another problem~ is the calculation of the frequency dependent effective dielectric constant, 68(0). Here, we will investigate the effect of disordered configurations on the frequency dependent effective dielectric Constant. Thus, we consider systems, including a randomized lattice, a vacancy lattice, random size-distributed fluids, and real silver in glass. Finally we will compare the effective dielectric constants, 88(0)) which result from both of the methods, 51 52 the minimum image method and the lattice sum approach, discussed in section 3.2. The results presented below demonstrate the use of these different models. 5.1 THE STATIC DIELECTRIC CONSTANT We consider a system of N non-overlapping, spherical inclusions of radius "a" and dielectric constant £1 embedded in a matrix phase of dielectric constant 82 all within a volume V. In the current problem, we attempt to calculate the effective (dielectric constants of composite materials accurately, discuss the effects of multipoles on the Ee’ and determine how the convergence of .the effective dielectric constant depends upon the various physical models. 5.1.1 CONDUCTING INCLUSIONS We first study perfectly conducting inclusions embedded in an insulating matrix. In other words, the dielectric constant of the inclusions is much larger than that of the matrix. A typical case, for example, is a cermet, a ceramic matrix with metallic inclusions. Analytic results of a similar multipole expansion method for 6e have been obtained by Dukhin‘51 and McPhedran19 for lattices. The computer simulation of the conductivity of the simple cubic lattice by McPhedran19h has a good convergence. Their calculation is restricted for m = 0 component among the ( 2£+1 ) moments for order 4. McPhedran 53 shows by direct comparision that the m¢0 components do not contribute to the dielectric constant for cubic lattices, though no proof of this observation is available. For a random structure we find that the m at 0 components must be included. Since our interests are focused on random systems with liquid-like volume fractions we must consider many multipoles with all the m values included. We choose the dielectric constant of the inclusions, e to be 1. 108 and set the dielectric constant of the insulating matrix, 62, equal to 1. The volume fractions are chosen to produce fluid structure. The range of the fluid structure is defined from a extremely dilute gas to a liquid-solid compatible state, where the volume fraction is about 0.45. We have implemented the calculations for multipole moments € from 0 up to 4. For convenience, we set C + 1 = n so that, for example, I) = 1 corresponds to the dipole term. We carried out the simulations at seven volume fractions: 4) = 0.001, 0.01, 0.1, 0.2, 0.3, 0.4, and 0.45. The results for 5e as a function of the number, n, of multipole moments versus the volume fraction 4 are shown in Figure 5.1. These have been compared with the convergent results, corresponding to n —) 0° for a simple cubic lattice, calculated by McPhedran and McKenzie‘gb, which are indicated by the solid curve in Figure 5.1. As can be seen from Figure 5.1, for O below 0.3, the interaction among inclusions can be treated with dipole-dipole interactions only. Thus, the calculation at the dipole level is enough to obtain an GD Figure 5.1 The effective dielectric constant (Ea/e:2 for conducting inclusions (81:108) embedded in an insulating matrix (8231) as a function of the volume fraction O for n=1(o), n=2(A), n=3(D), n=4(0), and n=5(V). The error bar is indicated. These results have been compared with the exact results, indicated by the solid curve, for the simple cubic lattice (11:0) calculated by McKenzie and McPhedran‘”. 55 accurate effective dielectric constant. As O increases above 0.3, the higher order interactions must be included. The results show that multipole effects are important at high volume fractions. At O higher than 0.45, the n = 5 results do not converge. This conductor in insulator case is the most difficult to converge as the electric field in the gap between a pair of nearly touching inclusions has a very rapid spatial variation. 5.1.2 AN INVERSE CASE Figure 5.2 shows the inverse case of the previous uniform case, where the system consists of poorly conducting inclusions in a conducting medium. The ratio of the effective dielectric constant to the dielectric constant of metallic matrix, eels:2 versus the volume fraction O from 0.05 to 0.7 has been presented. The reason we can calculate the Ee at higher volume fractions in this case is that these insulating inclusions avoid a divergence of 8e' Physical models include impurities in metals. Comparison of the variation of 66/52 with increasing n indicates that for n = 5, even when 82/51 -—> 0) , quite good convergence is obtained even for O near to the hexagonal close packed volume fraction OC = 0.74. In our calculation the highest O is chosen to be 0.7, which is between the random close packed OC = 0.63 and the hexagonal close 1 Off y . . - . . - - , J 4 J I : O.8(l . J - 4 N J 1 co 1 . \ 0.6+ . . (D 4 ‘ (O l e . ‘ e . 0044 ‘ T J J l o . J 4 02 i-e--eerrxxrggg 0.0 0 2 0.4 0.6 0 a Figure 5.2 The effective dielectric constant 88/62 of a system with poorly-conducting inclusions (81:1) surrounded by conducting medium (827-108) as a function of the volume fraction O from 0.05 to 0.7 for n=l(o), n=2(A), n=3(D), n=4(o), and n=5(V). 57 packed OC = 0.74. The difference between the results with n = 4 and n = 5 is very small, and cannot be resolved in Figure 5.2. From Figure 5.1 and Figure 5.2, explicitly, the effective dielectric constants of the inverse case converge faster than those of the conducting case. The rate of convergence of 5e with n exhibits an even-odd pattern at high O. It is important to realize that, in principle, the high .density liquid structure is quite symmetric, and therefore most of the contribution to the polarization is from the odd n multipoles. For example, in a simple cubic lattice, due to the symmetrical structure, the dipole moment at a given inclusion will have equal in magnitude but opposite in sign contributions from symmetrically located pairs of inclusions when n is even. 5.1.3 A SPECIAL CASE From the results of the uniform case, when 81/622 —) 00, the effective dielectric constant will diverge at high O. It is interesting to study how the closeness of. the inclusions affects the effective dielectric constant. We generate two sets of configurations: one set at volume fraction O = 0.4 and 0.6, and the other initially at O = 0.5 and 0.7, but then the inclusions are shrunk sufficiently to adjust the volume fraction to be O = 0.4 and 0.6. A characteristic feature of these two sets of configurations is the different value of the closest distance between inclusions. 58 We compare the results from these two different sets of configurations in Table l with 51 : 100 and £2 = 1. Clearly, by directly generating configurations (the former method), the inclusions can get closer than in the latter method. Thus, the effective dielectric constant Ce evaluated from the former configurations should be higher than that from the latter configurations, as is borne out by the data. Also, note that at the higher O value, the difference between n and n+1 is much larger than low O illustrating the non-convergence of the multipole method at very high volume fractions. This is similar to the case of oxide-coated metallic inclusion in insulating matrix which we consider in the next section. The conclusion there is that a thin layer (i.e. oxide-coated or insulating material) is capable of preventing the divergence of the effective dielectric constant arising from touching inclusions. 5. 1.4 COATED INCLUSIONS An interesting extension of our methodology is to consider a problem of composite composites. In this section, we investigate the effect on the 8e of a coated-metal inclusion with an arbitrary thickness of coating material. That is, we consider spherical inclusions which are composed of spherical layers. For example we consider the modification of the effective dielectric constant by coating an inclusion with a thin layer of material. The conductivity of the coating material is relatively Table l 0.6/0.6 set sets are generated at O 56 e/ 2 = 0.5 and 0.7 and the inclusions are then two sets of configurations: is generated as for Figure 5.1. shrunk to correspond to O = 0.4 and 0.6, respectively. The 0.4/0.4 and The 0.4/0.5 and 0.6/0.7 n=l n=2 n: n=4 nzo 0.4/0.4 2.9848 3.0705 3.1527 3.2069 3.2333 0.4/0.5 2.9562 3.0107 3.0631 3.0969 3.1121 0.6/0.6 5.4109 5.4764 5.6676 5.7692 6.1042 0.6/0.7 5.4121 5.4350 5.5822 5.5899 5.8884 60 small compared with that of the medium, and the interior of the inclusion is a good conductor. For parameter values we choose 51 : 108, e = 1, and Em : 103, where £1, 8 refer, respectively, to the 2 2’ 8m dielectric constant of the metal inclusion, coating and medium. Now, if the multipole polarizability for this composite inclusion is known, then we can solve the problem. Let us define that the inner radius of the composite as a, corresponding to a metal inclusion, and the outer is R, for the metal inclusion and coating. The expression for the multipolar polarizabilities of a composite inclusion has been given by Maxwells as 2441+ 2104-1 ((1)-INN {+1 )+ (40R a p—u)(£+ u( (+1))a - {u( (1+1 ) + “CHM {+1 ) + 120+ a 01H wu)(1—u)(a/R) “4 2t+1 (5.1) Here )1 : (El/eIn and u : ez/em. In Table 2 we list the data from the simulations, the parameter values, and the volume fraction Oi and O0 corresponding to the inner and outer sphere volume fractions. The first three data sets are for a/R ~ 0.9999, a thin coating. Comparison of Figure 5.1 (uniformly-conducting inclusion) and Figure 5.3 (this case) shows that a thin, protective layer rapidly reduces the dielectric constant relative to the uniform inclusion case. These simulations also describe the effective thermal conductivity of composites. This simulation method can be used to predict the degeneration by oxidative coating of metal-insulating matrix composites. Meanwhile, note that the multipole expansion for the composite-composites converges at much higher volume fraction than those for the uniform composite. 61 Table 2 86/8“) for composite-composites. The first three data sets are for a thin insulating coating; the final two are for a thick -). insulating coating, as discussed in the text ()1 = 105, u = 10 Oi/Oo n=1 n=2 n: n=4 n=5 0.3999/0.4 2.4004 2.4263 2.4446 2.4538 2.4569 0.4499/ 0.45 2.7490 2.7838 2.8138 . 2.8315 2.8384 0.4999/0.5 3.1615 3.2067 3.2531 3.2844 3.2998 0.3/0.45 0.4567 0.4545 0.4510 0.4487 0.4475 0.4/0.45 0.4695 0.4666 0.4633 0.4611 0.4602 62 00m Figure 5.3 The effective dielectric constant as”? of metal inclusions coated with an insulator as a function of the volume fraction O for n=l(o), n=2(A), n=3(D), n=4(O), and n=5(V). See Table 2 for the chosen parameters. 63 The last two data sets in Table 2 correspond to thicker coatings a/R = 0.9615 and 0.8735, respectively, where the composite-composites are now clearly acting as insulators in the sense that ee/em < 1. Conditions on the values of u, u, and a/R were previously obtained by Sridharan and Cukier15 which indicate which regime (se/ern < 1 or > 1) composite-composites would fall in for the second virial coefficient. The results presented here obey those conditions. The results displayed in Figure 5.3 show that a thinly coated inclusion permits the multipole expansion to converge quickly in comparison with a uniform inclusion. The thickness of the coating controls the composite-composites behavior as insulator or conductor-like. 5.1.5 CONCLUSION For the static case, if we know the expression of the n rank tensor and the analytic form of the multipole polarizabilities, then the analytic simulation of the higher order multipole contributions to the effective dielectric constants of a composite material can be carried out straightforwardly. By adding as many multipoles as we can, we can obtain more accurate results for the effective dielectric constant. However, the application of the multipole simulation method runs into the following problem. In Appendix B we have represented an n rank tensor as a power series in the reciprocal distance between inclusion pairs and summarized it into a general expression to arbitrary 11 value (n z 0). 64 1.OE+04-.-.-....I- %‘ 8000.0“ 4 C 0 4 8 (D 6000.04 .1 C C/ 1 . (D g 4000.04 } *J < 8. O 2000.0( 1 1 .0« r—T.rrew—+Tr O 0.0 10.0 20.0 30.0 40.0 50.0 60.0 Dimension of array Figure 5.4 Plot of the dimension of the matrix versus the total CPU time. . I V I V I ‘ I v ' v ( B . 1 1.0 , x , ¢ 4 4 J J 4 a 4 1 1 v f V ‘r 0.5 007 0'8 O Figure 5.5 The effective dielectric constant 86/82 versus a function of the volume fraction O for the different types of configurations. The solid line is generated by the multipole method of Ref. 19b for the simple cubic lattice. For the Metropolis (randomized lattice) configurations the random walk method is denoted by A (x). For the face-centered cubic lattice the random walk method is denoted by CL The o’s are generated for the Metropolis configurations by the multipole simulation method (see section 5.1.1) for a conducting case. 66 Thus, it is not difficult to figure out. But, because the computation involves complicated matrix algebra, the computing time and computing space of the computer system must be of concern. Figure 5.4 shows the total CPU (Central Proceding Unit) time versus the dimension of the matrix of these calculations with 25 inclusions and 100 configurations. It demonstrates that if the higher multipoles are included, the CPU time increases parabolically. Due to the non-convergence of 8e for the conducting case, the comparison of Figure 5.1 and the 8e resulting from a random walk methodsz has been shown in Figure 5.5. The benefit of the random walk method is that it can be applied to higher O than a multipole simulation. As can be seen from the above, the multipole simulation of the conducting case is difficult to converge at O higher than 0.45. From Figure 5.5, at low O (~ 0.3) the multipole simulation (n = 5) and random walk method results are in quite good agreement with each other. It also shows that a multipole simulation is sufficient to obtain an accurate result. Thus, at high O the use of the random walk can supply the data not available from the multipole simulation. 5.2 THE FREQUENCY DEPENDENT 86(0)) We have studied cam) of a composite material consisting of inclusions of dielectric constant 81(0)), which have a Drude model resonant form, embedded in a background medium of frequency-independent dielectric constant 62. 67 Our interest is to find how the optical properties of such materials are affected due to the electrostatic interactions among the inclusions. The frequency dependence of the effective dielectric constant as a function of the inclusions’ geometry, size-distribution, spatial distribution and concentration is explored. Our method is the analytic simulation method described in Chapter 4. We assume that the quasi-static Maxwell equations are valid. Because the wavelength of the applied field ( ~ 3600 X ) that we choose is long compared with the inclusions’ size ( ~ 100 I: ) and the mean distance between inclusions, the electrostatic approximation with 61(0)) can be used here. The multipolar effects on 56(0)) at high volume fractions of inclusions have been inVestigated in several hypothetical cases. From the previous static multipole simulation, at low volume fraction the dipole approximation is enough to obtain 66 quite accurately. Thus, a comparision with the Clausius-Mossotti formula, which applies to dilute inclusions, is shown in the Figures. In the following, we have set the plasma frequency w = P 9.4x1015s‘1, the damping constant 1" = 10143-1, £100: 1 and £2: 1 for a hypothetical case. Only the resonant, imaginary part of 68(0)) is displayed in the Figures. All of the parameters are defined in Chapter 2. 68 5.2. 1 RANDOMIZED LATTICE In Figure 5.6, sew) of the randomized lattice is compared with 88(0)) of the fluid structure at volume fraction O = 0.03. The parameter DS is chosen to be 1 or 3/4. DS = 1 corresponds to the maximum randomness for this method of randomization. Each curve corresponds to an aVerage over 200 configurations. They demonstrate that the fluid structure broadens the spectrum. That is, each inclusion is presented with a different local environment. Since the Clausius-Mossotti theory predicts that the plasma resonant peak position (3 depends on O (see Chapter 2), the local environments associated with slightly different local densities lead to a broadening of the resonance. An important feature of the definition of (3 is that as O increases, the resonant peak has a redshift. Similarly, this is the case for a randomized lattice. DS = 0 for the randomized lattice corresponds to the periodic case. We then certainly obtain the Clausius-Mossotti result. With some randomization of the configurations the resonant absorption peak is broadened. The reason has been explained previously. The broadening increases with increasing randomization, i.e., increasing DS value of the configurations. From Figure 5.6, for the maximum possible excursion, DS = 1, the broadening is still considerably less than that of the fluid. The broadening for the randomized lattice is very symmetric, which indicates that on average the local fluctuations in (3 are themselves symmetric. 69 6.04. . ,fiTT, ”Irflrrrfij . - —-CM , 5'07 s 05 -1 . j ‘ 401 g (on) 1 4 3.01 .E j 1 204 1 i 104 1 1 0.0 are. we.“ WW If“ Frfr 0.545 0.555 0.565 0.575 0.585 0.595 w/wp Figure 5.6 Im 56(0) for fluid .(A) and randomized lattice (op). The solid line is the Clausius-Mossotti expression. See text for the definition of DS. 5.2.2 VACANCY LATTICE In the following, we apply the same parameters to a vacancy lattice. Figure 5.7a and 5.7b presents the 68(0)) for a vacancy lattice. There are four cases we will discuss. First, we start with 125 inclusions on a simple cubic lattice and randomly delete inclusions until there are the number remaining indicated on the figures. The remaining inclusions in the cubic box are expanded to obtain the original O. A second run starts with 343 inclusions with 63 remaining. As the system becomes more dilute in terms of the number of inclusions, the spectra develop a bimodal distribution, with the weight of the higher frequency part growing at the expense of the lower frequency part. In principle, at high O the Clausius-Mossotti formula for the effective dielectric constant of composite material predicts a redshift of the resonant peak. In other words, the lower (a peak corresponds to more dense regions. Hence, it appears that such systems consist of low and high density regions, and this is responsible for the development of the bimodal distribution. The notation in the Figure, 40/125, indicates that 125 inclusions were present initially in a cubic box and then inclusions were randomly removed until 40 were left. All results have been compared with the Clausius-Mossotti formula with the same parameters. Each curve corresponds to an average over 200 configurations. 71 6.0 r—yy—rvfi‘r‘rr‘vrfvrv'rri —- CM , 95 = 0-03 0 75/125 .. 5-0‘1 A 0 55/125 4.0- a’ - q, i 2 3.0—1 ' " O _. ‘ . e 2.0-1 "3\ ‘ .l \l' 3 1.0-1 I. "‘ .l 'I' . 040 fff rfi V ' VW—f Fr r ‘— 0.1 0.545 0.555 0.565 0.575 .5 w/wp 0.595 Figure 5.7a Im 56(0)) for vacancy lattices. The notation 75/125 indicates that of 125 inclusions on a cubic lattice, all but 75' were randomly removed. 72 60 r l‘rjlr' Irvlrvr 1 . ""' CM .( ‘2 = 0-03 0 40/125 5'0? 63/343 404 4 “0) J 4 3.0- .. .E 4 4 2.0-1 .1 i ‘ i _, ' .. . .1 10 , s\"' q 0.545 0.555 0.565 0.575 0.585 0.595 w/wp Figure 5.7b Im 60(6) for vacancy lattices. 73 5.2.3 RANDOM SIZE-DISTRIBUTED EFFECT . 53 . . . . . . Experimental results on the Size-distribution of inclusmns has been fitted to a log normal distribution 2 1 [ In, _x_ f(x) = exp - (5-2) {21:}:an 2(1n0)2 where f(x) is a normalized distribution in the inclusion radius, a(x) :: xa, u is the mean inclusion size, and o is the standard deviation of the distribution. The expression for the volume fraction due to the size-distribution of inclusions is given as _9_ 4). ) (5.3) ¢:(£3!-)na3u3exp( where n is the number density of the inclusions and 1 is defined as .1. 2 ( In a )2. (5.4) l :: The log normal distributions with different standard deviation values a and with u = 3.3 have been plotted in Figure 5.8. The a = 0 corresponds to equal size inclusions. We consider fluid configurations which are generated by throwing inclusions into a cubic box, while avoiding overlap, and choosing each new inclusion’s size from the log normal distribution until the desired volume fraction (b is obtained by adjusting the parameter, fac, defined as 74 LOG-NORMAL DISTRIBUTION 1000 trrrV—Ij—V—‘TIVUV'j—U—r'r' T . — a- 1.40 --- c- 1.65 n BOO-I -- 0' =- 2.02 --I _3 1 -- v=2.72 5 BOO-I ”I 1 °' ‘. 95 4 Jf\‘ fl = 3.30 q L. 0 .0 E a z Figure 5.8 The log-normal distribution for different standard deviation values a. ~I 0! 1/3 In practice, it is difficult to generate fluid configurations with large 0, because overlap occurs frequently with increasing a. The maximum 0 in our present work is 2.72. In Figure 5.9 the results for equal size and three sets of the different size distributions of inclusions are compared. For the equal size and the a = 1.4 fluids, the results are not significantly different. This a is chosen because the size distributed data were fitted with this a in the experimental work”. If we increase the width to a = 2.02, there occurs a low frequency shoulder relative to the equal size case. It appears that high density regions occur because of large size inclusions. All the results correspond to an average over 200 configurations. In order to investigate the effects of the spatial distribution of these inclusions, we exhibit two different runs at a = 2.02 in Figure 5.10. Each run has 200 configurations averaged. The difference between the random-size distributed and equal-size results is not an artifact of the choice of configurations. There is no background contribution to the imaginary part of 88(0) from the matrix absorption. There are several mechanisms which will lead to such additional lineshift and line broadening provided by the size distribution: (1) The local density fluctuates and results in line-broadening. (2) The 76 1nmtmfit'flrpmmfivmwfi-na I -— W a 4.0 1.65 3.5 3:38 I 3.0 2.5 2.0 1.5 1.0 0.5 0.0 ' 0.45 0.50 0.55 0.60 0.65 0.70 w/wp Illllll D X A qqq Mice "Dr «no Figure 5.9 1111 86(0) for equal-sized (0:0.00) and size-distributed fluids. 77 nmfinvflmmm'm'gamm 4.5 ¢ = 0.03 o a a 2.02“) 40 o a - 2.02(2 35 3o 25 20 15 1.0 3‘ 0.5 83 ‘ 0.0 0.45 0.50 0.55 0.60 0.65 0.70 m/wp LLLLLLLJ Ll.“ Muse Figure 5.10 In com) for size-distributed fluids with the same width, but different sets of configurations. 78 large size of inclusions causes local high density compared with the average density and a red lineshift occurs. Under the condition of such a low volume fraction as we have chosen, the red lineshift is not obvious. 5.2.4 SILVER IN GLASS Now, let us consider a system with a suspension of silver spheres in glass, which is a more realistic model for a experimental investigation”’“. Because bulk silver has a low damping constant, it is used for studying optical spectra. According to the values used in Felderhof’s work”, we choose the plasma frequency for silver as up : 1.46x10163-1, the bulk damping constant as 1“, = 0.24x1014s-1, the interband contribution as rib = 0.87x1014s-1, the Fermi velocity as V F = 1.44x108cm/s, E' = 4.5, and the dielectric constant of glass as 52 : 100 2.25. Then the damping constant 7a of Drude form for a sphere of ° 14 -1 14 -1 radius, 100A, is equal to 1.68x10 8 , the width is I‘ = 2.55x10 8 , the resonance frequency is (as : 4.87310158-1, and dim : 0.16. The results are presented in Figure 5.11. The width parameter of the Drude form we use is larger relative to us than in the preceding hypothetical case (section 5.2). Thus, the peak of the Clausius-Mossotti formula is sharper. Even at 0 = 0.1, the broadening is less than in a hypothetical case at 0 = 0.03. In such a circumstance when the spheres are size distributed with the narrower distribution, a : 1.40, the results are indistinguishable from those with equal size. Even for the wider distribution, a : 2.02, the results are not resolved. _‘ 79 "Ag In gloss T 'VV‘TYUTV “T'TUVVVIUVUV‘V'UVIUUTU'UY'II'V'U'U'U‘ —— CM A a = 0.00 a a =- 1.40 o 0 = 2.02 9“ o ¢=0.l h 9’ o llLlLLLLLl4llllJlLllLlALlLLL Im 5e f O l" o ALLALLLILUJLILIAAIIJLLllLLllJJlLl lllillll P o P -.* 0.2 0.3 0.4 0.5 0.6 Figure 5.11 In 66(0) for parameters chosen to represent Ag in a glass matrix. Equal-sized and two size-distributed configurations are compared. Note that 1‘ depends on the inclusion's size. 80 We found that for the volume fractions considered here the dipole approximation is accurate. We have probed this in send) by including the quadrupolar effects, where they are expected to be greatest, and have found minimal changes with respect to the dipole simulations. 5.2.5 COMPARISON OF MINIMUM IMAGE AND LATTICE SUM METHODS In this section we compare €8le evaluated by using the minimum image and lattice sum approaches. Our computational technique used in the previous calculations of the effective dielectric constant 68(0)) is based on the minimum image method. In addition, we have extended the calculation of the frequency dependent effective dielectric constant by using the lattice sum approach. It is worthwhile to note that the speed of calculation using the minimum image method is considerably faster than a lattice sum approach. We have performed our calculations of 66(0) for silver in glass case at volume fraction 0 = 0.2 using both techniques, shown in Figure 5.12, in the dipole approximation. Here we exclude the interband contribution to the width parameter. Because we want to compare with the simulation results by Felderhof and Cichocki“, we choose the same parameters as their values. Felderhof and Cichocki used the lattice sum approach with 500 inclusions in the primary box, but we use only 75 inclusions. The Cole-Cole plot of the imaginary part of 58(0), Im(€e(u)), versus the real part Re(€e(u)) clearly shows differences between the results from the Clausius—Mossotti expression 81 and those of the various simulations. In Figure 5.12 the large circle follows from the Clausius-Mosstti expression. The small solid curve is a result of the analytic-simulation technique using 75 inclusions; solid dots resulted from the lattice sum approach using 75 inclusions; pluses are points excerpted from Figure 5 in reference 48. The Cole-Cole plot shows that the results from both techniques nearly agree with each other. The spectral shape of the imaginary part [m(£e(w)) using both techniques has also been plotted in Figure 5.13. We find that the analytic-simulation technique yields results that agree quantitatively with those of the lattice sum calculations. Hence, we have demonstrated that the effective dielectric constant can be accurately, directly, and quickly obtained by using the analytic-simulation technique. Similar conclusions have been reached, in a somewhat different context; by Gales et al.5‘, who studied Coulomb and dipole effects in disordered solids. Finally, we note that at 0 = 0.2, dipole level approximations are not sufficient to describe correctly 66(0), or ce' However, the sums of higher multipole terms do not exhibit conditional convergence due to their much shorter range interaction. Hence the special techniques associated with the analytic-simulation method or lattice sum approach need not be used to obtain accurate results for higher multipole interactions. 82 r ,— r r . he ff -5.0 -3.0-1.0 1.0 3.0 5.0 7.0 Re 83/82 Figure 5.12 Cole-Cole plot of the dielectric constant of silver spheres in glass (neglecting interband transitions), at volume fraction “0.2, in the dipole approximation. The large circle follows from the Clausius-Mossotti expression. The small solid curve is a result of the analytic-simulation technique using 75 inclusions; solid dots result from the lattice-sum approach using 75 inclusions; pluses are points excerpted from Figure 5 of Ref. 48. 83 12.0-4 ¢ = 0.2 10.04 8.0-1 6.04 Im 88/82 4.04 Figure 5.13 In com) for monodisperse silver spheres in glass (neglecting interband transitions). Solid curve, Clausius-Mossotti expression; dashed curve, analytic-simulation technique, 75 inclusions; dots, lattice-sum approach, 75 inclusions. 6. CONCLUSION The use of computer simulations provides the ability to look at detailed microscopic properties which are difficult to probe theoretically and experimentally, and permits the use of more realistic models in order to make comparison with experiments. Therefore, as long as the analytic formalism is known, then the computer simulation can give practical results. In the present study we have shown that the effective dielectric and optical properties of composite materials or nonpolar polarizable systems can be obtained by the analytic-simulation method. As long as the multipolar polarizabilities at of the inclusions are known, the analytic-simulation method can be employed routinely to obtain the effective dielectric constant and optical properties of a macroscopic sample. In the static case, we were able to deal with any ratio of 51/62, as well as composite-composites. The drawbacks of the analytic-simulation method are the difficulty of convergence at high volume fraction 0, and the requirement of obtaining c:c analytically. For the most difficult case of infinitely conductive inclusions in an insulating matrix, higher-order multipoles than we have used are required for convergence at volume fractions beyond $0.45. In other words, when the average electric field varies rapidly in space it becomes necessary to include more multipoles to describe this rapid 84 spatial variation. For the frequency-dependent dielectric constant, the interactions among the inclusions lead to the following strong deviations from the Clausius-Mossotti result. The lineshape broadens and shifts to lower frequency. The broadening arises from the distribution of peak frequencies in the slightly different local environments. In our work, the applied field is limited to long wavelengths relative to the size of inclusions, so the electromagnetic interactions among the inclusions can be replaced by electrostatic interactions. It is important to investigate systems with the wavelength of the applied field comparable to the size of inclusions, as this is a common experimental situation. Therefore, in addition to electric effects, magnetic effects should be consideredZ—I’ZS’“. About the calculation of optical properties under this circumstance, one may go back to Mieas. He first found the color variation of colloidal metal inclusions and he derived a theory to calculate the optical properties when the inclusion size is comparable to the wavelength of the applied field. Mie provided the solution to the problem of the interaction of a plane electromagnetic wave with a sphere of arbitrary size. A depolarization effect associated with the inclusion boundary usually causes an apparent size-dependent optical constant because electrons are confined by the surfaces of inclusion and produce large depolarization effects; these effects cause the electron cloud to resonate at a quite different frequency than in the bulk. 86 The problem of the conditional convergence of the dipolar tensors for an infinite system was handled via the Ewald summation method. In our work, the long-ranged interaction of the dipolar field has been overcome by using a minimum image method to evaluate the induced multipole moments on each inclusion. We found that this latter method gave the same results as the former, and it is much faster computationally. Another application of simulation methods is to non-spherical inclusions, which could be shaped like ellipsoids, slabs, needles, etc. For ellipsoids the generation of Metropolis configurations can be developed via a contact functionss. Although the generation of any geometrical distribution of inclusions is not hard, the problem is the difficulty in deriving the multipolar (polarizabilities for these non-spherically shaped inclusions. In these cases, the random walk simulation method is more powerful because it does .not require simple inclusion geometries. We have evaluated the static and frequency-dependent dielectric constant in our present work. For wavevector dependent applied fields and the resulting wavevector dependent dielectric constant, Felderhof .56’57'58 used a diagrammatic expansion to analyze the optical et a1 properties of a random medium. The direct use of the analytic simulation method to study this problem is quite promising. A similar problem is to investigate the localization of an electromagnetic wave in a dielectric medium. Arya et al.59 have discussed the localization of classical waves in a dielectric medium of randomly distributed metal 87 particles. The localization of one electron in a random medium can also be addressed. Finally, chemical reaction problems60 can be equally well solved by the analytic-simulation method. APPENDICES APPENDIX A A. DERIVATION OF THE LATTICE SUM FOR INDUCED DIPOLAR SYSTEMS 3(r+n)(r+n) _:2 - ~ ~ 'SInIZ (Al) I 5 e . |r+n|3 |r+n| where the prime on the sum indicates that if r = 0, the term with n = 0 should be omitted. Let us separate 1‘ into two parts, A and B. I 1 - -s|n|2 , Azi 3 e (A-ZI .. |£+n| Applying Eq. (3.30), the integral representation, gives II :3 M IF 3 r—h—x r——-\ h-fi 8 N + ‘——3 h—J O. H C" 7s— + :3 ('9 g._v__.a I m 5 N (A1) + (A2) (A.3) There are two situations considered: (1) If r at Q, then the sum over n includes 0 and the prime is removed. (2) If r = 0, then n = Q is 88 89 excluded. (Here n = (n , n , n I with n , n , n from -00 to 00 in x z x y z 3-dimensions). IN CASE (1) : m 2 2 2 (A1) :21 {J dt t“2 e-|r+n| t } e-SInl (AA) ’ {n 2 n a . . 112 We substitute u = t lr+nl, 2 '9 uz _ 2 -| '2 (A1):21 J du2-———3-eu e8n (A.5) n ' {n |r+n| g|r+n| and then integrate by parts, subsequently taking the limit 3 -+ 0, 2 l 2 (A1) :21 ___3_ “IE‘I'DI e-(dIB-nl) ' 4' u |r+nl n ~ I n + erfc(a|r+n|) . (A.6) a 2 2 2 (13(2):.21 {J dt t‘” e-lgml ‘ } e‘sllrll (A.7) ' {T n 0 We combine the exponential terms and rearrange as shown tr stlgl2 - r+n2t-sln|2=-(s+t)|n+ . (A.8) s+t|- s+t 90 We can directly apply Jacobi’s transformation (in k space) because 0 is included in 2 . Thus, we obtain n 2 2 a 3/2 (A2):Z 1 JdttUZI “ ) - n s + t n 0 n2 2 i2n(n-r)t stlrl2 n ~ - e(-s+t s+t )e(-s+t ). (A.9) For n = 0 and s —) 0, there is divergence in part (82). We extract n = 0 from 2 (in k space) and add it on again. n Now, (A2) : (A2.l) + (82.2) 2 012 Iran2 (A2.1):: 2 1 I dt t-i n3/2 e( - 4-12n n.§ ) s —-) O 0*? J n O (A.10) We substitute u = t.1 to get on '21! n°r -1 -n2n2u (A2.1) : 2 1 2 n e1 ~ du u e . (A.11) mtg Us‘2 “2 2 1,2 x 3/2 (A22) : 1 dt t ( —_t_ ) - 11—1:— 0 s + nznz i2n(n-r)t st|r~|2 e(_s+t + s+t )e(-s+t :3 II (C 1/2 1! (-———) Jdtt (——_s+t)e s+t (A.12) ‘l/Z I! J dt t I _s—T-t—- I + 0(8) (A.l3) . 1/2 /2 . . . Setting t = s tanG and usmg the trigonometric method leads to 20:2 s s (1122):;21! 1n( )+ln[ 2+1+ 1+ 2] 8 2a 20: a2 1/2 - 2 ( 2 ) + 0(5) 8 + (X 3—9 0 . = ( - 2 n 1n(s/2a2) + 2 1t 1n(2) - 4 n ) 1. (A.14) Now, we pull out the divergent part, ln(s/2a2), and keep it separate. IN CASE (2) : The term (A1) is equal to (A1) from case (1), evaluated with r = 0. Due to the application of Jacobi's transformation, we need to include the n = 0 term in the summation over n ( in r space ) and then subtract it separately. (A2) : (112.1”r : o + (A2.2)|r : o - (112.3) (A.15) ~ ~ a 2 2 2 (A23) :; {I dt tuz e‘IEmI t } e‘SInI (A.16) We follow an analogous derivation for B : ,3(r+n)(r+n) 2 ” " -sIn| e w H ”M |r+n|5 Q 4 . 2 2 3 I r + n )( r + n ){ I dt (,3,2 e‘IEInI t } e’SInI I 3F“— 0 (A.17) Applying the identity, R R : - V6 v5 618-5 5:0, (A.18) let us write m 4 ’ . 2 2 B = _ 2 v5 v5 e-1§.(£+n) J dt t3/2 e-|£+n| t e-slnl I n n 0 5 = 9 B : (Bl) + (82) (A.19) IN CASE (1) : (D 4 2 2 (Bl) : E ( r + n )( r + n )J dt t3"? e'IE+nI t e’SInI n " ” 2 n a 2 S _) O 40:3 e'IaI€+nII = 2 (r + n )( r + n) 2 In n ~ I r + n I 60: e-I“|£+n|)2 + 2 ( r + n )( r + n) 4 ’I n ~ ~ I r + n I n ‘8 erfc( aIr+nI ) +32(r+n)(r+n) 5 (A.ZOI ~ ~ I r + n I n ~ 2 4 ‘1 2 (32) = _ E V V e_i€'(£+n) J dt t3/2 e—Ir+nI t. (T ~5 ~15 n 0 -sInI2 5 =9 Rearranging the exponential parts gives : -IE+nI2t-sInI2-iE-(r+n) tr i5 2 "‘“t” I“ s+t)+-2—(s+—t)——I 52 stIrI2 isE-r 4mm * s + t + s . t (3.22) Taking Jacobi's transformation for (82) yields O(2 4 3/2 1: 3/2 (82):- ZYEYEIdtt (5+t) I 1: n 0 Iran2 i211 n-rt ng-n -P e(-s+t+s+t - s+t)e 5 =9 (A.23) (82) : (82.1) + (82.2) (12 ‘ 4 3/2 x 3/2 (82.1)=- EVEVEIdtt (t) s —> O :I n n¢9 0 Hana nE-n E2 8‘ “2“"? t ’eI‘Tt—l 5 =9 2 2 n n 2 “n" ( ) (12:: n-r) - - -— e 2 e ~ 2 a n¢0 InI «2 nznz + 2 n; 2 I dt t“ e" t "2" “'5’ (A.24) mtg 0 a2 4 n (82.2) : _ v5 VEJ dt t3/2 ( a“ )3/2 -P 1! ~ ~ 0 E : g, n :9 a2 1 stIrI2 _ 3/2 5/2 (- I -21($Idtt (3+t) e s+t 0 a2 1 StIEIZ 2 3/2 7/2 (- ) +4nsr£Idtt (st-t) e s+t 0 s -’ 0 4 l 1 2 =4n1[-—+—ln(2)-—ln(s/a)] (A.25) - 3 2 2 IN CASE(2): (81) : (81) from case (1), evaluated at r = 0. (A.26) (82.1) = (82.1) from case (1), evaluated at r = g. (A.27) (182.2) : (82.2) from case (1), evaluated at r = 0. (A.28) (82)|£:g,n:9 = 0 (1129) Summarizing, T ( case (1) ) : A - B 2 1 In 1 til- I 1: Irft'nI:3 (r+n)(r+n){ 4G3 _ 2 z n - ~ I Z {aIr+nI 940:)an) + 2 erfc(aI r+n|) } e-(o:I1:1-n|)2 6a -(aI'g+nI) erfc( aIan ) + 2 + 3 - 3 I It I r + n I I g + n I nznz «lunn (- ) . - . 4n + ___3_ e “2 e(121! n r) + __ 2 . (A.30) n¢0 Inl 3 T ( case (2) ) : A - B 2 1 .2 n : 1 {aInI e—(aInI) + erfc(aInI) } - [—— 3 n¢0 1: Ir” 96 2 3 -(a|nI) _ _2.’i_{ 4“ g(«lnhz. 6“ e .3 "“‘“'“"I 2 2 3 MO Inl (T (T Inl InI ~ 2 2 nn 3 4nnn (- 2 I 41! 40: + 2 ———e a +( - )é. (131.31) I2 3 3In mtg In The divergent terms in (82) and (Q) cancel each other; thus, the divergence is removed. APPENDIX B B. IRREDUCIBLE CARTESIAN TENSORS The study of tensors is important in electrostatics, electromagnetics, and hydrodynamics. The analysis and properties of tensors are detailes in textbooks61 and the literature”. Here, in order to determine gadr) defined in. Eq. (3.20), we use a potential-theory approach‘sz3 to the construction of irreducible cartesian tensors and focus on the technical development of these tensors. The definition of an irreducible tensor I n I , of rank n , is 2n+1 r r—n r“ =(-1)n ( ) <—). (3.1) (Zn-1) !! 6r r The I I means irreducible. An irreducible tensor of rank n is characterized by being: (1) traceless and (2) symmetric. A tensor is traceless if it vanishes on contraction on any pair of indices, i.e. Z 1 VW- --( V2 T ) = 0, r ¢ 0. For a symmetric tensor it is arbitrary as ~~~ to which indices are used to make the contraction. As we differentiate in Eq. (8.1) to order n, we can represent the 97 98 . I n I . expanswn of r in terms of the sum of Irr(n-2m)...1_;(m)I factors. ~ Here n is the rank of the tensor, and m is the number of unit tensors, _1_. The definition of Irr(n-2m)...,_L.;(m)j is, for example, for n = 1, 2 andm:1,2 I—fi - r; — raafiy + r3507 +r7503 I—'_I - r5; - r0385“, + rarychb + rar¢581 + r Iii-I : 5a351¢ + 31.75041 + r8r¢5a7+ r1r¢5a8 50:15” I 5375a¢ and so on. Where 60:3 : Va r8 — . a, 8, 1, 4) each run on the cartesmn indices x, y, z. Due to permutations of indices, Irr(n-2m)..u(m)I has n n(n-l)....(n—2m+1) 02m (2m-l)!! = (2m-1)!! terms. (2m) ! The coefficient of each term is determined by the condition that I n I . . r must be traceless. For example, the coeff1c1ent for ~ m 2m (-1) r Irr(n-—2m)..11(mW is . ~~ -- (211-1) (2D‘3).....(2n-2m+1) 99 A formal expression then is: 2 ("“1 “ r————1 rn : r r r....r --——— rr(n-2)l + . ~l -2 .3 ~n w - (Zn-l) J. r [FFIU-4ILL] _ ............ + IZn-l)(2n-3) H (_1)m I.2m 'rr(n-2m)...11(6')" (8.2) (2n-l)(2n-3)....(2n-2m+1) ”7 -- . ' I . CC . . . . Since C. (r) 18 3 {+15 +2 rank tensor, for convenience, we let C+C'+2 : n, and thus: . l (2n-1)!! [—1 CC 9 (g) = ( ) I-l)n 2n+1 rn . (8.3) - 41:8 r " 2 . I n I. . . CC' . . . . Since r 18 an irreduc1ble tensor, G (r) is symmetric in all pairs ~ of its indices. For instance, for the 3-rank tensor, Q01“), there are the equalities: GMBIEI = 'GBaaII-‘I = GaBaIEI' Generally speaking, for the tensor of rank n, G (1:), there are 3n ( n 2 1 ) components. By using this symmetric point of view, we can reduce the 3n components to a smaller number. In the following table, we have reduced the first five n. Here, No is the reduced number. 100 3“.) 9 2781243’ Nc 3 6 10 15 21 Each reduced set of components must be accompanied by a suitable coefficient indicating the number of times it occurs in the its expression. This coefficient is determined by the number of permutations. For example, with n = 5, there are the equalities: GaaBnIEI 3 GaBanIEI = GnfldaIEI = --------- . The total number of permutations for indices as, B, and 17 is 30. Also, we can reduce the dimension of the (£+1)th rank tensor, 156(5), by using the same reduced numbers. The dimension of the (366 (r) matrix array not only depends on the number of components, it also depends on the number of inclusions in the primary cell. The use of the above described reduction method permits us to incorporate a relatively large number of multipole moments in the simulation. APPENDIX C C. LIST OF COMPUTER PROGRAMS This section presents the programs used for the calculation of the static and frequency-dependent effective dielectric constants in the simulations of Chapter 5. 101 102 Ct8*831*¥*¥*¥¥¥¥t**¥¥*¥¥¥*¥*¥¥t**¥**¥¥¥¥¥*¥**iiittttltittttt PROGRAM TITLE : STDIEL x PURPOSE : THE MULTIPOLE MOMENTS EFFECTS ON THE * STATIC DIELECTRIC CONSTANT OF COMPOSITE * MATERIALS FOR UNIFORM OR COATED INCLUSIONS.: THE CONTRIBUTIONS OF MULTIPOLES ARE UP TO x L=5. 1 METHOD : BY THE MINIMUM IMAGE METHOD 1 CONTAINS : ONE MAIN PROGRAM, THREE SUBROUTINES, AND * ONE CALLING MATH-LIBILINPACK OR FPS-164) * FOR SOLVING LINEAR EQUATIONS. t DATE : 1989 x xxxxxxxthtxxx:txxxxxtxxxxxxxttxtxitxxxtxtxxxxxxxthxxxxxxx DEFINITION : 3 DI : DIELECTRIC CONSTANT CF INCLUSION 1: D2 : DIELECTRIC CONSTANT OF MEDIUM 1 E1 : III-COMPONENT OF APPLIED FIELD 1: E2 : Y-COMPONENT OF APPLIED FIELD t E3 : Z-COMPONENT OF APPLIED FIELD t NA(MAX) AND NS : NUMBER OF INCLUSIONS x NCELL : INTEGER NUMBER FOR IMAGE CELL t NOON : NUMBER OF CONFIGURATIONS * NIVL : FINAL VOLUME FRACTION TO REACH * NPIMAX) AND NPOLE : DIMENSIONS OF MULTIPOLES * NPL : FINAL NPOLE TO REACH * NSINK : NUMBER OF INCLUSIONS * RADSINK : RADIUS OF INCLUSION x RADSYS : RADIUS OF SYSTEM * - VFRAC : VOLUME FRACTION t txtxtxxttxtxtxtxxttxxxxxxxxxxxxxxxtxxxxxxtxxtxtxxttxxxxxxxx INPUT : TININP,TININ2 ¥ OUTPUT : TINOUT * cxxxxxXtttxtxxxt*txt*xxttxxtxxxxxxxxxtxxxxxxxxxxxxxtxxtxxtxx OOOOOOOOOOOOOOOOO00000000000000 PROGRAM STDIEL IMPLICIT REAL*4 (A-H,O-Z) IMPLICIT INTEGER¥4 (I-N) PARAMETER (NA = 75, NP = 3) REAL PX(NA).PY(NA).PZ(NA) REAL OC(NP*NA),DP(NP*NAJ,TTWNP*NA,NP*NA),B(NP*NA) REAL WV(NP*NA+1),AINV(NP*NA,NP*NA) INTEGER IPVT(NA¥NP) (XXQIEI/ST/ TIIO) (Irtffil/SPV NCELL,PRD,PRD2,RADCUT2 OG‘TDN /MAIN/ A7,A9,A11,A13,A15,A17,A19 OPEN(1,FILE = ’TININP’) OPEN(2,FILE = ’TININZ') OPEN(7,FILE = ’TINOUT’) 35 103 PT : 4.03 ATANI1.0) P12 = 2.0*PI ANCUT=1.0 A7=105. A9=-A7¥9. A11=-A9t11. A13z-Alltl3. A15z-A13‘15. Al7=-A15*17. A19=-A17!l9. GET INFORMATION FROM CONFIGURATION DATA SET READ(1,*) NSINK, NCONG, VFRAC, RADSYS, RADSINK GET VARIABLES FROM INPUT DATA READIZ,” NSINK, NOON, RADSYS, NPL, NIVL READIZ,‘) D2, D1, El, E2, E3, NCELL SCAN N-POLE; FROM 3,9,19,34,55 DO 200 IPL = l, NPL IF(IPL.EQ.1) NPOLE:3 IF(IPL.EQ.2) NPOLE=9 IF(IPL.EQ.3) NPOLE=19 IF(IPL.EQ.4) NPOLE:34 IF(IPL.EQ.5) NPOLE=55 SCAN VOLUME FRACTION DO 100 IVOL = l, NIVL GET VOLUME FRACTION FROMCINPUT DATA READ(2,¥) VFRAC RADSINK=(3.0*VFRAC/(4.0¥PI¥NSINK))‘*(1.0/3.0) RADCUT = ANCUT t RADSYS RADCUT2 = RADCUT t RADCUT TDIST = RADCUT + RADSYS PRD : 2.08RADSYS PRD2 = PRD*PRD RADSYSZ = RADSYStRADSYS FACT = 4.03PItD2 CREATE T(MP) VECTOR 0102:01-02 MP=7 DO 30 L=0,MP-1 A=1. 88:1. DO 35 I=1,L BB=FLCAT(I)*BB CONTINUE IF(BB.EQ.0.) 88:1. 33 3O 20 888 DC 33 I:1,L 104 A=FLOATIZFI+1I*A CONTINUE T(L+1)= (D1D2)*(RADSINKt¥(2*L+3I)/(A*BB*((L+l)*Dl+(L+2)*D2)) CONTINUE WRITEI6,20) FORMAT(2X.//,4OX,’ RANDOM SINK CONFIGURATION ’,//) WRITE(6,t) WRITE(6,*) WRITEI6,*) WRITE(6,t) WRITE(6,#) WRITEI6,¥) WRITE(6,¥) WRITE(6,¥) WRITE(6,:) WRITE(6,x) WRITE(6,1) WRITE(6,3) U.” U‘ h-‘-- NUMBER OF SINKS = ’,NSINK RADIUS OF SINK = ’,RADSINK RADIUS OF SYSTEM 2 ’,RADSYS VOLUME FRACTION = ',VFRAC NUMBER OF CONFIGURATIONS : ’,NCON NPCLE=’,NPCLE DIELECTRIC CONSTANT OF SINK: D1 : ’,Dl DIELECTRIC CONSTANT OF BACKGROUND: D2 = ’,D2 T-OPERATOR FOR SINK: TZZ&t11&T22 = ',TZZ,T11,T22 X-COMPONENT OF APPLIED FIELD = ’,El Y-COMPONENT OF APPLIED FIELD : ',E2 Z-COMPONENT OF APPLIED FIELD = ',E3 SET INITIAL VALUE SSMPX : 0.0 SSMPXZ = 0.0 BEGINNING CONFIGURATION AVERAGE DO 1000 ICON = l, NOON GET NEW CONFIGURATION READ(1,¥) ICI NS = 0 DO 888 I = l, NSINK NS = N8 + 1 READ(1,*) PX(I), PYII), PZ(I) CONTINUE NDIM=NS*NPOLE CREATE B VECTOR (CONSTANT PARTS) DO 51 I = 1, NS IPl =(I-1)tNPOLE+ 1 IP2 :(I-1)tNPOLE+ 2 DO 52 LA = 1, NPOLE ID = (I-1)tNPOLE+ LA IF(ID. EQ. 1P1) THEN B(ID) = FACTxEl ELSE IF(ID.EQ.IP2) THEN B(ID) = FACTIEZ ELSE B(ID) = FACT3E3 ENDIF 00 333 1000 105 CONTINUE CONTINUE FOLLOWING CALL FOR CREATING THE MATRIX ELEMENTS CALL POLEALL(PX,PY,PZ,NS,NDIM,'IT,NPOLE) FOLLOWING CALL FOR SOLVING THE LINEAR EQUATIONS FROM FPS-164 MATH-LIB CALL PFINVINDIM,’I'T,WV,AINV,IERR) DO 81 I=l,NDIM SUM=0. DO 82 J=1,NDIM SLM:SLM+AINV(I,J)*B(J) CONTINLE DP(I):SUM CONTINUE FOLLOWING CALL FOR SOLVING THE LINEAR EQUATIONS FRCM LINPACK MATH-LIB CALL SGEFA('I'I‘,NDIM,NDIM,IPVT,IERR) CALL SGESL(I'T,NDIM,NDIM,IPVT,B,O) DO 81 I=l,NDIM DP(I)=B(I) CONTINUE THIS SECTION EVALUATES THE EFFECTIVE DIELECTRIC CONSTANT NR = 0 II) 333 I =1, NDIM, NPOLE NK : NK + 1 CC(NK) = DP(I)/D2 CONTINUE SMPX D015 1, ns SMPX - SMPX + CC(I) CONTINUE 0.0 I : IHH SSMPX=SSMPX+SMPX SSMPX2=SSMPX2+SMPX*SMPX' CONTINUE SSMPX = SSMPX/FLOATINCDN) IF(NCXDN.GT.1) THEN TEMP = ABS((SSMPX2 - FLOATING)N)*(SSMPX**2)I/ FLOATINCDNHNCDN-IHI SDSMPX = SQR'I‘ITFMP) 106 WRITE(6,501) SSMPX, SDSMPx ENDIF 501 FORMAT(2X,’ X-COMPONENT OF DIPOLE MOMENT = ' *.F14.8.' ( ’.F14-8.’ )'./) c P'TSH = SSMPX/I(4.D0¥PI/3.DO)*RADSYS**3) P'TSH:SS§’1PX/l.0 'IOPDE : 1.0 + 2.0XPTSH/(3.0¥E1) BOTDE = 1.0 - FISH/(3.0tE1) DEDZ = TOPDE/BOI'DE BO'TZ = BOTDEHZ SDDE : ABSI (SDSMPX/I (2.0*PI*El)*RADSYS¥*3)) t t (BO'I'DE+O.503’IOPDE)/BOT2) WRITE(6,66)DED2,SDDE c WRITE(7,43)VFRAC,DED2 66 FORMAT(2X,/,3X,’ EFFECTIVE DIELECTRIC CONSTANT: DE = ’, tF20.13,’ ( ’,F20.13,’ ) ’,/) 043 format(2x,2f15.8) 100 CONTINUE 200 CONTINUE STOP END C3ttt38333888838833881883883183833833$3888383333383883888888 c THIS SUBROUTINE CREATE THE MATRIX ELEMENTS x C¥titt83*3838888838881888388883381838888888¥8338838383833388 SUBRCUTINE POLEALL(PX,PY,PZ,NS,NDIM,TT,NN) IMPLICIT REALx4 (AFH,o-2) IMPLICIT INTEGERt4(I-N) DIMENSION AK(55),G(55,55) DIMENSION PX(NS),PY(NS),PZ(NS),TT(NDIM,NDIM) Demon /ST/ T110) Damon /sp/ NCELL,PRD,PRDZ,RADCUT2 m /MR/ RIlI) DATA (AK(I),I=1,55)/ 1.,1.,1.,-1.,-1.,-1. ,-2.,-2.,-2.,1.,1.,1.,3.,3.,3.,3.,3.,3. ,6.,-1.,-1.,-1.,-4.,-4.,-4.,-4.,-4.,-4. ,-6.,-6.,-6.,-12.,-12.,-12.,1.,1.,1.,5. ,5.,5.,5.,5.,5.,1o.,10.,1o.,1o.,1o.,1o. ,2o.,20.,20.,3o..30.,3o./ ”06*.“ C ROWING CALL FOR CREATING t COEFFICIENTS CALL CREATETINDlMgNNMS/I'T) 0061 =1,NS DO3J =1,NS 00100K=-NCELL,NCELL 00110L=-NCELL,NCELL m120Mz-NCELI..NCELL 21 13 16 130 120 110 100 107 DISTZ = (K*K+L*L+M*M)*PRD*PRD IF(I.EQ.J.AND.DIST2.EQ.0.DO) GO TO 130 XJ : PX(J)+K*PRD YJ = PY(J)+L*PRD ZJ : PZ(J)+M*PRD R2 : (PX(I)-XJ)xx2+(PY(I)-YJ)**2+(PZ(I)-ZJ)¥*2 IF(R2.EQ.0.) GO TO 130 IF(R2.GT.RADCU'I'2) GO TO 130 R1 = SQRT(R2) R(1)=R1 R(2)=R2 DO 21 II = 3,11 R(II) = 1.0/(R(1)**II) CONTINUE RX = (PX(I) - XJ)/R(1) RY : (PHI) - YJ)/R(1) R2 = (PZ(I) - ZJ)/R(1) DO 16 ID=1,NN DO 13 JD=1,NN KN=NN*(I-1)+ID JN=NN*(J-l)+JD IF(ID.LE.JD) THEN CALL ELEMENT(ID,JD,H,RX,RY,RZ) C(ID,JD)=H ENDIF IF(ID.GT.JD) THEN G(ID,JD)=G(JD,ID) ENDIF TI‘(IQJ,JN)=’IT(IQJ,JN)-ABS(AK(ID) )tG( ID,JD)*AK(JD) CONTINUE CONTINUE CONTINUE CONTINUE CONTINUE CONTINUE CONTINUE CONTINUE RETURN END Ctttitttttitttttttxtttittttittitttttttttttittitttttltxttlt C THIS SUBIDUTINE CREATE THE t mEFFICIENTS * C***3*¥ittttttitttittttttitttitttt¥#¥**¥‘*¥*t8¥t¥¥¥¥¥¥¥t¥t SUBROUTINE CREATET(NDIM,NN,NS{TT) IMPLICIT REAL*4(A2H,O-Z) IMPLICIT INTEGER*4(I-N) 108 DIMENSION IT(NDIM,NDIM) CG‘Q'DN /ST/ THO) I = 1,NS IS = 1,NN (I-1)*NN+IS J = 1,NS JS = 1,NN (J-1)*N'N+JS IF(ID.EQ.JD) THEN IF(IS.LE.3) IT(ID,JD) = 1./T(1) 888888 IF(IS.GT.3.A1\D.IS.LE.6) 'I'I‘(ID,.JD) = 1. /T(2) IF(IS.GT.6.AND.IS.LE.9) TT(ID,JD) = 2. /T(2) IF(IS.GT.9.AND. IS. LE. 12) T'NID, JD) = 1. /T(3) IF(IS.GT.12.AND. IS. LE. 18) 'I'I‘(ID, JD) : 3. /T(3) IF(IS.EQ.19) 'I'I‘(ID,JD)= 6. /T(3) IF(IS.GT.19.AND.IS.LE.22) 'I'I'(ID,JD) = l. /T(4) IF(IS.GT.22.AND.IS. LE. 28) 'I'I'(ID,.TD) : 4. /T(4) IF(IS.GT.28.AND.IS.LE.31) TT(ID,JD) = 6. /T(4) IF(IS.C}T.31.AND. IS.LE.34) T'I‘(ID,JD) : 12./T(4) IF(IS.GT.34.AND.IS. LE. 37) 'I'I‘(ID,JD) : 1. /T(5) IF(IS.GT.37.AND.IS. LE. 43) T'I‘(ID,JD) : 5. /T(5) IF(IS.GP.43.AM).IS.LE.49) I'I'(ID,JD) = 10./T(5) IF(IS.GT.49.AND.IS.LE.52) 'IT(ID,JD) = 20./T(5) IF(IS.GT.52.AND.IS. LE.55) '1'1‘_(ID,JD) = 30./T(5) ELSE TT(ID,JD):0. ENDIF 4 (DNTINUE 2 CDNTINUE 3 mNTINUE 1 CDNTINUE RETURN END Ctttiti333*13331**¥**¥**3¥¥t**!#titttitttttittitltttttttttt 0 THIS SUEmUTINE CREATE THE MATRIX W t C***i¥¥33$333"It*ttiiittittttittittittiitttx3t3¥3¥3¥*****i SUEROUTINE W(I,J,H,RX,RY,RZ) IMPLICIT REALM (A-H,O-Z) IMPLICIT INTEGER“ (I-N) comm Mil Rm) comm /MA.IN/ A7,A9,A11,A13,A15,A17,A19 KCX,KCY,KCZ REFER TO (XXJNT THE NUMBER OF X,Y,Z RX,RY,RZ ARE THE DISTANCE BETWEEN INCLUSIONS IN X,Y,Z DIRECTIONS COO KCX=0 KCY =0 KCZ=0 109 DC 1 II=1,2 IF(II.EQ.1) THEN =I ELSE N=J ENDIF IF(N.EQ.1) THEN KCX:KCX+1 KCY:KCY KCZ=KCZ ELSE IF(N.EQ.2) THEN KCX=KCX KCY=KCY+1 KCZ=KCZ ELSE IF(N.EQ.3) THEN KCX=KCX KCYzKCY KCZ:KCZ+I ELSE IF(N.EQ.4) THEN KCX:KCX+2 KCYzKCY KCZ=KCZ ELSE IF(N.EQ.5) THEN KCX=KCX KCY=KCY+2 KCZ=KCZ ELSE IF(N.EQ.6) THEN =KCX KCY=KCY KCZ=KCZ+2 ELSE IF(N.EQ.7) THEN KCX=KCX+1 KCY=KCY+1 KCZ=KCZ ELSE IF(N.EQ.8) THEN KCX=KCX+1 KCY=KCY KCZ=KCZ+1 ELSE IF(N.EQ.9) THEN KCX=KCX KCYzKCY+1 KCZ=KCZ+1 ELSE IF(N.EQ.10) THEN KCX=KCX+3 ' KCY=KCY KCZ=KCZ ELSE IF(N.EQ.11) THEN KCX=KCX KCY=KCY+3 KCZ=KCZ ELSE IF(N.EQ.12) KCX=KCX KCYzKCY KCZ:KCZ+3 ELSE IF(N.EQ.13) KCX=KCX+2 KCY:KCY+1 KCZ:KCZ ELSE IF(N.EQo14) KCX=KCX+2 KCYzKCY KCZ:KCZ+1 ELSE IF(N.EQ.15) KCX=KCX+1 KCY=KCY+2 KCZ:KCZ ELSE IF(N.EQ.16) KCX=KCX KCY=KCY+2 KCZ=KCZ+1 EISE IF(N.EQ.17) KCX=KCX+1 KCY=KCY KCZ=KCZ+Z ELSE IF(N.EQ.18) KCX=KCX KCY:KCY+1 KCZ=KCZ+2 ELSE IF(N.EQ.19) KCX=KCX+1 KCY=KCY+ 1 KCZ=KCZ+1 EISE IF(N.EQ.20) KCX=KCX+4 KCY=KCY KCZ=KCZ ELSE IF(N.EQ.21) KCXzKCX KCY=KCY+4 KCZ=KCZ ELSE IF(N.EQ.22) =KCX KCY=KCY KCZ=KCZ+4 ELSE IF(N.EQ.23) KCX:KCX+3 110 KCY=KCY+1 KCZzKCZ ELSE IF(N.EQ.24) KCX:KCX+3 KCY:KCY KCZ=KCZ+1 ELSE IF(N.EQ.25) KCX=KCX+1 KCY=KCY+3 KCZ:KCZ ELSE IF(N.EQ.26) KCX:KCX KCY :KCY +3 KCZ=KCZ+1 ElSE IF(N.EQ.27) KCX:KCX+ 1 KCYzKCY KCZ=KCZ+3 ELSE IF(N.EQ.28) KCX=KCX KCY=KCY+1 KCZ=KCZ+3 ELSE IF(N.EQ.29) KCX=KCX+2 KCY=KCY+2 KCZzKCZ ELSE IF(N.EQ.30) KCX=KCX+2 =KCY KCZ=KCZ+2 ELSE IF(N.EQ.31) chzKCX KCY:KCY+2 KCZ=KCZ+2 ELSE IF(N.EQ.32) KCX=KCX+2 KCY=KCY+1 KCZ:KCZ+ l ELSE IF(N.EQ.33) KCX=KCX+1 =KCY+2 KCZ:KCZ+1 ELSE IF(N.EQ.34) ' KCX=KCX+1 KCY=KCY+1 KCZ=KCZ+2 ELSE IF(N.EQ.35) KCX=KCX+5 =KCY 111 KCZ=KCZ ELSE IF(N.EQ.36) KCX:KCX KCY:KCY+5 KCZ=KCZ ELSE IF(N.EQ.37) KCX=KCX KCYzKCY KCZ=KCZ+5 ELSE IF(N.EQ.38) KCX=KCX+4 KCY=KCY+1 KCZ=KCZ ELSE IF(N.EQ.39) KCX=KCX+4 KCY=KCY KCZzKCZ+l ELSE IF(N.EIQ.40) KCX=KCX+1 KCY=KCY+4 KCZ=KCZ ELSE IF(N.EQ.41) KCX=KCX KCY=KCY+4 KCZ=KCZ+1 ELSE IF(NJ'HAZ) KCX=KCX+1 KCY=KCY KCZ=KCZ+4 ELSE IF(N.ED.43) KCX=KCX KCY=KCY+1 KCZ=KCZ+4 ELSE IF(N.m.44) =KCX+3 KCYzKCY+2 KCZ=KCZ ELSE IF(N.m.45) KCX=KCX+3 KCY=KCY KCZ:KCZ+2 ELSE IF(N.m.46) :KCX+2 KCY=KCY+3 KCZ=KCZ ELSE IF(N.EQ.47) KCX=KCX KCY:KCY+3 KCZ=KCZ+2 ELSE IF(N.ED.48) 112 113 KCX=KCX+2 KCY:KCY KCZ:KCZ+3 ELSE IF(N.EQ.49) THEN KCX=KCX KCY:KCY+2 KCZ=KCZ+3 ELSE IF(N.EQ.50) THEN KCX=KCX+3 KCY=KCY+1 KCZ:KCZ+1 ELSE IF(N.EQ.51) THEN KCX:KCX+1 KCYzKCY+3 KCZ:KCZ+1 ELSE IF(N.EQ.52) THEN KCX=KCX+1 KCYzKCY+1 KCZ=KCZ+3 ELSE IF(N.EQ.53) THEN KCX=KCX+2 KCY=KCY+2 KCZ=KCZ+1 ELSE IF(N.EQ.54) THEN KCX=KCX+2 KCY=KCY+1 KCZ=KCZ+2 ELSE IF(N.EQ.55) THEN KCX=KCX+1 KCY=KCY+2 KCZ=KCZ+2 ENDIF CONTINUE NT=KCX+KCY+KCZ KX=MAX(kCX,KCY,KCZ) KZ=MIN(KCX,KCY,KCZ) KY=NT-KX-KZ IF(KCX.EQ.KX.AND.KCY.EQ.KY.AND.KCZ.EQ.KZ) THEN =RX =RY Z=RZ ELSE IF(KCX.EQ.KY.AND.KCY.EQ.KX.AND.KCZ.EQ.KZ) THEN X=RY Y=RX Z=RZ ‘ ELSE IF(KCX.EQ.KZ.AND.KCY.EQ.KY.AND.KCZ.EQ.KX) THEN X=RZ =RY 114 Z=RX ELSE IF(KCX.EQ.KX.AND.KCY.EQ.KZ.AND.KCZ.EQ.KY) THEN X=RX Y=RZ Z=RY ELSE IF(KCX.EQ.KY.AND.KCY.EQ.KZ.AND.KCZ.EQ.KX) THEN X=RY Y=RZ Z=RX ELSE IF(KCX.EQ.KZ.AND.KCY.EQ.KX.AND.KCZ.EQ.KY) THEN XzRZ =RX Z=RY ENDIF IF(NT.EQ.2) GO TO 2 IF(NT.EQ.3) GO TO 3 IF(NT.EQ.4) GO TO 4 IF(NT.EQ.5) GO TO 5 IF(NT.EQ.6) GO TO 6 IF(NT.EQ.7) GO TO 7 IF(NT.EQ.8) GO TO 8 IF(NT.EQ.9) GO TO 9 IF(NT.EQ.10) GO TO 10 IF(KX.EQ.2) THEN H = R(3)*(3.0*X¥*2 - 1.0) ELSE IF(KX.EQ.1) THEN H = R(3)*(3.03X*Y) ENDIF GO TO 999 IF(KX.EQ.3) THEN H = -R(4)*(15.0*X*‘3 -9.0*X) ELSE IF(KX.EQ.2) THEN H = -R(4)*(15.0tY¥X**2-3.0¥Y) ELSE IF(KX.EQ.1) THEN H = -R(4)¥(15.0*X*Y*Z) ENDIF . GO TO 999 IF(KX.EQ.4) THEN H = —R(5)t(90.08Xtt2-105108Xtt4-9.0) ELSE IF(KX.EQ.3) THEN H = -R(5)*(45.o*x*Y-105.*thtta) ELSE IF(KX.EQ.2.AND.KY.EQ.2) THEN H = -R(5)*(15.0¥(X**2+Y**2)-105.0#(X*Y)**2-3.0) ELSE IF(KX.EQ.2.AND.KY.EQ.1) THEN H = -R(5)*(15.0*Y*Z-105.0¥Y*Z*X**2) ENDIF X t 6 t t t t t t t 115 GO TO 999 IF(KX.EQ.5) THEN H: ELSE ELSE ELSE H: R(6)¥(-225.¥X+ 1050.*X**3 —945.¥X¥*5) IF(KX.EQ.4) THEN R(6)*(-45.tY+630.*(X**2)*Y-945.*(X**4)*Y) IF(KX.EQ.3.AND.KY.EQ.2) THEN R(6)*(-45.*X+315.*X*Y*¥2 -945.¥(X¥¥3)*(Y*¥2)+105.¥X¥*3) IF(KX.EQ.3.AND.KY.EQ.1) THEN R(6)*(315.*X*Y*Z-945.¥(X*‘3)*Y*Z) IF(KX.EQ.2.AND.KY.EQ.2) THEN R(6)*(-15.*Z+105.‘(X**2)*Z+105.*(Y**2)*Z -945.*(X**2)‘(Y‘*2)*Z) ENDIF GO TO 999 IF(KX.EQ.6) THEN H = R(7)*(-225.+ 4725.*Xt*2 -14175.*X**4+ 10395.tX**6) ELSE IF(KX.EQ.5) THEN H = R(7)‘(1575.*X*Y -9450.¥X*¥3¥Y+ 10395.*X*¥5#Y) ELSE IF(KX.EQ.4.AND.KY.EQ.2) THEN H = R(7)t(-45.-5670.¥X**2‘Y‘*2+10395.*X*¥4tY*‘2 + 630.*X‘¥2 -945.¥X¥*4+ 315.*Y**2) ELSE IF(KX.EQ.4.AND.KY.EQ.1) THEN H = R(7)*(315.*Y*Z -5670.‘Xt*2‘YtZ+ 10395.‘X*¥4*Y*Z) ELSE IF(KX.EQ.3.AND.KY.EQ.3) THEN H = R(7)*(945.‘X‘Y -2835.‘X¥Y¥‘3-2835.*X¥t3*Y + 10395.¥X**3¥Yt33) ELSE IF(KX.EQ.3.AND.KY.EQ.2) THEN H = R(7)*(315.*XtZ-945.*X*33*Z-2835.¥ XtY‘*2‘Z+ 10395.¥X**3*Y**2*Z) ELSE IF(KX.EQ.2.AND.KY.EQ.2) THEN H = R(7)*(-15. -945.¥Xt¥2¥Y*t2 -945.1X332*Z¥t2 -945.*Y*¥2*ZX*2+ 10395.*X**2*Y**2*Z**2 + 105.*X**2+ 105.¥Y**2+ 105.*Z**2) ENDIF GO TO 999 IF(KX.EQ.7) THEN H = R(8)¥(A13*X¥*7+A11321.*X**5+A9*105.*Xt33+A7*105.*X) ELSE IF(KX.EQ.6) THEN H = R(8)*(A13*X¥‘6*Y+A11*15.*X*¥4*Y+A9$45.*X**2*Y ELSE H: +A7¥15.¥Y) IF(KX.EQ.5.AND.KY.EQ.2) THEN R(8)¥(A13*X¥*53Y¥*2+A11*(10.*X*¥3*Y*t2+xt¥5)+ A9*(15.¥X¥Y**2+10.*X*¥3)+A7*15.*X) 116 ELSE IF(KX.EQ.5.AND.KY.EQ.1) THEN H : R(8)*(A13*X**5*Y*Z+A11*10.*X**3*Y*Z+A9*15.¥XtY*Z) ELSE IF(KX.EQ.4.AND.KY.EQ.3) THEN H = R(8)t(A13tX**4*Ytt3+A11t(6.*X**2*Y**3+3.*X**4¥Y)+ * A9*(3.*Yt*3+18.*Xttth)+A7*9.*Y) ELSE IF(H.EQ.4.AND.KY.EQ.2) THEN H = R(8)t(A13*X**4*Y**2tZ+A11*(6.tx**2*Y**2*Z+x**4*Z)+ x A9*(3.*Yt*2*Z+6.*Xtt2*Z)+A7*3.*Z) ELSE IF(KX.EQ.3.AND.KY.BQ.3) THEN H = R(8)¥(A13¥X¥*3*Y**3¥Z+A11*(3.tX¥Y*¥3*Z+3.*X**3*Y*Z)+ * A9*9.*XtY*Z) ELSE IF(KX.EQ.3.AND.KY.EQ.2) THEN H : R(8)3(A13¥X¥*3¥Y*¥2t2*¥2+A11t(3.*X¥Y**2*Z*¥2+ : thatth2+th3th32)+A9t(3.tXtZt32+3.EX¥Y¥*2+ x xxt3)+A7t3.tx) . ENDIF GO TO 999 8 IF(KX.EQ.8) THEN H = R(9)*(A15*X**8+A13*28.¥X*¥6+A11*210.*Xtt4+ t A93420.*X**2+A7*A7) ELSE IF(KX.EQ.7) THEN H = R(9)t(A15*X*¥7tY+A13t21.*X**5*Y+A11*105.*X¥*3*Y+ t A93105.tX¥Y) ELSE IF(KX.EQ.6.AND.KY.EQ.2Y THEN H = R(9)3(A15tX336tht2+A133(15.*X**41Y**2+X*16)+ * A11¥(45.xxtt2*Y*12+15.¥X**4)+A9*(Y**2+45.tXtt2)+ * A7t15.) ELSE IF(KX.EQ.6.AND.KY.EQ.1) THEN H = R(9)‘(A15¥X¥*6*Y¥Z+A13*15.txtt4tY¥Z+Allt45.xxxx23Ytz+ * A9*15.thZ) ELSE IF(KX.EQ.5.AND.KY.EQ.3) THEN H = R(9)¥(A15*X**5tY333+A13¥(10.¥X**3*Y*¥3+3.tx*t5tY)+ * A11*(15.tXtht3+30.tx**3*Y)+A9¥45.*X*Y) ELSE IF(KX.EQ.5.AND.KY.EQ.2) THEN H = 8(9)¥(A15*X**5*Y**2*Z+A13t(10.¥X**3*Yt*2*Z+X¥*5*Z)+ * A11t(15.xX*Y*¥2*Z+10.*X¥*3*Z)+A9tl5.*X*Z) ELSE IF(KX.EQ.4.AND.KY.EQ.3) THEN H = 8(9):(A15xX*t4*Y**3tZ+A13X(6.¥X¥*21Y3t3*z+3.*X**4*Y*Z) * +A11t(3*Y*¥3*Z+18.*X**2¥Y*Z)+A9*9.¥Y*Z) ELSE IF(KX.EQ.4.AND.KY.EQ.2) THEN H = R(9)i(Alstxtt4thtthtt2+A13t(6.*X**2*Y**2¥Z*t2+ t xxx4tzxt2+xtt4thxz)+A11t(3.thxzxztxz+ t 6.xxtxzxztt2+6.*XthxYxtz+xtx4)+A9*(3.xzxxz+ t 3. ¥Yt¥2+6. *X#¥Z)+A7*3.) ELSE IF(KX.EQ.4.AND.KY.EQ.4) THEN H = R(9)*(A15tXt*4¥Y**4+A13*(6.*X**2*Y*t4+6.¥X¥*4*Y**2)+ * A11*(3.*Y**4+3.txtt4+36.*th2*Y**2)+ * A9¥(18.*Y**2+18.*X**2)+A7*9.) ELSE IF(KX.EQ.3.AND.KY.EQ.3) THEN t t t 9 t H : 117 R(9)t(Alstxxxsxyxx3t2tx2+A13x(3.xxxYxx3xzxt2 +3.:x:x3xszxt2 +th3tht3)+A11x(9.xx*Yx2tx2+3.xXthx3+3.*xxx3xY)+ A9*9.*X*Y) ENDIF GO TO 999 IF(KX.EQ.9) THEN H : ELSE R(10)¥(A17*Xt*9+A15*36.¥X**7+A13*378.*X**5+ A11¥1260.¥X*‘3+A9*9.*A7*X) IF(KX.EQ.8) THEN R(10)¥(A17*X**8tY+A15*28.*X**6*Y+A13*210.*X¥*4*Y+ A11*420.*X*¥2¥Y+A9*A7*Y) IF(KX.EQ.7.AND.KY.EQ.2) THEN R(10)*(A17tX**7*Y**2+A15*(21.*X‘*5*Y**2+X**7)+ A13*(105.*X**3*Y*¥2+21.*X*¥5)+ A11*(105.*X*Y*¥2+105.‘X**3)+A9*105.*X) IF(KX.EQ.7.AND.KY.EQ.1) THEN R(10)‘(A17*X*¥7*Y‘Z+A15*21.*X**5*Y*Z+ u A13*105.*X**3¥Y*Z+A11*105.*X*Y¥Z) IF(KX.EQ.6.AND.KY.EQ.3) THEN R(10)*(A17¥X**6*Y*33+A15*(15.*X**4*Yt*3+3.*X**6*Y)+ A133(45.*X**2*Y*¥3+45.*X**4*Y)+A11¥(15.*Y**3+ 135.3Xt*2*Y)+A9*45.*Y) IF(KX.EQ.6.AND.KY.EQ.2) THEN R(10)$(A17*X**6‘Y‘*2*Z+A15t(15.¥X¥*4*Y**2*Z+X**6*Z)+ A13¥(45.*X**2*Y**2*Z+15.*X**4*Z)+ A11*(15.*YI*Z*Z+45.*X**2*Z)+A9*15.*Z) IF(KX.EQ.5.AND.KY.EQ.3) THEN R(10)*(A17¥X‘t53Y**3*Z+A15*(10.*X**3*Y3*3*Z+3.*Xtt53Y32) +A13*(15.*XtY3‘3¥Z+30.‘X**3*Y*Z)+A11*45.*X*Y‘Z) IF(KX.EQ.5.AND.KY.EQ.2) THEN R(10)‘(A17‘X*¥5‘Y**2*Z*HZ+A15t(10.*X**3*Y**2*Z**2+ X33532132+Xt15¥Y¥32)+A13*(15.¥X*Y*¥2XZ**2+ 10.tx*‘332**2+10.¥X**3iY*32+X**5)+ A11*(15.*X*Z*32+15.*XiYi*2+10.*X¥¥3)+A9*15.*X) IF(KX.EQ.5.AND.KY.EQ.4) THEN R(10)*(Al7*X**5‘Y**4+A15*(10.*X**3*Y‘*4+6.*X**5*Y*i2)+ A133(15.*X¥Y**4+3.¥X335+60.*X!t3*Y¥*2)+ A11*(90.*X*Y**2+30.*X¥*3)+A9*45.*X) IF(KX.EQ.4.AND.KY.EQ.3) THEN R(10)¥(A17¥X3*4$Y*#3*Z*#2+ A15*(6.*X¥¥2*Y**3*Z**2+3.*X**4*Y*Z**2+ X134¥Y*¥3)+A13*(3.*Y**3323*2+18.¥X**2*Y*Z**2+ 6.*X¥*2*Y**3+3.*X**4¥Y)+A113(9.*Y*Zt*2+ 3.‘Y**3+18.¥X**2¥Y)+A9*9.*Y) IF(KX.EQ.4.AND.KY.EQ.4) THEN R(10)*(A173X334¥Y¥¥4¥Z+A15¥(6.*X*¥2*Y**4*Z+ 6.¥X*¥4*Y**2¥Z)+A13*(3.¥Y**4*Z+3.*X**4*Z+ 36.*X**2*Y**2*Z)+A11*(18.¥Y**2*Z+18.*X**2*Z)+ 10 t t l * ELSE H: 118 A9¥9.‘Z) IF(KX.EQ.3.AND.KY.EQ.3) THEN R(10)*(A17*X**3*Y**3*Z**3+A15*(3.*X*Y**3*Z**3+ 3.*X**3*Yt*3¥2+3.*X‘*3*Y*Z**3)+A13*(9.¥X*Y*Z**3 +9.¥X*Z*Y**3+9.*X**3*Y*Z)+A11*27.*X*Y*Z) ENDIF GO TO 999 IF(KX.EQ.10) THEN H: ElSE R(11)*(A19*X**10+A17*45.*Xtt8+A15*630.*X**6+ A13t3150.*xtt4+A11t45.¥A7*X**2-A9*A9) IF(KX.EQ.9) THEN R(11)*(A19*X**9¥Y+A17*36.*X**7*Y+A15*378.*X**5*Y+ A13*1260.*X**3*Y+A11*9.*A7*X*Y) IF(KX.EQ.8.AND.KY.EQ.2)ITHEN R(11)*(A19¥X*¥8*Y#*Z+A17*(28.¥X**6*Y**2+X**8)+ A15*(210.*X**4*Y**2+28.*X**6)+ A13*(420.tX¥t2¥Yttz+210.*Xtt4)+ A11¥(A7¥Yt¥2+420.tX**2)+A9¥A7) IF(KX.EQ.8.AND.KY.EQ.1) THEN R(11)t(A19txtt8tY*Z+A17*28.¥X**6*Y*Z+ A151210.¥X**4*Y*Z+A13*420.txx*2*Y*Z+All*A7*Y*Z) IF(KX.EQ.7.AND.KY.BQ.3) THEN R(11)¥(A19txt*7tY¥*3+A17*(21.*X**5*Y**3+3.tXtt7tY)+ A153(105.*Xt*3¥Y**3+63.*X¥*5*Y)+A13*(105.*X*Y**3+ 315.*X**3¥Y)+A11*315.*XtY) IF(KX.EQ.7.AND.KY.EQ.2) THEN R(11)t(A19*X*¥7tY*32*Z+A17*(21.*X¥*5*Yt*2*Z+X**7*Z)+ A15*(105.*Xt*3*Y¥*2*Z+21.*X*¥5*Z)+ A13¥(105.:xtht232+105.*X**3*Z)+A11*105.*X*Z) IF(KX.EQ.6.AND.KY.EQ.4) THEN R(11)*(A19*X¥*6*Y$*4+A17*(15.*X**4*Y**4+6.*X**6*Y**2) +A15¥(45.tx**2*Y**4+3.*X**6+90*X*t4*Y**2)+ A13*(15.tY**4+45.*X**4+270.*X**2*Y*t2)+ A11*(90.tht2+135.*X¥*2)+A9*45.) IF(KX.EQ.6.AND.KY.EQ.3) THEN R(11)t(A193X*:6tY*¥3*Z+A17*(15.*X**4tY**3*Z+ 3.*Xt*6¥Y*Z)+A153(45.*X**2*Y**3*Z +45.*X¥*4¥Y*Z)+A13*(15.*Y**3*Z+135.*X*t2*Y*Z)+ A11t45.*Y*Z) IF(EX.EQ.6.AND.KY.EQ.2) THEN R(11)t(A19¥th63Y**2tZ**2+A17*(15.*X**4*Y**2¥Z**2+ xxxexzxx2+xtxstxxz)+A158(45.xxttzxvxtztzxxz+ 15.*Xtt4tzt*2+15.*xtt4th*2+X**6)+ A13¥(45.*X**2¥Z**2+45.*Xt82*Y**2+15.*X**4+ 15.xyxx2:zx:2)+A11x(15.xzxx2+15.xyxtz+45.xxxx2) +A9¥15.) IF(KX.EQ.5.AND.KY.EQ.5) THEN R(11)*(A19tX**5*Y**5+A17*(10.*X*¥3*Y*15+10.*xtt5*Y**3)+ A15*(15.*Xth*5+15.*x**5*Y+100.*X**3*Y*¥3)+ 999 *fififl *N‘I‘I’N 119 A13¥(150.*X*Y**3+150.*X**3*Y)+A11*225.*X*Y) ELSE IF(KX.EQ.5.AND.KY.EQ.3) THEN H = R(11)*(A19*X**5*Y*‘3*Z$*2+ A17*(10.*X¥*3*Y**3*Z**2+3.*X**5*Y*Z**2+ X**5*Y**3)+A15#(15.*X*Y**3*Z**2+30.*X¥:3*Y*Zt¥2+ 10.¥X**3*Y**3+3.*X**5*Y)+A13*(45.¥X¥Y*Z¥*2+ 15.*X‘Yt‘3+30.*X**3*Y)+A11*45.*X*Y) ELSE IF(KX.EQ.5.AND.KY.EQ.4) THEN H = R(11)*(A19*X**5*Y**4*Z+A17*(10.*X**3*Y**4*Z+ 6.*X¥*5*Yt*2*2)+A15*(15.*X*Y**4*Z+3.*X‘*5*Z+ 60.¥X**3*Y**2*Z)+A13*(90.*X*Y**2*Z+30.¥X**3*Z)+ A11*45.‘X¥Z) ELSE IF(KX.EQ.4.AND.KY.EQ.3) THEN H = R(11)3(A19*X**4*Y**3*Z**3+A17‘(6.*X**2*Y**3*Z**3+ 3.‘X**4iY*Z**3+3.tX*¥4*Y**3*Z)+A15*(3.*Y**3 tZ¥¥3+18. tx3t2¥YtZ¥¥3+18. ¥X¥¥2*Y¥*3*Z+ 9.*X**4*Y*Z)+A13¥(9.*Y*Z**3+9.*Y**3*Z+ 54.tX¥*2¥Y‘Z)+A11*27.‘Y‘Z) ELSE IF(KX.EQ.4.AND.KY.EQ.4) THEN H = R(11)*(A19*X‘*4tY**4¥Z*‘2+A17*(6.¥X*¥2*Y**4*Z**2 +6 . ¥X¥¥4¥Yt3232¥¥2+X1343Y3#4 ) +A15t ( 3 . tY‘X4tZit2'} 3.*X**4*Z‘*2+6.*X**2*Y“4+6.*X**4*Y**2+36.t X¥¥23Y¥¥212332 ) +A13*( 18 . ‘Y*¥2¥Z*¥2+3 . *Y**4+ 3.*X**4+36.*X*‘2‘Y**2+18.‘X**2*Z*32)+ A11*(9.*Z¥*2+18.*Y*‘2+18.¥X**2)+A9¥9.) ENDIF CONTINUE RETURN 120 Ctttti¥*******¥**¥****t*********¥**¥***¥*¥****¥**¥XX*¥***fit! 0000000000OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO PW TITLE : FmDIEL PURPOSE : SIMULATION OF THE FREQUENCY -DEPENDENT DIELECTRIC CONSTANT OF OCMPOSITE MATERIALS FOR UNIFORM OR COATED INCLUSION. METHOD : CHOOSE ONE (1) THE MINIMUM IMAGE METHOD OR (2) THE LATTICE SUM METHOD CONTAINS : ONE MAIN PROGRAM ,THREE SUBROUTINES, AND ONE CALLING MATH-LIB FOR SOLVING THE COMPLEX LINEAR EQUATIONS DATE : 1989 **************¥¥*****¥*##3##ttitt1**********t**¥tt*****¥t* DEFINITION : ANCUT = RADCUT/RADSYS APHA : CONVERGE PARAMETER BF : BEGINNING FREQUENCY: 0.2D0 CD : THE CONDUCTIVITY D1 : DIELECTRIC CONSTANT OF PARTICLE D2 : DIELECTRIC CONSTANT OF MEDIUM DE : THE EFFECTIVE DIELECTRIC CONSTANT DREl : RE(BULK DIELECTRIC CONSTANT FOR INCLUSION) DIMI : IM(BULK DIELECTRIC CONSTANT FOR INCLUSION) E1 : X-CO‘IPONENT OF APPLIED FIELD E2 : Y-COMPONENT OF APPLIED FIELD E3 : Z-COMPONENT OF APPLIED FIELD EF : ENDING FRmUENCY: 0.55D0 II(MAX) AND NUMW : THE NUMBER OF FREQUENCY POINTS N2(MAX) : THE NUMBER OF INCLUSIONS IN PRIMARY CELL NCELL : INTEGER IN k-SPACE NCONG : NUMBER OF CONFIGURATIONS NN : THE DIMENSIONS OF DIPOLES; NN=3 FOR L=1 NSINK : NUMBER OF INCLUSIONS MCELL : INTEGER IN r-SPACE PRD : UNIT LENGTH OF CELL RADSYS : RADIUS OF SYSTEM RADSINK : RADIUS OF INCLUSION TAOINV : DAMPING CONSTANT VFRAC : VOLUME FRACTION VF : ELECTRON FERMI VELOCITY W(I) : FREQUENCY WP : PLASMA EREQUENCY WS : RESONANCE FREQUENCY INPUT : TININP,TININ2 OUTPUT : TINQJT , TININ2 *‘I’itfl’tfiflfifififl’fififififlflfi'bt’tfi******************** 3*3****¥*#*¥**¥*¥**¥*¥*****¥**¥**¥***¥¥****tttttttitttttttt t * Cit****¥***¥**¥*¥i¥t*¥*¥¥¥t¥*********¥***t**tt*¥¥***¥*3*¥t¥¥ PROGRAM FRQDIEL IMPLICIT REAL*8 (A-H,O—Z) IMPLICIT INTEGER34 (I-N) 20 121 121 PARAMETER (N2 = 25,II : 100,NN : 3) REAL PX(N2),PY(N2),PZ(N2),XC(N2),YC(NZ),ZC(N2),R(N2) REAL W(II),SSREAL(II),SSIMAG(II),SSREZ(II),SSIM2(II) REAL SSREIM(II),SMPXRE(II),SMPXIM(II) COMPLEX*16 TMAT(NN*N2,NN*N2),B(NN*N2),C(NN*N2),CC(NN*N2) COMPLEX*16'THHAT(N2¥NN,N2*NN){HMT(NN*N2,NN*N2),SSMPX(II) COMPLEX*16 SMPX(II),Dl,DIN1,DDINO,SSMPXT,TOPDE,BOTDE CQPLEX*16 DE,TZZ,T(10),U,V,UV,D1D2,'I‘OPI‘,BOTT,V1 INTEGER IP'I'I‘(N'N¥N2) WN /PP/ NCEL,."CELL,PRD,APHA,RADCUT2,PI,PI2 OPEN(1,FILE = ’TININP’) OPEN(2,FILE = ’TININZ’) OPEN(7,FILE = ’TINOUI") OPEN(8,FILE = ’TINOUl’) PI = 4.DO* DATAN(1.DO) PIZ = 2.DO*PI ANCUT = 1.DO (DNPIGURATION DATA INPUT m TININP READ(1,¥)NSIM{,NOONG,VFRAC,RADSYS,RADSINK GET VARIABLES FRO! TININZ READ(2,*)NSINK,N(I)N,VFRAC,RADSYS READ(2,¥)D2,m,E1,E2,E3,NLW,NCELL,PCELL,APHA,BF,EF RADSINK : (3.DO*VFRAC/(4.D0*PI*NSINK))tx(1.D0/3.DO) RAD=(3.¥PHI/(4.*PI¥NSINK)NHL/3.) AR=RADSINIURAD RAIXZ‘UT : ANCUT t RAIBYS RAIX'JUTZ : RAWI‘ t RADCUT PRD = 2.D0*RADSYS FACT = 4.DO*PIilh WRITE(6,20) EDWANZXJ/AOX,’ RANIXM SINK CDNFIGURATION ' ,//) WRITE(6,¥) ’ NUMBER OF SINKS " ',NSINK WRITE(6,*) ’ RADIUS OF SINK = ' ,RADSINK WRITE(6,*) ' RADIUS OF SINK = ’,RAD WRITE(6,*) ’ RADIUS OF SYSTEM = ’,RADSYS WRITE(6,*) ’ VOILME FRACTION = ’,VFRAC WRITE(6,t) ’ VOLUME FRACTION = ',PHI WRITE(6,*) ' NUMBER OF WEIGURATIONS = ’,NmN WRITE(6,121) WT(2X,//,4OX,’ PARAMETERS OF THE SYSTEM ’,//) WRITE(6,H ’ DIELECTRIC mNSTANT OF MEDIUM:D2 = ’,DZ WRITE(6,*) ’ X-CG'IPONENT OF APPLIED FIELD = ’,El 49 000 888 U] 0| 0| 122 WRITE(6,¥) ’ Y-COMPONENT OF APPLIED FIELD WRITE(6,¥) ’ Z-COMPONENT OF APPLIED FIELD ’,E2 ’,E3 SET INITIAL VALUES FOR EACH FREQUENCY POINT DO 49 I = l, NUMW SSMPX(I) = CMPLX(0.D0,0.D0) SSRE2(I) : 0.00 SSIM2(I) = 0.D0 SSREIM(I) : O.D0 CONTINUE FOLLOWING IS CONFIGURATION AVERAGE DO 1000 ICON : 1, NOON GET NEW CONFIGURATION; RANDOMLY GENERATE THE CENTER OF NON-OVERLAPPING INCLUSIONS, XC,YC,ZC BETWEEN 0 AND 1. PERIODIC BOUNDARY CONDITION USED READ(1,*) ICI DO 888 I : 1,NSINK READ(1,¥) XC(I),YC(I),ZC(I) CONTINUE CHECK PARTICLES ONLY INSIDE THE SPHERE(RADCUT) OR PRIMARY CELL WRITE(6,¥)’ICON=',ICON PRDH = PRD/2.00 NS = 0 _ DO 555 M = 1, NSINK IF( ABS(XC(M)-RADSYS) .GT. PRDH ) GO TO 555 IF( ABS(YC(M)-RADSYS) .GT. PRDH ) GO TO 555 IF( ABS(ZC(M)-RADSYS) .GT. PRDH ) GO TO 555 NS = NS + 1 PX(NS) XC(M) PY(NS) YC(M) PZ(NS) ZC(M) CONTINUE NDIM IS TOTAL DIMENSION OF MATRIX NDIM : NNSNS CONSTRUCTION OF B VECTOR(CONSTANT PARTS) DO 51 I = 1, NS IPl = (I-1)*NN+1 1P2 = (I-1)tNN+2 DO 52 IA = 1, NN ID = (I-1)*NN+IA IF(ID .EQ. 1P1) THEN B(ID) = FACT t E1 ELSE IF(ID .EQ. IP2) THEN 00000000 35 123 B(ID) FACT ¥ E2 ELSE B(ID) END IF CONTINUE CONTINUE FACT # E3 MIMATRIX AND LSMATRIX ONLY DEPEND UPON THE POSITIONS OF INCLUSIONS. IF USE LATTICE SUM METHOD, CALL ISMATRIX CALL LSMATRIX(TMAT,PX,PY,PZ,NN,NS,NDIM) IF USE MINIMUM IMAGE METHOD, CALL MIMATRIX CALL MIMATRIX(TMAT,PX,PY,PZ,NN,NS,NDIM) SCAN FREQUENCY DO 50 IN = 1, NUMW FWP : FLOAT(IN)*(EF-BF)/FDOAT(NUMW)+BF WP=1.46D+16 GAMAI=O.24E+14 TAOIB=0.87E+14 VF:1.44E+14 GAMAAzGAMAI+VF/(RADSINK*(10.E-5)) TAO=GAMAA+TAOIB WS=WP/SG¢T(DRE1+2.3D2) WSWP=WS**3/(WP**2) DIM1=(TAO-GAMAA)/WSWP TAOINV: 1 . 68D+ 1 4 DRE1:4.5DO DIMI=0.16DO W(IN) = FWP*WP WWP : W(IN)/WP TAOWP =TAOINV/WP DINl = CMPLX(DRE1,DIM1) DDINO = CMPLX(WWP,TAOWP) D1 = DINl - I.D0/(WWP*DDINO) =D2/Dm V:Dl/Dm V1=V—1. =UPV D1D2=D1~D2 NP=2 IX) 30 L=O,NP-1 A=1. 83:1. DO 35 I=1,L BB=FLOAT(I+1)*BB CONTINUE IF(BB.EQ.0.) BB=1. D0 33 I=1,L 33 30 47 333 151 1000 124 A=FDOAT(2*I+1)*A CONTINUE L1:L+1 L2:L+2 L23=23L+3 TOPT: L1*V1*(V*L2+U*L1)*RAD**L23+L1*UV*(L1+V*L2)*RADSINK**L23 BOTT ((V*L2+U*L1)‘(L2+V*L1)+L1¥L2*UV*V1*AR**L23)*A*BB T(L1)=TOPT/K7I'T CONTINUE TZZ=T(1) MATRIX ADDS DIAGONAL PARTS(SELF TERM IN PRIMARY CELL) ELECTRIC DO 47 I=l,NDIM IX) 47 J:1,NDIM IF(I.EQ.J) TTMAT(I,J)=TMAT(I,J)+1.DO/TZZ IF(I.NE.J)'PHHAT(I,J)=TMAT(I,J) CONTINUE SOLVING LINEAR EQUATION (TTMATXC:B ) CALL SOLVEC(TTMAT,C,B,NN,NDIM,NS,TZZ,IPTTfHflT) ONLY PICK DIPOLE MOMENTS AT X-OOMPONENT NK = 0 DO 333 I = 1, NDIM, NN NK = NK + I COME) = C(I)/Dm CONTINUE SUM ALL DIPOLE MOMENTS SMPX(IN) = CMPLX(O.D0,0.D0) DO 151 I = 1, NS SMPX(IN) = SMPX(IN) + CC(I) CONTINUE SMPXRE(IN) = REAL(SMPX(IN)) SMPXIM(IN) = AIMAG(SMPX(IN)) SSMPX(IN) : SSMPX(IN) + SMPX(IN) SSREZ(IN) = SSRE2(IN) + SMPXRE(IN)*SMPXRE(IN) SSIMB(IN) = SSIM2(IN) + SMPXIM(IN)*SMPXIM(IN) SSREIM(IN) = SSREIM(IN) + SMPXRE(IN)*SMPXIM(IN) CONTINUE CONTINUE RADNE3 (3.D0/(4.DOiPI))*(2.D0*RADSYS)¥t3 RADNEW RADNE3**(1.DO/3.D0) 43 60 125 AL = 1.DO/(4.DO*PI*E1*RADNEW**3) ALZ = AL¥AL ANALYZE FREQUENCY IX) 60 I = 1, NUMW SSMPXT : SSMPX(I)/N(X)N TOPDE = 1.D0 + 2.D0*AL*SSMPXT BOTDE = 1.D0 - AL*SSMPXT DE DmHTOPDE/WTDE) CD W(I)*AIMAG(DE)/(4.DO*PI) SECTION FOR ERROR ANALYSIS IF(NOON.GT.1) THEN SRET = REAL(SSMPXT) SIMI‘ = AIMANSSMPXT) SREIMI‘ = SREI'XSIMI‘ SIGA ABS(SSRE2(I) - N®N*SRET*SRET)/(NCON-1) SIGB ABS(SSIM2(I) - NmN*SIMT*SIMI')/(NCON-1) SIGAB = ABS(SSREIMU) - NCONIISRETIRSIMU/(NmN-l) TOP = 1.DO+ AL*SRET- 2.D0*AL2HSRE'T*SRE'T+ SIMHSIMT) BOT = 1.DO- 2.DO*AL*SRET+ ALZHSRE’NSREN SIMI‘tSIMl‘) 8012 = 30er DlDA (Bort(AL-4.DoxAszSREr)-mpt(4.DOtAL2xSREr t -2.DO*AL))/Kfl‘2 ' DlDB = (BOT*(-4.DO*ALZtSIMI‘)- TOP¥(2.DO*AL2*SIMT)) * IBO'I‘2 TZOP : 3.DO¥AL*SIMI‘ DZDA = T20P*(-2.DO*ALZtSRET + 2.*AL)/BOI‘2 DZDB = (Borx(3.DoxAL) - T20P*(2.xAL2*SIMl‘))/BUI’2 SD11 = SIGAtDlDAtDIDA SD12 = SIGBxDIDBxDIDB 8013 = SIGABtDlDAxDlDB SD21 = SIGAtDZDAtDZDA SD22 = SIGBxDZDBxDZDB 3023 = SIGAB*D2DA*D2DB SIGDl : SGZTHSDII + SD12 + SD13)/NCDN) SICEZ = SQQTHSDZI + SD22 + SD23)/NCX)N) ENDIF PRINT OUT DATA WRITE(7,43) W(I)/WP,REAL(DE),SIGD1 WRI'TE(8,43) W(I)/WP.AIMAG(DE),SIGD2 MAT(ZX,f10.6,1X,F10.6,1X,F10.6) CDNTINUE STOP END C133*!*¥******¥***Sittlttttiittittit333*31‘33**¥**¥***¥¥*t***** C THIS SUBFDUTINE CREATEB A WEEK-MATRIX BY LATTICE SUM * Cit3333*!¥¥‘¥¥**¥¥***¥**¥*¥****t*¥*¥t*¥******¥¥¥¥¥¥33113338*ttt «I 00 “‘1 126 SUBWIINE LSMATRIX(WAT,PX,PY,PZ,NN,NS,NDIM) IMPLICIT REAL*8 (a-h,o-z) IMPLICIT INTEGER‘4(I-N) C(I'IPLEX‘IS IMAT(NDIM,NDIM) REAL PX(NS),PY(NS),PZ(NS) COWION /PP/ NCEL,1‘CELL,PRD,APHA,RADCUI2,PI,PIZ EXTERNAL ERIC SPI = DSQR'HPI) TREL = 4.D0¥PI/3.D0-4.DO¥APHAt*3/(3.DO*SPI) CREATES THE INITIAL ELDIEN'I‘S OF THE MATRIX DO 777 I = 1, NDIM DO 777 J : 1, NDIM IF(I.EQ.J) TMAT(I,J)=CMPI..X(TRH,0.DO) IF(I.NE.J) 'IMAT(I,J)=CMPLX(0.D0,0.DO) CONTINUE START TO CREATE THE ELEMENTS OF MATRIX THE W IS THE SUM OF k-SPACE AND r-SPACE 1, NS 1,NN )*NN + IS 1, NS 1,NN )tNN + J3 HUJ HIIIIHIIH 888888 THIS CHECKS A SYWIE'TRIC MATRIX. IF(ID.LT.JD) GO TO 3 THIS PART IDES THE k-SPACE SUM m100K=-NCELL,NCELL m110L=-NCELL,NCELL m120M=-NCELL,NCELL DIST2 = (K3K + L‘L + MEMHPRDHRD IF(DIST2.LE.REAL(NCELL¥NCELL).AND.DIST2.NE.0.) THEN XX 2 REAL(K) YK : REAL(L) ZK =. REAL(M) m : (PX(I) - PX(J)) m = (PY(I) - PY(J)) m ‘-' (PZUZ) - PZ(J)) RDK = (M*M(K)+M*M(L)+M*M(M))¥m THE'TA = 2.D0tPI$RDK D = 4 .DO‘PI3EXP( -PI*PI*DIST2/(APHA¥APHA) )‘(XBVTI-IE'TA) x /DIST2 IF(IS.EH.1) THEN IF(JS .m. 1) THEN 'IMAT(ID,JD) = 'IMAT(ID,JD) + D*XK3XK ELSE IF(JS .m. 2) THEN 120 110 100 127 'IMAT(ID,J'D) = IMAT(ID,JD) ELSE IF(JS .EQ. 3) THEN "IMAT(ID,JD) = TMAT(ID,JD) ENDIF ELSE IF(IS .EQ. 2) THEN IF(JS .EQ. 1) TTEN "IMAT(ID,JD) = 'IMAT(ID,JD) ELSE IF(JS .EQ. 2) THEN TMAT(ID,JD) = 'IMAT(ID,JD) ELSE IF(JS .EQ. 3) THEN 'IMAT(ID,JD) = TMAT(ID,JD) ENDIF ELSE IF(IS .EH. 3) THEN IF(JS .m. 1) THEN IMAT(ID,JD) : IMAT(ID,JD) ELSE IF(JS .EQ. 2) THEN IMAT(ID,JD) = 'IMAT(ID,JD) ELSE IF(JS .m. 3) THEN TMAT(ID,JD) = '1MAT(ID,JD) ENDIF ENDIF ENDIF CDNTINUE CDNTINUE WINUE THIS PART 113-THE r-SPACE SUM II) 140 K II) 150 L TX) 160 M m,m m,m -MZELL,MIELL D*XK*YK D3XK1‘ZK D*YK*XK D‘YIUYK D¥YK*ZK D*ZK*XK DtZK*YK DXZIGZK DIST2 = REAL ( K¥K+L$L+M¥M) *PRDiPRD IF(DISTZ.GT.REAL(1‘CELL»CELL)) (I) TO 130 XJ YJ ZJ mu RRY m R2 = RRX3FRX+RRYiRRY+RRZ¥RRZ PX(J)+KtI=RD PY(J)+L¥PRD PZ(J)+M:PM) (PX(I) - XJ) (PY(I) - Y3) (P'Z(I) - ZJ') IF(R2.m.0.DO) (I) TO 130 IF(RZ.GT.RAIXDUT2) CD '10 130 R1 = DSQI'NRZ) R3 = 1.DO/(R1H3) RX = REL/RI RY : FRY/RI R2 = RRZ/Rl ARl = APHAle 128 EARL = ERFC(AR1) EARZ = ECP(-APHA’APHA*R2) B EAR1¥R3+2 .DOtAPEIA*EAR2/(R2*SPI) C 4 .DOtAmAn3tEAR2/SPH6 . DOtAPHAtEARZ/ ( R2*SPI) * +3.DO¥EAR1tR3 IF(IS.EQ.1) THEN IF(JS .EQ. 1) then IF(I.EH.J) THEN MAT(ID,JD) = 'IMAT(ID,JD) - CtRxxRX+B ESE IF(I.NE.J) THEN IF(DIST2.NE.O.) THEN 'IMAT(ID,JD) = TMAT(ID,JD) - CtRXIKRX+B ELSE "IMAT(ID,JD) = 'I‘MAT(ID,JD) — CtRXtRX+B+4.DO*PI/3.DO ENDIF ENDIF ELSE IF(JS .m. 2) THEN 'IMAT(ID,JD) = 'IMAT(ID,JD) - C*RX*RY ESE IF(JS .m. 3) THEN 'I‘MA'I‘(ID,JD) = 'IMAT(ID,JD) - CtRIURZ ENDIF ENDIF IF(IS .FR. 2) THEN IF(JS .m. 1) THEN ’1‘MAT(ID,JD) = 'IMAT(ID,JD) - CtRYflIX ESE IF(JS .m. 2) then IF(IJ‘DJ) THEN 'I‘MAT(ID,JD) = 'IMAT(ID,JD) — CtRYtRY+B ESE IF(I.NE.J) mm IF(DIS’I‘2.NE.O.) THEN 'IMAT(ID,JD) = 'IMAT(ID,JD) - CXRYtRY+B ESE 'IMAT(ID,JD) = 'I‘MAT(ID,JD) - CtRY¥RY+B+4.D0*PI/3.DO ENDIF ENDIF ESE IF(JS .m. 3) THEN 'IMAT(ID,JD) : 'IMAT(ID,JD) - CXRth ENDIF ENDIF IF(IS .133. 3) THEN IF(JS .m. 1) then 'IMAT(ID,JD) = TMAT(ID,JD) - CtRZflUI ESE IF(JS .m. 2) then TMAT(ID,JD) = 'I‘MAT(ID,JD) - C*RZ*RY ESE IF(JS .m. 3) THEN IF(IJHJ) THEN TMAT(ID,JD) = TMAT of the atoms is defined as l (2.11) The following algorithm was written to facilitate obtaining (R3 in a system when periodic boundary conditions are used. For example, let us consider only the x direction. rx(i) = O. DOSOi:1,N rx(i) = rx(i) + ( x(i)-xp(i) ) 151 .O 2? Hp‘:09 4 Q) 3:: E 01.04 A t, C71 0'0 ffi‘f'r * r f r v i r f v 1 Figure 2.1 The helium-helium radial distribution function at reduced number density p‘:0.9 with T‘:30.23. The number of helium atoms is 107. 50 continue Then : rx(i)z + ry(i)2 + rz(i)2 where x(i) is the forward position of the ith molecule after one molecular dynamics step, and xp(i) is the previous position of the ith molecule located inside the simulation box. By using the algorithm we cumulated the difference of displacement between every step to obtain . After we check the boundary condition of the molecules, we then renamed xPU) = 3(0) The diffusion coefficient is given in terms of the dimensionless 2 - (Ra ) and tan as 02 at: > D 2 am —— : ftm ( ) (2.12) t 6t km 6 cum“ # of MDR Figure 2.2 shows a plot of the mean square displacement (RD of helium versus time calculated at pt = 0.9 and T‘ = 30.23. The number 4 of particles is 864, A? is 9.441110- (2.397x10-1ssec), the side MD 153 12.0 Y T I I Y T' I T V T r J 10.0- _ 8.0- ~ /\ N -4 t 6.0: L. V 4.0- - 2.0~ - 0-0 T r * T ' i V i ' I v 6‘ 0.0 2000.0 4000.0 6000.0 8000.0 1.0E+O4 1.2E+ t Figure 2.2 Plot of the mean square displacement (R2) of helium versus time at 920.901 and T‘:30.23. The estimated diffusion constant is 4.86xi0". 154 Table 2.1 The list of the diffusion coefficient of helium at given reduced number density and reduced temperature. p‘ i“ D 0.5 30.23 1.09:10' 2 0.701 30.23 7.27x10" 0.901 30.23 4.86x10“ 0.5 - 0.75 4.65x10‘5 0.7 0.75 2.70xio'5 0.3 0.75 1.55x10‘5 0.9 0.7 5 — 155 length of the box, L, is 9.840, and the total number of steps is 10 ‘, 2.397 picoseconds (ps). Table 2.1 lists the simulation results of the diffusion coefficients with different reduced number density for helium. In the argon case, at T* = 0.75 and p‘ :: 0.805, the number of particles is 864, A?” is 0.032 (10"‘sec), the side length of box is 10.230, and the number of steps is 5000 (50 picoseconds). The slope of (R2) versus time yields a diffusion coefficient 2.35x10n‘5 cmz/sec, which has 3% relative error compared with the data obtained by Rahmanso 52,53 and the experimental value , 2.431(10.5 cmZ/sec. 2.2 THE scaao'nmcsa EQUATION The method we use to study the localization of an excess electron in fluids is based on the symmetrically split operator technique presented by Feit, Fleck, and Steiger‘:3 for solving the one-electron Schrodinger equation. The time-dependent Schrodinger equation is . 6w h 2 V.p(§) I at 2 -- 71-11: V (P + ——I\_— ID (2.13) The split operator Fast Fourier Transform (FFT) method43 provides a very efficient technique for propagating the solution. At each small time step At the propagator is split into kinetic and potential parts. According to the second-order accuracy of the symmetrically split 156 operator solution to the propagation equation, we write the electronic wavefunction as iAth V (r) iAth 2 . ep ~ 2 —-— v - —_ —-— lp([“,t+At) - exp( I ) exp( 1At h ) exp( '1 V ) (“5,0 +0013) (2.14) w(r,t) is a given initial electronic wavefunction. An error proportional to (At ) is due to the noncommutativ1ty of the kinetic and . 3 . . potential energy operators. The symbol 0(At ) means the solution is 3 54 accurate to order (At ) . The FFT method yields the electronic wavefunction on a grid defined by integers N1, N2, and N3 which span on the simulation box with side lengths Lx, Ly, L2. The grid spacing in the direction of x, y, z is szLx/Nl, AyzLy/NZ, and Az=Lz/N3. The grid representation for the wavefunction implies a spatial periodicity. The split operator method involves performing the kinetic part of propagation in momentum space where the kinetic energy is diagonal and the potential part of propagation in coordinate space. The coordinate space Mgt) = Willi“) and the momentum space WW”) are related by iii/z 112/2 N312 (ix my w.. (t) : Z 2 2 (p (t) exp[ 2m ( -— + _— "" Cali/2+1,mz-uz/zn,n:-u3/2+1 Om L" 1" 112 + L2 ) ] (2.15) When the wavefunction propagates for a half time step At/2, the iAth solution of exp(-4—a— V ) w(_r,t) can be obtained by the use of the 9 . . . . . 43 band-limited Fourier series representation , iAth 2 c 2 .71 2 '1 2 wémn(t+m/2) : worm“) exp - (7T) (21!) [( Lx ) +( Ly ) +( L2 ) ] (2.16) Eq. (2.13) can be rewritten as the dimensionless Schrodinger equation, 1 all! h 1 2 Vepuz) i 3 - = ' V (Pi ---—-w. (23]) CF at 2m‘a 02 0 h where 'i' : t/CF, (2.18) x(J : x/O, (2.19) 2 2me0 and CF 3 —'h—. (2.20) 1.1286x104s (sec) for the electron in the The scale factor is CF helium system. The algorithm we use in our calculation is summarized in the following flow chart (see program Feit), 158 A general flow chart of Feit’s Algorithm: Input initial wavefunction, wijk”) (Normalized) (1) ‘- FFT ‘v Wnrmlt) (2) * Factorl ownlt-TAt/Z) J .- IFFT (3) W. . (t+At/2) 111: (4) l * Factor2 wUk(t+At/2) Factor2 (5) «- FFT WW’1ti'At/2) (6) t Factorl wfmn' (t+At) (7) .- IFFT wUkm-At) (8) Normalization DO LOOP, T = ntime * At 159 wijk(t+T) At :3). Time correlation function, . Spectral analysis, Eigenvalues. iAt J) l). The Hamiltonian operator: kinetic and potential energies 2). Ground state: eigenfunction and eigenvalue 3). Projection operator method: excited states 4). Physical quantities: diffusion coefficients In this flow chart, the following symbols and abbreviations are used : FFT : Fast Fourier Transform. IFFT : Inverse Fast Fourier Transform. «- : Perform 1:»: Multiply iAth 2 f 2 m 2 n 2 Factorl 7- exp - (TI-E...) (27!) [(v) 4%?) +(—I:;—) ] (2.21) V (g) Factor2 :: exp(-iAt —3-%— ) (2.22) A full-step propagation, At, is from step (1) to step (8). The computation proceeds as a succession ntime full-step propagations. This flow chart is suitable for both real and imaginary time computation. In the 3-dimensional case, for the initial wavefunction we choose a Gaussian distribution function by adjusting the half-width 160 C according to the boundary where the wavefunction must go to zero. In our calculation, we place the central peak of the Gaussian distribution function on the position (xnpick, ynpick, znpick), which is the position of the removed molecule from the system with position nearest to the center point of the simulation box. The non-normalized initial wavefunction is . 2 . 2 . 2 z wijku) = exp { -[ (x-xnpick) +(y-ynpick) +(z-znpick) ]/ 2C } (2.23) where (x,y,z) is the position of grid point (i,j,k). The boundary condition W(O,t) = w(L,t) -—9 0 must be obeyed. Here L is the side length of the simulation box in units of 0 and C is the standard deviation of Gaussian distribution in. units of 0. We renormalize the wavefunction 1p” i‘(t+At) for every full-step propagation. Most of the computing time is taken in doing the FFT and IFFT. One quench step is defined as a full-step propagation with iAt. In the following we are going to summarize the applications of Feit’s algorithm: 1). The Hamiltonian operator if: The Hamiltonian operator 7! of the electron is 3.2 where the first term on the right hand side is the kinetic energy operator and the second term on the right hand side is the potential energy operator. For a given Mr) the eigenvalues can be obtained from 161 2 0(5) 2 E w(g)“'55. Since the wavefunction is represented on a set of discrete grid points, the expectation value (V) of the potential energy 7(r) is evaluated by multiplying the electron density p = w‘w at grid point (i,j,k) in coordinate space by 7(r) at grid point (1.1.1:) and summing over (i,j,k). The expectation values (K) of the kinetic energy is evaluated by Fourier transforming to momentum space, and multiplying the electron density at grid point (£,m,n) by TZ/Zme and summing over (é,m,n.). The momentum 3" is equal to (ZRC/Lx, Zara/Ly, 21m/Lz). 2). Ground state: eigenfunction and eigenvalue The procedure involves following 'the evolution of the wavefunction in imaginary time t = “301.12.21.55. Therefore, the solution of Schrodinger equation is given in terms of a sum of exponentially decaying terms w(r,t) :ngoa.n ¢n(r) exp(-Enr/h) (2.25) where 0.n is the amplitude, ¢n(r) is the eigenfunction and En is the eigenvalue for the nth state. The principle is that as t is sufficiently large (low temperature), Feit’s algorithm quenches the electronic state down to minimum E o’ the lowest eigenvalue (ground state). Meanwhile, the component of the nth state eigenfunction is reduced relative to the ground state by exp(oAEt), AE : En-Eo. If AE I >> 1 is true, then the propagation time 1: is enough. On the other 162 hand, we do procedure 1). until the energy converges within some tolerance value, and then E0 and wohz) are the ground state energy and the ground state wavefunction, respectively. 3). Projection operator method: . . 21,55 . . The pr0jection operator method permits us to obtain the excited states of the electron. The operation of the projection operator, 7’, on an arbitrary wavefunction (ptr(r) is defined as flip”) : eXp(-t4€)|wtr> (2.26) As 2’ can be represented by a spectral series, the above equation can be rewritten as letr> =ngoexp(-rEn> lwn> (2.2.) From the above procedures 1) and 2), we can obtain- the ground state. Then, we use this projection operator method to project out the ground state to obtain the next excited state. When a new wavefunction is produced after one quench step, we project out those known eigenstates with lower eigenvalues during quenching in order to obtain the excited state correctly. In our calculation, we reorthogonalized ’4’ after every quench step. In order to decrease the computation time, we can reorthogonalize after every two or more quench steps. The formal expression for reorthogonalization at nth eigenstate is 163 n-l w = wn -lgowliwllwn > (2.28) n Every step produces a new wn', which is to be used as the initial wavefunction of the next quench step. If quenched enough, the final wavefunction will be the nth eigenstate. The new trial function is € II n In our test, it is efficient to use a decreasing At during quenching. 2 1 for the A trial test shows that we can use a large step At- (0.2) initial few quench steps until the energy converges within tolerance (0.001). Then, we out At in half and run the same procedure of quenching and energy-convergent testing to obtain the ground state. The procedure is iterated until At is less than 0.01. For the excited state we could not use such a large initial At; for the first excited state we used At = 0.01 as the initial value. (See program Projector) 4). The spectral analysis: For a real time evolution the wavefunction (pi |‘(t+T) consists of J a superposition of oscillating wavefunctions, which is useful for the spectral analysis of eigen-energiesn'12’43’54. During a real time evolution we can calculate the time correlation function C(t+At) = at every full-step time propagation, and after ntime iterations we Fourier transform C(t), (I(E)=IC(t)exp(-itEn/h)dt), to obtain the spectrum corresponding to the eigen-energies. The time step size At determines the maximum absolute spectral energy that can be 164 obtained. The relation is fmax: 1/2At, where fmx is the maximum 56 . . . frequency . The number of the t1me increments, ntime, controls the resolution of the spectrum, 1 1 Af : ntime At : —1-_-. (2.30) A large ntime permits the resolution of fine spectral structure. 2.3 THE ADIABATIC SIMULATION METHOD: The procedure12 used to study the localization of an excess electron in classical fluids is discussed in this section. Let us consider a system consisting of an electron and several hundred classical molecules. The dynamics of the electron and molecules are characterized by widely different time scales, with the motion of the quantum electron much faster than the motion of the molecules. The motion of molecules can be treated classically. Therefore, the adiabatic or Born-Oppenheimer approximation is valid, and its use leads to the terminology adiabatic simulation method. For a given nuclear configuration we obtain the ground state of the electron by performing an imaginary time evolution. with Feit’s algorithm. The adiabatic simulation method is summarized as follows: 1). We freeze the electronic wavefunction and integrate the classical equations of motion for a time step ATHD. 2). A new potential is produced after the molecular motion. The 165 electronic wavefunction certainly will be no longer in the ground state with the new potential. Therefore, we freeze the molecular configuration and quench the electron to the new ground state by use of the imaginary time evolution method. We repeat the procedures: step 1) and step 2). The interaction of the electron and molecules is given by a psuedopotential“ A B V (r) : —- ( -—-— - 1 ) ( in a.u.l (2.31) ep ~ 4 a 6 r c + r where for electron-helium case A = 0.655, B = 89099, and C = 12608 (in a.u.). 5 A B V (E) : 3.157x10 -: ( -——-—6-— - 1 ) ( in K ) _ (2.32) ep r C + r a (K) here r = r00“, and 0“ : . 0.529 (K/e.u.) The units of V.p(g) and distance r between electron and molecules are in K and in a.u.. Every one or several molecular dynamics steps is followed by ntime quench steps in order to quench the electronic quantum state to the ground state or excited-state for a given frozen nuclear configuration. The number of quenching steps is determined by the 166 21 energy convergence test . 6-7,ii-12 . , we must include When we consider the molecular motion the contribution of the quantum force arising from the interaction between the electron and the molecules. The classical equations of motion are then written as 2 d R. N 6V (R..) EV (r,(R..}) m--———”1 : - ————-M 1 - d3r )tp(r,t)|2 ep Q (2.33) 2 . . a R. ~ ~ 5 R. dt jati i i 2 m0 The factor -- is extracted from the above equation in order to C MD produce the dimensionless equation, 2 2 dR . N_ C 6V (r,(R..}) -—:£L‘— : X F(R..) - --1°— d3r lw(r,t))2 “P ‘1 (2.34) -2 . .~ ~ij 2 ~ ~ 6 R. dt J¢1 m0 1 - Cm,Z avuml ) where, F(R..) : - 4 (2.35) ~ ~ij m02 6 R1 The second term on the right hand side is the quantum force 2 .— IIID 3 2 ~ ~1 3 2 OP ' "TJ d E ““5"" '3 a. J = ' 1 d 28 Wart" TEE—T m0 1 0,1 (2.36) - 01402 where V9p : ---—2- V.p(ro,(Ra’ij}). (2.37) m0 167 There are two places we need to include this quantum force: 1) In 1 the molecular equation of motion, V : Vep, and 2) In the C ep 48E .. . . —. - F Schrodinger equation, Vep - h Vep. — 3 2 6 Van The Vep and [— d 50 [MEGAN WI] terms are calculated by 2 47,48 the minimum image convention method . Since we quench the electron to the ground state with the lowest energy in a closed system, the decreasing energy will lower the temperature of the classical system and result in a deviation of temperature of the system from the original value. Therefore, the temperature fluctuations of 15% were allowed during the simulation. Otherwise, we renormalized velocities of molecules to restore the original temperature. 2.4 THE MOVING GRID TECHNIQUE: (SEE PROGRAM MGT) This technique” was developed to monitor the behavior of the electron in a finite system extended to an infinite system via the application of periodic boundary conditionsn’m. The ground-state electronic wavefunction produced after a quench has very small values on the boundary of a grid box which is small compared to the simulation box. Therefore, a smaller grid box can be used relative to the simulation box. The advantage of this technique is that it allows us to use a smaller grid spacing and number of grid points which decreases 168 the computation time and increases the spatial resolution. The grid box is moved to follow the center of mass of the electronic wavefunction density distribution, - X m" (t) (2 38) r ' j {ljk wijk wijk ' where {in denotes the position of the grid point (1.1.1:), and wUku) is the electronic wavefunction on the grid point (i,j,k). For a given frozen nuclear configuration and a fixed grid box we quench the electronic state to the ground state and calculate the center of mass of the electron. Whenever we obtain the center of mass of the electron, we move the center point of the grid box close to the new center of mass of the electron and re-express the electronic wavefunction relative to this new origin of the grid box. When the center of mass of the electron doesn’t move, the grid box won’t move. 3. RESULTS AND DISCUSSIONS In this Chapter, the various equilibrium and dynamical properties for both the electron and the helium, obtained by using the molecular dynamics adiabatic simulation method described in Section 2.2, will be presented and discussed. Figure 3.1 shows that information concerning the spatial extent of the electron can be obtained from a correlation function such as the radial distribution function of the electronic probability density measured from the mean position of the electron, 1;” = I I: Mr) IZ d1: — the center of mass of the electron. The results reported here are indicative of the various sizes of the localized electron versus the 0.1, 0.5, 0.7, and 0.9 in a helium system reduced number density p‘ at reduced temperature T‘ = 30.23. In the present calculation, the electronic wavefunction was represented on 163 grid points and the number of the helium atoms is 107. Atypical imaginary time step used in this process was At}, = 0.01063. We have used 200 quench steps to obtain the ground state of electron for each different p‘. At p‘=0.1 the electron is not localizedaib, because the probability does not go to zero of large distance from the electron center of mass. Comparison of the results as the molecular density 9* is increased shows that the size of the localized electron decreases and that it is well localized in a given equilibrated helium configuration. 169 170 1.2 7 o p'- 0.1 q a p'- 0.5 1.0-1 - - x p‘- 0.7 .1 A p'- 0.9 0.84 1’ g 0.64 3 m 0.4-1 0.2“ ‘- ~o . 4 ‘. O O . . . 0,0 , I W 0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 r/A Figure 3.1 Plots of the radial distribution function. of the electron’s probability density vs. the distance from the electron center of mass at various densities in helium at T‘=30.23. 171 5.00 1 e , ; , . , . , . J 4.004 4 4 1 300-4 .4 >. 4 4 2.004 . 1 . 1.00-1 J 4 0.00 . r T r T , , T , 0.00 1.00 2.00 3.00 4.00 5.00 X Figure 3.2 Contours of the excess electron density corresponding to the localized state wavefunction in the x-y plane. 172 More details of the probability density distribution, lflrllz, of the localized electron can be obtained from the contour plot. Figure 3.2 shows p(x,y), a contour plot on the x-y plane, corresponding to the simulation of Figure 3.1, at the reduced density 0.9. Here p(x,y):f|¢(r))2dz. The values indicated on the contour plot are in arbitrary units. As can be seen, the ground state of the localized electron is quite symmetrical and has an s-like structure. (For data contour plotting, see subroutine CONTOUR). Information about the local fluid environment around the electron was obtained from the helium radial distribution measured from the center of mass of the electron. In Figure 3:3, the number of the helium atoms is 107, the reduced number density is 0.7, and the temperature is 309 K. In this run, the first 1.1><10"s molecular dynamics steps were used to equilibrate the helium and the radial distribution function was calculated in the following 10‘ steps. For the ground state quenches we have used 20 quench steps after each molecular dynamics step. The imaginary time step and the molecular dynamics time step are, respectively, 0.01063 and 9.4x10'4. In this run, the total real evolution time is 2.3 ps. The results show that there is zero probability of finding .the helium atoms in the region of the localized electron. As mentioned in Chapter 1, as the helium density is increased to liquid density, the electron becomes localized in a small region. A bubble-like structure occurs because the repulsive electron-helium interaction produces a cage effect which excludes the helium atoms from the region of the electron. 173 U .s... .. ..___J N L g(r) cm—He m... 0 N I 6 r/A Figure 3.3 The radial distribution function of helium measured from the center of mass of the electron. 174 4O ' I V I V— I 1'— T fl I fir r T T i J 2: 304 /'—1‘ ‘ V l 8 1 3‘ ( ,r\ (.0. 5. E 20-1 .1 ......... q x L—‘H >~ 1". C7) 4 \l. 4 a V 1‘. .......... 15 10- %.\ , J A K T“ _----."."_'::::; ’---'-’ """"""""""" g 0 j T— T r I fir : r Y r 1f : f r I— 0 20 4O 60 80 100 120 140 150 number of quenching steps Figure 3.4 Plots of the electronic energy versus the number of quench steps. 175 In Figure 3.4 we show that the energy convergence of the ground state electron varies with the number of quench steps and the quench time step from 0.2 to 0.01. From this plot we know that, in this case, 150 quench steps are enough to obtain the lowest energy state. We show the typical electronic configurations in Figure 3.5, which exhibits the localized ground state 00 and first excited state 01 obtained via the projection operator method, as mentioned in Section 2.2. Here, the number of helium atoms was 107, the reduced density was 0.9 and temperature was 309 K. In order to obtain the ground state, the initial dt-F was chosen to be 0.5 and the tolerance value of the energy convergence was 0.001. The total number of quench steps was '155 when the final dTF was less than 0.01. The corresponding kinetic, potential, and total energy of the ground state electron are, respectively, 3.79, 1.78, and 5.57 (3.25 eV) in rescaled units (one rescaled unit is equal to 6,768 K). The parameters used to obtain the first excited state by projecting out the ground state were the initial d-t'F= 0.5, energy convergence tolerance 0.01, and final di} 5 0.001. The total number of quench steps was 550. The corresponding kinetic, potential, and total energy of the first excited state are, respectively, 5.44, 2.87, and 8.31 (4.85 eV). Figure 3.5 shows that the ground and first excited state have, respectively, s-like and p-like structure. During the simulation process, 00 and 01 fluctuate leading to electronic configurations with considerable occasional distortion from s-like and p—like structure. It is worth comparing the present results for the ratio of the Assam “5:88 ass: .3 can 033m «5:30am .3 15:53 a .5“ 5503233200 015.3020 of mo 3:03:00 md 9235 17 177 30 f I I r r x 4 «'3 204 . 1" CD :4 1 1 X 8'10.) W Ed 5 L 1 or ‘T r f I r f 0 10m 2000 3000 30. - , . , - I 1 . X 8 204 .1 [x Q X 1} 4 5‘- 3 10“ K4 C o ‘ . ‘ W 0 a ' E r - r 10” 2000 3000 ac r I T I x 4 3'3 ”1 1 a X I 4 3 5.. 104 VJ) 5 4 , ‘ c 'i I fi r ‘E I )0” 2000 3000 number of time steps Figure 3.6 Plots of the fluctuation of the total electronic energy (E). kinetic energy (K). and potential energy (V) in the molecular dynamics adiabatic simulation. 178 total energies of the ground and first excited states with those , 31b . . obtained by Coker and Berne for an electron in a helium configuration taken from the path-integral Monte Carlo simulationau. The ratio is 1.66 in their result and our ratio is 1.49 at 920.9. In addition, our results for the energy are roughly 1.7 times larger than theirs. Our electron is systematically smaller and, as a result, our kinetic and potential energies are larger. These differences arise because Coker and Bernemb quenched the electron without allowing the solvent configuration to relax. Since the Path—Integral Monte Carlo produces a thermally distributed electron, which is larger than the ground state electron we create, our solvent equilibrated electron will be smaller than their, leading to a higher kinetic energy. Also, our solvent molecules will be closer to the electron leading to a higher potential energy. Figure 3.6 shows that the excess electron’s ground state energy E, kinetic energy K, and potential energy V fluctuate as the helium system evolves in time. In this run, the number of particles and reduced density are, respectively, 255 and 0.9 at I‘ = 309. The molecular dynamics time step is 2.61940"3 and the quench time step is chosen from 0.2 down to 0.01 followed by the energy convergence. Quenching to the ground state followed after every 10 molecular dynamics steps. The real time span is 3000 steps (Zps). The computation time on the Titan was 32 hours. In this context we note that the time variation of the electronic potential energy exhibits a stronger fluctuation compared with that of the kinetic energy. In other words, the electron is very sensitive to the local fluid 179 environment. The confined electron usually has higher kinetic energy than the potential energy. In accord with the model of a particle in a spherical boxa, the simulation demonstrates a high kinetic energy for a well localized electron, which agrees with expectations. The averages of the electron kinetic and potential energy over the 3000 steps are, respectively, 3.49 and 3.45. Examination of the ground state energies shown in Figure 3.6, shows that the kinetic energy does not change significantly, but the potential energy does. It appears that the electron must go over a potential barrier during the electron transition from one localized state to another localized state. There is an intermediate state of high potential. The diffusion constants of an excess electron in molten salt6’7, polar water20 and ammonia“ have been obtained by simulation, and the results are in agreement with experiment57'59. The diffusion constants of the electron simulated in water‘20 and in KCl6 are 3.31043 and 2.0x10'3 cm2 seed, respectively. It is of great interest to simulate the diffusion constant of the excess electron in dense helium and we hope to understand the electron diffusion mechanism. We assume that the electron is well localized at all times. We monitored the electron dynamics by following the motion of the center of 'mass of the electron, as if it were a classical particle. Figure 3.7 shows the mean-square displacement of the center of mass of the electron as a function of time for an electron in helium. The number of helium atoms is 255 in the simulation box at p‘: 0.9 and T‘: 30.23. The length of the grid box was chosen to be about 75% of the simulation box in order to obtain the same resolution as in the system of 107 particles. The grid points are fixed electron 4.48x10‘3 - .6 bellom , in water, 180 in number, at 163. The result given is an average over 12 trajectories. The diffusion constant obtained is about cm2 sec”. Based on the results obtained by Rosskyzo and the diffusion of the electron in helium is faster than that and it is the same order as in KCl. The mechanism is not clear. .‘vlore numerical experiments are essential. 181 U‘ 0') O O usudssssfl l J 1 .1 .1 J 1 1 p. O (A O 10 3 3 3 '3 3 3 '3 3 20 3 O l T r j I I f f' I ' I I r T I r 0 200 400 600 800 1000 1200 1400 1600 number of time steps Figure 3.7 The time dependent mean-square displacement of the center-of-mass of an electron in helium. 4. CONCLUSIONS We have described and programmed computational algorithms for a time—dependent quantum molecular dynamics simulation with application to both the molecular and electronic configurations in fluids. It is evident that the technique combines the quantum and classical aspects, and efficiently and quickly solves the problem. In the present context, we applied the method to the particular system, an excess electron in helium. We were able to obtain the ground state, the excited state, the corresponding energies of the electron, and the correlation functions of the electron and the local environment molecules. We have studied the diffusive behavior of an excess ground state electron in helium. In order to progress further, there is the essential need for electron-molecular pseudopotentials, which are difficult to determine, and molecule-molecular interaction potentials, which are better known. An optimal molecular dynamics time step is necessary. With the electron dynamics process, in order to decrease the integration error, the final quench time step must be small enough and there must be a sufficient number of quench steps. This technique will play an important role in future simulations because of its accuracy, flexibility, and efficiency, leading to 182 ' 183 greater insight into dynamical molecular processes. The problems for future investigation include electron transfer dynamics in solution and in biological systems, and the electron transport rate. The transition from an excited state to the ground state, and the resulting optical absorption spectra can also be obtained by this approach. APPENDICES APPENDIX A Reference Table A.l lists the information on the simulation parameters to facilitate comparison. Table .—\.l The list of the reference parameters: 6‘ 0.1 0.3 0.5 0.7 0.9 N 107 107 107 107 107 ,; 10.23 7.09 5.98 5.35 4.92 (a) L 26.14 18.13 15.29 13.67 12.56 (91), He 0N 5.9911021 1.8x10922 3.0::1022 4.2211022 5.39111032 (cm'a) pm 0.04 0.12 0.20 0.28 0.36 (g cm'31 0 0.05 0.16 0.26 0.37 0.47 The simulation box side length is L = C a. p‘K : The reduced number of density, N/éa. N : The number of particles. C : The box side length in units of the diameter of the particle, a. 0N : The number density of particles, N/La. pm : The mass density, x m. 9N m : The mass of particles, mm: 4.0026xl.6747x10'2‘gram. (b : The volume fraction, pat/6. 184 B. APPENDIXB LIST OF C‘Q‘IPUTER W Ct*ttt!¥¥¥¥*¥****tt*¥****¥¥¥¥t¥tt¥t¥¥¥¥¥¥¥t¥t¥t¥$¥¥t¥titi‘tt! 000000000000000OOOOOOOOOOOOOOOOOOOO(300000000 LIST OF NOIATIONS ABC : ENERGY CONVERGENCE I‘OlmANCE, 0.001 ABCD : MINIMUN QUENCH TIME STEP, 0.01 Ci‘lD : THE RESCALE FACTOR FOR IDLECULAR DYNAMICS CMX,CMY,CMZ,C§DCP,CMYP,CMZP : THE CENTER OF MASS OF ELECTRON CX(N),CY(N),CZ(N) : THE CUMULATED DISPLACEMENT OF PARTICLES CXE,CYE,CZE : THE CLMUIATED DISPMCEMENT OF ELECTRON DIS : MEAN SQUARE DISPLACBIENT OF PARTICLES DIFF : DIFFUSION CONSTANT OF PARTICLES DISPMX : LARGEST DISPLACEMENT DX : TIE GRID SPACING IN x DIRECTION DY : TIE GRID SPACING IN Y DIRECTION DZ : TIE GRID SPACING IN 2 DIRECTION DMD, II : THETIMESTEPFORMOLECULARDYNAMICS DT, DTFT : TIE TIME STEP FOR QLENCII DIAM : TIE DIAMETER OF PARTICLE IN ANGSTI'KJII ’ N : KINETIC ENERGY OF PARTICLES F : FORCE FX(N),FY(N),FZ(N) : TIE FORCES OF PARTICLES GIIE) : RADIAL DIS. FUNC OF ELECTRON GN * Cit!t***t*1t********¥**************¥***¥¥*¥*******************¥ 21 SUBROUTINE NORM(W1,W2,CSUMJ IMPLICIT INTEGER¥4(I-N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER ( N1=16, N2:16, N3=16 ) COMPLEX316 W1(N1,N2,N3) * ,W2(N1,N2,N3) ’ ,CP(N1,N2),CH(N1),CSUM MN /Cl/ DX,DY,DZ DO 21 IX=1,N1 DO 21 IY:1,N2 CP(IX,IY)=(0.0,0.0) DO 21 IZ=1,N3-1 CP(IX,IY)=CP(IX,IY)+DZ*(CONJG(W1(IX,IY,IZ))*W2(IX,IY,IZ)+ * OONJG(W1(IX,IY,IZ+1))*W2(IX,IY,IZ+1))*0;5 CONTINUE DO 31 IX=1,N1 CH(IX)=(0.0,0.0) DO 31 IY=1,N2-1 - CH(IX)=CH(IX)+DY*(CP(IX,IY)+CP(IX,IY+1))*O.5 199 31 CONTINUE CSUM:(0.0,0.0) DO 41 I=l,N1-1 CSCM=CSUM+(CH(I)+CH(I+1))*DX*O.5 41 CONTINUE RETURN END cx:xxt:xttxtttxxtxtxtxxttxxxtxtxxxxtxxxtxxxxttxxx:x:xxxxxxx:xtx C SCEROUTINE INTEGRAL : INTEGRATION IN 3-DIMENSION * C USAGE : wz CORRESPONDS TO POT, KIN AND TOTAL ENERGY * Ctittxtttttttxtttttttt*ttttttttxxttxttt¥¥***¥t*¥*¥***tttxttxttx SUBROUTINE INTEGRAL(W2,CSUM) IMPLICIT INTEGER¥4(I-N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER ( N1=16, N2=16, N3=16 ) COMPLEX*16 W2(N1,N2,N3) * ,CP(N1,N2),CH(N1),CSUM comm /C1/ DX,DY,DZ DO 21 IX=1,N1 DO 21 IY=1,N2 CP(IX,IY)=(0.0,0.0) DO 21 IZ=1,N3-1 CP(IX,IY)=CP(IX,IY)+DZ*(W2(IX,IY,IZ)+ t W2(IX,IY,IZ+1))‘0.5 21 CONTINUE DO 31 IX=1,N1 CH(IX):(0.0,0.0) DO 31 IY=1,N2-1 , CH(IX)=CH(IX)+DY‘(CP(IX,IY)+CP(IX,IY+1))*O.5 31 CONTINUE CSUM:(0.0,0.0) DO 41 I=l,N1-1 CSUM=CSUM+(CH(I)+CH(I+1))tDXtO.5 41 CONTINUE RETURN END Ctttiitttittittttititt¥!¥¥t*¥**tittttttttttttittXttttttttttt C SUEROUTINE OONTIUR : PLOT OF CONTOURS * C USAGE : ISIGN DECIDEE ON THE X-Y OR Y-Z OR X-Z PLANES¥ C ISIGNzl INTEGRATE IN X DIRECTION * C ISICN=2 INTEGRATE IN Y DIRECTION * C ISIGN=3 INTEGRATE IN 2 DIRECTION * Ct*¥**1*tttt*tt!3it¥33*3tit!Iittittttttttttttt¥¥***¥¥tt¥¥**t SUBROUTINE CONTOUR(ISIGN) IMPLICIT INTEGER34(I-N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(N1=16,N2=16,N3=16) COMPLEthS W(N1,N2,N3) "L 200 COMMON /A/ W COMMON /C1/ DX,DY,DZ ———WRITE OUT DATA IN UNIT:NFILE-—- IF(ISIGN.EQ.1) NFILE=51 IF(ISIGN.EQ.2) NFILEz52 IF(ISIGN.EQ.3) NFILE=53 DO 42 I=l,N1+1 DO 42 J=1,N2+1 IX=I IY=J IF(IX.EQ.N1+1) IX=1 IF(IY.EQ.N2+1) IY=1 SUM=0. DO 56 K=1,N3 IZ:K+1 IF(IZ.EQ.N3+1) IZ=1 IF(ISIGN.EQ.1) THEN SUM:SUM+(CONJG(W(K,IX,IY))xW(K,IX.IY)+ ¥ CONJG(W(IZ,IX,IY))*W(IZ,IX,IY))*DX*0.5 ELSE IF(ISIGN.EQ.2) THEN SUM=SUM+(CONJG(W(IX,K,IY))¥W(IX,K,IY)+ * CONJG(W(IX,IZ,IY))*W(IX,IZ,IY))*DX¥O.5 ElSE IF(ISIGN.EQ.3) THEN SUM=SUM+(OONJG(W(IX,IY,K))IWKIX,IY5K)+ t CONJG(W(IX,IY,IZ))*W(IX,IY,IZ))tDX*0.5 ENDIF CONTINUE NNN=MOD(IX,2) PT!k¢KEHIY,2) IF(NNNLEQ.1MMELQQELEQ.1)THfifld IF(SUM.LT.1OE-6) SUM=0. WRITE(NFILE,¥)REAL((I—1)¥DX),REAL((J-1)*DY),REAL(SUM) ENDIF CONTINUE RETURN END 201 ext:xxxtx::xxxxxtxtxxxxxxxtxxxxxttxxxxxxxxxxtxtxxxxxxxxttxxx C SLEROUI‘INE : PEIT, PROJECT, HGT t C USAGE : SOLVE THE SCHRODINGER EQUATION * C PROJECT OUT THE GROUND STATE t C BY MOVING GRID TECHNIQUE . x cxx::xzxxxxxt:xt:xxtxxxxxxtxxxxxxxxxxx:xx*xxxxtxxxxx:xtxxxtx SUBROUTINE FEIT1N,RX,RY,RZ,VEIX,VEIY,VEIZ,DT,NTIME) IMPLICIT INTEGER*4(I-N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER ( NDIM=3, N1=16, N2=16 ,N3=16, 4 TEI=46.625,TEIMD=643.55) DIMENSION NN(NDIM),VEI(N1,N2,N3) .RX1N).RY(N),RZ(N) ,VEIX(N),VEIY(N),VEIZ(N) yH(N1).P(N1.N2) ,G(50),RAD(50+1),GN(50),RADS(50+1) COMPLEthS W(N1,N2,N3),W0(N1,N2,N3) ,THETA1(N1,N2,N3) ,FACT2(N1,N2,N3),FACT1(N1,N2,N3) ,CP(N1,N2),CH(N1),CSUM,TSUM,PSUM ,POT(N1,N2,N3),SW(N1,N2,N3) ,CW(N1,N2,N3),TW(N1,N2,N3) ,TOT(N1,N2,N3),TE(N1,N2,N3) ,PP(N1,N2,N3),WN(N1,N2,N3) ”*‘I’fi *flflfififlfl INTEGER CLOCK, TIMEI COMMON /A/ w W /B/ SIDE, RSIDE, SICMA, CLOCK, QISUM COMMON lBI/ XNPICK, YNPICK, ZNPICK COMMON /C1/ DX,DY,DZ,CXE,CYE,CZE COMMON /C2/ ITXS,ITYS,ITZS,CMXP,CMYP,CMZP,PREF MN /C3/ ABC,ABCD,IE,IH,DIAM,NW,TIMEI NN(1):N1 NN(2):N2 NN(3)=N3 NT=N1*N2*N3 NT2=2*NT SIDEH=SIDE/2. GSIDE=(107./0.7)*¥(1./3.) GSIDEH=GSIDE/2. XL=GSIDE SIDE4=(SIDE-GSIDE)/2. SIDEF=SIDE4+GSIDE DX=XL/FLOAT(N1) DY=XL/FDOAT(N2) DZ=XL/FDOAT(N3) PI=4.*ATAN(1.) PIZ:2.*PI CSUMP=0. C ---SET INITIAL WAVE FUNCTION S POTENTIAL --- RR=4.83DO 202 AA:O.655DO BB=89OQQ.DO CC=12608.DO STANDV:SIGMA*0.5DO ----CLOCK COND MUST AGREE WITH PSL2.F--- IF(CLOCK.EQ.TIMEI) THEN O \"IMGzo . DO 71 IX=1,N1 DO 71 IY=1,N2 DO 71 IZ=1,N3 C --- SHIFT POINTS TO SYMMETRIC INTERVAL FOR FFT--- C --- INITAL GRID BOX IN CENTER OF SIMULATION BOX--- SX:SIDE4+FDOAT(IX-1)¥DX SY=SIDE4+FLOAT(IY-l)*DY SZ=SIDE4+FLOAT(IZ-l)*DZ VREAL=((EXP(-((SX-XNPICK)*(SX-XNPICK)+ x (SY-YNPICK) T: (SY-YNPICK) + x . (SZ-ZNPICK)*(SZ-ZNPICK)) I /(2. *STANDVIISTANDV) ) )) W0(IX,IY,IZ)=CMPLX(VREAL,VIMG) 71 CONTINUE C --- NORMALIZE INITIAL WAVEFUNCTIO ----I CALL NORM(WO,WO,SUM) RSQRTSUM=1./SQRT(SUM) DO 232 IX=1,N1 DO 232 IY=1,N2 DO 232 IZ=1,N3 WO(IX,IY,IZ)=WO(IX,IY,IZ)*RSQRTSUM W(IX,IY,IZ)=WO(IX,IY,IZ) 232 CONTINUE ENDIF --- POTENTIAL MEASURED AT GRID POINT--- ----- SET UP SIDES OF GRID BO --- —---FIRST TIME THRU ITXS=0 FROM MAIN, THEN OUTPUT OF FEIT--- ---FIRST TIME cmcp=o FROM MAIN, THEN OUTPUT OF FEET-~- -----DECIDE WHICH BOXES TO USE ----- NXP=ANIN'I‘(CMXP/SIDE) NYP=ANINT(G‘1YP/SIDE) NZP=ANINI‘(CMZP/SIDE) 00000 C ------ LOOP 1 OVER GRID---- DO 1 IX=1,N1 DO 1 IY=1,N2 DO 1 IZ=1,N3 SX=SIDE4+FLOAT(IX-1+ITXS)¥DX SY:SIDE4+FLOAT(IY-1+ITYS)*DY SZ=SIDE4+FIDAT( IZ-1+ITZS ) *DZ VEI(IX,IY,IZ):0.0 C ----- IF PARTICLE FROM ONE BOX IS IN GRID KEEP IT, OTHERWISE-— 000 133 203 ---CHECK IF IMAGE PARTICLE IS IN GRID -------- ---FOR PARTICLE IN GRID EVALUATE XX, PARTICLE -GRID DISTANCE-- ------- LOOP 21 OVER PARTICLES ——— -— DO 21 I:1,N XP=RX(I)+SIDE*NXP XX=SX-XP IF(ABS(XX).GT.GSIDEH.AND.ABS(XX+SIDE).GT.GSIDEH) GO TO 21 IF(XX.LT.-SIDEH) XX=XX+SIDE IF(XX.GT. SIDEH) XX=XX-SIDE YP=RY(I)+SIDE¥NYP YY=SY-YP IF(ABS(YY).GT.GSIDEH.AND.ABS(YY+SIDE).GT.GSIDEH) GO TO 21 IF(YY.LT.-SIDEH) YY=YY+SIDE IF(YY.GT. SIDEH) YY=YY~SIDE ZP=RZ( I ) +SIDE¥NZP ZZ=SZ—ZP IF(ABS(ZZ).GT.GSIDEH.AND.ABS(ZZ+SIDE).GT.GSIDEH) GO TO 21 IF(ZZ. LT. -SIDEH) ZZ=ZZ+SIDE IF(ZZ.GT. SIDEH) ZZ=ZZ-SIDE RIJSQ=1XX¥XXH (YYtYYH (221122) VEI(IX,IY,IZ)=VEI(IX,IY,IZ)+TEI*AA/((RRtRSIDEHIMH (BB/(CC+RIJSQ*RIJSQ*RIJSQ*(RSIDExRRHtS)-1. ) / (RIJSQtRIJSQ) ----- CLCBE ILDP OVER PARTICLES------- mNTINUE CIDNTINUE ------ BEGIN PEIT STEP -------- DO 2 IT=1,NTIME ------ CALCULATE ENEHIY -------- ------- POTENTIAL ENERGY--------- m 133 LY:1,N1 DO 133 IY=1,N2 IX) 133 IZ=1,N3 CW(IX,IY,IZ)=(X)NJG(W(IX,IY,IZ)) PUT(IX,IY,IZ)=CDNJG(W(D(,IY,IZ))#W(IX,IY,IZ)*VEI(IX,IY,IZ) wNTINUE --IN USING FFI‘ ONLY W IS TRANSFORMED PASSED IN DMDN A--- ----- FFI‘ «47R!!! R SPACE TO K SPACE-- CALL FUJRN AND DIFF. CONSTANT IN UNIT=330-- WRITE(66,*)CIDCK~TLVIEI,C§D(,CMY,CMZ WRITE(6,¥) ’NO OF QUENCH’ ,IT-l kRITE(44,*)CImK,REAL(PSUM),REAL(TSUM) ,REAL(CSLM) DIFFzFREFtRSQ/FLDANCIQIK) kRITE(330,¥)CILXIK,RSQ,DIFF ENDIF ----RELABEL CM TO ’PREVIOUS' CM ------ CMXP:CMX ‘ CMY'P=CMY CMZP=CMZ IF( (CLOCK/NW)*NW.m.CLOCK) THEN ---PAIR DISTRIBUTION FL’NCI'ION--- RAD(l):0.0 DR:GSIDEH/FIDAT( IE-l ) CONST:4.0*3.141590/3.0 DO 13 I=1,IE G(I):0. MD(I+1):RAD(I)+DR CONTINUE ------- -G(R) - CM - E---- SUMPon. DO 84 I=l,N1 DO 84 J:1,N2 DO 84 K:1,N3 SX=SIDE4+FLOAT( I-1+ITXS ) #DX-CMXP SY :SIDE4+FLOAT( J - 1+ITYS ) tDY-Q‘IYP SZ=SIDE4+FLOAT(K-1+ITZS) *DZ-QWZP DIS=SG?I‘( SX**2+SY"2+SZ"2) PPP=OONJG(W(I,J,K))¥W(I,J,K) DO 51 IR=1,IE IF(DIS.GI‘.RAD(IR).AND.DIS.LE.RAD(IR+1)) THEN G(IR)=G(IR)+PPP*(1.0/((RAD(IR+1)H3-RAD(IR)H3)¥CX)NST)) SLMPP=SUMPP+PPP GO TO 84 ELSE ENDIF MIME CONTINUE RHOW=SUMP/(GSIDEH3) WRI'I'E(10,*)CI£XZK DO 7 K=1,IE r:(ra.d(k)+rad(k+1))*DIAM/2. write(10,*) r,G(K)/RHOW continue WRITE(51,¥)CIOCK CALL WU) ENDIF 000 33 00000000000 0000000000 208 ---- CALCULATE PARTICLE-ELECTRON FORCE ------- ---LOOP 169 OVER PARTICLES, 170 OVER GRID ------- ----FOR FIXED PARTICLE SUM OVER GRID ------ FIND PARTICLE LOCATIONS IN GRID BOX---- ---IF PARTICLE NOT IN GRID BOX, SET VEIX,..=O---- IF((CLOCK/NW)¥NW.EQ.CLOCK) THEN RAA=GSIDEH RADS(1):0.0 DRS:RAA/FIDAT( IH-l ) CONST=4.0*3.141590/3.0 DO 33 I:1,IH RADS(I+1):RADS(I)+DRS GN(I):O. CONTINUE ENDIF —-SET THE RANGE OF GRID BOX--— DISI=SIDEF+ITXS*DX DISZ:SIDE4+ITXS*DX DISB=SIDEF+ITYS*DY DIS4=SIDE4+ITYS*DY DISS=SIDEF+ITZS*DZ DI86=SIDE4+ITZS*DZ NXP=NINT(CMXP/SIDE) NXM=NINT(CMXP/SIDE)-1 NYP=NINT(CMYP/SIDE) NYM=NINT(CMYP/SIDE)-l NZP=NINT ( M/ S IDE ) NZM:NINT(CMZP/SIDE)-1 DO 169 J=1,N XM:RX(I)+SIDEtNXM XP=RX(I)+SIDEtNXP IF(XM.LT.DISI.AND.XM.GT.DISZ) THEN XX=XM ELSE IF(XP.LT.DISl.AND.XP.GT.DISZ) THEN XX=XP ELSE GO TO 666 ENDIF IF(XX.LT.-SIDEH) XX=XX+SIDE IF(XX.GT. SIDEH) XX=XX-SIDE YP=RY(I)+SIDEtNYP YM:RY(I)+SIDE*NYM IF(YM.LT.DIS3.AND.YM.GT.DIS4) THEN YY=YM ELSE IF(YP.LT.DIS3.AND.YP.GT.DIS4) THEN YY=YP ELSE GO TO 666 ENDIF IF(YY.LT.-SIDEH) YY=YY+SIDE 00000000000 0 0 0 112 209 IF(YY.GT. SIDEH) YYzYY-SIDE ZP=RZ(I)+SIDE*NZP ZM=R2(I)+SIDE*NZM IF(ZM.LT.DISS.AND.ZM.GT.DISS) THEN ZZ=ZM ELSE IF(ZP.LT.DISS.AND.ZP.GT.DISS) THEN ZZ=ZP ELSE GO TO 666 ENDIF ' IF(ZZ.LT.-SIDEH) ZZ=ZZ+SIDE IF(ZZ.GT. SIDEH) ZZ=ZZ-SIDE ----- —CALCULATE RADIAL DISTRIBUTION OF PARTICLES--- IF( (CIDCK/NWHNW.EQ.CLOCK) THEN ---G(R) - CM — PARTICLES---- DO 112 IR:1,IH XX:RX(J)+SIDE*NXP-CMXP IF(XX.LT.-SIDEH) XX=XX+SIDE IF(XX.GT. SIDEH) XX=XX-SIDE YY:RY(J)+SIDE*NYP-CMYP IF(YY.LT.-SIDEH) YY=YY+SIDE IF(YY.GT. SIDEH) YY=YY-SIDE ZZ:RZ(J)+SIDE*NZP-CMZP IF(ZZ.LT.-SIDEH) ZZ=ZZ+SIDE IF(ZZ.GT. SIDEH) ZZ:ZZ-SIDE ----- -CALCULATE RADIAL DISTRIBUTION OF PARTICLES—-- DIS=sqrt((XX)¥¥2+(YY)**2+(ZZ)¥*2) IF(DIS.GT.RADS(IR).AND.DIS.LE.RADS(IR+1)) THEN GN(IR)=GN(IR)+1./((RADS(IR+1)“3-RADS(IR)¥¥3)*CONST) ENDIF CONTINUE ENDIF ----- EVALUATE FORCE COMPONENTS ----SUM’OVER GRID POINTS, TO DO INTEGRAL ------- DO 170 I:1,3 DO 1172 IX=1,N1 DO 1172 IY=1,N2 P(IX,IY):0.0 DO 1172 IZ=1,N3-1 SX=SIDE4+FIDAT( IX- 1 +ITXS ) *DX SY:SIDE4+FDOAT(IY-1+ITYS)¥DY SZ=SIDE4+FLDAT( IZ-1+ITZS ) *DZ SZl=SIDE4+FIOAT(IZ+ITZS)*DZ ----- PARTICLE-CRIB DISTANCES-------- --- EVALUATE FORCES BY MINIMJN IMAGE METHOD-- XP=RX(J)+SIDE*NXP XXX=SX-XP IF(XXX.LT.-SIDEH) XXX=XXX+SIDE 61 31 41 1172 1173 170 210 IF(XXX.GT. SIDEH) XXX=XXX-SIDE YP=RY(J)+SIDE*NYP YYY=SY-YP IE(YYY.LT.-SIDEH) YYY=YYY+SIDE IF(YYY.GT. SIDEH) YYY=YYY-SIDE ZP=RZ(J)+SIDE*NZP ZZZ=SZ-ZP IF(ZZZ.LT.-SIDEH) ZZZ:ZZZ+SIDE IF(ZZZ.GT. SIDEH) ZZZ:ZZZ-SIDE ZZl=SZl-ZP IF(ZZl.LT.-SIDEH) ZZl=ZZl+SIDE IE(ZZl.GT. SIDEH) ZZl=ZZl-SIDE RIJSQ: (XXX*XXX+ YYY*YYY+ ZZZ*ZZZ)*RR*RR RIJ6=RIJSQ*RIJSQ*RIJSQ WIJzTEIMDtAA¥((4.*(1./RIJ6))¥(BB*(1./(CC+RIJ6))-1.) * +6.*BB*(1./((CC+RIJ6)*(CC+RIJ6)))) RIJSQ1= (XXX*XXX+ YYY*YYY+ ZZl*ZZl)*RR*RR RIJ61=RIJSQ1*RIJSQ1*RIJSQ1 WIJ1=TEIMD*AA*((4./RIJSI)*(BB/(CC+RIJ61)-1.) * +6.*BB/((CC+RIJ61)*(CC+RIJ61))) ----EVALUATE THE DIFFERENT FORCE COMPONENTS---- IF(I.EQ.1) GO TO 61 IF(I.ER.2) GO TO 31 IF(I.EQ.3) GO TO 41 P(IX,IY)=P(IX,IY)+DZ*(XXXtRR*WIJ*CONJG(W(IX,IY,IZ))*W(IX,IY,IZ) ¥ +XXX*RR*WIJ1*CONJG(W(IX,IY,IZ+1))*W(IX,IY,IZ+1))*0.5 GO TO 1172 P(IX,IY)=P(IX,IY)+DZ*(YYY*RR*WIJ*CONJG(W(IX,IY,IZ))*W(IX,IY,IZ) * +YYY*RR*WIJ1*CONJG(W(IX,IY,IZ+1))*W(IX,IY,IZ+1))*0.5 GO TO 1172 P(IX,IY)=P(IX,IY)+DZ*(ZZZ*RR*WIJ*CONJG(W(IX,IY,IZ))*W(IX,IY,IZ) t +ZZl*RR*WIJ1*CONJG(W(IX,IY,IZ+1))*W(IX,IY,IZ+1))*O.5 GO TO 1172 CONTINUE DO 1173 IX=1,N1 H(IX):0.0 D0 1173 IY=1,N2-1 H(IX)=H(IX)+DY*(P(IX,IY)+P(IX,IY+1))tO.5 CONTINUE SUM=0.0 D0 1174 K=1,N1-1 SUM=SUM+(H(K)+H(K+1) )xDXx0.5 CONTINUE IF(I.EQ.1) VEIX(J)=-SUM IF(I.EQ.2) VEIY(J):-SUM IF(I.EQ.3) VEIZ(J)=-SUM CONTINUE 169 17 211 CONTINUE IF((CIDCK/NW)*NW.EH.CIDCK) THEN RHO=FIDAT(N)/(SIDE**3) WRITE(11,*)CLOCK DO 17 K:1,IH R:(RADS(K)+RADS(K+1))tDIAM/(Z.) WRITE(11,*) R, G(K)/RHOW, GN(K)/RHO CONTINUE ENDIF RETURN END 212 cxxtxttxxxxtttxxtxxtxxxxxxxtxxtxxxtxttxxxxxxtxx*xxxx:x*xxxxxx PROGRAM TITLE : PROJ3D x PURPOSE : OBTAIN THE EIGENSTATES * METHOD : USE THE PROJECTION OPERATOR METHOD t .CONTAIN : CALL FOURN, CONTOUR * DOUBLE PRECISION 3 DATE : MAY, 1990 * INPUT : W(N1,N2.N3), VEI(N1,N2,N3), TINP * ttxxxxttittxttttttt¥*t***¥¥¥¥***¥t***¥**t**¥¥******ttttxittx 00000000 SUBROUTINE PROJ3D(W,VEI) IMPLICIT INTEGER*4(I-N) EKPLICIT REAL*8(A-H,O-Z) PARAMETER (NDIM=3, N1=16, N2=16, N3=16) DIMENSION NN(NDIM),VEI(N1,N2,N3),THETA1(N1,N2,N3) COMPLEXtIG W(N1,N2,N3),WO(N1,N2,N3),WNO(N1,N2,N3) ,WN1(N1,N2,N3),WTR(N1,N2,N3) ,FACT2(N1,N2,N3),EACT1(N1,N2,N3) ,CSUM,CSUM1 ,POT(N1,N2,N3),SW(N1,N2,N3) ,CW(N1,N2,N3),TW(N1,N2,N3) ,‘I‘O'I‘(N1,N2,N3) ,TE(N1 ,N2,N3) COMMON /A/ W ' COMMON /B/ SIDE,SIDEH,SIGMA,NCLOCK COMDN /C1/ DX,DY,DZ ****** OPEN(1,FILE = ’TINP’) READ(1,#)UI‘,NTIME,NCUX?K,NSTATE NN(1)=N1 NN(2)=N2 NN(3)=N3 NT=N1*N2*N3 NT2=2*NT XL=SIDE DX=SIDE/FDOAT( N 1 ) DY=S IDE/ FLOAT ( N2 ) DZ=SIDE/FDOAT(N3) PI:4.*DATAN(1.DO) PIZ=2.*PI DTH:DT/2.0 C ————— NORMALIZE INITIAL WAVEFUNCTION W(O)---- C ----- CREATE FACTI AND FACTZ CALL NORM(W,W,CSUM) DO 232 I=l,N1 DO 232 J=1,N2 DO 232 K=1,N3 W(I,J,K)=W(I,J,K)/SQRT(CSUM) WO(I,J,K)=W(I,J,K) 232 0 133 213 CONTINUE ------ BEGIN QCENCHING STEP ------------ DO 2 IT:1,NTIME ------ CALCUIATE ENERGY -------- DO 133 121,511 DO 133 J=1,N2 DO 133 K=1,N3 CW(I,J,K)=CONJG(W(I,J,K)) RIT(I,J,K)=CDNJG(W(I,J,K))*W(I,J,K)*VEI(I,J,K) CONTINUE ------ EFT --FRO'1 R SPACE TO K SPACE--------- CALL FOCRN(NN,NDIM,NT2,1) DO 3 I=l,N1 DO 3 J=1,N2 DO 3 K=1,N3 IX:N1/2-ABS(I-1-N1/2) IY=N2/2-ABS(J-l-N2/2) IZ=N3/2-ABS(K-1-N3/2) THETA1(I,J,K)=DT*4.¥PI*PI*FLOAT(IX*IX+IY*IY+IZ*IZ) /(2.*XL*XL) FACT1(I,J,K)=EXP(-THETA1(I,J,K)) TW(I,J,K)=W(I,J,K)¥THETA1(I,J,K)/DTH W(I,J,K)=W(I,J,K)*FACT1(I,J,K) CONTINUE ------ IFF'I‘ "FRO“! K SPCE TO R SPACE--------— CALL F®RN(NN,NDIM,N’1‘2,-l) DO 4 I=l,N1 DO 4 J=1,N2 DO 4 K=I,N3 'IHEI‘AzzDTWENLJm) FACI‘2(I,J,K)=EXP(-TI-IE1‘A2) SW(I,J,K)=W(I,J,K)iFACT2(I,J,K) W(I,J,K)='I‘W(I,J,K) CONTINUE ----- CALCULATE KINETIC ENERGY IF((IT/NCIDCK)¥NCUXZK.m.NCIOCK) THEN ------ IFFI‘ «PM K SPCE '10 R SPACE-------- CALL FOURN