A M O N T E C A R L O S T U D Y O F FLUIDS W I T H D E G R E E S O F O R I E N T A T I O N A L F R E E D O M By Mark James Blair B . Sc. (Chemical Physics) U . of Toronto, 1989 M . Sc. (Chemistry) U . of Toronto, 1991 A THESIS S U B M I T T E D IN PARTIAL F U L F I L L M E N T O F THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in T H E F A C U L T Y O F G R A D U A T E STUDIES (DEPARTMENT OF CHEMISTRY) We accept this thesis as conforming to the required standard T H E UNIVERSITY O F BRITISH C O L U M B I A July 1996 © Mark James Blair, 1996 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of m y department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Department of Chemistry) T h e University of British Columbia 2075 Wesbrook Place Vancouver, Canada V 6 T 1Z1 Date: Abstract We have modeled an electrorheological ( E R ) fluid as hard-sphere particles, each with smaller hard-sphere ions constrained to roll on the sphere's inside surface. W h e n the model E R fluid is placed in an electric field, each particle becomes polarized (due to rearrangement of the ions confined within the particle) with a dipole moment depending on both the field and interactions with its neighbors. Using N V T Monte Carlo simulations, we have shown that our model can display chain formation as seen in real E R fluids. Chaining occurred at a field where the net moment no longer varied linearly with the field. This model was extended to include particles shaped as ellipsoids of revolution. In the prolate case, slightly non-spherical particles were readily ordered by the field. In the oblate case, the induced dipole is roughly perpendicular to the symmetry axis. Oblate particles may then form a biaxial phase in an applied field. N P T and Gibbs ensemble Monte Carlo simulations were performed for spherical particles modified by an anisotropic potential of the form — 4 A £ ( < r / r ) P ( c o s 7 ) . 6 2 We have investigated both a soft-core model of Lennard-Jones particles and a hard-core model of hard spheres. For the soft-core model at anisotropy parameter A = 0 . 3 , we have con- structed the isotropic-nematic (IN) coexistence curve using Gibbs ensemble Monte Carlo simulations (the Gibbs method). For IN coexistence, we modified the Gibbs method so that particle exchanges were rotationally biased. For the hard-core model, we have determined the isotropic fluid-nematic coexistence curve using the Gibbs method. No gas-isotropic liquid transition was found for this hard-core model. Gibbs ensemble Monte Carlo simulations were performed for a mixture of neutral and dipolar hard spheres. A t fixed pressure, the coexistence curve for the demixing 11 transition was constructed for several values of the diameter of the neutral hard spheres. We extrapolate the critical point temperature to a vanishing diameter for the neutral hard spheres. In the limit of vanishing neutral hard spheres, the demixing transition of the mixture resembles, if it exists, the gas-isotropic liquid transition for dipolar hard spheres. O u r extrapolation suggests that the gas-isotropic liquid transition would occur at a very low temperature. 111 Table of Contents Abstract ii Table of Contents vi List of Tables vii List of Figures viii Table of Symbols x Acknowledgement xv 1 Introduction 1 2 C o m p u t e r Simulation 8 3 2.1 Statistical Mechanics and the Monte Carlo Method 10 2.2 Phase Transitions and the Gibbs Method 15 2.3 Pair Potentials and Boundary Conditions 22 2.4 Chemical Potential 30 2.5 Pressure 35 2.6 H a r d Spheres 40 2.7 Orientational Degrees of Freedom and Order Parameters 41 M o d e l Electrorheological F l u i d 46 3.1 46 Introduction iv 4 3.2 The Model 49 3.3 Computational Details 51 3.4 Determination of Ordering and Structure 53 3.5 M o d e l Building 55 3.6 Computational Tests 57 3.7 Results and Discussion 62 3.7.1 E R Behavior as a Function of Electric Field 62 3.7.2 E R Behavior as a Function of Particle Shape 75 3.7.3 Ordering of an E R F l u i d of Oblate Ellipsoids as a Function of Density 82 3.7.4 Summary 85 Liquid Crystal Models 87 4.1 Introduction 87 4.2 M o d e l and Computational Details 89 4.3 Results and Discussion 95 4.3.1 Results at A = 0.15 95 4.3.2 Results at A = 0.3 109 4.3.3 Results for the Hard-core Model 112 4.4 5 Summary 115 Mixtures of Neutral and Dipolar Hard Spheres 118 5.1 Introduction 118 5.2 M o d e l and Computational Details 120 5.3 Computational Tests 123 5.4 T h e Critical Dipole Moment as a Function of Pressure 126 5.5 T h e Critical Dipole Moment as a Function of Diameter of the Neutral H a r d Spheres 128 v 6 Conclusions 134 Bibliography 138 Appendices 143 A R o t a t i o n of a Particle 143 B C h o o s i n g a N e w Ion Position Inside an Ellipsoid 145 C R o t a t i o n of Charge Coordinates 149 vi L i s t o f Tables 3.1 Results for b/a = 1, p* = 0.75, q* = 0.2, and E* = 50 55 3.2 Results for b/a = 1, p* = 0.5, and E* = 10 56 3.3 Results for b/a = 1, p* = 0.5, and E* = 100 56 3.4 Results for b/a = 1, p* = 0.75, q* = 0.35, a* on = 0.05, and E* = 20. . . . 58 3.5 Results for b/a = 1, p* = 0.75, q* = 0.35, a* on = 0.05, and E* = 100 . . . 59 3.6 Results for b/a = 1, />* = 0.5, 3.7 Results for p* = 0.75, q* = 0.35, a* 3.8 Results for 6/a = 1, p* = 0.75, q* = 0.35, and <r* = 0.05 63 3.9 Results for b/a = 3, p* = 0.75, q* = 0.35, and a* 75 = 0.5, < on n = 0.1, and E* = 100 = 0.05, and E* = 40 on on 3.10 Results for p* = 0.75, q* - 0.35, a* on = 0.05 - 0.05, and E* = 40 62 77 3.11 Results for b/a = 0.333, q> = 0.35, a* 4.1 Previous results for soft model with A = 0.15 107 4.2 Gibbs results for hard-core model with A = 0.15 114 5.1 Gibbs results for a Lennard-Jones mixture at T* = 0.928, P* = 0.076. . . 124 5.2 Results obtained from starting with two different initial configurations for on = 0.05, and E* = 100 60 85 dipole mixture 125 5.3 Gibbs results at P* = 2.0, p* = 2.2 for dipole mixture 125 5.4 p* versus P* for dipole mixture. 127 5.5 Coexistence results at P* = 2.0 for dipole mixture 128 5.6 p* versus a* for dipole mixture 131 vii List of Figures 2.1 Phase diagram for a typical, one component (pure) system 3.1 Sample configurations for N = 128 and N = 120with b/a = 1, p* = 0.5, q* = 0.5, and a* ion 3.2 = 0.1 at E* = 100 61 T h e mean square displacement for b/a = 1, p* — 0.75, q* = 0.35, and < n 3.3 = 0-05 64 (a) (/j.*) and (b) (Pj) versus E* for b/a = 1, p* = 0.75, q* = 0.35, and < = 0-05 n 66 3.4 (a) g (z) 3.5 Dipole moment distributions for b/a = n C 3.6 and (b) g(r) for b/a = 1, p* = 0.75, q* = 0.35, and cr* = 0.05. . on 69 Sample configurations for b/a — 1, p* = 0.75, ( P ) versus 3.8 0H(z; u) for 6/a = 3, p* = 0.75, 3.9 Sample configurations for l/D a* = 0.05 on 70 E* for 6/a = 3, p* = 0.75, * = 0.35, and a* = 0.05 3.7 at (a) = 0.35, and E* = 50 and (b) E* = 60 2 72 on ? = 0.35, and a* = 0.05 73 on b/a = 3, p* = 0.75, 9 * = 0.35, and a* = 0.05 on E* = 0, (b) E* = 10, and (c) E* = 100 vs. £ * for (a) 2 74 b/a = 1 and (b) b/a = 3 at p* = 0.75, 4 * = 0.35, and = 0.05 < n 3.11 (a) p \ 68 1, p* = 0.75, q* = 0.35, and = 0-05 at (a) 3.10 16 76 (b) ( P ) , and (c) (cos 7) versus b/a for / 2 = 0.75, q* = 0.35, a* = 0.05 at £ * = 40 78 on 3.12 g\\(z] u) for p* = 0.75, q* = 0.35, < o r i viii = 0.05 at £ * = 40 80 3.13 Sample configuration for on E* = 40 at 3.14 p/pmax 3.15 b/a = 0.333, p* = 0.75, q* = 0.35, and a* = 0.05 (a) versus 81 E* for a single particle for q* = 0.35, and a* = 0.05. . . . on g\\(z;n) and (b) g\\(z;u) for b/a = 0.333, q* = 0.35, and a* = 0.05 at on E* = 100 84 4.1 (a) p* and (b) (P ) versus T * at P * = 0.1 for soft model with A = 0.15. 4.2 (a) 4.3 M e a n square displacements at P * = 0.1 for soft model with A = 0.15. 4.4 p* versus T* for N = (a) 32, (b) 64, (c) 108, and (d) 150 for soft model 2 g(r) and (b) g {r) 220 . at P * = 0.1 for soft model with A = 0.15 96 98 . . with A = 0.15 4.5 83 99 101 (P ) versus T* for N = (a) 32, (b) 64, (c) 108, and (d) 150 for soft model 2 with A = 0.15 102 4.6 u* and u* versus 1/N at T* — 0.7 and 0.9 for soft model with A = 0.15. . 4.7 Gas-liquid coexistence curve for soft model with A = 0.15 108 4.8 (a) p* and (b) (P ) versus T* for soft model with A = 0.3 110 4.9 Phase diagram for soft model with A = 0.3 Ill 2 104 4.10 Phase diagram for hard-core model with A = 0.15 113 5.1 Coexistence curves at P * = 2.0 129 5.2 Sample configurations from each Gibbs box for TV = 500 particles, a* = 0.8, x n = 0.8, P * = 2.0, and p* = 2.4 IX 132 Table of S y m b o l s Symbol Description a Ellipsoid of revolution degenerate axis. A Helmholtz free energy E — TS. b Ellipsoid of revolution unique axis. d Unit vector nematic director. E Total energy K + U. E External applied electric field for E R fluid. g( ) Radial distribution function. g (r) T e r m in expansion of full angle dependent distribution func- r 220 tion. 9\\( ) z Distribution function measuring ordering alonj1 induced dipole of spherical E R particle. 9\\( l ) z u Distribution function measuring ordering along symmetry axis ellipsoidal E R particle. 9\\{z\n) Distribution function measuring dipole of ellipsoidal E R particle. G Gibbs free energy h Planck's constant. k Boltzmann's constant. K Total kinetic energy. L Box length. m Particle mass. fiN. X ordering alonjI induced M Net moment, or total polarization, of the system. n Integer component lattice vector. n M a x i m u m length of vector N Number of particles. max n included in lattice vector sum. Ni Number of ions confined in each E R particle. P Pressure. P Polarization order parameter. P Nematic order parameter. on 1 2 P N 2 Nematic order parameter estimated from largest eigenvalue of ordering matrix Q. P (a;) Second order Legendre polynomial. q Ion charge. Q Classical partition function. Q Ordering matrix. Q Biaxial order parameter. r Particle position. 2 22 r Spherical cut-off separation for pair potential. nj Pair separation (under m i n i m u m image convention). n Ion position relative to center of E R particle. s Position or separation in terms of box length L. S Entropy. T Absolute temperature (in Kelvins). u Potential (one-body, two body, etc.). cut on u (12) Anisotropic contribution to pair potential. UDD Dipole-dipole pair potential. a xi UHS Hard sphere pair potential. U{ Coulombic (ion-ion) pair potential. ULJ Lennard-Jones pair potential. U Total potential energy. U Total field energy for E R fluid ( - M • E). Ui Total interaction energy for E R fluid (U — Us). ON F U cor Correction due to small system size to total potential energy for Lennard-Jones system. Us Total self energy for E R fluid. V Volume. Z Classical configuration integral. Z Classical configuration integral in scaled coordinates. a Convergence parameter in Ewald sum. s f3 1/kT 7 Shear strain. 7 Strain rate. 8 Shell size for calculating hard core contribution to the pressure. Ar Ar M a x i m u m translational step size along each coordinate axis. i o n M a x i m u m step size for hard-sphere ion contained within an E R particle. AU Difference in potential energy between a trial configuration and the previous configuration in the Markov chain. AU + Difference in potential energy due to trial insertion of a par- tide. . xn AUA->B Difference in potential energy due to converting a particle of type A to one of type B. AV M a x i m u m volume change. AO M a x i m u m angular step size. Aip M a x i m u m spin size about symmetry axis or induced dipole for E R particle, e Dielectric constant of the system. e Dielectric constant surrounding spherical assembly of periodic images. e 0 Permittivity of free space, e Energy parameter. TJ Packing fraction |/?*. TJ Shear viscosity. 7/0 Plastic viscosity. A Anisotropy strength parameter for liquid crystal model. A T h e r m a l wavelength. p Chemical potential. /j, Dipole (scalar quantity is the dipole moment). p Excess chemical potential p — pid- ex p Ideal chemical potential kT In pA . 3 id p M a x i m u m possible induced dipole in E R particle. p Residual chemical potential nd Limiting distribution for Markov chain. p Number density a Particle diameter. max res N/V. Xlll p— kTln(A/a) . 3 Gi Diameter of ion confined in E R particle, r Shear stress. on r 0 f2 Y i e l d stress. Particle orientation either in spherical polar coordinates or Euler angles (9, <f), ip). (...) Indicates an ensemble average. * Superscript indicates a reduced quantity. xiv (j>) Acknowledgement I thank Professor Gren Patey for providing both guidance and freedom throughout m y time at U B C . M a n y people have passed through our lab, all of whom have contributed to an enjoyable and stimulating work environment. In particular, I am grateful to John Shelley for the generous use of his scientific visualization program. I am also grateful to Gary A y t o n , who has provided computer code and assistance with the word processing of this thesis. I must also thank K u r t Ding, who undertook the unpleasant task of teaching me the art of science and scientific programming. Lastly, I thank my wife, Golnar Rastar, although her sole contribution has been to act as a distraction to the completion of this thesis. xv Chapter 1 Introduction T h i s thesis covers three loosely related projects. T h e first project was the modeling of, and investigation of the orientational ordering in, an electrorheological ( E R ) fluid. T h e second project was an investigation of the isotropic-nematic (IN) phase transition for a simple liquid crystal model. T h e final project was an investigation of the demixing phase transition for a mixture of neutral and dipolar hard spheres. Although the three types of system studied were all quite different, there is a unifying theme. Each project was an investigation of the phase transitions in systems with orientational degrees of freedom. Monte Carlo ( M C ) computer simulations were used to study each of the three types of system. Liquid crystals and E R fluids are both examples of systems with long-range orientational order. As such, they have many properties in common. Some interesting properties of E R fluids are closely related to the behavior of liquid crystals and we have explored this connection. T o date, there has been little interaction between the two fields of E R fluids and liquid crystals. In studying liquid crystals and E R fluids, we are dealing with complicated, many particle systems. Thus, we rely heavily on computer simulation for calculating their properties. It should be rioted that, for the most part, we do not discuss dynamical properties. Although the dynamical behavior of these systems, particularly the dynamical behavior of E R fluids, is interesting, we focus on the equilibrium behavior. T h e M C method is simply a way of estimating the multi-dimensional integrals that arise in classical statistical mechanics. Thus we may evaluate all the equilibrium properties of 1 Chapter 1. Introduction 2 our system of interest using the M C method. O f course, we cannot get any dynamical information from M C calculations. T h e truth of the matter is that liquid crystals and E R fluids are intrinsically interesting. Nevertheless, perhaps something should be said about the potential applications of these two systems. E R fluids could find application where materials with variable and controlled viscosities are required. E R fluids have possible engineering applications ranging from automotive devices to the stabilization of buildings during an earthquake [1]. Liquid crystals are used in flat screen display devices [2]. In the practical use of both E R fluids and liquid crystals, it is the switching times, and the required energy for switching, that are of crucial interest. Given the current trend towards engineered and advanced materials, the importance of understanding the properties of E R fluids and liquid crystals is clear. Electrorheological fluids consist of colloidal particles suspended in liquids, such as oils, of low dielectric constant. More precisely, an E R fluid is simply a colloidal suspension of polarizable particles in a non-polarizable solvent. Although the preparation of an E R fluid is simple, an E R fluid is a very complicated system. E R fluids display large changes in their rheological properties upon the application of an electric field [3]. Rheology relates the way a fluid deforms to the various forces which may be applied to that fluid. T h e fluid rheology, which may be Newtonian in absence of a field, can become distinctly non-Newtonian as the field strength is increased. For example, one often finds a field induced Bingham plastic [4,5]. Undoubtedly, the most dramatic behavior of an E R fluid is that its viscosity, or its resistance to flow, can be increased by orders of magnitude when placed in an electric field. Solidification can occur at fields above a critical value. It is as if water has turned to molasses in the presence of an electric field. This behavior is complicated because it involves an electric field, coupled with particle flow in a reasonably concentrated colloidal suspension. T h e reason for studying these systems is that we Chapter 1. Introduction 3 believe the essential features of E R fluids can be incorporated into a simple model. E R fluids were first discovered by Willis M . Winslow in the 1939. two patents based on E R fluids before publishing a paper on the Induced Suspensions He took out Fibration of in 1949 [6]. Winslow himself first realized the dependence of shear viscosity on the square of the electric field. W i t h such a fluid, one could, for example, control the spinning of a shaft surrounded by an E R fluid by the application of an electric field. Such a device, with a m i n i m u m of moving parts would prove very useful in the automotive and other industries. It is most likely that the patents did not lead to riches for Winslow, in that the E R fluids at the time did not have the structural strength to be of practical use. There has been resurgence of both practical [7] and theoretical [8-17] interest in E R fluids. F r o m experimental, theoretical, and computer simulation studies, it is now more or less accepted that the interesting behavior of E R fluids stems from the spatial ordering of field induced dipoles. A.system of permanent dipoles will spontaneously form a ferroelectric fluid phase at sufficiently low temperature [18,19]. A n E R fluid resembles such a system except that the E R particle has an induced dipole and an electric field induces the ordering to a ferroelectric fluid phase. It is generally agreed that in an E R fluid, the field polarizes the suspended particles. Thus each suspended particle acts as a dipole directed along the electric field. T h e dipoledipole interaction tends to align the particles in chains or fibers, again in the direction of the field. For such a simple model, the origin of the dramatic, and field dependent, viscosity increase is apparent. Particle flow involves either distorting or breaking the chains of suspended particles. Once flow has been established, it is not clear how much chain structure remains. One chapter of this thesis deals with the mechanism for orientational ordering in E R fluids. Our objective was to understand the structural properties of model systems by applying M C computer simulation methods. T h e polarization of an E R particle is due Chapter 1. Introduction 4 to the rearrangement of ions confined to the particle's surface. We have modeled the E R fluid as a collection of hard spheres or hard ellipsoids of revolution with small, mobile hard-sphere ions confined to their inside surface. It is possible that there are preferred sites for the ions on the surface of a real E R particle, however, any specific binding is left out of the model. Each hard E R particle is electrically neutral and moves in a continuum background. Each E R particle is readily polarized, and will build up a dipole moment to a m a x i m u m value of all positive charges on one end and all negative charges on the opposite end. As a function of applied electric field, the average dipole moment per E R particle obviously increases. Furthermore, the dipole in each E R particle aligns with the field, giving an orientationally ordered system. Unquestioningly, the orientation is due to the applied field. Nevertheless, the dipole-dipole interaction between the E R particles also enhances the orientation. component to the E R fluid. Non-spherical E R particles introduce a liquid crystal T h e dependence of the orientation of the E R particle's induced dipole on particle shape was also studied. T h e dipole ordering may be enhanced due to a larger induced dipole along the long axis of an ellipsoidal E R particle. This is a single particle, geometrical effect. T h e dipole ordering may also be enhanced by an orientational ordering due to the liquid crystal properties of the ellipsoidal E R particles. Another interesting property to consider is the total moment of the system. This is a collective response of the system to the applied field. Once each E R particle develops a sufficient dipole moment, the E R particles will spontaneously start to chain, acting to further enhance the total dipole moment of the entire system. This effect will be quite noticeable provided the dipole of the E R particles is still growing as a function of field when chaining occurs. In other words, the effect will be large if we are in the linear response regime. If chaining does not occur until the E R particle is almost completely polarized, chaining will not have much effect on the net moment. We now turn our attention to liquid crystals, a liquid-like phase of matter with Chapter 1. Introduction 5 anisotropic properties. As one might guess, some kind of anisotropy at the particle level is required for a liquid crystal phase. A l l real liquid crystals are composed of elongated or disk-like particles which interact via both short-ranged repulsions and longer-ranged attractions. Computer simulations have shown that sufficiently non-spherical hard particles, with no interactions other than that the particles cannot penetrate one another, can form liquid crystal phases [20-25]. Simulations have also shown that hard spheres (each embedded with a unit vector) interacting via a sufficiently anisotropic term in the pair potential may also form a liquid crystal phase [18,26]. W h y is there a liquid crystal phase in the first place? Consider heating a solid. T w o things will eventually happen. One is that the particles will start to spin and tumble and the other is that the lattice structure breaks down. There is no reason to assume that these two events happen simultaneously or even that one will always occur before the other. In the case where the lattice structure breaks first, we have a liquid crystal. A s we continue to heat the liquid crystal, the particles will start to spin in all directions. A t this temperature, we have a transition from a liquid crystal phase to an isotropic phase. In general, the liquid crystal may have various types of orientational order, leading to many different types of liquid crystals. T h e simplest type of liquid crystal is a nematic, in which there is no long-range translational order but orientational order in one direction. T h e direction of orientational order is indicated by a unit vector director d. There may also be translational order in one direction, in which case we have a smectic liquid crystal. A smectic phase is composed of two dimensional layers of liquids stacked on top of one another. plane. T h e orientational order in the liquid layer need not be perpendicular to the T h e angle of ordering leads to different classifications of smectic phases. The particles may be chiral, leading to chiral or twisted liquid crystal phases (also called cholesteric liquid crystals). T h e orientational order may be determined by the particle's dipole, and if the system has a net moment, we have a ferroelectric liquid crystal phase 6 Chapter 1. Introduction [18,19]. T h e ferroelectric liquid crystal phase is unusual in that the liquid crystal, in the absence of an electric field, generally adopts a structure minimizing the net moment. There are also discotic liquid crystal phases, formed from disk-like particles. Another topic discussed here deals with the IN equilibrium for simple models displaying liquid crystal behavior [27,28]. Methods for simulating liquids are readily extended to simulating liquid crystals. Thus the phase diagram for liquid crystal models may be mapped out in a straightforward manner. A n efficient and elegant way is to use the relatively new Gibbs ensemble Monte Carlo method (Gibbs method). T h e Gibbs method is an extension of the conventional M C method but designed to maintain thermal, mechanical, and material equilibrium between two systems. Although very successful when one of the systems is a gas and the other is a liquid, to date, no one has been able to use the Gibbs method when one of the systems is an isotropic liquid and the other is a nematic liquid crystal. This is accomplished in the present thesis. In the Gibbs method, material equilibrium is maintained by exchanging particles between the two systems. We have implemented a modification such that if one of the systems is orientationally ordered (as in a liquid crystal) then the exchanges are performed such as to preserve this order while still maintaining material equilibrium between the two systems. M u c h of our work has been employing the Gibbs method to simulate the IN coexistence. We have constructed much of the phase diagram for some simple models. Recent computer simulations suggest that no gas-isotropic liquid transition exists for dipolar hard spheres [29,30]. In other words, dipolar interactions alone are not sufficient to stabilize the liquid phase. These results are consistent with our results for a liquid crystal model of hard spheres embedded with an anisotropic interaction term. For this liquid crystal model, no gas-isotropic liquid transition was found. T h i s system, like dipolar hard spheres, only has a single isotropic phase. We were motivated by this finding to consider a mixture of neutral and dipolar hard spheres. T h e mixture undergoes a Chapter 1. 7 Introduction demixing phase transition, in which it separates into a neutral-rich phase and an isotropic, dipole-rich phase. T h e demixing transition provides another route by which to investigate the existence of a gas-isotropic phase for dipolar hard spheres. T h e demixing transition is similar to a gas-isotropic transition for dipolar hard spheres in a background of neutral hard spheres. In the limit that the diameter of the neutral hard spheres vanishes, the demixing transition, if it existed, would correspond to the gas-isotropic liquid transition. A mixture of neutral and dipolar hard spheres is a simple model for a polar substance in a non-polar solvent. A t a fixed pressure, the mixture undergoes a demixing phase transition as the temperature is lowered. Such a demixing transition was predicted by integral equation theory [31,32]. We have confirmed this prediction with our computer simulations. It is interesting that theory and simulation agree on the existence of the demixing transition because integral equation theory also predicts a gas-isotropic liquid transition for dipolar hard spheres [33,34]. Such a transition for dipolar hard spheres has not been found by computer simulation, despite several thorough searches [29,30]. There is the possibility that either the transition temperature is too low to achieve converged simulation results, or else that another phase transition pre-empts the gas-isotropic liquid transition. We have determined the critical temperature for the demixing transition as a function of the diameter of the neutral hard spheres. T h e demixing transition is expected to approach the gas-isotropic liquid transition for dipolar hard spheres as the diameter of the neutral hard-spheres vanishes. It is this possibility that we have explored here. Chapter 2 C o m p u t e r Simulation Computer simulation has become a powerful tool in the study of systems of interacting particles [35-38]. Two general approaches to the simulation of such systems have been developed. Molecular dynamics ( M D ) simulation is conceptually the simplest. One integrates the classical equations of motion for all the particles and physical quantities of interest are calculated as the time averages over a single trajectory through phase space. Monte Carlo simulation involves generating a distribution of systems which correspond to a desired statistical mechanical ensemble. Quantities are then calculated as averages over the ensemble of systems. M D and M C methods have developed in parallel, with advances in one often spurring advances in the other. For example, because of its close connection to statistical mechanics, methods for working in different ensembles (eg. N P T , uVT) were first developed for M C simulations. Soon after, pseudo-Hamiltonians (pseudo because Hamilton's equations are not satisfied) for M D simulations were proposed which gave the equations of motion corresponding to systems other than at constant N V E . Many of the classic papers on computer simulation are collected in the book of Liquids and Solids Simulation [36]. T h e two approaches, M D and M C , should be seen as complimentary to one another. M D gives information on the time dependence of the system. For example, one can calculate time correlation functions, from which one gets the transport coefficients, from M D simulation. In M C simulation, there is no real time, only a series of (highly correlated) configurations of the system. Since the M C method is not tied to the equations 8 Chapter 2. Computer Simulation 9 of motion, discontinuous moves may be made, which would be impossible using M D . As examples, particles may be created and destroyed, or entire groups of particles may be moved in a single step. This aspect of being able to destroy and create particles makes the M C method particularly suited to the study of phase transitions. T h e reason is as follows. T h e number of particles in a system is determined by its chemical potential. W h e n two phases are in coexistence, one condition is that each phase has the same chemical potential. Creating and destroying particles allows one to fix the chemical potential, facilitating the location of a phase transition. It should be noted that M D methods have been developed which continuously grow and shrink particles, thus varying the number of particles at fixed chemical potential [39,40]. Computer simulations are becoming ever more sophisticated. T h e rapid growth in the calculating speed of computers encourages longer, larger, and more complicated simulations. Developments in the methodology allows one to select the most convenient ensemble with which to work [36,41,42]. W i t h the development of non-equilibrium M D , one is not even restricted to equilibrium simulations [43]. Nor is one restricted to bulk properties. Surface properties [44,45], interfaces [46], and finite size effects [47] have all been studied with computer simulations. Q u a n t u m mechanical M C and M D calcula- tions have also been done [36]. T h e interaction potentials between particles can be made realistic enough to both reproduce and predict experimental results [38,48]. Special techniques have been developed so that many different systems, from polymers [48,49] to molten electrolytes [36] to water [36], are susceptible to computer simulation. This thesis is restricted to M C simulations in the N V T , N P T , and Gibbs ensembles of classical particles which interact via a pairwise additive potential. T w o of the systems under consideration are pure, one component systems and the third system is a mixture of two components. It is only the bulk, thermodynamic and structural properties that are of interest for the purposes of this thesis. Chapter 2. Computer Simulation 10 B y restricting ourselves to M C simulations, we have excluded the study of the dynamical properties of the systems under consideration. This may seem to be an unnecessary sacrifice, considering that one can determine both the equilibrium and dynamical properties from M D simulations. There is the question of what type of dynamics to perform, for example, molecular dynamics or Brownian dynamics [50], but the appropriate choice should be clear from the specific system of interest. Although one can calculate equilibrium properties from such simulations, the investigation of equilibrium properties need not be so straightforward. Using the assumption that time averages are equal to ensemble averages, we can calculate ensemble averages instead of integrating the equations of motion. This is the M C approach, and it also gives the equilibrium, but not dynamical, properties. However, for calculating ensemble averages, there are many tricks that one can use [51]. This is not the case for calculating time averages. 2.1 S t a t i s t i c a l M e c h a n i c s and the M o n t e C a r l o M e t h o d Classical equilibrium statistical mechanics can be expressed entirely in terms of multidimensional integrals and ratios of multidimensional integrals [51,52]. T h e Monte Carlo method is the fanciful name for the method of estimating integrals using pseudo-random numbers [53]. It is essentially the only method of estimating non-trivial high-dimensional integrals. T h e pseudo-random numbers are used to sample values in the space of the integral (the region over which the integration is performed). T h e sampling is then used to give estimates of the integral. It is easiest to first consider the Monte Carlo method for fixed number of particles TV, volume V, and temperature T. Thus it is the canonical ensemble that concerns us. Ideally it is the classical configuration integral Z that we wish to evaluate. T h e integral Chapter 2. Computer Simulation 11 is defined as (2.1.1) Z(N,V,T) = J...jexp(-/3U(r ))dr , N where r N N is an abbreviation for the positions 7*1, r , ...f;v of the TV particles, /? = 1/fcT 2 where k is Boltzmann's constant and U is the total potential (or configurational) energy. This is an integration over all TV particle coordinates which forms a 37V dimensional space called the configuration space. Each point in this space corresponds to a single configuration of the TV particles. Sampling from points in the integration space corresponds to sampling from random configurations. In practice, almost all configurations will make a negligible contribution to the integral, thus the estimate will be too uncertain to be of any use. Due to particle overlaps in a random configuration, the potential energy would be so high that the exponential term gives such a configuration a zero contribution to the integral. Only realistic configurations, those without significant overlap, contribute to the integral. Thus Z is a measure of the accessible configuration space. Due to the extreme sensitivity of U to the coordinates, the accessible configuration space is a complicated, high-dimensional shape. It is not possible to find the shape's hypervolume, thus directly evaluating Z is hopeless. Nevertheless, it is possible to estimate ensemble averages involving the ratio of integrals of the form f---fO(r )eM-PU(r ))dr f---fexp(-f3U{r ))dr N ( n s { = ' N N N ( 2 N ' 1 2 ] ' ' ' where the brackets denote an ensemble average and 0 is any configurational quantity. A configurational quantity is one that is explicitly a function of the particle positions. This excludes interesting statistical quantities such as the chemical potential. T h e ratio of two integrals is a more promising quantity to evaluate since we no longer need to know the total accessible configuration space. Instead, we need only to know the ratio of accessible configuration space for the two integrals. Estimating (O) using the Monte Carlo method with randomly chosen configurations Chapter 2. Computer Simulation 12 still fails due to bad sampling of the space of the integrals. Instead, we use importance sampling Monte Carlo which estimates (0) by sampling from a distribution of configu- rations that span the regions important to both the numerator and denominator. Such configurations are generated by a Markov chain. E a c h state in the Markov chain is a possible configuration. T h e chain is constructed using a transition probability matrix such that the chain approaches a unique limit distribution TTd- Metropolis and co-workers [54] chose TTd to correspond to the Boltzmann distribution such that 7r (r ) N d oc exp(-(3U(r )) N . (2.1.3) Other choices for ird are possible. Beyond leading to configurations which contribute to both integrals in the expression of (0), the unique limit distribution is arbitrary. One also has a choice of transition probability matrices which lead to the B o l t z m a n n distribution. Metropolis and co-workers take the probability of going from one configuration to the next as probability 1 if U decreases, exp(—/SAU) if U increases. This prescription has become known as the Metropolis Monte Carlo method. Indeed, it is also our method of choice. How does the Metropolis choice allow us to sample from the accessible volume of phase space? W h a t we want is a sampling distribution Td of configurations which is a small subset of all possible, random configurations. If we select a configuration from all possible configurations with a probability proportional to e _ / ? c / , then we will only allow realistic configurations into our sampling distribution. We are now estimating the ensemble average (O) on a distribution TTd which contains only those configurations which make a significant contribution to the integral. This is a far, far smaller distribution than that of all random configurations. Now the problem lies in generating this distribution of configurations. This is where the power of computers is exploited. In practice, the Markov chain is started from some initial configuration of N particles Chapter 2. Computer Simulation 13 in a cubic box of volume V. T h e initial configuration may be either a lattice configuration, random configuration, or a configuration from a previous simulation. A particle is then chosen at random and we attempt to move it to a new position randomly chosen within a cube with length 2Ar (where Ar is the m a x i m u m translational step size along each coordinate axis) centered on the old position. If the new position is outside of the central cell, it is translated by the cell length back inside the cell. T h e change in potential energy AU is then calculated. Following the Metropolis choice of the transition probability matrix, the move is accepted if the energy is decreased. If the energy is increased, the move is accepted with probability exp(—/3AU), otherwise it is rejected. Whether the move is accepted or rejected, the configuration becomes the next step in the Markov chain. This is repeated until long after the system reaches equilibrium. Conventional wisdom says that accepting approximately half the trial moves after reaching equilibrium is a reasonable compromise between adequate sampling and computational time. We claim equilibrium is reached when the energy appears to be fluctuating about a constant value. There is always the possibility that the system has become trapped in a metastable region of configuration space. We report when we suspect that this has happened. We have described the Metropolis Monte Carlo method for the canonical ( N V T ) ensemble. T h e description for the isothermal-isobaric ( N P T ) ensemble differs only slightly. During an N P T simulation, one allows volume fluctuations such that the ensemble average of the pressure equals the desired pressure P. In the N P T ensemble, the limiting distribution is ir {r ,V) N d oc exp [-B(U+ oc V e-W ^ N +pv PV - NkTInV)] (2.1.4) Chapter 2. Computer Simulation 14 and the probability in going from one configuration to the next is m i n ( l , exp(Aiu)), where (2.1.5) and A V is a value randomly chosen from the interval [—AV , A V max m a a ; ] . Recall that for N V T simulations, Aw = —/3AU. In practice, we randomly decide for each step whether to attempt a particle translation or a volume change. T h e probability is set such that a reasonable fraction (1/10 to 1/4) of the total number of steps are attempted volume changes. accepted. AV max is adjusted so that roughly half of the attempted volume changes are W h e n the volume change is accepted, the volume changes by AV. Often, the particle coordinates are scaled by the box volume and the scaling properties of the potential energy are exploited so that attempted volume changes are not computationally expensive. After equilibrium, we begin to collect averages of various mechanical quantities using each configuration of the Markov chain. O f course, each step is highly correlated with the previous step so that we must average over a large number of steps to get a single independent value. Before a simulation is started, we anticipate a reasonable number of steps, or block size, required to give independent averages. During the simulation, the Markov chain is typically broken up into n& equal length blocks with averages collected over each block. T h e standard error is estimated in the usual way by considering the block average values {Oi} as independent observations from a parent distribution of width (standard deviation) UQ- T h e estimate of the mean is (0) (2.1.6) T h e standard deviation of the mean is a. (2.1.7a) Chapter 2. Computer Simulation 15 where a = J^-rKO ) 2 0 V n - 1 - (0?\ • (2.1.7b) b T h e standard deviation o~o is somewhat independent of rib, but the standard deviation of the mean a m becomes smaller for longer simulations (rib increasing with the length of each block remaining fixed). According to statistics, we may be confident that the mean lies within a m from our estimate roughly 65% of the time. We choose to report (0) ± 3cr . • We may be confident that the mean lies within 3<7 roughly 99% of the time. m m Reporting the larger error has two advantages. T h e first is that if values disagree by more than the reported error, we may be reasonably certain that the difference is meaningful. T h e second advantage is that in practice, all simulations were for rib — 10 blocks, so that ao ~ 3<r . m One should check that the block size is larger than the correlation length of the Markov chain. Usually, correlations are readily seen as a systematic drift about the equilibrium value in the block average values. T h e Metropolis method is constantly driving the energy to its equilibrium value. Thus one could imagine that correlations in other quantities could persist well after the energies had become uncorrelated. Thus one should check for correlations in all quantities of interest, not only in the energy. 2.2 Phase Transitions and the Gibbs Method T h e equilibrium phases for a pure substance may be represented on a P vs. T diagram or a T vs. p diagram. A simple and typical example of each type is shown in F i g . 2.1. T h e phase diagram may be more complicated, with many solid phases (such as the many solid phases of water) and unusual liquid crystal phases or superfluid phases. vs. T diagram is made up of a set of coexistence curves. The P A n y point in the P vs. T diagram not on a coexistence curve represents a single phase of the system. In other Figure 2.1: Phase diagram for a typical, one component (pure) system represented in a (a) P vs. T diagram and a (b) T vs. p diagram. Chapter 2. Computer Simulation 17 words, there is only one thermodynamically stable phase of the system at that particular pressure and temperature. Along the coexistence curve, there are two phases which are thermodynamically stable. Three phases may exist at a point called the triple point, but certainly no more than three phases may coexist for a pure substance. A point on the coexistence curve locates the pressure and temperature at which the phase transition occurs. For a first order phase transition, there is a density difference between the two phases. T h e density difference implies a discontinuity in the derivative of the Gibbs free energy G with respect to pressure, since (2.2.1) A second order transition is one in which the discontinuity occurs in the next higher order derivative of the Gibbs free energy. A density difference implies that no single phase exists at an intermediate density. T h i s corresponds to the two phase region in the T vs. p phase diagram. We now consider how to determine the phase diagram by computer simulation. T h i s requires the location of phase transitions. If an N V T simulation is started at a temperature and density in the two phase region of the T vs. p phase diagram, one might expect the system to reach the thermodynamically stable state of two phases in coexistence. T h i s does not happen [55]. Instead, the system remains in a single phase. T h e reason is that the energy for setting up the interface region between two phases in a small system is prohibitively high. Relative to Avogadro's number, all computer simulations are of small systems. A more promising approach to locating a phase transition would be to perform N P T simulations. in section 2.1. Recall that the method for performing N P T simulations was discussed This corresponds more closely to a search through the P vs. T diagram, rather than through the T vs. p phase diagram. phase As the coexistence curve Chapter 2. Computer Simulation 18 is crossed by varying either T or P, one would expect the system to undergo a phase transition. Indeed, this does occur, but not necessarily at the true coexisting temperature and pressure. Consider varying the pressure at fixed temperature. remain mechanically stable, (dV/dP)x < 0 , beyond T h e system may where the phase transition should have occured. In this case, the system is in a metastable state. Examples of metastable states are super-cooled and super-heated liquids. There are other, indirect ways of locating phase transitions from computer simulations. One brute force method is to perform thermodynamic integration using the equation of state (pressure as a function of temperature and density) generated from simulation [51]. For example, by integrating the thermodynamic relation (2.2.2) one can determine the free energy for a system in a particular phase. One then determines where the free energies of two phases are equal, thus locating the phase transition. One clever way of locating a phase transition is using Monte Carlo scaling techniques, which give results over a range of temperatures, densities, etc. from a single simulation [51]. There is a simulation method which allows for the simulation of two phases in equilibrium with one another. Thus it is a wonderful way for mapping out phase diagrams. T h e three conditions for equilibrium between two phases are equal temperatures, equal pressures, and equal chemical potentials. T h e Gibbs ensemble Monte Carlo method is designed to establish these three conditions between the two simulation boxes. T h e Gibbs method, developed in 1987, is a relatively new simulation technique [56]. T h e Gibbs method is essentially a canonical ensemble ( N V T ) Monte Carlo simulation of a system divided into two subsystems, or boxes, with volume and particle exchanges allowed between the boxes. T h e advantage of the Gibbs method is that, starting with a total density (total number of particles divided by the total volume of the two boxes) corresponding Chapter 2. Computer Simulation 19 to a two phase system, one should eventually get a different phase is each of the two boxes. T h e phases will be in equilibrium with one another. T h e Gibbs method reduces to a simulation of two equivalent N V T systems for an initial density corresponding to a single phase system. T h e Gibbs method is now considered in more detail. Equilibrium is established by the following three types of moves. Particles are translated, just as in an N V T Monte Carlo simulation. Volume exchanges are also performed between the two boxes such that the total volume remains fixed. This is like a conventional N P T Monte Carlo move, except that the boxes are coupled. This establishes mechanical equilibrium, thus equating the pressures. Finally, particle exchanges are performed between the boxes. Again, this is like a conventional ^ V T Monte Carlo move except that the boxes are coupled. This establishes material equilibrium, thus equating the chemical potentials. T h e original version of the Gibbs method used macroscopic, thermodynamic arguments to derive the probability of accepting each type of move. T h e Gibbs method was put on a firmer footing using statistical mechanical arguments [57]. T h e partition function, appropriate to the ensemble from which the Gibbs method samples, was determined. Smit et al. [58] then showed that the Gibbs and canonical ensembles are equivalent in the thermodynamic limit. T h e y derive the acceptance probabilities starting from the partition function for the Gibbs ensemble. For the canonical ensemble, the Boltzmann factor is exp(—/3U). For the Gibbs ensemble, Smit and Frenkel derive the pseudo-Boltzmann factor [59] exp(zw) = exp In { WIN"] where the superscripts / and / / label each of the two simulation boxes. T h e probability for the acceptance of any type of move is m i n ( l , exp(Aiw)). For example, for a volume Chapter 2. Computer Simulation exchange of AV A „ = N * l n 20 between boxes / and 77, the expression for Aw is (£+£V) + n N l n ( Y ! L ^ _ ? A V . _ . For particle exchanges (say for exchanging a particle from box / to box II), for (2.2.4) the expression Aw is / Aw = In i yll\ N [j^yr) ~ P&U 1 - (3AU» . (2.2.5) Note that unlikely particle removals may be made overall favorable if insertion into the other box is favorable. In Panagiotopoulos's original formulation of the Gibbs method [56], he proposed a sequential choice of moves. First iV.translations, then one volume change and n h exc par- ticle exchanges were attempted. It was found that if more than about 3% if the particles were actually exchanged, the system was pushed so far out of equilibrium that the Gibbs method did not give reliable results. Smit and Frenkel [59] recommend performing the moves in random order. T h e y claim that this makes the Gibbs method more reliable as well as giving lower uncertainties in the chemical potential. Thus, we have used this random order selection of move type. Panagiotopoulos has subsequently introduced a fourth type of move, applicable for simulation of mixtures [60]. It involves changing an A particle to a B particle in one box and simultaneously performing the reverse move in the other box. T h e resulting Aw for the identity exchange move (converting an A to a B particle in box I and a B into an A particle in box II) A w is = (/vJTl) ln One of the two boxes (/ or II) + ln (TVFTt) ~ ^ U *- B ~ • (2 - " 2 6) is chosen with equal probability. T h e type (A or B) of particle is then chosen with equal probability. A particle of the chosen type in the chosen box is then randomly selected. A particle of the other type in the other box is also Chapter 2. Computer Simulation 21 randomly selected. There are other algorithms for selecting two particles to participate in the identity exchange. T h e important point is that the selection algorithm must satisfy the condition of microscopic reversibility [61]. T h e effect of the identity exchange move is to ensure that the difference in chemical potentials between the two components is the same in each box. As such, along with identity exchanges, it is still necessary to perform particle exchanges for one component. T h e advantage is that one is free to choose the component for which to perform particle exchanges. T h e sensible choice is to pick the component for which particle exchanges are easiest. T h e Gibbs method for a mixture differs from that for a pure system in another way. In principle, there is no reason why a Gibbs ensemble simulation cannot continue even if one of the boxes has emptied [58]. One includes the contribution of the empty box to the accumulating average quantities and continues the simulation. In practice, one usually restarts the simulation with a different total volume. For mixtures, it may be quite common for one of the boxes to be devoid of one component. In this case, the simulation really is continued. For example, we may attempt to remove a particle of a type which is not in the selected box. In this trial move as an attempted, but rejected, and the averages are updated accordingly. It is necessary to know how many quantities may be fixed in a computer simulation. T h e answer is supplied by the Gibbs phase rule (which actually was, unlike the Gibbs ensemble, developed by Josiah W i l l a r d Gibbs). T h e phase rule says that the number of thermodynamic degrees of freedom / (independent intensive thermodynamic variables) is given by / = 2+ n where n c c - n is the number of components and n p p , (2.2.7) is the number of coexisting phases at Chapter 2. Computer Simulation 22 equilibrium. For a pure system with two phases in coexistence, there is only one thermodynamic degree of freedom. In practice, this is almost always taken to be the temperature. Thus, when we map out a coexistence curve using the Gibbs method, the coexistence values are measured at a series of temperatures. For simulating a single phase of a pure system, we are free to specify both temperature T and pressure P. ulation is perfectly valid, even for a pure system. Thus an N P T sim- However, one may not fix both T and P for a Gibbs ensemble simulation for coexistence, as this violates the phase rule. Nor can one perform a pPT simulation for a pure system for the same reason. If one attempted to do so, the particle number TV would not have an equilibrium value, and would drift aimlessly. Note that the phase rule does permit a Gibbs ensemble simulation of coexistence of a binary mixture at fixed T and P. Instead of performing particle exchanges for both components, we may perform particle exchanges for one component along with identity exchanges. For an identity exchange, we convert a neutral into a dipole in one box and at the same time convert a dipole into a neutral in the other box. It is easier to exchange a neutral particle, since one must only find room for it in the receiving box. Exchanging a dipole usually involves an energy penalty for removing a dipole from an energetically favorable position and inserting it into an energetically unfavorable position. T h e benefit of exchanging neutral spheres only becomes better as the diameter of the neutral spheres shrinks. 2.3 Pair Potentials and B o u n d a r y Conditions Consider a system of TV interacting particles. T h e potential energy may be written as the sum over individual particles, pair interactions, triplet interactions etc. to give U(r ) = J2n(r ) N l i + J2<r ,r )+ l i<j J £ i<j<k u(r , r r ) . . . % i? k . (2.3.1) Chapter 2. Computer Simulation 23 For the purposes of our computer simulations, this expansion was truncated after the second sum. There are cases where the expansion of the potential energy in such a form is not useful because the many-body interactions make a large contribution [62]. A n extreme example is superfluid liquid helium. Even for less exotic systems, the threebody sum may make a contribution of a few percent to the potential energy. Nevertheless, we felt justified in making this truncation since much of the behavior of the real systems is captured by a pairwise interaction. T h e first sum accounts for the effects of applied fields. T h e second sum is a sum of the pair potential for the interacting particles. the pair potential may be written as For a translationally invariant system, where = Tj — V{. T h e subscript ij is sometimes dropped for convenience where it is understood that we mean the pair separation. T h e pair potential may be spherically symmetric, in which case, it depends only on the magnitude r. Perhaps the two most famous examples are the hard-sphere potential, t UHs(r) = i 0 r<a oo r >a , (2.3.2) where a is the hard-sphere diameter, and the Lennard-Jones potential, CT\ 12 ULj(r) = 4e la (2.3.3) . r, where — e is the m i n i m u m potential energy and a is the separation at which the potential energy vanishes. As for the hard-sphere system, a may be considered to be the particle "diameter". T h e pair potential may depend on the orientation of the particles. In this case, we use the notation u(12) where 12 = r, J? and fl represents the Euler angles 2 describing the orientation of a particle. W h e n a translational move takes a particle out of the simulation cell, the particle is translated back by the cell length so that it remains within the cell. This assures that our Chapter 2. Computer Simulation 24 system has a well-defined number of particles and volume for any particular configuration. However, it is not sufficient to remove surface effects. A particle in the corner of a simulation cell is clearly not in a bulk phase environment. T h e simplest way remove surface effects is to impose m i n i m u m image (MI) boundary conditions. T o impose M I boundary conditions, one calculates the pair separations using the M I convention. Under the M I convention, the separation between particle i and j is the separation between particle i and the nearest image of particle j . T h e nearest image of particle j may be the particle itself within the simulation cell, or its image in the periodically replicated cubic simulation cell. Thus, the maximum value of r,j is \J\L where L = V / is the 1 3 simulation box length. For short-range potentials, those that decay faster than r - 3 as r —» oo, simulation has proved remarkably successful. For such systems, the bulk properties of the infinite system are readily inferred from simulations of a small number of particles with M I boundary conditions. This is successful because the potentials can be truncated beyond either a sphere or box of reasonable size (typically several particle diameters). Corrections may be applied to account for the small system size. For the Lennard-Jones potential with spherical cut-off at r where p = N/V c u t , the correction to the total potential energy is [37] is the number density. T h e corrections to the potential energy are applied during the simulation, so that quantities calculated from the potential energy do not need to be adjusted later on. For example, both the pressure and chemical potential are related to the potential energy. T h e correction assumes that particles are uncorrelated beyond a separation of rcut. Thus, the correction is not valid for a solid phase or for too small. O f course, in practice, r cut must not exceed half the box length. ensembles, the box size varies, so it is most convenient to always set r cut r cut For some to half the box 25 Chapter 2. Computer Simulation length. In this case, r cut varies with box size. Determining the bulk properties of an infinite system involving long-range potentials from a small number of particles is not so easy. It is the potential energy, the sum of all the pair potential energies, which appears in the statistical mechanical integrals. Its correct evaluation is crucial to produce results appropriate to the thermodynamic limit. T h e evaluation may be carried out as done for short-ranged potentials. For short-range potentials, m i n i m u m image (MI) boundary conditions are imposed. T h e problem with long-range potentials is that the m i n i m u m image sum for the potential energy may not well represent the potential energy for the small system immersed in the bulk system. Various boundary conditions have been proposed to address this problem. T h e long-range potentials used in this thesis are the coulombic (ion-ion) potential u (r) = ^ im , (2.3.5) where q is the ion charge, and the dipolar (dipole-dipole) potential u (12) = DD ^ , (2.3.6) where / / is the dipole. T h e scalar quantity is the dipole moment p. In rationalized M K S units (SI units), these two potentials should include the factor l / 4 7 r e , where eo is the 0 permittivity of free space. T h e two potentials should also include a factor 1/e, where e is the dielectric constant (relative permittivity of the medium). These two factors are readily absorbed into the definition of the reduced charge or reduced dipole moment, so they are dropped for clarity. First consider a net charge neutral system of ions; T h e m i n i m u m image (MI) expression for the total potential energy is U = Y,— , (2-3.7) Chapter 2. Computer Simulation 26 where r;j is the separation between particle i and the nearest image, not necessarily within the simulation cell, of particle j . It is possible to take into account the interaction of a particle with all its neighbors, including their images as well as with its own images. T h e system is contained in a cube which is periodically replicated in a space-filling way to yield U = IJ2 E E , m ' , • (2-3.8) T h e prime on the sum indicates that terms for i = j are not included when n = o. Terms for i = j correspond to a particle interacting with its own image. T h e lattice sum is a sum over integer component vectors n — n, n) y z , where n, n, n x y z = 0, ± 1 , ± 2 , . . . . (2.3.9) T h e lattice sum is only conditionally convergent, meaning that the result depends on the order of summation. T h e technique for evaluating the lattice sum was originally developed by Ewald [63]. Unfortunately, the problems associated with the conditional convergence were not fully appreciated. In fact, there are several ways [64] of getting the wrong answer! T o obtain a unique limit for a conditionally convergent sum, the method of summation (the order in which the sum is performed) must be defined. T h e total potential energy of the central cube may be written as a sum over approximately spherical shells of cubic images. T h e limit as the radius of the spherical array becomes arbitrarily large is then taken. This introduces a positive term proportional to the square of the net dipole moment which depends on the dielectric constant e' of the continuum surrounding the arbitrarily large array. One must define a value of e' in order to apply periodic boundary conditions ( P B C ) . For a conducting surrounding continuum (e' —> oo) the extra term vanishes. T h e e' —> oo case is sometimes referred to as "tin foil" boundary conditions. Simulations employing the Chapter 2. Computer Simulation 27 original Ewald method implicitly set e' to infinity. T h e opposite limit is a surrounding vacuum (e' = 1), where the extra term makes its m a x i m u m contribution. conducting and vacuum limits of P B C are used. B o t h the Although avoiding the problems of m i n i m u m image, periodic boundary conditions are also not ideal as they replace the true, isotropic potential with an anisotropic, effective potential. Nevertheless, the full Ewald expression is ^ + ^ + rijto e n j -v n /a 2 « J2 / i 2 \ J i ' \ 2 5 — l5Zftexp(27rm-s,-)| 2 i 7 r n /(OyEfc-.-] } 2 , (2-3.10) where positions and separations are in units of the box length (s,- = r , - / L and S y = rij/L). T h e complimentary error function is defined as erfc(x) = - = / e~* dt . (2.3.11) T h e function f(e') depends on the dielectric constant of the continuum surrounding the system of periodic images. For this case, it is f(e') = — ^ — . (2.3.12) T h e potential energy has been expressed as the sum of two converging series. T h e convergence rate of each series is controlled by a. In practice, a is chosen so that the contributions from the complementary error function terms vanish (within the m a x i m u m allowed error) beyond half the box length. T h e m a x i m u m length n max of the lattice vectors determines the number of terms in the reciprocal space sum. Effectively, this determines the radius of the spherical assembly of images. n max is taken sufficiently large Chapter 2. Computer Simulation 28 so that the truncation error for both sums match each other. A s a final simplification, consider the sum over n. equals the — n th T h e terms have the symmetry property that the n term th term. As a consequence of this symmetry property, we may write where the star on the sum indicates that only vectors pointing up are included. T h e final term in Equation (2.3.10) is proportional to the net moment, M or total polarization, of the system. = ( - 2 t 3 1 4 ) A computationally efficient form for the potential energy with periodic boundary conditions is Tit >\ 1 I eiic(asij) rtmax * + + 2 2 T , — | V] qi e x p ( 2 7 r m • S j ) £ a Typically, a value of a ~ -7r n /a 2 e P-1 f ( ^ M \ (2.3.15) 5 is sufficient so that only the first term i n the real space (complementary error function) sum need be included. O f course this implies including many more terms in the reciprocal space sum. A value of n 2 max = 14, implying a sum over 124 lattice vectors, allows for the calculation of the potential energy to sufficient accuracy that any calculated system property remains unchanged (within the statistical uncertainty) with the inclusion of more lattice vectors. T h e Ewald expression is readily generalized to a rectangular box. shortest side is taken as the box length L. In this case, the T h e box dimensions L , L , L x y z are given in Chapter 2. Computer Simulation terms of L so that L\ = L{/L 29 where i = x,y,z. A t least one of L' ,L', X L' z will be 1, and all three will be one for a cubic simulation cell. T h e lattice sum is no longer a sum over integer component vectors but a sum over vectors v defined as The Ewald expression for a rectangular box is then TT, a j 1 erfc(o;3 ) tj [i<j h J Si Y + v -k g — m a x ir v /a 2 2 2 —T~ TT7 7 I2>exp(27rii/-a,-)| — 2 27T + m ^ M . 2 (2.3.17) A few comments are in order about the scaled coordinates. B o t h in computer simulation and theoretical derivations (see, for example, section 2.4), it is convenient to work with scaled coordinates s = r/V / . 1 3 For the most part, we employ a cubic simulation cell, in which case, L = V / and the scaled coordinates s — rjL 1 3 as before. A rectangular box has V / / 1 3 rectangular box, we have chosen to use have the same meaning L, where L is the shortest box length. For a s = r/L, as opposed to s = r/V / . 1 3 O f course, this makes no difference to the final results. The periodic summation may also be performed for the dipolar pair potential. T h e computationally efficient Ewald expression for dipoles (in a cubic box) is [37] [i<j L Umax * ^ + iJ e — K J2 n#o ^ U S 2 n S 2 2 n / a 2 I - ) exp(27rm • s ) \ n 2 { i E ^ + f ^ M 2 , (2.3.18a) Chapter 2. Computer Simulation 30 where B(s) = erfc(a.s) + 2as (2.3.18b) and 4(as) C{s) = B(s) + — 3 (2.3.18c) —6 T h e remaining quantities all have the same meaning as in Equation (2.3.15) for ions. 2.4 Chemical Potential Recall that three conditions must be satisfied for equilibrium between two phases of a system. T h e three conditions are equal temperatures, equal pressures and equal chem- ical potentials (for each component in a mixture) in each phase. T h e Gibbs method establishes these three conditions without the need for knowing the actual values of the pressure or chemical potential. Nevertheless, it is both useful and instructive to calculate the chemical potential. Most people would claim to have an intuitive sense of what is meant by temperature and pressure. W h a t about the chemical potential? Computer simulations, by forcing us to see how quantities are actually calculated, lead to physical insight into the actual meaning of these quantities. T o prove this point, there is no better example than the chemical potential. T h e computational method of particle insertion leads to insights into the meaning of the chemical potential that would be difficult, if not impossible, to get from any current textbook. T h e chemical potential is a less familiar quantity than the temperature or pressure, so it deserves some discussion. Consider the partition function Q(N, V, T) for the canonical ensemble for a system with only translational degrees of freedom. In the classical limit, the partition function is (2.4.1) Chapter 2. Computer Simulation 31 T h e total energy is E(r ,p ) N = K( ) + U(r ), N N N (2.4.2a) P where K(P ) (2-4.2b) = 7^EP1 N is the total kinetic energy and m is the particle mass. T h e integrations over the canonical momenta p N are readily performed to yield We define the thermal wavelength as and the classical configuration integral as Z(N, V,T) = J...J We may work with scaled coordinates s = V e-WUr" _ 1 . / r , such that dr = Vds, 3 (2.4.5) so that the limits of integration no longer depend on volume to obtain yN ^ ~ N\A 3NZS = QZ id where Q u , a (2.4.6) is the ideal gas partition function and Z = J ...je~^ ^ds sN N s We note that Z s . (2.4.7) still depends on the volume through the total potential energy. For an ideal gas, the potential energy U is zero, thus Z = 1. s Chapter 2. Computer Simulation 32 T h e appropriate thermodynamic potential for the canonical ensemble is the Helmholtz free energy A — E — TS. T h e connection to statistical mechanics is through A(N,V,T) = -kTlnQ . (2.4.8) F r o m Equation (2.4.6), we may write the free energy as A = A + A ld , ex (2.4.9a) where — kT = N ln A /V + ln NI . 3 (2.4.9b) Using Stirling's approximation (TV! m N ln N — AC), we find - ^ - = l NkT 9A -i . 3 n / y y (2.4.10) ! F r o m the first and second laws, we have dE = TdS - PdV + fidN . (2.4.11) It follows that dA = -SdT - PdV + udN , (2.4.12) anc Thus the ideal gas contribution to the chemical potential is given by k T = l n p A We now look at the excess contribution. Since , u v / N,T 3 . (2.4.14) Chapter 2. Computer Simulation 33 converting to a partial derivative with respect to density gives, It follows that dp ( A^T")N,T {pkr) (2.4.17) p' As a function of density, the excess contribution to the free energy may be expressed as A {p) A(p) - ex NkT ld NkT A (po) NkT ex In the limit po —* 0, A (p ) ex A (p) j" (JP_ _ J [p'kT + /' P0 ) d£ ) p' 1I ^ • (2-4.18) also goes to 0 and we assume the integral is well-behaved to 0 obtain Differentiating with respect to p (which is equivalent to differentiating with respect to N since V is held fixed) gives H^-HtsHf '- (2 4 20) We now derive Widom's [65] expression for the chemical potential. T h e chemical potential is = A{N + 1,V,T)-A(N,V,T) = _jfcrin%±I, (2.4.21) QN where we assume N is sufficiently large that we may treat it as a continuous variable. The ratio is QN (iV + l ) A 3 /.../(-«'"«.» ' 1 ' Chapter 2. Computer Simulation 34 thus giving T h e subscript TV indicates that it is an ensemble average over the N particle system. AU+ is the interaction potential energy of the (N + l) th particle with the rest. T h e integral in brackets may be evaluated using brute force (no importance sampling) Monte Carlo. Thus we choose positions in the unit volume at random and calculate AU+. For sufficiently large TV, we may replace N + 1 with JV, so that the first term on the righthand-side is pid/kT. T h e ensemble average in the second term on the right-hand-side is an average of the brute force M C evaluations of the integral. One may simply think of it as an ensemble average of the quantity e x=u l i d _ / 3 A L / + . Thus we get the expression -kT\n(e-^) , N (2.4.24) due to W i d o m [65], for the chemical potential. This is the particle insertion method. u We may define the residual chemical potential H = kTln(^j res +u res through [59] , (2.4.25) giving M r „ = - f c r i n ( | ^ I e - ^ ) . (2.4.26) We have introduced the residual chemical potential because it is this quantity (not p ) ex which must be equal for two phases in equilibrium. T h e particle diameter cr has been included to keep the arguments of the logarithms dimensionless. V* = V/cr is the 3 reduced volume. Note that V* and N are included within the ensemble average. In the Gibbs ensemble, both V and N fluctuate. In this case, the ensemble average is no longer for fixed TV, although AU particle. is of course still the energy for inserting an (TV + l) th + ghost Chapter 2. Computer Simulation 35 We should make one final note about calculating the chemical potential from computer simulations. For systems of moderate densities and beyond, the vast majority of the computational time is spent determining whether or not the inserted particle overlaps with particles in the system. One must write the algorithm for determining overlaps with this in mind. For example, we check for overlaps by checking the components of the separation vector between a particle and the inserted particle. If any component is larger than the largest particle diameter, the particles do not overlap, and we do not need to test the remaining components. This way, we only explicitly test for overlaps with the immediate neighbors of the inserted particle. 2.5 Pressure In computer simulation, one typically fixes the temperature through the Boltzmann factor. T h e calculation of the chemical potential was discussed in section 2.4. We now turn our attention to the calculation of the pressure. In computer simulation, one may fix the pressure, in which case it is unnecessary to explicitly calculate the pressure. Recall that for the Gibbs method, one does not need the actual value of the pressure. Nevertheless, it is an excellent check on the simulation results to compute the pressure. Like the chemical potential, the pressure may be related to the Helmholtz free energy, (2.5.1) Recall that A = E — TS, so that the pressure is Also recall that E = K + U, and that the kinetic energy is independent of volume. We also make use of the Maxwell relation (2.5.3) Chapter 2. Computer Simulation 36 such that This equation is known as the thermodynamic equation of state. T h e pressure may be broken into P = Pid + P , ex where the ideal pressure is PM = NkT/V and the excess pressure is '- (£L-(£Lr T h e pair potential may include both a continuous part and a discontinuous, hard-core term. T h e excess pressure may be broken into P = P ,hard ex ex + Pex,soft these two contributions to the total pressure. We first consider P ,soft, ex to account for which is identified as Pex,sofi ex,so ft = - [ % ) • (2-5-6) T h e potential energy often has the form f/oc^n • (2-5.7) We refer to this as a soft potential since it is continuous, unlike a hard-core potential. In a cubic box of length L with volume V = L , one may write 3 The quantity being summed no longer depends on the volume. A t equilibrium, the ensemble average of the potential energy may be identified with the thermodynamic potential energy which appears in Equation (2.5.6). Noting that d _ L d dV ~ WdL ' ( 2 - 5 9 ) the excess pressure for the soft potential is given by Pex,soft = ^ • (2.5.10) Chapter 2. Computer Simulation 37 It may happen that the potential energy is a sum of terms, each which scales with a different power n. For example, the Lennard-Jones potential contains a repulsive term with n = 12 and an attractive term with n = 6. Provided one calculates the contribution to the energy from each term in the potential separately, the contribution to the pressure is readily calculated. A further advantage to calculating the excess pressure from the potential energy is that any corrections due to small system size need only be applied to the potential energy. Corrections to the pressure are then automatically included through Equation (2.5.10). T h e expression for the excess pressure was derived assuming a cubic box shape. In fact, the expression is valid for an arbitrarily shaped volume. As mentioned, the potential may also include a discontinuous, hard-core term. T h e hard-core contribution to the pressure may be identified with the first term on the righthand-side of Equation (2.5.5), however, this is of little use for explicitly evaluating P ,hardex T h e hard-core contribution to the pressure is most easily calculated from the virial equation (which in turn is derived from the thermodynamic equation of state) — = l-- -[g(r)r(^)dV, (2.5.11) £ NkT 6kT J y y ' \dr) K ! where g(r) is the radial distribution function. T h e radial distribution function is defined by 9(r) > ( - 2 5 1 2 ) where 6(r) is the Dirac delta function. T h e radial distribution function is the ratio of the angle-averaged local density about a particle to the bulk density. T h e first term on the right-hand-side of Equation (2.5.11) is the ideal gas contribution. We are only interested in the excess contribution to the pressure from the hard core. Assuming a spherically symmetric, hard-sphere potential, the excess contribution to the potential is Pe*Mrd = ~ P 2 J9(r) r ' ^ d r . (2.5.13) Chapter 2. Computer Simulation 38 T h e derivative in the integral is awkward, so one typically converts it to a derivative of a step function, giving Pe*Mrd = Y^^I 9{r) {jr ~ ePUHS e )* PUHS r dr • -- (2 5 14) T h e derivative of the step function gives a delta function, thus PesMrd = Y P 9(°)° kT 2 • 3 (2-5-15) To calculate the hard-core contribution to the pressure, one accumulates g{r) during the simulation. Once the simulation is finished, g(r) is extrapolated back to the contact value at r = cr. It would be convenient if the pressure could be calculated without waiting until the end of the simulation. Indeed, one can define a contact function F which permits the calculation of the pressure during the simulation. If we count the number of separations N s which are within a shell of width ACT about each hard sphere, we get an estimate for g(cr). T h e estimate is g(a) « T h e factor 2/(N ° N 2 A s . — 1) accounts for the fact that we are using all N(N of the N particle system. N s (2.5.16) — l ) / 2 separations is estimated using the contact function F where N s = J2F(i,j). (2.5.17) T h e prime on the sum indicates that only terms for pair separation r < o + A<r are included. For a hard-sphere system, the contact function is trivially given by 2 F(z,j)=[-) . (2.5.18) T h e m a x i m u m value for F is then Fmax ~ 1 + — • (2-5.19) 39 Chapter 2. Computer Simulation In general, the contact function has the property that t WW < 1 if 1 and 2 overlap = 1 if 1 and 2 are in contact > 1 if 1 and 2 do not overlap (PV/NkT) T h e expression for the hard-sphere compressibility PV NkT = 1 + 7^7 ^ Some workers define 8 = 2Aa/a 3N l i m A<T-O Ao- and sum over 1 < A typical value is Aa/a is (2.5.21) -z-Y.F(hJ) f£ y ' J ' F <1+ T h e problem has been reduced from extrapolating g(r) reasonable value for Aa/a. (2.5.20a) 8. to contact to choosing a 0.01. T h e contact function method for calculating the pressure is surprisingly insensitive to the value up to moderate reduced densities (up to p* ?s 0.5). A t higher reduced density, one must reach a compromise between the number of counted separations and the uncertainty. T h e contact function method has the advantage that there is no need to accumulate g(r) if one only wants the pressure. However, the true power of the contact function method becomes apparent when calculating the pressure for hard, non-spherical particles. In this thesis, we consider hard ellipsoids of revolution. T h e expression for the compressibility is still given by Equation (2.5.21), however, the contact function F is now much more complicated [66]. W i t h no hard-core interaction, one typically reports the reduced pressure W i t h a hard-core interaction, one typically reports the reduced pressure P* = Pa /e. 3 P* — Pa /kT. 3 Note that the reduced pressure depends on a. Thus if we arbitrarily measured the particle diameter as 2a, the reduced pressure would change. T h e dependence of P* on the choice of length scale a means that one must be careful in defining P * for a mixture. In this thesis, we consider binary mixtures with spherical hard cores. Labeling the diameters of Chapter 2. Computer Simulation 40 the components as a A and O \ B , we arbitrarily select a A as the unit of length. Thus the hard-core contribution to the reduced pressure P * = Po~A/kT is z= P* ex,hard where /9* = —O* o > /0<r^, 2 O ^ B x _|_ (-^L) V crA A = (cr^ + / xx A B 9AB{O-AB) + (—'] W A / 2 B Ni/N (i = A, B) and £ ; = 0\B)/2 x g B(o-B) B (2.5.22) are the mole fractions. As for the pure case, the radial distribution functions at contact may be estimated using the contact function method. 2.6 H a r d Spheres A system of pure hard spheres provides a useful reference or test system. There has been a great deal of computational and theoretical work on hard spheres. One of the earli- est important computer simulation results showed that hard spheres undergo a freezing transition [36]. There is no gas-liquid transition for hard spheres. M u c h of the theoretical work on hard spheres can be encompassed in the Carnahan-Starling equation of state for the hard-sphere fluid [67]. This equation accurately describes the hard-sphere pressure up to the density of the freezing transition at p* ~ 0.95, although the equation does not predict the freezing transition. T h e virial (or pressure) equation is op oo ^ = l + £ i W T > " . P (2.6.1) n=l Carnahan and Starling observed a pattern in the virial coefficients for a hard-sphere system. T h e y were then able to perform the summation, giving BP 1+ P 7] +n - n 2 3 (i-v) (2.6.2) ' 3 where n — \p*• From Equation (2.4.19), one can show that A ex NkT~ r/(4 - 3r/) (1-n) 2 • -- (2 6 3) Chapter 2. Computer Simulation To this, we add 41 — 1 to get the excess chemical potential. After some algebra, the result is Equations (2.6.2) and (2.6.4) provide useful checks for the computer simulation of hard spheres. 2.7 Orientational Degrees of Freedom and Order Parameters The expressions for the chemical potential may be readily generalized to include rotational degrees of freedom. A particle may have two rotational degrees of freedom, determined by the polar angles 9 and <j> or three, determined by the Euler angles </>, and ip. In either case, we denote the angular degrees of freedom by Q. Integration over the angular degrees of freedom gives the angular volumes W 47T = < for 2 angles 8TT for 3 angles . 2 As was the case for the translational conjugate momenta, we can readily integrate over the rotational conjugate momenta. T h e classical partition function for a particle with two rotational degrees of freedom is «=m^(£)"^/-""'"^ - --> e , oW 2( 7al where Q = h /87r Ik . 2 rot 2 (2.7.1b) A particle with two angular degrees of freedom has two degenerate moments of inertia /. A particle with three angular degrees of freedom may have three different moments of inertia. This leads to three different 0 <, labeled QA, 0 B , and 0 c . For the partition ro Chapter 2. Computer Simulation 42 function for particles with three angular degrees of freedom, one replaces (T/Q ) ROT N in Equation (2.7.1a) with \e e &cJ A (2.7.2) ' B In either case, Widom's expression for the chemical potential may be written as p = Hid,r ~ *(^ kTl 0t / e-^ds dn ^j N N , where the effects of the rotational momenta are absorbed into (2.7.3) pid,rot- T h e expression for the residual chemical potential is identical to that for a system with no orientational degrees of freedom, Pr„ = - * T l n ( ^ ^ c - ^ ) . (2.7.4) T h e factor of W is canceled out when the integral in Equation (2.7.3) is estimated by inserting randomly oriented ghost particles at random positions in the system. Order parameters are quantities which change from zero to a non-zero value for a system passing through a phase transition. In this thesis, orientational transitions are of interest. Such transitions are characterized by simple, scalar order parameters. general, the order parameter may not be so simple. In For exotic phase transitions, the order parameter may be a more complicated quantity, such as a tensor. For a dipolar system, a P i order parameter is defined as A = II>,-|/I>. i (2-7-5) i T h e quantity calculated in simulations is the ensemble average of P (Pi) = (M)lN(p) . l 5 (2.7.6) P i measures the degree to which the dipoles are pointing in the same direction. Essentially, P i is the normalized polarization per particle. Chapter 2. Computer Simulation 43 For a system of uniaxial particles (two rotational degrees of freedom), the P 2 order parameter is defined as ^ = ^E^(3cos ^-l), (2.7.7) 2 where # is the angle between the unit vector iii along the symmetry axis of the 8 i th particle and the director d (defined below). In practice, however, the order parameter is calculated using the ordering matrix [20]. T h e ordering matrix is Q = ^E^(3AA-J). (2.7.8) It has three eigenvalues, in decreasing order, A , A , A _ . + 0 T h e director d is defined as the eigenvector associated with the largest eigenvalue A+. Various combinations of the eigenvalues may be used as an estimate for the order parameter P . In the isotropic 2 phase, P is zero. Due to finite size effects, the estimate of P 2 2 computer simulation is not quite zero. in the isotropic phase from Rather, the eigenvalues in the isotropic phase have the N dependence = Bv) • A+ ° = (if) ' 0 A 0 M D - = (TF) ' x 0 '' <2 7 9) In the nematic phase, the eigenvalues have the N dependence K = P 2 + (Ly 0 A o = _lp 2 + 0 ( - ^ ) , and A_=-I f t + 0 (-i=) . (2.7.10) The best estimate of P 2 in the isotropic phase is —2Ao, while the best estimate in the nematic phase is A . One must be consistent in the choice of estimate. We have elected + to use P 2 = —2A , as it is easiest to distinguish between an isotropic and weakly nematic 0 phase when the estimate for P 2 in the isotropic phase is as close to zero as possible. Nevertheless, we have calculated both quantities to facilitate comparison with published results since both choices appear in the literature [37]. T h e difference between the two Chapter 2. Computer Simulation 44 estimates is small, but large enough to be annoying for the purpose of quantitative comparisons. T h e notation P N 2 results using A + is used to indicate when we compare with published as an estimate of the order parameter. Again, we note that it is the ensemble average (P2), calculated from the instantaneous values of P2, which is of interest. Particles may be orientationally ordered in more than one direction. we must measure the biaxial order parameter Ql [23]. 2 In this case, In general, the orientation of a biaxial particle is described by the Euler angles 0,<f>,ip. T h e biaxial order parameter is then defined as [23] '(1 + cos di) 2 Ql\2 — 2 cos 2<j>i cos 2ipi — cos 9{ sin 2(j){ sin 2ipi (2.7.11) This is not a computationally useful form because the Euler angles are not measured with respect to a lab-fixed frame, but measured in a frame given by unit vectors X, Z (as measured in the lab-fixed frame). Consider the unit vectors ii , ii , ii x y z Y, pointing along the axes of the biaxial particle. We may then construct three ordering matrices Q X \ Q , YY Q Z \ where Q XX with similar expressions for Q VV ^E^(3«f«f = and Q . - I) , (2-7-12) T h e eigenvectors corresponding to the largest ZZ eigenvalues of each of the three matrices then define X, Y, Z. O u r situation of a uniaxial particle in a field allows for considerable simplification. Since the field is in the z direction, Z is then the lab-frame z axis. T h e unit vector ii = (u , u , u ) x y z for the particle is then rotated into the 2 = 0 plane, giving a vector if = . \l l u = T h e vector u v (U ,Uy,0) 1 X +y u « , < , 0 ) - is readily constructed by the requirement that u ii = «,-<,o) y (2.7.13) x • ii y = 0. T h e result is Chapter 2. Computer Simulation 45 0) • T h e matrices Q , XX Q and the corresponding eigenvectors X, Y ave readily constructed, YY from which one can calculate Q\ 2 Q\ = | ( X • Q XX 2 (2.7.14) using the expression •X + Y •Q YY •Y - X •Q YY •X - Y •Q XX • Y) . (2.7.15) Chapter 3 M o d e l Electrorheological F l u i d 3.1 Introduction Rheology relates the way a fluid deforms or flows under the various forces which may be applied to that fluid. Electrorheological ( E R ) fluids are those whose rheology varies dramatically with applied electric field. Although the field dependence of the fluid properties can be complicated, the most dramatic effect is the rapid increase in viscosity with applied field [7]. In some cases, a transition from a fluid to a solid is observed. A n E R fluid is prepared from a suspension of polarizable particles in a solvent of low dielectric constant. A n E R particle is mesoscopic in size and the polarization is due to the rear- rangement of ions confined to the particle's surface. In an electric field, the particles form chains which may then aggregate into columns. T h e viscosity depends on the degree of chain formation which, in turn, depends on the electric field. Thus the viscosity of the material may be varied and controlled with an electric field. The microscopic description of the chain formation is as follows. As a function of applied electric field, the average induced dipole moment per particle will increase. T h e dipole in each particle tends to align with the field, giving an orientationally ordered system. T h e dipole-dipole interactions between the particles further enhances the orientational order. T h e dipolar forces between particles then act to form chains along the field direction. These chains can extend across the entire fluid. This chaining process happens on a time scale of a few milliseconds. Chain formation alone is enough to cause 46 47 Chapter 3. Model Electrorheological Fluid the transition to a solid. O n a much longer time scale, the chains aggregate into columns of staggered chains. Rheology relates the shear stress r to the shear strain 7. For example, for a Newtonian fluid, r = 7 7 7 , where 7 is the time derivative of the shear strain and is called the strain rate. In fact, for a Newtonian fluid, this defines the shear viscosity 77 (usually just called the viscosity). For a non-Newtonian fluid, r = 77(7)7, which defines the apparent viscosity 77(7). One class of non-Newtonian fluids is a plastic (also called viscoplastic or B i n g h a m plastic). T h e plastic is described by the equation T = T + 77o7 , (3.1.1) 0 for T greater than the yield stress T , where 770 is the plastic viscosity. 0 A n E R fluid is observed to behave as a B i n g h a m plastic in an electric field such that T ( E , 7 ) = T ( E ) + 77O7, (3.1.2) 0 for shearing perpendicular to the field direction. It is typically found that 770 does not depend on field, T (0) = 0 (Newtonian fluid in zero field) and T (E) O 0 OC E N where n ^ 2 ( E R effect). T h e enhanced apparent viscosity is defined as 7 7 ( E , 7 ) - 77(0,7) = ^ Showing that this quantity behaves as E 2 7 • (3.1.3) is a demonstration of the E R effect. Previous models for an E R fluid begin at the level of point polarizable particles and did not explicitly include the ions which give rise to the induced dipole. For example, Halsey and Toor have modeled an E R fluid using polarizable spheres [8]. Tao and Sun have also used this model to determine the ground state of an E R fluid [9,10,68]. T h e y found the lowest energy to be a tetragonal-I lattice. This result has been verified experimentally by laser diffraction [69]. Tao has also simulated a low density lattice and, as a function of Chapter 3. Model Electrorheological Fluid 48 the electric field, found a transition to a chained fluid, and then to an unchained fluid as the field was decreased [70]. T h e liquid phase of an E R fluid has been studied within the mean spherical approximation, and the enhancement of the dipole moment due to manybody effects was investigated [71]. T h e dynamics have also been studied using Stokesian dynamics with electrostatics [11], Brownian dynamics [12], and molecular dynamics with hydrodynamic resistance [13]. Others have investigated the particle-particle interactions of E R fluids using various expansions of the potential [11,14,15]. It was found that the behavior was sensitive to the charge distribution within a particle at small particle separations. In particular, the short-ranged interaction was significantly underestimated by including only point dipoles. Zhang and W i d o m have studied the force between particles for the magnetic analogue of an E R fluid [72]. T h e y compared the force between isolated particles with the force between particles within a chain. It was noted that the chain structure can roughly double the force between dipoles. Some consideration has been given to the distortion of particle shape due to the applied field [7,16,72]. Particles in a typical E R fluid act as spherical droplets. W h e n polarized by a uniform applied field, such particles will distort to a prolate ellipsoidal shape. T h e distortion is understood by the competition between surface tension and field energy. T h e surface tension is sufficiently strong that the elongation is slight. We consider a model E R fluid in which the ions are treated explicitly. A simple model of an E R fluid is that of a system of large spheres, each with a number of small hard-sphere ions free to roll along the sphere's inner surface. T h e net charge on each large sphere is zero but each sphere will have its own field-induced dipole moment. T h e advantage of treating the ions explicitly is that the details of the particle-particle interaction are accounted for. For example, polarization and multipolar effects are naturally included in a realistic way. We discuss the parameters for the model for which we expect to see the E R effect. Although we cannot simulate as many ions per E R particle as in a Chapter 3. Model Electrorheological Fluid 49 real E R fluid, much of the behavior of real E R fluids is reproduced. We also examine the influence of particle shape on the electrorheological effect. Instead of considering small deformations in shape due to the field, we consider particles initially constructed to be non-spherical. B o t h oblate and prolate ellipsoids of significant eccentricity are investigated. For prolate particles, the symmetry axis and the induced dipole are strongly coupled such that the field orders even slightly prolate particles. Also, the field required to order prolate particles decreases as the particles become more and more elongated. In contrast, oblate particles order with their symmetry axis perpendicular to the field. We show that, unlike prolate particles, oblate particles may then form a biaxial phase. 3.2 The Model We consider a general ellipsoidal model. A n ellipsoidof revolution is characterized by its length to breadth ratio b/a, where 2b is the length of the ellipsoid along its symmetry axis and 2a is the length of an axis perpendicular to the symmetry axis. length, a, is defined through a 3 T h e unit of = 8a b. Physically, this is equivalent to taking the unit 2 of length as the diameter of a sphere of equal volume to the ellipsoid. T h e unit of energy, ekT, is the product of the dielectric constant e of the solvent, Boltzmann's constant and temperature T. T h e dielectric constant e of the solvent is taken to be unity. k, As mentioned, an E R fluid is a colloidal suspension in which the suspended particles have a dielectric constant sufficiently higher than the suspending medium. O u r model does not explicitly contain this feature. In fact, the dielectric constants of both materials are taken to be the same. However, the main effect of the difference in dielectric constant is to keep the ions confined to the suspended particles. Thus our model implicitly accounts for the difference in dielectric constant by confining the ions inside the E R particles. There are N E R particles and each particle contains Ni on ions. There are an equal number Chapter 3. Model Electrorheological Fluid 50 of positive and negative ions so that each E R particle is neutral. To characterize the system for fixed 7Y and i V , - , the reduced density p* = Ncr /V, 3 the reduced ion charge the reduced ion diameter o*on = <Ti /a, and the reduced electric on q* = q/^4:Tre ekTa, 0 on field E* = ^4:Tre ea /kTE 3 0 We must be set. wish to make a connection between the reduced parameters employed here and typical experimental values. Typical E R particle diameters range from lpm [11]. Since this covers quite a large range, for now we take a — dpm, to 100/mi where d is a variable parameter. A t room temperature ( T = 300K), the reduced fundamental charge is q* ~ 0 . 2 3 6 J / . A typical electric field at which the E R effect is observed is l k V / m m . -1 2 In reduced units, this field gives E* w 164d / . T h e reduced field is sensitive to the particular value of d and can range from E* « 164 for d = 1 to E* « 3 2 164,000 for J = 100. We typically report reduced particle energies, i.e., in units of ekT. particle energy u* for an infinite chain of point dipoles is [10] u* « pIV'4irtoekT'a 3 — T h e reduced per 2Ap* , where p* = 2 is the reduced dipole moment. For chain formation to be energetically favorable, the electrostatic interactions must overcome the thermal motions. This occurs for \u*\ > 1, which implies p* > 0.65. Although we are not dealing with point dipoles, there must be enough ions per E R particle to form reduced dipole moments of this order of magnitude. T h e m a x i m u m reduced dipole moment, p m a x , is approximately 9 * o t o ,/*/2, where q*otai = Nionq* is the total of the absolute value of all the charges contained by an E R particle and /* is the maximum reduced separation of the positive and negative ions within the particle. Thus we have bY 1/3 max U ~ 2^ i* ion( where the quantity in square brackets is /*. p* max as Nlonq*/2. (3.2.1) max For a spherical particle with a*on <C 1, This implies that we would need Nion « 54 for q* = 0.024 ( E R Chapter 3. Model Electrorheological Fluid 51 particles with diameter of lpm) for the dipolar energy to just match the thermal energy. T h e dimensionless parameter p* has been used to characterize E R fluids [8,71]. For 2 p* 1, chaining behavior is expected. 2 Unfortunately, with present day state-of-the-art practical limit to values of / V , o n computer workstations,, a realistic is about ten. T o make chain formation energetically favorable, we should have q* > 2p* /Ni . max For 7V on ton = 8, this requires q* > 0.16. Consider splitting a pair of ions using an electric field. Assuming cr* <C 1, the total on gain in reduced energy is q*E* whereas the loss due to the lower coulombic interaction is 1* / *on2 a Thus the field required to break up the pair is roughly given by equating these two energies, which implies E* « q*/u* . cr* is a particularly awkward parameter to on on choose. It must be small enough that the ions are relatively mobile. In other words, most of the ion moves should not be rejected due to overlap with other ions i n the same E R particle. Nevertheless, a* on should not be so small that the unlike ions collapse into pairs. In such a situation, the E R effect would be severely limited. Some further restrictions on <r* on 3.3 are given in Appendix B . Computational Details Periodic boundary conditions were imposed on the simulation cell and the Ewald summation method [63,64] was employed to calculate the total potential energy. Recall that the Ewald summation method was discussed in section 2.3. T h e total energy U was broken up into a self energy Us and an interaction energy Ui, as defined below, and the field energy UF- T h e total potential energy of the system of ions confined within the E R particles was calculated using the Ewald summation method. then defined as k i<3 i,j€k T h e self energy Us was Chapter 3. Model Electrorheological Fluid where u,j = qiqj/rij 52 is the usual coulombic pair potential. labels an E R particle, and i and j label ions within the k th T h e summation index k E R particle. We then use Ui — U — Us — UF to determine the interaction energy. T h e dipole of the k E R particle is th /•** = £ ?fr,-. (3.3.2) T h e field energy is then U =- £ F k E-ii k = -E • £ n k = -E • M . (3.3.3) k Note that for a polarizable system, the net lowering in energy due to the field is [73] (assuming the net moment M grows linearly with E). the result is — E • M. —\E-M For a non-polarizable system, We have a polarizable system, thus our field energy UF is not the net lowering in energy due to the field. Rather, half of the energy UF goes into pulling the charges apart, thus increasing the energy by approximately \E in a net lowering in energy of approximately — \E • M. •M and resulting T h e energy changes due to the charge separation are accounted for in Us and U\. In the present simulations, three types of Monte Carlo moves were performed. For each move, we either translate an E R particle, rotate an E R particle or roll an ion within an E R particle. Details of how to rotate a particle are given in Appendix A . Details of how to select a trial ion position within an ellipsoidal E R particle are given in Appendix B . In principle, rotational moves were not necessary for spherical E R particles as rolling the ions has the effect of reorienting the induced dipole direction. Nevertheless, rotational moves, which collectively move all the ions in a single E R particle, speed convergence to equilibrium. Details on how to rotate the ions along with the E R particle are given in Appendix C . T o assist the approach to equilibrium, attempts were made to transfer an ion to a position on the opposite side of the E R particle as well as rolling the ion. T h e Chapter 3. Model Electrorheological Fluid 53 trial position of the ion was reflected through the center of the E R particle. Such moves were attempted with a probability somewhere between 0 to 0.1 every time an attempt was made to roll a single ion. O n average, a M C sweep was composed of TV translations, 2JV rotations (TV rotations for spheres) and JViV,- on ion rolls. For non-spherical particles, rotations were performed about both the symmetry axis and the induced dipole. E a c h time, the type of move to perform was randomly selected with the appropriate probability. T h e step sizes were adjusted to give an acceptance ratio of ~ 0.5 for each of the three types of moves. R u n lengths were from 2,000 to 10,000 M C sweeps depending on how long it took for quantities to fluctuate about an equilibrium value. 3.4 D e t e r m i n a t i o n of O r d e r i n g and S t r u c t u r e As the E R fluid becomes polarized, each E R particle acquires an induced dipole p. T h e vector sum of all TV dipoles gives the net moment M. In section 2.7, we defined the (M)/N(p), normalized polarization per particle as (Pi) = where (p) is the ensemble average dipole moment of a particle. (Pi) acts as an order parameter for the E R fluid. A value of (Pi) close to 1 indicates that the system is polarized, however, it does not necessarily indicate the formation of chains. T h e orientational order of the symmetry axes is described by the (P ) order parameter [20]. Recall from section 2.7 that it is 2 defined as (P ) = 2 (i(3cos 0-1)} 2 , (3.4.1) where cos 9 is the dot product of the director d with the unit vector it along the symmetry axis of the ellipsoid. Except for finite size effects, (P ) is zero in an isotropic phase but 2 increases towards unity for an orientationally ordered phase. T h e distribution function, <7||(.z) was used to check for chain formation [74]. g\\(z) is the probability of finding a particle a distance \z\ away from a particle oriented along the Chapter 3. Model Electrorheological Fluid 54 z axis at the origin and within a cylinder of radius a/2 about the z axis. In a chained system g\\(z) shows peaks roughly where z is a multiple of the particle axis forming the chains. g\\(z) is given by 9 " { Z ) = 2 {N-l) ^a/2Y ' P ( 3 A 2 a ) where ij = z 6(x) •r\ {j , r\ 3 = \r i3 - z^e \ { , (3.4.2b) is the Dirac delta function, 9(x) is the Heaviside step function, and e; is the unit vector along some direction associated with the i th particle. To see chain formation along the direction of the field, the obvious direction for e; is along the induced dipole, although for non-spherical particles g\\(z) is also calculated using the symmetry axis. To distinguish between the two for ellipsoidal particles, we use the notation g^(z;u) and g\\(z; p) to indicate that the function has been calculated using the symmetry axis or the induced dipole direction, respectively. In addition to these characterizations of ordering, the usual quantities such as the radial distribution function g{r) and the mean square displacement (\r(t) — r(0)\ ) (where 2 time t is not a real time, but the number of M C sweeps in the simulation) were also accumulated. Such information gives some indication as to the liquid or solid nature of the system. Such plots are usually reported for molecular dynamics simulations where M C sweeps corresponds to time. One then calculates the diffusion coefficient from the slope of such a plot. Although a M C sweep does not correspond to a physically meaningful interval of time, such plots are useful for identifying a solid phase. Note that if the particles are executing a random walk, then the mean square displacement gives a straight line of finite slope. For a random walk, the quantity (\r(s) — r(0)j ) is proportional to the number of 2 steps s in the walk. Although the average reduced dipole moment (p*) is reported, this does not give any indication about the distribution of dipole moments. A distribution Chapter 3. Model Electrorheological Fluid 55 sharply peaked about the average would indicate that the E R fluid is well-represented by particles all of the same dipole moment. For this reason, the dipole moment distribution f(fJ^) was also accumulated. 3.5 Model Building T h e present E R models are characterized by a large number of parameters. T h e de- pendence of the models on these parameters has been investigated. We have tested the dependence of the results on the ion diameter while keeping all other values fixed (see Table 3.1). T h e results are weakly dependent on the ion diameter. However, the dependence is weak only when the ions within each E R particle are relatively mobile. W h e n the system is fully polarized, the ions are forced close together, so that the results are sensitive to the exact size of the ions. Table 3.1: Results for b/a = 1, p* = 0.75, q* = 0.2, and E* = 50. * ion uncertainty 0.01 ± </0 («*/> 0.010 ± 0.04 0.88 0.639 ± 0.03 ± 0.6 ± (Pi) 0.02 -0.67 -31.0 0.97 0.97 0.96 0.05 0.604 0.85 -0.60 -29.2 0.1 0.558 0.67 -0.49 -26.9 A set of simulations was performed varying the number of ions per E R particle while keeping the total charge (q*Ni ) on and the reduced density of the ions (Ni ) on {Ni a* ) 3 on n fixed. T h e choice of quantities to keep fixed is somewhat arbitrary. For example, instead of fixing the charge density, the ion size could be held constant. T h e results are only expected to be qualitatively similar as the number of ions is varied. A set of simulations was performed both at a low reduced field of E* = 10 and a reduced field of E* — 100 which was sufficiently strong to cause chain formation. T h e results at E* = 10 and Chapter 3. Model Electrorheological Fluid 56 E* = 100 are shown in Tables 3.2 and Tables 3.3, respectively. A t low field, the results depend quite strongly on the number of ions. dramatic. A t high field, the dependence is not as In the high field case, the E R particles were highly polarized for all three systems, and the properties were mainly determined by the induced dipole. Since the total charge was the same for each of the three systems, the induced dipole in each system was also roughly the same. This was not the case for the lower field. A t high fields, it is reassuring that chain formation is a robust feature of the model and is reasonably insensitive to the number of ions per particle. Nevertheless, the particular value of field at which chaining occurs will depend, to some extent, on i V j . on A l l further simulations were performed with 8 ions per particle. Table 3.2: Results for 6/a q* N a* lon (p*) on uncertainty = !,/>* = 0.5, an d E* = 10. ± 0.010 ± K) 0.12 (Pi) <«}> ± 0.11 ± 0.1 ± 0.02 1.0 4 0.126 1.075 -4.22 -2.12 -9.9 0.92 0.5 8 0.1 0.829 -1.57 -0.85 -7.4 0.89 0.25 16 0.0794 0.585 -0.71 -0.33 -5.0 0.85 Table 3.3: Results for b/a = 1, p* = 0.5, and E* = 100 as a function of number of ions q* and a* were adjusted so that charge density and packing fraction per E R particle, on within E R particle remained fixed. a* 1 Nion ly a* ton (Pi) ( « 5 > u uncertainty ± 0.001 ± 0.01 ± 0.3 ± 0.1 ± 0.001 -168.8 0.996 1.0 4 0.126 1.696 6.30 -14.0 0.5 8 0.1 1.688 10.25 -11.2 -167.9 0.995 0.25 16 0.0794 1.628 10.30 -8.1 -161.9 0.994 Chapter 3. Model Electrorheological Fluid 3.6 57 C o m p u t a t i o n a l Tests We first checked how our results depend on the choice of boundary conditions (see T a ble 3.4). T h e Ewald expression depends on the dielectric constant e' of the continuum surrounding the assembly of periodic images. This affects the interaction energies which, in turn, affect all the other quantities. We compared the limiting cases of e' = 1 (vacuum) and e' —> oo (a conductor). T h e difference is that in the vacuum case, there is a positive contribution to the energy proportional to M . 2 In other words, the e' = 1 case suppresses the net moment of the system. There is a distinct difference between the e' — 1 and e' —> oo results. W h i c h value of d should be used? Ideally, a value of e' should be chosen consistent with the dielectric constant of the system. composed of materials of low dielectric constant. A typical E R fluid is However, even a moderate value of e' significantly damps out the contribution to the energy from the net moment (recall from Equation (2.3.12) that f(e') ~ 1/e')- Also, it is not clear that we would want to pick the boundary condition of e' = 1 which suppresses the net moment for a system which is very polarized. Thus, we have assumed that the dielectric constant of an E R fluid is sufficiently high that it is justified to use e' —> oo. T h e qualitative features displayed by the model remain unchanged by the value of e', although the particular choice of e' affects the numerical results. Results given in Table 3.4 also test whether the parameters a and n m a x in the Ewald sum were giving sufficiently accurate results. Larger values of a. and n m a x imply greater accuracy at the expense of more computations. It turned out that the values of a = 5 and n m a x = 14 were a reasonable compromise between speed and accuracy. These values are typical for simulations of ions and dipoles [37]. We have also tested for convergence problems by using a zero field and a high field configuration as the initial configuration for a simulation at an intermediate field. T h e results were then checked that both initial configurations lead to the same equilibrium Chapter 3. Model Electrorheological Fluid Table 3.4: Results for b/a = 1, p* = 0.75, 58 q* = 0.35, a* on = 0.05, and E* = n indicates the boundary condition used. For the Ewald method, listed are a, BC uncertainty values. ± 0.004 ± 0.03 ± 0.03 ± 0.1 (Pi) 0.848 0.63 0.941 0.900 0.98 -0.23 -1.22 -16.0 3,14,oo -17.1 0.950 5,14,oo 0.905 1.02 -1.42 -17.2 0.950 -17.2 0.951 -17.3 0.950 7,14,oo 0.906 1.03 5,25,oo 0.907 1.04 -1.44 BC and e'. ± 0.008 5,14,1 -1.44 max 20. For certain sets of parameters, this was the case. A t more extreme parameter values, the two systems would not converge to the same results over the length of a typical simulation. For example, convergence problems were readily apparent for small values of cr* . For the smaller ion diameter, the ions form bound pairs within the E R particle at on zero field. A t high field, the bound pairs of ions are formed across E R particles, locking a pair of E R particles together. In other words, a positive ion at the top end of one E R particle binds with a negative ion at the bottom of another E R particle (above the first E R particle). A n intermediate value of the field is not sufficient to break up the ion pairs in either case. Such effects most likely occur in real E R fluids as well. For example, a real, highly polarized E R fluid may take a long time to re-equilibrate once the applied field is removed. Parameters are avoided where such convergence problems are expected or encountered. T h e dependence of the results on the shape of the simulation cell was investigated. F r o m the length of the simulations, it was difficult to determine whether a chained system was i n a highly viscous liquid phase or a solid phase. A solid phase could be frustrated both by the fact that a chain could not completely span the box length and that the chains were not optimally spaced apart. T h e ground state lattice was assumed to be the Chapter 3. Model Electrorheological Fluid 59 same as that found for the model E R fluid used by Tao and Sun [9,10]. Simulations were started from a random configuration and from a lattice in both a cubic and rectangular cell. T h e lattice configuration in the rectangular cell was the ground state configuration (the cell dimensions were chosen to exactly fit the ground state) at the chosen reduced density of p* = 0.75. T h e results are given in Table 3.5. T h e parameters were chosen such that a chained system was found by starting from an initial random configuration. It is reassuring that starting from an initial lattice configuration in a cubic cell gave the same equilibrium results. In the rectangular cell, the main difference was that the reduced interaction energy dropped by almost 2 units. Chaining would most likely occur at a lower field for an optimally shaped simulation cell. We are more interested in the qualitative behavior of the model E R fluid, so no further consideration was given to the optimal cell shape. Since chaining is due to the interaction energy and the interaction energy was most affected by the cell shape, a more detailed study should account for the optimal cell dimensions. Table 3.5: Results for b/a = l, p* = 0.75, q* = 0.35, a* on box ratio initial configuration uncertainty < « 5 > </*•> ± 0.001 ± 0.05 = 0.05, and E* = 100 (w/> ± 0.17 ± 0.1 ± (PI) 0.001 random 1.242 5.37 -5.08 -123.3 0.993 1:1:1 lattice 1.242 5.38 -5.10 -123.3 0.993 6.53:3.27:4 lattice 1.248 5.71 -6.86 -123.9 0.993 1:1:1 As a check of the ./V dependence for the spherical E R particle model, we have compared 64 particles in a cubic box, 120 particles in a 6 : 5 : 4 rectangular box, and 128 particles in a 1 : 1 : 2 rectangular box (see Table 3.6). T h e results for all three system sizes did not quite agree with one another within the statistical uncertainties. Other tests of iV dependence for unchained systems showed little difference between 64 and 128 particle Chapter 3. Model Electrorheological Fluid results. 60 T h e slight differences for this chained system were attributed to the fitting of the chains within the simulation cell. Sample configurations for the 128 and 120 particle systems are shown in F i g . 3.1. As can be seen from the sample configurations, the cell shapes were quite different for each system. Table 3.6: Results for b/a = !,/>* = 0.5, q* = 0.5, a* = 0.1, and E* = 100. on TV uncertainty ± 0.001 ± 0.04 K> ± 0.1 ± 0.1 (Pi) ± 0.0008 64 1.688 10.24 -11.2 -167.9 0.9951 120 1.686 10.13 -10.6 -167.7 0.9949 128 1.683 9.99 -10.0 -167.4 0.9946 Results have been measured as a function of particle shape. As a particle axis becomes longer at fixed density, it takes fewer particles to span the entire simulation cell. Note that in terms of the unit of length, the particle axes are given by 2a* = ( f r / a ) / - 1 2b* 3 and — (b/a) / . Some tests for N dependence were performed at the extreme limits of 2 3 particle shape considered in this thesis. T h e results comparing 64 particles in a cubic box and 128 particles in a 1 : 1 : 2 rectangular box are shown in Table 3.7. Indeed, there is a noticeable system size dependence, especially in the interaction energy, for the b/a = 3 system. A little more than two particles end-to-end will span the box length for b/a = 3 ellipsoids in a cubic box. In contrast, it takes about three particles at b/a = 2 or b/a = 0.333. No N dependence was observed at these two elongations. Further E R simulations were performed with either 64 or 128 particles. T h e larger system size was used when finite size effects were expected for a 64 particle system. Chapter 3. Model Electrorheological Fluid 6 1 Figure 3.1: Sample configuration showing chains for 128 and 120 particles with b/a = 1, p* = 0.5, q* = 0.5, and a* = 0.1 at E* = 100. Field direction is indicated by arrow in on middle of box. Chapter 3. Model Electrorheological Fluid Table 3.7: Results for 62 p* —0.75, q* = 0.35, a* = 0.05, and E*= 40 on </**> («s> ( *i) <«*•> (Pi) (P*) (cos 7) ± 0.003 ± 0.04 ± 0.11 ± 0.1 ± 0.003 ± 0.001 ± 0.0008 3 2.501 6.41 -3.68 -99.1 0.991 0.970 0.9985 128 3 2.513 6.56 -4.74 -99.7 0.992 0.972 0.9986 64 2 1.856 5.15 -4.02 -73.4 0.988 0.950 0.9948 128 2 1.859 5.18 -3.99 -73.5 0.988 0.951 0.9949 64 0.333 1.687 4.15 -3.43 -66.6 0.988 0.0348 128 0.333 1.683 4.09 -3.44 -66.4 0.987 0.0354 N b/a uncertainty 64 3.7 u R e s u l t s a n d Discussion T h e configurations in F i g . 3.1 already show our model displaying E R fluid behavior of chain formation in an applied field. Note that the chains span the entire length of the box in each case, indicating that it is unlikely that the chain formation is an artifact of the simulation. For example, it may have been that a chain of 4 particles could span the box length, but the entropy for a chain of 8 particles would prevent such an occurrence. T w o points should be noted on this matter. T h e periodic boundary conditions imply that a chain stretching across the simulation cell represents an infinitely long chain. Also, while chains i n a real E R fluid may span the entire fluid, there may well be a natural chain length not detected by our simulations. In other words, for a sufficiently large system size, the chains may not be infinitely long, but rather strands of some finite length. T h e behavior of spherical and ellipsoidal E R fluid models is now investigated more systematically. 3.7.1 E R B e h a v i o r as a F u n c t i o n of E l e c t r i c F i e l d We now proceed to investigate the behavior of the model of spherical E R particles as a function of applied field. Some results for 64 spherical E R particles at a reduced density Chapter 3. Model Electrorheological Fluid 63 of p* = 0.75 are presented in Table 3.8. T h e system remains in a fluid phase for all but the highest field studied. A t E* = 100, the system was in a chained phase. T h e motion of the chains was slow, and they appeared to wriggle about an equilibrium position. The mean square displacements, shown in F i g . 3.2, confirm this picture. T h e mean square displacement at E* = 100 appears to level off, indicating a solid phase. T h e mean square displacements at all other reduced fields increased with the number of M C sweeps, indicating fluid behavior. Nevertheless, after a larger number of M C sweeps, the chains at E* = 100 may slip past one another, in which case, the mean square displacement would continue to slowly increase. Such behavior would indicate a viscous fluid. Ideally, one would like to measure the mean square displacement at E* = 100 for many more M C sweeps to determine whether the curve levels off, or drifts upwards. Table 3.8: Results for b/a = 1, p* = 0.75, q* = 0.35, and a* = 0.05 as a function of on E*. T h e number in brackets is the uncertainty (three standard deviations of the mean) in the final digit. E* </0 < « 5 > ± 0.04 ± <«J> 0.13 uncertainty ± 0.002 0 0.351 -2.00 -0.08 ± 0.1 (Pi) 0.13(1) 10 0.606 -0.82 -0.55 -5.2 0.858(9) 20 0.905 1.02 -1.42 -17.2 0.952(2) 50 1.158 3.37 -3.04 -57.0 0.985(3) 60 1.187 3.88 -3.59 -70.3 0.988(4) 75 1.215 4.52 -4.25 -90.2 0.9904(8) -5.08 -123.3 0.9930(6) 100 1.242 5.37 T h e average dipole moment (the magnitude of the dipole) and the order parameter (Pi) are plotted as a function of the reduced field in F i g . 3.3(a) and F i g . 3.3(b), respectively. T h e dipole moment measures the response of an individual particle to the field. (Pi) is proportional to the polarization of the whole system and is a measure of the Chapter 3. Model Electrorheological Fluid MC 64 sweeps Figure 3.2: T h e mean square displacement for b/a = 1, p* — 0.75, q* = 0.35, and a* on 100. — 0.05. Solid, dotted, dashed, and long-dashed curves are for E* — 0,20,50, and Chapter 3. Model Electrorheological Fluid 65 collective response of the system to the field. (Pi) grows faster than the dipole moment as a function of field. This happens because once each sphere develops a sufficient dipole moment, the dipole-dipole interactions act to further enhance the total moment of the entire system. A t zero field, the dipoles are randomly oriented so that (Pi) is roughly zero. However, due to the influence of surrounding E R particles, the induced dipole moment is quite large even at zero field. Consider a single E R particle. T h e small number of ions within the E R particle will give it a dipole at any instant in time. A t each instant, the dipole moment may be large, even though the vector sum of the dipole over time is zero. T h e neighboring E R particles will also have instantaneous dipoles and the dipole-dipole interactions will further enhance the dipole moment. This is very similar to the origin of van der Waals forces, which are due to induced-dipole-induced-dipole interactions. As the field is increased, the dipole moment approaches its m a x i m u m value where all the positive ions are pulled to one end of the E R particle with the negative ions at the opposite end. Furthermore, the system is essentially fully polarized. It is only under such conditions, far from the linear response regime (polarization M field E), proportional to the that chain formation occurs. From F i g . 3.3(b), the linear response regime does not extend beyond E* = 20, yet chaining occurs at a reduced field well above 50. Note that all the other simulations of E R fluids of which we are aware assume linear response [10,11,13,70]. A t high field, the total energy is dominated by the field energy with the self and interaction energies roughly cancelling each other out. T h e repulsion from an ion's likecharged neighbors within the E R particle are roughly cancelled by the attraction to the counter ions in the next E R particle in the chain. T h e reduced field energy is roughly given by (p*)E* (becoming more accurate the higher the field). In a strong field, chain formation is indicated by the peaks in g\\(z) in F i g . 3.4(a) at Chapter 3. Model Electrorheological Fluid 66 0 20 40 60 80 100 Figure 3.3: (a) Reduced dipole moment and (b) order parameter (Pi) as a function of reduced electric field for same system as in F i g . 3.2. Chapter 3. Model Electrorheological Fluid 67 1 diameter and 2 diameters. g(r) is given in F i g . 3.4(b) to show that, although it may indicate a structural change at E* = 100, it gives little information on the nature of the new structure. T h e dipole moment distributions are plotted in F i g . 3.5. T h e dipole moment distributions provide further evidence that system at E* — 100 is quite different from the same system at a lower field. In the high field case, all the E R particles essentially have the same dipole moment. However, this is not necessarily an indication that the E R particles are chained. For example, at E* = 50, all the E R particles also have roughly the same dipole moment, but the particles are not chained. It is a narrowing of the dipole moment distribution, rather than a shift to higher dipole moments, that accompanies the chaining. T h e system is chained at E* = 100, but not at E* = 50. Sample configurations show that the E R fluid at intermediate fields of E* = 60 and E* = 75 is chained, although there does not seem to be a sharp transition to a chained system. Sample configurations at E* = 50 and E* = 60 are shown in F i g . 3.6. T h e configurations indicate that chaining occurs at a reduced field somewhere between 50 and 60. Simulations were performed to show the difference between a model where the ions are mobile and one in which the E R particles have a fixed dipole. T h e final configuration for the system at E* — 100 was used as the starting configuration for a series of simulations with the ions held fixed. In effect, the system was composed of real (as opposed to point) dipoles with a reduced dipole moment of p* « 1.24. As the field was lowered, the chains persisted even to a low reduced field of E* = 5. A t zero field, g\\(z) still showed a large peak at 1 particle diameter, although the peak at 2 diameters was essentially flattened. Clearly, a model of fixed dipoles would predict a much lower field at which to observe chaining than our model. Simulations were performed starting with a single chain of 16 E R particles, where the chain exactly spanned the length of the cell. This single, virtually isolated, chain was Chapter 3. Model Electrorheological Fluid Figure 3.4: (a) 68 g\\(z) and (b) g(r) and for same system as in F i g . 3.2. Solid, dotted, dashed, and long-dashed curves are for E* = 0, 20, 50, and 100. Figure 3.5: Dipole moment distributions for same system as in F i g . 3.2. dashed, and long-dashed curves are for E* — 0,20,50 and 100. Solid, dotted, Chapter 3. Model Electrorheological Fluid 7 0 (A) (B) Figure 3.6: Sample configurations for same system as in F i g . 3.2 at (a) E* = 50 and (b) E* = 60. Chapter 3. Model Electrorheological Fluid 71 placed in a very high field and then the field was lowered until the chain became unstable. T h e chain was unstable at a field at which chain formation was observed at a reduced density of p* = 0.75. In fact, fields several times higher were required to stabilize the single chain. This means that there must be excluded volume and electrostatic effects which force particles to remain chained at high densities when the isolated chain would fall apart. Results as a function of field for a fixed prolate shape with b/a = 3 are shown in Table 3.9. For a fixed particle shape of b/a = 3, the field was varied and the (P2) order parameter was measured, (see F i g . 3.7). Even a small field is very effective at ordering the symmetry axes of a system of prolate ellipsoids. We show g\\(z; u) for a few values of the field (see F i g . 3.8). Since the results are for 128 particles in a 1:1:2 box, the function may be extended twice as far as for 64 particles in a cubic box at the same density. This allows us to see the second peak even at an elongation of b/a = 3. A large second peak is a good indication of chaining. Sample configurations are shown for reduced fields of 0, 10, and 100 (see F i g . 3.9). A t E* — 0, the system is disordered. A t a small reduced field of E* = 10, the system is strongly ordered, but not chained. A t E* — 100, the system is chained. T h e interesting feature for an E R fluid is that the viscosity generally grows as the square of the field. Although one cannot get true dynamical information from M C simulations, a pseudo-diffusion coefficient D may still be measured from the slope of the mean square displacement. We plot the inverse of D as a function of the square of the reduced field E* for spheres [Fig. 3.10(a)] and b/a = 3 ellipsoids [Fig. 3.10(b)]. We feel justified in comparing the slopes from different simulations because the step sizes were adjusted to give an acceptance rate of roughly a half for all the runs. Also, from Stokes' law [35], the product of the viscosity and diffusion is roughly constant for fixed particle mass, temperature, and density. Making the additional assumption that Stokes' law also Figure 3.7: (P ) as a function of field for 2 b/a = 3, p* = 0.75, q* = 0.35, and a* = 0.05. on Chapter 3. Model Electrorheological Fluid Figure 3.8: g\\(z\u) for same system as in F i g . 3.7. T h e solid, dotted, and dashed are for E* = 20, 50, and 100. Chapter 3. Model Electrorheological Fluid Figure 3.9: Sample configurations for b/a = 3, p* (a) E* = 0, (b) E* = 10, and (c) E* = 100. 74 0.75, q* = 0.35, and cr* on = 0.05 at Chapter 3. Model Electrorheological Fluid 75 Table 3.9: Results for b/a = 3, p* = 0.75, q* = 0.35, and cr* = 0.05 as a function of on reduced field. Results are for 128 particles in a 1:1:2 box. E* </0 (Pi) < « 5 > uncertainty ± 0.013 ± 0.04 ± 0.16 0 0.420 -1.84 -0.11 5 0.838 -0.90 -0.51 ± 0.1 -3.6 (P*) 0.10(2) 0.15(4) 0.80(1) 0.59(3) 10 1.638 1.59 -1.82 -15.8 0.962(2) 0.87(1) 20 2.226 4.14 -3.35 -43.8 0.984(3) 0.944(2) 40 2.513 6.56 -4.74 -99.7 0.9917(9) 0.9715(8) 50 2.573 7.48 -5.27 -127.8 0.9933(6) 0.9768(4) 60 2.618 8.37 -5.94 -156.2 0.9944(8) 0.9807(3) 75 2.661 9.46 -6.93 -198.6 0.9955(6) 0.9845(2) 100 2.702 10.88 -7.99 -269.2 0.9964(7) 0.9878(5) applies to our pseudo-coefficient, then the inverse of the pseudo-diffusion coefficient gives some indication of the viscosity. It should be noted that the pseudo-diffusion coefficient perpendicular to the field (D±) differs from the pseudo-diffusion coefficient along the field (D||) . W h a t we have measured is the trace of the diffusion tensor (D = | D | | + f - D j J Initially, the curves appear quite linear. A s chains are formed, the viscosity may grow faster than E . 2 However, the deviation from linearity is slight, and there is no accurate estimate of the uncertainty in the data points. Eventually, the system solidifies, in which case, the viscosity becomes roughly independent of the field in the high field limit. 3.7.2 E R B e h a v i o r as a F u n c t i o n of P a r t i c l e Shape T h e effect of E R particle shape on the electrorheological behavior was investigated. Results as a function of b/a at E* = 40 are presented in Table 3.10. T h e field was strong enough to orient the E R particles, but not to cause chain formation. T h e mean square displacements indicated that each system was in a fluid phase under the chosen conditions. T h e dipole moment as a function of b/a is shown in F i g . 3.11(a). T h e shape Chapter 3. Model Electrorheological Fluid 111 1 11 1 11 76 1 111 1 11 1 (a); 2 h1 /D 1h 0 1 1 1 1 1 1 11 1 1 1 I I I 1 I I I 0 20 40 60 80 100 E* / 100 2 3 o5 1111111111111111111 1 /D 0 20 40 60 80 100 E* / 100 2 Figure 3.10: Inverse of pseudo-diffusion coefficient as a function of square of reduced field for (a) spheres (b/a = 1) and (b) prolates (b/a — 3) at p* = 0.75, q* = 0.35, and < n = 0-05. Chapter 3. Model Electrorheological Fluid 77 of the curve is mainly due to the ellipsoidal shape of the E R particles. Recall from Equation (3.2.1) that 1 Umax ~ fb\ 2/S 2 \a) * Nionq 1 fCax ~ 2 i m q (prelates) (3.7.1a) ~ ° L* * (oblates). (3.7.1b) ion /6\~ N ~ N q*b* 1 / 3 \a) Ni n( a For prolates, the reduced dipole moment grows as (b/a) / 2 (b/a)~ l . l 3 3 and for oblates it grows as This is due to the choice of length scale <7, where a is the diameter of a sphere of equal volume to the ellipsoids. W h e n comparing results as a function of particle size, we are always comparing E R particles of equal volume. Table 3.10: Results for p* = 0.75, q* = 0.35, a\ on = 0.05, and E* = 40 as a function of particle shape. Results are for 128 particles in a 1:1:2 box at b/a — 3. b/a </0 («*/> ± 0.05 ± 0.13 (Pi) <«F> ± 0.2 ± 0.003 (cos 7) uncertainty ± 0.004 3 2.513 6.56 -4.74 -99.7 0.992 0.99858(4) 2 1.856 5.15 -4.02 -73.4 0.988 0.9948(1) 1.25 1.290 3.45 -2.96 -50.7 0.982 0.954(4) 1.1 1.165 2.99 -2.70 -45.7 0.981 0.83(1) 1 1.116 2.81 -2.59 -43.7 0.980 0.49(1) 0.9 1.149 2.94 -2.64 -45.0 0.981 0.25(1) 0.8 1.202 3.13 -2.81 -47.2 0.981 0.162(5) 0.5 1.437 3.69 -3.00 -56.6 0.984 0.064(1) 0.333 1.687 4.15 -3.43 -66.6 0.988 0.0348(6) As indicated by the order parameter (P ) in F i g . 3.11(b), even slightly prolate particles 2 were readily ordered by the field. T h e oblate case is quite different. For oblate ellipsoids, there is little ordering with respect to the symmetry axes. However, the induced dipole within the oblate particle orders quite readily in a sufficiently strong electric field. T h e particles tend to maximize their dipole moment by orienting themselves with the longest Chapter 3. Model Electrorheological Fluid 3 1 II I I I | I I I I | I I I I | I M ' - _l 111 11 1 1 11 1 | 1 l-ljl 1 • (a) I I 1 1 1 1 T1 -1 1 1 I — I ' —n ' ' ' ' ' ' ' ' ' I ' 0 . 2 F_l -3 1 1 1 1 1 11 1 2.5 78 "3 1 2 3 b/a Figure 3.11: (a) T h e reduced dipole moment, (b) order parameter (P2), and (c) average of cosine of angle between the symmetry axis and induced dipole as a function of particle shape for f^max' p* = 0.75, q* = 0.35, a* = 0.05, and E* = 40. T h e dashed curve in (a) shows on Chapter 3. Model Electrorheological Fluid 79 axis along the field direction. Thus the symmetry axis for an oblate particle is roughly perpendicular to the field. This effect is most dramatically shown from the cosine of the angle ((cos7)) between the induced dipole and the symmetry axis. This quantity as a function of particle shape is shown in F i g . 3.11(c). It is almost a step function. T h e field orders both oblate and prolate particles along their longest axis. It is surprising that even slightly non-spherical particles respond so dramatically to the field. g\\(z;u) is shown in F i g . 3.12 for various particle shapes. indicates that g\\(z) Recall that the notation was measured with respect to the symmetry axes. For prolate parti- cles, the symmetry axis and induced dipole lie along the same direction. As noted, this is not the case for oblate ellipsoids, where the two directions are perpendicular to one another. Thus, a large first neighbor peak in g\\(z;u) is only seen for prolate particles. T h e situation for oblate particles is shown by a sample configuration for 128 particles with elongation b/a = 0.333. F i g . 3.13 shows how the particles are aligned with the field, yet the symmetry axis of each particle is free to rotate in the plane perpendicular to the field. One interesting question to investigate is whether the dipole moment grows faster with the field for non-spherical particles. Thus we have simulated a single E R particle as a function of field for various particle shapes. T h e dipole moment increased from the sphere to oblate ellipsoid to prolate ellipsoid. However, when scaled by the m a x i m u m dipole moment for each shape, the behavior of the dipole moment as a function of field was similar for each shape (see F i g . 3.14). T h e dipole moment shows a rapid increase over the same range of the electric field for each particle shape. T h e actual value of the dipole moment of an isolated E R particle was only a few percent lower than the value for an E R particle in a system at reduced density of p* = 0.75. This was true even when the dense system was chained. T h e reason is that the particle is almost fully polarized at the high fields at which chaining occurs. Neighboring particles can only slightly enhance Figure 3.12: g\\(z;u) for same system as in F i g . 3.11. T h e solid, dotted, dashed, and long-dashed curves are for b/a = 3,2,1,0.5, and 0.333. Curve at b/a = 3 for 128 particles. Chapter 3. Model Electrorheological Fluid Figure 3.13: Sample configuration for b/a = 0.333, p* = 0.75, q* = 0.35, and a\ 81 on at E* = 40. = 0.05 Chapter 3. Model Electrorheological Fluid 82 the dipole moment. 3.7.3 O r d e r i n g of an E R F l u i d of Oblate E l l i p s o i d s as a F u n c t i o n of D e n s i t y Recall that a system of oblate ellipsoids ordered in the field such that the symmetry axis was roughly perpendicular to the field direction. It is possible for the symmetry axes of oblate particles to align along a single direction perpendicular to the field. Indeed, we have attempted to form a biaxial phase by increasing the density of a system of oblate ellipsoids in an orienting field. T h e results for b/a = 0.333 as a function of density are presented i n Table 3.11. Higher density configurations were generated by using N P T simulations with a very high pressure. This way, we avoided the problem of starting high density simulations from a lattice configuration. T h e biaxial order parameter (Q\ ) w a s 2 calculated for the oblate ellipsoids. Recall that Q\ 2 was discussed i n section 2.7. (Q\ ) 2 ' s zero (except for finite size effects) in both an isotropic and uniaxial phase but increases towards unity in a biaxial phase. (Q22) — 0-79 indicates that the system is clearly in a biaxial phase at p* = 1.05 (at p* — 0.9, the system is in a weakly biaxial phase). g\\(z; p) and g\\(z;u) are shown in Figs. 3.15(a) and 3.15(b), respectively. g\\(z; p) is peaked around 2a* RS 1.44, indicating ordering in the field along the degenerate oblate ellipsoid axis. T h e only feature of note i n g\\(z; p) is that the peak becomes broadened as the density is increased. T h e density was increased by shrinking the simulation cell, resulting in a tilting or sliding of the oblate particles within the chain. This causes the broadening of the peak. T h e broadening is a result of the box length being incommensurate with an integer number of particles in a chain. As such, no real significance should be attributed to the broadening of the peak. T h e g\\(z\u) curves prove more interesting. g\\(z; u) at p* = 1.05 shows peaks at roughly integer multiples of 2b* « 0.48. These peaks confirm that the oblate ellipsoids are ordering along their symmetry axes. T h e mean square displacements indicate that Figure 3.14: Scaled dipole moment as a function of reduced electric field for a single particle for q* = 0.35, and a* = 0.05. Solid, dashed, and dot-dashed curves for b/a = 3,1, and 0.333. on Chapter 3. Model Electrorheological Fluid 84 10 8 6 gn(z;M) A 4 2 0 0 0.5 1 1.5 2 i i i i i i i i i i i i i_ -" " l i i z/cr 1 gn(z;u) II II II = r (b) j l;i III i\ w \ I \ 11 0 11 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 0.5 1 1.5 z/a 2 Figure 3.15: (a) g\\(z; p) and (b) g\\(z;u) for 6/a = 0.333, q* = 0.35, and a* E* = 100. Solid, dotted, and dashed curves are for p* = 0.75,0.9, and 1.05. on = 0.05 at Chapter 3. Model Electrorheological Fluid 85 the system was in a highly viscous or solid phase at all densities considered. Table 3.11: Results for b/a = 0.333, q* = 0.35, a* = 0.05, and E* = 100 as a function on of density. p* uncertainty 0.75 3.7.4 </**> ± 0.001 1.848 («5> ± 0.04 7.23 ± 0.09 -6.96 <«F> ± 0.1 -184.0 (Qh) ± (Pi) (cos 7) ± 0.0009 ± 0.0003 0.22 0.9958 0.0190 0.0201 0.0205 0.10 0.9 1.840 6.93 -6.55 -183.1 0.48 0.9952 1.05 1.836 6.79 -6.53 -182.7 0.79 0.9953 Summary A model for an E R fluid has been presented. Unlike previous models, our model explicitly considered the ions confined to an E R particle's surface. It is these ions which give rise to the large induced dipole of an E R particle in an electric field. T h e interesting properties of an E R fluid are due to the chaining of the E R particles along the field direction. O u r model reproduces this structural change in the E R fluid. Previous studies of the E R fluid have assumed that chaining occurs in the linear response regime. We find that chaining occurred well past the field for which linear response was valid. T h e dependence of E R fluid behavior on E R particle shape was investigated. The most dramatic effect was due to the increase in the induced dipole, which formed along the longest axis in the ellipsoidal E R particle. This meant that a prolate particle aligned with its symmetry axis along the field whereas an oblate particle aligned with its symmetry axis perpendicular to the field. Thus it was possible to form a biaxial phase for a system of oblate E R particles in a field. Densities and length to breadth ratios were avoided for which the ellipsoidal particles orientationally order in the absence of a field. T h e electric field readily ordered the ellipsoidal E R particles well before chaining occurred. Chapter 3. Model Electrorheological Fluid 86 This allows for the formation of liquid crystal phases using.an electric field for systems which would otherwise be isotropic. Chapter 4 Liquid Crystal Models 4.1 Introduction Some of the earliest attempts to investigate the isotropic-nematic (IN) transition used a lattice model where nearest-neighbor particles interact via the Lebwohl-Lasher potential [75,76] u(12) AP (cos7) , 2 where A is a constant and the angle dependence enters through P 2 ( c o s 7 ) , which is defined in section 4.2. This potential was chosen to allow comparisons between simulation results and predictions of Maier-Saupe theory [77]. to simulate a true nematic liquid on a lattice. O f course, one could not expect Nevertheless, the IN transition may be somewhat independent of the translational degrees of freedom [37]. Luckhurst and Romano [78] included translational degrees of freedom using a Lennard-Jones potential plus an anisotropic term similar in form to the Lebwohl-Lasher potential. Luckhurst and Romano provided evidence for an IN transition in this system. Others soon used their location of the IN transition temperature to compare with integral equation [79] and density functional calculations [80]. We have more completely studied the liquid crystal model of Luckhurst and Romano. T h e model is not the best one for a real liquid crystal. key features required for an IN transition. Nevertheless, it contains the O u r interest in this model was to use the Gibbs ensemble Monte Carlo method to simulate the IN coexistence. T h e conventional 87 Chapter 4. Liquid Crystal Models 88 Gibbs method readily gives the gas-isotropic and gas-nematic coexistence. Despite several attempts [81,82], the conventional Gibbs method has not been previously used to simulate IN coexistence. B y biasing the orientation of the inserted particles during the particle exchange step of the Gibbs method, we were able to achieve this. This is the first example of the Gibbs method being successfully used to simulate IN coexistence. A t first, we attempted to reproduce the N P T Monte Carlo results of Luckhurst and Romano [78]. We found a lower IN transition temperature than they report. Further- more, there was a large hysteresis in the cyclic process of melting the solid and then freezing the resulting liquid at constant pressure. During the freezing branch of the hysteresis loop, an IN transition occurred in the liquid. Thus, our IN transition temperature was within the range over which the hysteresis occurred. In other words, the nematic phase may be metastable with respect to the solid. T o avoid this uncertainty about the thermodynamic stability of the nematic phase, we have studied the same liquid crystal model with increased anisotropic interactions. IN coexistence was determined for this more strongly interacting system using Gibbs ensemble simulations. We have also studied a liquid crystal model of hard spheres interacting via an anisotropic potential. From a combination of N P T and Gibbs simulations, a portion of the phase diagram for this hard-core model was constructed. T h e phase behavior for this model is interesting since no gas-isotropic liquid coexistence was found. Rather, a coexistence between an isotropic fluid and a nematic was found. This is interesting because no gasisotropic liquid transition is found for dipolar hard spheres [29,30]. T h e hard-core liquid crystal model is another example where the attractive, anisotropic forces cannot stabilize a distinct isotropic liquid phase. This is explored further in the next chapter. In section 4.2 we define the liquid crystal models and give some computational details. In section 4.3 we report and discuss the N P T and Gibbs ensemble Monte Carlo results for the various liquid crystal models. In section 4.4 the results are summarized in a Chapter 4. Liquid Crystal Models 89 discussion of the phase diagrams for the liquid crystal models. 4.2 M o d e l and C o m p u t a t i o n a l D e t a i l s We consider axially symmetric particles which interact via a pair potential of the form « ( 1 2 ) = u ( r ) + ix (12) , 0 (4.2.1) 0 where uo(r) is a spherically symmetric interaction, u (12) is an anisotropic term and r is a the distance between particles 1 and 2. We consider two forms for u (r). 0 model, we take u (r) 0 to be the usual Lennard-Jones potential ULJ{T). also consider a hard-core model where u (r) 0 For a soft-core In addition, we is the usual hard-sphere potential u#s(r). T h e anisotropic contribution to the potential is u„(12) = -4Ae ( ^ ) P ( C O S ) , (4.2.2) 6 2 7 where P 2 ( c o s 7 ) is the usual second order Legendre polynomial [P (cc) = (3x — l)/2] and 2 2 7 is the angle between the symmetry axes of the particles. Although the particles are spherical, each particle contains an embedded unit vector giving it an arbitrary symmetry axis. In Equation (4.2.2), e and a are the Lennard-Jones parameters and A is a variable determining the strength of the anisotropic interaction. We note that for the hard- core model, the anisotropy parameter A and the reduced temperature T* = kT/e may be combined to give an effective temperature T * / A . In other words, for the hard-core model, A and T* are not independent variables. T h e potential was spherically truncated for separations greater than half the box length. T h e usual corrections were applied to the Lennard-Jones contribution to the energy [37]. T h e correction for the Lennard-Jones contribution to the total potential energy was given by Equation (2.3.4). We now discuss the correction due to the anisotropic term in the pair potential. T h e pair distribution function g(l2) plays an important part Chapter 4. Liquid Crystal Models 90 in the theory of liquids and is defined in several of the references (eg. see [35]). In the isotropic phase, the angle-dependent pair distribution function g(l2) of separation r for large r. It follows that g(l2) is a function only — g(r) for large r in the isotropic phase. Assuming g(r) ~ 1 beyond the cut-off separation r , the anisotropic contribution to the cut correction for the potential integrates to zero in the isotropic phase. T h e same argument cannot be made in the nematic phase. T h e projection g (r) 220 is one term in the general expansion of the pair distribution function [18,83]. In the nematic phase, the projection g (r) 220 must satisfy [18] (r) 220 9 where (P ) 2 - 5(P ) 2 2 , is the equilibrium order parameter [20]. energy is then X(P ) U -6, 2 2 cor r - o o , (4.2.3) T h e correction for the anisotropic where (4.2.4) is the correction to the energy from the attractive part of the Lennard-Jones potential. Including the correction to the anisotropic part of the potential is inconvenient since it requires knowing (P ) 2 consistently. beforehand. In practice, this requires determining (P ) 2 self- We have tested the effect of including the correction in an orientationally ordered phase. T h e difference made by including the correction should be most apparent in the anisotropic contribution to the energy for a small system. T h e difference was greater than the statistical uncertainty only for a 32 particle system. For larger systems, the anisotropic contribution to the correction is not significant. Thus, we do not include any correction to the anisotropic part of the potential. It should be noted that any correction which assumes g(r) m 1 beyond the cut-off is not valid in the solid phase. There are several ways of determining whether the system is in an isotropic or nematic phase. Typically, one calculates the P order parameter which is zero in an isotropic phase 2 Chapter 4. Liquid Crystal Models 91 and has a m a x i m u m value of one in the nematic phase. Recall that the calculation of P 2 was discussed in section 2.7. We have also calculated the radial distribution function g(r) and the projection g (r). 220 g(r) is of little help i n distinguishing between the isotropic and nematic phases. As previously noted, however, g {r) 220 decays to zero in the isotropic phase but has a finite value for large r i n the nematic phase. Technically, large r should be less than the scale of spatial director fluctuations, although this is not a problem i n practice due to the small system size. For an N P T simulation, we define a M C sweep as N translational attempted moves, N rotational attempted moves and a fraction of N attempted volume changes. N is the number of particles and the number of attempted volume changes typically ranges from Af/lO to N/2. Furthermore, we randomly chose which type of move to perform each time. Thus it is more accurate to say that, on average, we performed N translational attempted moves etc. each M C sweep. Step sizes were adjusted so that each of the three types of moves was accepted about half the time. A M C sweep for a Gibbs ensemble simulation is the same as that for an N P T simulation except that n h exc moves are also performed. n h exc particle exchange was adjusted so that between 1 and 3 percent of the particles were exchanged during each M C sweep. T h e reason for this target acceptance rate was discussed in section 2.2. In general, the Gibbs method may fail to simulate IN coexistence for two reasons. It may fail because the acceptance rate for particle exchanges between the two dense, liquid phases is too small [81,82,84]. W h y are particle exchanges so difficult between the isotropic and nematic phase? After all, we can exchange particles between a gas and liquid, so exchanging particles between two liquid phases cannot be much more difficult. T h e problem is that even if one finds room to insert a particle i n the nematic phase, unless it is oriented in approximately the same direction as the surrounding particles, the move will still be unfavorable. In other words, if you have a one in x chance of squeezing a Chapter 4. Liquid Crystal Models 92 particle into the nematic phase and a one in y chance of orienting it correctly, the chance of actually inserting the particle is the product of these two probabilities. Obviously, a helpful thing to do would be to orient the particle approximately correctly while trying to insert it. Indeed, this is precisely what we have attempted to do by biasing particle insertions in the Gibbs method for IN coexistence. T h e Gibbs method may also fail because the density difference between the isotropic and nematic phases is so small that random density fluctuations, due to particle ex- changes, pushes each box in and out of each phase. From our own N P T and Gibbs ensemble simulation results, we knew where to expect IN coexistence. Yet when we attempted to simulate IN coexistence using the conventional Gibbs method, each of the two simulation boxes would rapidly change identity back and forth between an isotropic and nematic phase. T h e problem was that the insertion of randomly oriented particles into the nematic phase constantly pushed it into the isotropic phase. T h e other types of moves would then cause the orientational order to build up again in one of the two boxes and the process would repeat itself. Again, a way to avoid this problem is to rotationally bias the inserted particles. T h e rotational biasing has the effect of stabilizing the nematic phase with respect to orientational fluctuations. W i t h rotational biasing, for the liquid crystal models under consideration, there was then a region along the IN coexistence curve where the density difference was sufficiently large that each box retained its identity long enough to collect useful statistical averages. A t any particular moment, a phase with orientational degrees of freedom has a director d, and the symmetry axis of each particle makes an angle 0 with that director. There is a slow, spatial variation in d for a large, nematic system. A typical simulation cell is sufficiently small that the spatial variation is not a problem in practice, and a single, coordinate-independent director adequately characterizes the system. We may Chapter 4. Liquid Crystal Models 93 accumulate an orientational distribution function /(cos 9) normalized such that (4.2.5) For an isotropic phase, /(cos 9) = 1, whereas /(cos 9) increases rapidly as cos 9 approaches unity in a nematic phase. During the course of the simulation, the distribution function /(cos 9) is accumulated for each of the two simulation boxes used in the Gibbs method. A t first, /(cos 9) is quite noisy, but it quickly settles down to a smooth curve. T h e orientation of the particle to be inserted is then sampled from the distribution function of the receiving box. T h e details on how this was performed are given below. Another technical detail is that a production run was typically started from an equilibrium configuration. T h e initial function /(cos 9) will then be close to the equilibrium distribution function. This is convenient, but not necessary. Eppenga and Frenkel [20] used a similar biasing method to calculate the chemical potential of hard platelets. T h e y performed N P T simulations and only the insertion of ghost particles was required to determine the chemical potential. A l l ghost particles were inserted with cos 9 = 1 and /(cos 9) was accumulated and extrapolated to cos 9 = 1. T h i s extrapolated value was needed to recover the unbiased value of the chemical potential [see Equation (4.2.8)]. Cracknell et al. [85] have also proposed a method for rotationally biasing insertions in the Gibbs method. T h e y used their biasing method in a simulation of liquid water, which is isotropic, with only local orientational ordering. Their biasing method is based on picking a random position in the receiving box and constructing a function based on the energy of insertion as a function of orientation. W i t h this method, they were able to increase the acceptance probability of an insertion by about a factor of 2. We now explain how the particle insertions are rotationally biased in the present calculations. A trial particle exchange from box I to box 77 is accepted with a probability Chapter 4. Liquid Crystal Models given by m i n ( l , 94 where P h,bias) exc _ / (cosfl) J i It.^.UI is the expression for an unbiased exchange [57]. Recall from section 2.2 that *exch,bias — .;;/ \±exch n r (cos6i) J and Pexc/i Aw is given by Equation (2.2.5). We must evaluate the two = exp(Au;) where Pexch quantities f (cos9) I and f (cosd). n Evaluating f (cos9) I is straightforward. T h e chosen particle to be removed (from box 7) has a unit vector it. One calculates cos 9 = \u-d\ and determines the corresponding value of f (cos9). I Determining f (cos9) n is slightly more complicated. T h e orientation of the particle to be inserted into box II is sampled from the orientational distribution function f (cos9) n this sampling is as follows. of box II. One method of performing A random number R is uniformly selected from the unit interval. This value is equated with the area under the curve R= C f 11 f (cos9) n from cos0 to 1, (cos 9')d cos 9'. (4.2.7) ./cos 8 T h e inserted particle then makes an angle 9 with the director of the receiving box (box II). Equation (4.2.7) is an implicit equation for determining cos0. since f 11 In other words, (cos 9) is accumulated throughout the simulation, the only unknown in E q u a - tion (4.2.7) is cos 9. For example, if the receiving box happened to be an isotropic phase, then f 11 (cos 9') = 1 and cos 9 = 1 — R. In this particular case, a value of cos 9 is uniformly selected from the unit interval, as is the case for the conventional Gibbs method. We use Widom's equation for calculating the residual chemical potential the trial particle insertions during the Gibbs simulation. p Tes from Since the particle insertions are rotationally biased, this bias must be removed to recover the correct estimate of the residual chemical potential. T h e equation with this correction is p res = -kT\n^ ^ exp-^ ^ u+ {N i)f , (4.2.8) Chapter 4. Liquid Crystal Models where AU + = U(N + 1) — U(N) 95 is the change in energy due to inserting the particle. / is the value of the orientational distribution function evaluated at cos where cos 9 is the angle the inserted particle makes with the director. 4.3 R e s u l t s and Discussion 4.3.1 R e s u l t s at A = 0.15 We have performed N P T simulations for systems of Lennard-Jones particles modified by an anisotropic interaction given by Equation (4.2.2). T h e reduced pressure P * = P a / e 3 was set at 0.1. A s a test, we confirmed that N V T results at densities corresponding to this pressure agreed with the N P T results within the statistical uncertainty. In this section, we consider the system at anisotropic parameter A = 0.15. We attempted to enter the nematic phase in two different ways. In one case, we started with a low temperature solid and slowly heated it. In the other case, we started with a high temperature isotropic phase and slowly cooled it. T h e system solidifies into an fee lattice. There are certain numbers of particles (4i , i = 1, 2, 3 ...) for which the fee lattice exactly fits the cubic simulation 3 cell. O f these numbers, we used 32, 108, and 256. To test the number dependence of the transition temperature, we also performed simulations using 64, 150, and 512 particles. We first discuss our results for 256 particles. U p o n heating the fee lattice, we found that it remained solid up to T* = 0.85. We performed long runs of 60, 000 M C sweeps at T* = 0.85, attempting to melt the solid. For T* > 0.85, the solid melts directly into the isotropic phase. We then cooled the isotropic phase. We found no dramatic j u m p in density signaling the transition to a nematic phase [see F i g . 4.1(a)]. the order parameter (P ) N 2 Nevertheless, quickly increased over the reduced temperature range from T* = 0.775 to T* = 0.725, indicating a transition from an isotropic to a nematic liquid phase at T* « 0.75 [see F i g . 4.1(b)]. We continued to cool this nematic phase and there Chapter 4. Liquid Crystal Models Figure 4.1: 96 Temperature dependence of (a) density and (b) order parameter for 256 particles. Results are from N P T M C simulations at P* = 0.1 of Lennard-Jones particles with an anisotropic interaction of strength A = 0.15. There is a low temperature solid branch with a nematic branch below and a high temperature isotropic branch. dashed squares are the M C results of Luckhurst and Romano [78]. The Chapter 4. Liquid Crystal Models 97 was a transition to a solid at T * = 0.6. T h e system spontaneously solidified into a (slightly imperfect) fee lattice at this reduced temperature. T h e radial distribution function and mean square displacement support this conclusion. Also, the order parameter of this system at T* = 0.6 is very close to that of the fee solid. It is interesting to note that such a relatively large system can show spontaneous freezing into a well-defined crystal structure. One way of viewing the N P T results is that the cyclic process of melting the solid and then freezing the resulting liquid was performed under constant pressure. This process generated a hysteresis loop composed of a melting branch with a melting to an isotropic liquid at T* ~ liquid at T* ~ 0.6. 0.85, and a freezing branch with a freezing from a nematic T h e hysteresis loop is reproducible and is due to metastability. It is tempting to presume that the solid phase remains metastable past its true melting temperature. However, without knowing the free energies, one cannot determine which phase is metastable with respect to another phase. is an IN transition. Along the freezing branch, there There is no hysteresis associated with the IN transition. If the IN transition is a true phase transition, in the sense that it is a transition between two thermodynamically stable phases, then it occurs at T* f» 0.75. We can only conclude that the solid is either metastable or the thermodynamically stable phase over the reduced temperature range from T* = 0.6 to T* = 0.85. We show g(r) and g (r) 220 for typical isotropic, nematic and solid phases in F i g . 4.2. We also include a plot of the mean square displacement versus the number of M C sweeps (see F i g . 4.3). Although particles in the nematic phase do not diffuse as quickly as in the isotropic phase, they are clearly diffusing. Recall from section 3.4 that a linear mean square displacement with non-zero slope indicates diffusion. T h e mean square displacement in the solid phase flattens out after a certain number of M C sweeps. Luckhurst and Romano have also performed 256 particle N P T simulations [78] for Chapter 4. Liquid Crystal Models 98 1 1.5 2 2.5 3 r/cr 1 1.5 2 2.5 3 r/<7 Figure 4.2: (a) Radial distribution functions and (b) projections g (r) 220 for the same system as in figure 4.1. T h e solid, dotted, and dashed curves are for T* = 0.5,0.7, and 0.825. 99 Chapter 4. Liquid Crystal Models 0.4 1 1 1 11 1 1 1 11 1 | 1 1 1 | 1 I <[r(s)-r(0)] > / a 2 1 1 2 / 0.3 - / — / / / / / / 0.2 / 1 / / 0.1 / f 1 ..-**" — - 11 ~ :// — / — V 0 I 0 I 1 I 1 1 1 1 1 I 1 I 200 400 600 MC sweeps 1 1 1 I I 1 800 1000 Figure 4.3: T h e mean square displacement as a function of the number of M C sweeps for the same system as in figure 4.1. T h e solid, dotted, and dashed curves are for T* = 0.5, 0.7, and 0.825. Each curve is an average of ten consecutive segments forming the entire run. Chapter 4. Liquid Crystal Models this model at P* = 0.1 and A = 0.15. 100 Their results are compared with ours in F i g . 4.1. Although their densities correspond to our I N branch, their order parameters more closely match our solid branch. Thus Luckhurst and Romano find a higher temperature (T* ~ 0.89) for the IN transition. Luckhurst and Romano determined the IN reduced transition temperature using the following procedure. A low temperature Lennard-Jones liquid configuration was used as the initial configuration. T h e embedded unit vectors were then aligned in a single direction and the anisotropic potential was turned on. was then heated in stages at constant pressure. T h e system T h e relatively sharp disappearance of orientational order marked the transition from the nematic to isotropic phase. We differ in several important ways from Luckhurst and Romano in how we performed our N P T simulations. Luckhurst and Romano rotate and translate a particle at the same time, in a single trial move. Thus the rotational step size and translational step size are combined to give a single acceptance probability for the trial move. Also, they only performed a single volume change per M C sweep. In principle, Luckhurst and Romano's method and our method should lead to the same results. However, by adjusting the orientational step size independently of the translational step size, as well as performing more volume changes, we speed up the convergence of the system to equilibrium. Their results are consistent with the conjecture that they did not allow enough moves for the orientational degrees of freedom to relax to equilibrium. Also, we could afford to perform much longer (at least 4 times as long) simulations to assure convergence. We now report results for other system sizes. T h e reduced density and (P ) 2 order parameter are plotted as a function of temperature for various system sizes in Figs. 4.4 and 4.5, respectively. First considered are results for N = 108 and N = 32, the other particle numbers which correspond to an fee lattice. For N — 108, as with N = 256, it was found that the solid remains stable or metastable up to T* = 0.85 [see F i g . 4.4(c)]. U p o n cooling the isotropic phase, the system shows some kind of transition at T* = 0.75 since it Chapter 4. Liquid Crystal Models 1.1 1 101 111111111111111111111111 = 0.9 - 0.9 j 111111111111111111111111 (b) I- 0.8 r 0.8 0.7 iiiI 'i i i i iI i i i i Ii i i i Im 0.5 0.6 0.7 0.8 0.9 1.1 1 0.9 1.1 1 P r f 1 0 7 (c) I III I II II I III III II M 0.5 0.6 0.7 0.8 0.9 r 1 (d) j 111111111111111111111111 111111111111111111111111 F- "l I I I 3 F- 0.8 h 0 7 hi 111111111111111111111 r 0.5 0.6 0.7 0.8 0.9 1 0.7 "i i i i Ii i i iIi i i i I i i iiI iir i 0.5 0.6 0.7 0.8 0 . 9 r 1 Figure 4.4: Temperature dependence of density for (a) 32, (b) 64, (c) 108, and (d) 150 particles. Results are for the same system as in figure 4.1, except that the system size was varied. A vertical line indicates that the system flipped between the two densities during the course of the run. Chapter 4. Liquid Crystal Models 0.5 0.6 0.7 0.8 0.9 102 1 0.5 0.6 0.7 0.8 0.9 1 Figure 4.5: Temperature dependence of order parameter for (a) 32, (b) 64, (c) 108, and (d) 150 particles. Results are for the same system as in figure 4.1, except that the system size was varied. Chapter 4. Liquid Crystal Models 103 occasionally slips from the density of an isotropic phase to a slightly higher density. T h i s may be evidence for a transition to either a nematic phase or a solid. Indeed, the N — 108 system spontaneously freezes into an imperfect fee lattice at T* = 0.7. Due to the small system size, the 32 particle system jumps back and forth between an isotropic phase and solid (with the energy being more negative in the solid) over a temperature range of 0.85 to 0.9 [see F i g . 4.4(a)]. T h e 32 particle system never formed a nematic phase. A n y simulations started below this temperature range spontaneously freeze except at the lowest temperature we considered of T* = 0.5. A glassy solid was found at T* = 0.5. A very long simulation of 5 million M C sweeps (10 times longer than for the other 32 particle systems) was performed to check that the system remained in this glassy phase. T h i s must be a quenched state, which is metastable with respect to the solid. Since the 64 particle system cannot fit exactly onto the fee lattice, we only cool the isotropic phase. We find that at T* — 0.75 the system flips between a liquid and solid phase [see F i g . 4.4(b)]. T h e liquid phase at T* = 0.75 has a small order parameter of ( P 2 ) ~ 0.28 [see F i g . 4.5(b)]. Furthermore, the order parameter increased smoothly to this value before the phase transition. If one considers this as a weakly nematic phase, then this implies that the IN transition occurs at a temperature just above this nematicsolid transition. However, this value of (P2) is barely above the N dependence in the estimate of the order parameter. For the 150 particle system, we again find a weakly nematic phase at T* = 0.75 [see Figs. 4.4(d) and 4.5(d)]. It is not surprising to find a nematic at T* = 0.7 since the solid is destabilized with respect to the nematic due to the incommensurability of the lattice with the number of particles. T h e system then solidified in going from T* = 0.7 to T* = 0.6. W i t h results at a variety of particle numbers, we looked for any systematic size dependence. In F i g . 4.6, both u* and u* are plotted as a function of 1/N. T h e energies for the isotropic phase at T* = 0.9 show no noticeable N dependence. As expected, a Chapter 4. Liquid Crystal Models 10 104 T—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—r 8 6 -<u> 4 - 2 LB 0 Figure 4.6: - B B - - ^ " B B IB--IS— BM-Bt—->•—H---P—+—!—>—+•••!—j-B I 0 0.01 0.02 1/N I 0.03 Energy (filled squares) and contribution of anisotropy term to the energy (open squares) at T* = 0.7 (solid line) and T* = 0.9 (dashed line) as a function of reciprocal particle number for the same system as in figure 4.1. A t T* = 0.7, the system was a solid and at T* = 0.9, it was in a nematic phase. Chapter 4. Liquid Crystal Models 105 the reduced anisotropic energy in this phase is close to zero. A t T* = 0.7, there is a significant difference between the small and large N results. A t larger N, the system is somewhat glassy and cannot freeze onto the lattice. A t smaller N, the system does freeze onto the fee lattice, thus the total energy is then more negative than for a glass or imperfect lattice. This difference is larger than any systematic N dependence. To summarize our N P T results at A = 0.15, we find that the fee solid melts to an isotropic liquid just above T* = 0.85. T h e solid phase may not have melted in the simulations until past the true melting temperature since we start with a simulation cell which exactly fits the fee lattice. Evidence against this is that the melting temperature was independent of the system size. U p o n cooling the liquid, we find a transition to a nematic phase around T * = 0.75. This transition is characterized by a smooth increase in the order parameter and with no dramatic j u m p in density. It is known that the I N transition must be first order [86]. However, the density difference between the coexisting phases is small so that the transition is only weakly first order. T h e discontinuity in density for the IN transition at A = 0.15 could not be resolved by the simulation results. A t or below T* — 0.75, the nematic phase spontaneously freezes for all but the smallest and largest systems studied. For the smallest size of 32 particles, the system froze directly from the isotropic liquid at T* — 0.85. phase whatsoever. For 32 particles, we found no nematic For 512 particles, the system was large and not commensurate with the fee lattice, so that freezing into lattice was not observed within the length of our simulations. T h e previous theoretical and computational results for the IN transition temperature are now discussed, with reference to our own M C results. Luckhurst and Romano found Tj N = 0.89 ± 0 . 0 1 and report that Maier-Saupe theory predicts Tf N = 0.811. Perera et al. [79] solved the integral equation theory using the reference hypernetted-chain ( R H N C ) closure for this potential at p* = 0.79. This is the nematic density at the transition Chapter 4. Liquid Crystal Models 106 temperature reported by Luckhurst and Romano, although it may also be taken as the isotropic density since the density change is ~ 1%. Integral equation theory determines the spinodal temperature, the temperature at which the system becomes orientationally unstable. T h e stability limit (SL) temperature was found to be T$ L = 0.786. Below this temperature, the system is in a phase with orientational order. This may be a nematic or solid. Density functional theory was then used to determine the transition, as opposed to the spinodal, temperature [80]. nematic phase at Tf N Density functional theory predicts a transition to a = 0.875 from an isotropic phase at p* = 0.79. Interestingly, the density functional theory reports the reduced density of the nematic phase as p* = 1.04. This amounts to a change in density of ~ 30%. T h e structure of the higher density phase, whether it is a nematic or solid, is an input for the density functional calculations. T h e reported results assumed a nematic phase. Such a large change in density suggests that the nematic phase may be physically unrealistic and that the high density phase is a solid. T h e density functional calculations would have to be repeated assuming a transition to a solid to test this. Luckhurst and Romano's M C result of Tj N the integral equation theory result of T§ L the density functional theory result of Tf N = 0.89 is consistent with = 0.786 and is in reasonable agreement with = 0.875. However, density functional theory tends to over-estimate the transition temperature. Furthermore, neither Luckhurst and Romano's results nor our results for the IN transition show the large change in density as calculated using density functional theory. T h e assorted results are given in Table 4.1. Due to the hysteresis in our N P T results, the existence of a thermodynamically stable nematic phase along the P* ='0.1 isobar is uncertain. Along the melting branch of the hysteresis loop, a solid-isotropic transition was found at T* ~ 0.85. It is tempting to claim that this agrees with density functional theory result. However, the density functional result is only valid for a transition to a nematic. In the process of freezing the liquid, we also found an orientational ordering of the isotropic at T* ~ 0.75. T h i s result is Chapter 4. Liquid Crystal Models 107 consistent with the integral equation result (recall that the integral equation result is for a slightly different density). We have not determined whether the solid or nematic is thermodynamically stable below T* ~ 0.75. Nevertheless, if there is a thermodynamically stable IN transition, then it occurs at T* ~ 0.75. This is in disagreement with Luckhurst and Romano's result of T* 0.89. Table 4.1: Previous results for Lennard-Jones particles with an anisotropic interaction of strength A = 0.15. M S is Maier-Saupe theory [78], L R is Luckhurst and Romano [78], S L is stability limit [79], and D F T is density functional theory [80]. T h e Monte Carlo results ( L R ) were at reduced pressure P* = 0.1. T h e number in brackets is the uncertainty (three standard deviations of the mean) in the final digit. (SL and D F T ) were for fixed isotropic reduced density p* IN IN 1 P*N Pi T h e theoretical predictions = 0.79. MS LR SL DFT 0.811 0.89(1) 0.786 0.875 0.79 0.79 1.04 0.790(2) 0.780(3) We now discuss the Gibbs ensemble simulation results. A t A = 0.15, we have determined the gas-liquid coexistence for this system using Gibbs ensemble simulations. We used a total number of 400 or 512 particles and adjusted the total volume of each run so that there were, very roughly, an equal number of particles in each box at equilibrium. T h e larger number of particles was used close to the critical point. T h e gas-liquid coexistence curve was constructed and is shown in F i g . 4.7. determined for the gas-isotropic transition at A = 0.15. T h e critical point values were Using the method of rectilinear diameters [56], we find T* fa 1.3, p* « 0.34, and P* w 0.1. As the temperature was lowered, an ordering transition occurred along the coexisting liquid curve. A t the lowest reduced temperature of T * = 0.7, the order parameter in the coexisting liquid was (P2) ~ 0.6. T h e transition occurred around T* = 0.75 at P* ss 0.005, which is consistent with our N P T results. This ordering transition is discussed further in section 4.4. Chapter 4. Liquid Crystal Models 108 i i i i i i I i i i I i ii 1.4 1.2 - o ho T* 1 0.8 • 0.6 0 i i i I i i l i 0.2 0.4 0.6 0.8 i i H 1 Figure 4.7: Gibbs ensemble results for the gas-liquid coexistence curve for Lennard-Jones particles with an anisotropic interaction of strength A = 0.15. Chapter 4. Liquid Crystal Models 4.3.2 109 Results at A = 0.3 A t the anisotropy parameter A = 0.15 for the soft-core liquid crystal model, we were unable to confirm that there was a thermodynamically stable nematic phase along the P* = 0.1 isobar. However, as the Gibbs ensemble results at A = 0.15 suggest, it may be much easier to find a stable nematic at lower pressure. We increased the anisotropy parameter to A = 0.3 and repeated the search for a nematic phase for the soft-core model. N P T simulations were performed for 256 particles at a reduced pressure of P* — 0.1. T h e system was started on an fee lattice and the temperature was slowly increased. T h e solid phase remained stable up to a reduced temperature of T * = 1.18. T h e solid then melted to a nematic liquid when the reduced temperature was increased to T* = 1.19. The IN transition then occurred at T* = 1.2. T h e isotropic phase was stable up to the gasisotropic reduced transition temperature T* = 1.3. We then attempted to freeze the liquid, however, the nematic phase did not solidify upon cooling from T* = 1.2. Rather, it became more viscous, reaching a low temperature glassy state. T h e N P T results are shown in F i g . 4.8. T o summarize, at P* = 0.1, we have found an IN transition at T* = 1.2. Furthermore, the reduced melting temperature is below T* = 1.19. Thus the nematic phase is certainly the thermodynamically stable phase, at least over the narrow reduced temperature range from 1.19 to 1.2. From F i g . 4.8(a), we see that the IN transition is clearly first order. Also, the jump in order parameter from the isotropic to nematic phase for A = 0.3 is much larger than for A = 0.15. Confident that a nematic phase exists at A = 0.3, we then performed Gibbs ensemble simulations to construct the phase diagram (see F i g . 4.9). These results are interesting because they give the entire phase diagram, apart from the nematic-solid and gas-solid coexistence, for our soft-core liquid crystal model. We used a total of 512 particles for the gas-isotropic and gas-nematic runs and 600 particles for the IN runs. Even Chapter 4. Liquid Crystal Models 110 1.2 1 0.8 p* 0.6 0.4 0.2 0 0.8 1 r 1.2 1.4 1.2 1.4 1 0.8 0.6 <P»>n „ 0.4 0.2 0 0.8 1 T* Figure 4.8: Temperature dependence of (a) density and (b) order parameter for 256 particles. Results are from N P T M C simulations at P* = 0.1 of Lennard-Jones particles with an anisotropic interaction of strength A = 0.3. There is a low temperature solid branch, below which there is a glassy nematic to isotropic branch. A t high temperature is the gas branch. Chapter 4. Liquid Crystal Models 111 |_l I I | I I I | I I I | I I I | I I l_| r-l l l I I I I II I I II I I I I I l-l 0 0.2 0.4 0.6 0.8 1 * P 1.1 1.2 r 1.3 Figure 4.9: Gibbs ensemble results giving part of the phase diagram for Lennard-Jones particles with anisotropic interaction of strength A = 0.3. Chapter 4. Liquid Crystal Models 112 with rotationally biased insertions, each box in the IN coexistence would occasionally exchange identity. T h e slightly larger system size of 600 particles helped to decrease the frequency of this occurrence. T h e frequency of boxes exchanging identity was low enough that useful statistical averages could be collected over much of the simulation. In practice, the frequency was roughly once every 10,000 M C sweeps. F r o m F i g . 4.9(a) and F i g . 4.9(b), one can readily identify the triple point and two critical points for this system. 4.3.3 R e s u l t s for the H a r d - c o r e M o d e l We may use hard spheres, instead of the Lennard-Jones particles, combined with the anisotropic potential in equation (4.2.2) to model a liquid crystal. This is our hard-core liquid crystal model. To compare with integral equation results [79], we set A = 0.15 and vary T*. Nevertheless, the thermodynamic state for this model is fully specified by P* and T*IX. Gibbs ensemble simulations were performed for this hard-core model using both conventional and rotationally biased particle exchanges. Production runs were from 5,000 to 10,000 M C sweeps for a total of 512 particles. T h e coexistence results are plotted in F i g . 4.10 and other values are given in Table 4.2. No gas-isotropic liquid transition was found for this model, so the coexistence is an isotropic fluid-nematic coexistence. One might expect three fluid phases for this system: an isotropic gas, an isotropic liquid, and a nematic liquid. Instead, only two fluid phases were found. T h e density difference over much of the coexistence branch was large enough to allow the Gibbs method to work with conventional particle exchanges. As previously noted, if the density difference had been too small then, with unbiased exchanges, the identity in each Gibbs box would flip between an ordered and disordered phase. Reassuringly, both methods of performing particle exchanges gave identical results within the statistical uncertainty. To show this, in Table 4.2, values are given for both methods at T* = 0.333. Chapter 4. Liquid Crystal Models 113 0.4 0.35 \0.3 h 0.25 0.2 1.2 J 0.4 I M 0.6 0.8 1 1 1 1 | I I | y(b) ] 1 x i 1 P* 0.8 ~ J / I 0.6 0.4 I I . - / "l 0.25 1 1 1 /l — , , , , 1, , , " 0.3 0.35 0.4 T* Figure 4.10: Gibbs ensemble results giving part of the phase diagram for hard spheres with an anisotropic interaction of strength A = 0.15. T h e dashed line in (b) gives the points at which a 108 particle system melts from a fee lattice. T h i s line was generated by melting the lattice at constant pressure using N P T M C simulations at a variety of pressures. Chapter 4. Liquid Crystal Models 114 As the temperature is lowered, the nematic phase must become metastable with respect to the solid. This may be the case for the lower temperature Gibbs results. T o test for this possibility, an fee lattice of 108 particles was heated under constant pressure until it melted. B y repeating the process at various pressures, we generated the dashed curve in F i g . 4.10(b). Loosely speaking, this is the stability limit of the solid (at least for 108 particles). We may then be confident that at least the 3 highest temperature Gibbs results give the thermodynamically stable coexisting phases. Table 4.2: Gibbs results for the isotropic fluid-nematic coexistence curve for the hard-core liquid crystal model with anisotropy parameter A = 0.15. Results are for 512 particles. p rea is the residual chemical potential. (/) is the value of the orientation distribution function /(cos 0) extrapolated to cos $ — 1. T h e final column is the percentage of the total number of particles exchanged per M C sweep. p* I Pi-es 1 kT II I II 0.353 1.05(4) 1.06(9) 1.76(7) 1.77(6) 0.333 0.93(2) 0.86(5) 1.39(4) 1.39(6) 0.333 0.95(3) 0.87(6) 1.39(6) 1.43(11) 0.316 0.80(2) 0.74(6) 0.99(4) 1.00(9) 0.300 0.68(1) 0.67(13) 0.69(5) 0.70(9) 0.22(5) 0.21(10) 0.286 0.53(2) 0.43(27) % exchanged (A) I II I II 6.8(14) 2.7(9) 0.13(9) 0.49(10) 1.9 11.2(8) 2.1(3) 0.07(5) 0.68(2) 1.6 - - 0.08(4) 0.68(3) 1.8 13.6(7) 2.1(3) 0.06(3) 0.74(1) 1.8 18.0(11) 1.8(1) 0.04(2) 0.810(8) 1.6 1.8(1) 0.030(9) 0.841(7) 1.9 21.1(10) T h e integral equation calculations were performed for this hard-core model at p* = 0.8 at A = 0.15 [79]. T h e temperature was lowered until the isotropic system became orientationally unstable. This gives the stability limit temperature for the I N transition. Chapter 4. Liquid Crystal Models 115 The reduced stability limit temperature using the R H N C closure was Tg L that integral equation theory with the R H N C closure predicted T§ L core model at p* = 0.79. T h i s value for T£ L = 0.85. Recall = 0.786 for the soft- agrees quite well with our M C results for the soft-core model. Replacing the soft-core (Lennard-Jones potential) by a hard-core particle was not expected to significantly perturb the IN transition. T h e integral equation calculations found an IN coexistence between two phases as typical liquid densities. Indeed, N V T simulations, performed for the hard-core model at p* = 0.8 at A = 0.15, show an orientational transition at a temperature in the neighborhood predicted by the integral equation calculations. However, our Gibbs results are for densities and temperatures well below the values used in the integral equation calculations. O u r Gibbs results most closely resemble results for a gas-liquid coexistence, except that the liquid phase is orientationally ordered. 4.4 Summary Let us now consider the three p vs. T phase diagrams. T h e p vs. T diagrams for the softcore model for A = 0.15 and A = 0.3 are shown in Figs. 4.7 and 4.9(a), respectively. T h e p vs. T diagram for the hard-core model is shown in F i g . 4.10(a). B o t h the diagrams for the soft-core model show a gas-isotropic liquid coexistence. This coexistence is due mainly to the Lennard-Jones contribution to the potential. Since both phases are orientationally disordered, the anisotropic contribution to the energy is minimal. T h u s the anisotropic potential does not significantly perturb the gas-isotropic liquid coexistence. Indeed, the transition is essentially the gas-isotropic transition for the pure Lennard-Jones system. In fact, the critical point values for this transition at both A = 0.15 and A = 0.3 are in good agreement with the critical point values for the Lennard-Jones system (A = 0) [56]. One feature that appears in the A = 0.3 phase diagram, and possibly in the A = 0.15 Chapter 4. Liquid Crystal Models case as well, is an IN coexistence. 116 This coexistence appears as a small hump in the coexisting liquid curve [right-hand-side of F i g . 4.9(a)]. It may be that such a hump does exist in the A = 0.15 case for the soft-core model, but the density difference between the isotropic and nematic phases may be too small to resolve by our computer simulations. Evidence for this is the fact that the coexisting liquid becomes nematic as the temperature was lowered. However, as suggested by the N P T results for A = 0.15, the IN transition may be pre-empted by an isotropic-solid transition. T h e phase diagram for the hard-core model is quite different from that for the softcore model. Replacing the Lennard-Jones particles by hard spheres, with no anistropic potential, results in the disappearance of the gas-liquid coexistence. hard spheres only have a fluid and solid phase. In other words, Replacing the Lennard-Jones particles by hard spheres, keeping the anistropic potential, results in a lowering in temperature of the gas-liquid critical point (the definition of the reduced temperature is consistent for all three of the phase diagrams). Furthermore, the liquid phase for the hard-core model is orientationally ordered. One may consider this as an IN coexistence, although the isotropic phase is a low density phase. For this IN transition, the lowering in energy due to orientational ordering goes into stabilizing the condensed phase relative to the gas. How is the IN critical point for the soft-core model at A = 0.3 affected by going to a hard-core model? In the soft-core model, the IN coexistence is between two phases at typical liquid densities. Thus the lowering in energy due to orientational ordering mainly compensates for the lowering of the entropy in the nematic liquid phase. A s such, the IN critical point found in the soft-core model should not be significantly perturbed by going to a hard-core model. This interpretation is consistent with integral equation calculations. Gibbs ensemble simulations were not performed to look for this higher density IN transtion in the hard-core model. However, N V T results indicate that the transition does exist in the neighborhood predicted by integral equation theory. Chapter 4. Liquid Crystal Models 117 As mentioned, no gas-isotropic liquid coexistence was found for the hard-core model. This situation is analogous to that for dipolar hard spheres, where the existence of a gas-isotropic liquid coexistence is uncertain [29]. T h e reason that a gas-isotropic liquid transition may not exist for dipolar hard spheres is that the entropy gain from forming a gas of dipolar clusters overwhelms the stabilizing energy from entering the isotropic liquid state. T h e energy difference between a gas of clusters and the isotropic liquid is small. For dipolar particles, the situation is complicated by the long-range nature of the dipole-dipole interaction. For the hard-core liquid crystal model, the situation is simpler. T h e reason that a gas-isotropic liquid coexistence does not exist for the hard-core model is as follows. T h e gas phase has essentially zero configurational energy. T h e energy of the isotropic liquid arises only from short-range orientational correlations between particles. This energy is not sufficient to stabilize the isotropic liquid phase with respect to the gas. Chapter 5 M i x t u r e s of N e u t r a l a n d D i p o l a r H a r d Spheres 5.1 Introduction T h e Gibbs method [56] has been used to investigate several types of binary fluid mixtures. Several binary mixtures of Lennard-Jones fluids were investigated by Panagiotopoulos and co-workers [57]. T h e fluid may undergo a gas-liquid transition as well as a demixing transition. T h e Lennard-Jones mixtures of Panagiotopoulos et al. provided examples of both types of transition. T h e y investigated the transitions as a function of pressure at fixed temperature. Experimentally, one typically varies the temperature at fixed pressure. This provides the phase diagram (temperature vs. mole fraction) of a binary mixture one usually finds in the literature. T h e Gibbs method is ideally suited to performing simulations in close analogy with the experimental situation. Recall from section 2.2 that for a pure system with two phases in coexistence, the Gibbs phase rule implies that there is only one independent variable. For a two component system, there is an additional independent variable (the chemical potential of the second component). Thus we may fix both the temperature and pressure (we are free to chose which two independent variables to fix) in the Gibbs ensemble Monte Carlo simulation of a binary mixture. A t constant pressure, the coexistence curve as a function of temperature may be constructed, as is done experimentally. Demixing phase transitions of non-additive soft disks [87] and nonadditive hard spheres [88] have been investigated using this method. A non-additive mixture of spheres of types A and B has the property that 2aAB ^ °~AA + &BB, where cr,j 118 Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 119 is the effective collision diameter for a pair of particles of types i and j . T h e demixing phase transition for a mixture of neutral (n) and dipolar (d) hard spheres has been investigated using integral equation theory [31,32]. T h e theoretical results indicated that such a demixing transition did exist. It is not obvious that computer simulations would confirm the existence of a demixing transition. Integral equation theory also predicts a gas-isotropic liquid transition for dipolar hard spheres [33,34], yet computer simulations have failed to find such a transition. Some workers have suggested that no such transition exists [29,30]. It is difficult to exclude entirely the possibility of a gas-isotropic transition. T h e transition may occur at a temperature too low to get reliable simulation results. T h e condensation of hard-sphere ions is not much affected by the presence of neutral hard spheres [89]. T h e condensation of hard-sphere ions in a background of neutral hard spheres would appear as a demixing phase transition. Our simulation results show that a mixture of neutral and dipolar hard spheres demixes. This suggests that the demixing of the neutral and dipolar hard-sphere mixture may be related to a gas-isotropic liquid transition in a background of neutral hard spheres. However, with no background of hard spheres, such a gas-isotropic liquid transition has not been found. Structural changes, such as a clustering of the vapor phase, or another transition may occur before the temperature at which the gas-isotropic liquid transition would occur. We have determined the critical temperature as a function of pressure for the demixing phase transition using Gibbs ensemble Monte Carlo simulations. For a system of dipolar hard spheres, the dipole moment is a variable equivalent to the temperature (see section 5.2). Thus we have determined the critical point dipole moment for the transition. For a two component system, the critical temperature varies with the pressure, thus it is more accurate to say that the critical temperature line was determined. T h e coexistence curve was then mapped out at a fixed pressure. Coexistence curves were constructed for Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 120 several values of the diameter of the neutral hard spheres. It was expected that the reduced critical dipole moment would increase as the diameter of the neutral hard spheres was decreased. T h e reason is that no gas-isotropic liquid transition was found for dipolar hard spheres at the parameters at which we found a demixing transition. In the limit of vanishing neutral hard spheres, the demixing transition resembles, if it exists, the gasisotropic liquid transition for dipolar hard spheres. A n extrapolation of the critical point dipole moment as the diameter of the hard spheres vanishes should give an indication of the temperature at which the gas-isotropic liquid transition might occur, at least at the fixed pressure. In section 5.2, we describe the model and give some computational details. Some computational tests are discussed in section 5.3. In section 5.4, we report our results for determining the critical dipole moment as a function of pressure. In section 5.5, we give the results for the critical dipole moment as a function of the diameter of the neutral hard spheres. 5.2 M o d e l a n d C o m p u t a t i o n a l Details A mixture consists of Nd dipolar hard spheres of diameter ad and N n of diameter a . n ( a d = 1). neutral hard spheres T h e diameter of the dipolar hard spheres is taken as the unit of length T h e ratio of component diameters is a* = a /ad n (reduced diameter of the neutral hard spheres). T h e total number of particles is N = Nd + N n density is p* = pa\, where p = N/V and the reduced is the total number density. T h e mole fractions are = N /N = 1 — Xd- Each dipolar hard sphere has a reduced dipole moment p* = p/'y4nCotkl'a . Equivalently, a reduced temperature may be defined Xd — Nd/N and x n n d through T * = l/p* 2 To avoid confusion, we have consistently chosen to characterize a system in terms of p*, as opposed to T*. Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 121 Due to the long-range nature of the dipole-dipole interaction, periodic boundary conditions were employed using the Ewald summation method [63,64]. T h e numerical parameters used for the Ewald method were a convergence parameter a = 5 with a m a x i m u m lattice vector length of e' = 1. n = 14 and a surrounding continuum of dielectric constant max Whenever using the Ewald summation method, one is faced with choosing an appropriate value of e'. For a mixture rich in neutral hard spheres, e' = 1 is a sensible choice. Systems of hard sphere dipoles can have a large dielectric constant [90]. Thus a value of c! —> oo may be a reasonable choice for a mixture rich in dipolar hard spheres. Either choice of e' has little effect on the thermodynamic properties, but does affect the net moment fluctuations. Such fluctuations are not of interest for the purposes of investigating the demixing phase transition. T h e Gibbs method involves two boxes, labeled I and II. of Nd, N , n and the reduced total volume V* = V/a d One may specify the values for the two boxes. Instead of initially specifying the total volume (thus the total density), we may specify the reduced pressure P* = Pa /kT. d u* and cr*, Along with this completely characterizes the system. Trial translations are performed for each particle within each box. T h e step size for each component was adjusted so that the acceptance rate for each component was roughly a half. For mixtures with different sized components, the step sizes for each component will be quite different. Trial rotations are performed for each dipole within each box. Trial volume exchanges are performed between the two boxes such that the total volume remains fixed. Recall that such moves ensure that the pressures in each box are equal to one another. Particle exchanges for each component are performed between the two boxes such that the total number of neutral and dipolar hard spheres also remains fixed. Recall that such moves ensure that the chemical potentials of each component in each box are equal to one another (p^ = p / 1 and p l = pi ). 1 n Typically, the acceptance probability for a particle exchange of one component is much different from that for the other component. Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 122 In this case, we pick the more difficult component to exchange with a fixed probability greater than 1/2. If the initial values are such that the uniform mixture is unstable, the system will phase separate into a dipole-rich phase in one box and a neutral-rich phase in the other box. Instead of performing particle exchanges for both components, we may perform particle exchanges for one component along with identity exchanges. For an identity exchange, we convert a neutral into a dipole in one box and at the same time convert a dipole into a neutral in the other box. It is easier to exchange a neutral particle, since one must only find room for it in the receiving box. Exchanging a dipole usually involves an energy penalty for removing a dipole from an energetically favorable position and inserting it into an energetically unfavorable position. T h e benefit of exchanging neutral spheres only becomes better as the diameter of the neutral spheres shrinks. Standard deviations were calculated in the usual way of dividing the run into ten equal length blocks and treating block sub-averages as independent measurements. A l l simulations were of 10,000 M C sweeps (N trial translations, Nd trial rotations, 50 to 100 trial volume changes and enough trial exchanges so that roughly 10 particles of each type are exchanged per M C sweep). Often, it was clear that block sub-averages were correlated, so simulations were continued for further sets of 10,000 M C sweeps. Numerically, this did not much affect the equilibrium values. T h e main reason for continuing simulations was to ensure convergence to equilibrium. For the higher reduced dipole moments, the acceptance rate for particle exchanges was quite low, and convergence required many M C sweeps. T h e identity exchange simulations gave identical results with a much higher acceptance rate for exchange moves. The identity exchanges allowed us to perform simulations as higher reduced dipole moments. T h e only drawback was that one could not calculate the chemical potential of the dipoles from the identity exchanges. Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 5.3 123 Computational Tests As a test of both the correctness of our computer code and as a check on the Ewald boundary condition parameters, we have reproduced published N V T results for a system of dipolar hard spheres [90,91]. We have also reproduced the demixing phase transition for previously published systems using the Gibbs method both with particle exchanges and with particle exchanges combined with identity exchanges. We have reproduced results for demixing non-additive hard spheres [88] and demixing non-additive soft disks [87]. We have also reproduced the Gibbs results of Panagiotopoulos et al. [57] for LennardJones mixtures. T h e Lennard-Jones mixture was composed of two components, labeled A and B. T h e total (i.e., in both boxes) mole fraction was then XA,totai = NA/{NA + NB)- The total mole fraction is constant throughout the Gibbs ensemble simulation, although the mole fraction in each box is free to vary. Provided we have a total mole fraction within the demixing region (within the two phase area of the T vs. XA phase diagram), the equilibrium results should be independent of this total mole fraction. As a test that the coexistence curve at fixed T and P is independent of total mole fraction, we have varied the total mole fraction for one of the Lennard-Jones mixtures. Our results, along with the published results are shown in Table 5.1. As expected, the results for XA and p* were independent (within the statistical uncertainty) of XA,totai- Of course, the number of particles in each box had to vary to maintain the same equilibrium values as XA,totai was varied. Recall that the usual mixing rules for the Lennard-Jones cross parameters are O~AB — {&A + C B ) / 2 XA,totai and EAB = \J^A^B- For a total mole fraction of = 0-2, two identical gas mixtures were found and at XA,totai = 0-7, two identical liquid mixtures were found. These total mole fractions lie outside the two phase region of the T vs. xA phase diagram. Provided that the total mole fraction was not too close to a coexistence mole fraction, the equilibrium results were independent of total mole Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 124 fraction. A system with a mole fraction near a coexistence mole fraction will remain in a metastable state instead of demixing. T h e drawback is that the only way to test if one is too close to a coexistence mole fraction, is to perform a series of simulations at a variety of total mole fractions. Table 5.1: Gibbs results for a Lennard-Jones mixture at T * = 0.928, P* = 0.076. Lennard-Jones parameters for the B component were <JB = 0.768 and SB = 0.597. T h e standard mixing rules were used to determine the Lennard-Jones cross parameters (see text). T h e first row shows the published results [57]. T h e number in brackets is the uncertainty (three standard deviations of the mean) in the final digit(s). X A,total x A I II I (P*.) (N) II I II 116(5) 484(5) 0.666(14) 0.195(17) 0.848(15) 0.098(4) 0.3 0.69(1) 0.207(5) 0.836(9) 0.099(1) 0.433 0.67(1) 0.20(1) 0.841(11) 0.098(6) 299(12) 301(12) 0.22(2) 0.836(6) 0.100(3) 431(30) 169(30) 0.567 0.70(3) A s a test for the demixing of neutral and dipolar hard-sphere mixtures, we performed simulations in several different ways. We checked that, when using identity exchanges, the results were independent of the component for which particle exchanges were performed. In other words, although performing particle exchanges using neutral hard spheres is more efficient, the same results were obtained exchanging dipolar hard spheres. T h e dependence on the initial configuration was also tested. Indeed, at sufficiently high pressure and/or dipole moment, the results did depend on the starting configuration. A n example of such a situation is shown in Table 5.2. T h e systems are identical except that the starting configuration was mixed for one system and entirely demixed for the other system. Each of the two systems was run for at least 50,000 M C sweeps before averages were collected for 10,000 M C sweeps. Note that the mole fractions in the dipolerich phase are quite different. It cannot be determined which, if either, system has reached Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 125 equilibrium. Converged, demixed results could only be found in a small region of P * x p* space. For too low values of either P * or //*, the system is super-critical. In other words, there is no demixing transition at the low values. For too high values, it is enormously difficult to get converged results. Through trial and error, values for P* and p* were determined at which converged, demixed results could be obtained. Table 5.2: Results obtained from starting with two different initial configurations for a P* = 3.7 Gibbs ensemble simulation for 500 particles at p* = 2.0. The / II I 0.712(4) 0.685(7) 0.745(6) 0.683(3) (x n) (AO II I II 0.674(34) 0.948(39) 270(4) 230(4) 0.394(28) 0.968(14) 146(5) 354(5) dependence of the demixing transition on the total mole fraction was tested (see Table 5.3). Varying the total mole fraction implies that the number of particles in each box must change to maintain the equilibrium mole fractions. T h e results were expected to be slightly different due to weak N dependence. Nevertheless, the results are reasonably independent of total mole fraction. Table 5.3: P* = 2.0 Gibbs results at p* = 2.2 for two different values of total mole fraction x n (x I (AO II / // 0.8 0.160(23) 0.983(9) 111(5) 389(5) 0.6 0.175(23) 0.985(11) 238(6) 262(6) We have tested the JV dependence for the demixing dipolar and neutral hard spheres. Indeed, the i V dependence was a frustrating problem. For a small system ( i V = 2 5 0 ) , it was possible for the dipoles to form a chain spanning the entire simulation cell. T h e Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 126 periodic boundary conditions stabilized this structure as an infinitely long chain. For a larger system (N = 500), the chained structure was not present. see clusters, strands, or rings of dipoles. Instead, one might We used the largest system size that was computationally practical to avoid, as much as possible, any unphysical dipole structures. A t high reduced dipole moment, the dipole-rich phase may become ferroelectric. Simulations were performed to test for a ferroelectric phase at the highest reduced dipole moment used in our calculations. No ferroelectric phase was found, confirming that our simulations were of two orientationally isotropic phases. 5.4 T h e C r i t i c a l D i p o l e M o m e n t as a F u n c t i o n of P r e s s u r e We first attempted to confirm the theoretical predictions that the demixing transition did, in fact, occur. Chen and Forstmann [32] have determined the stability of the mixture at p* = y/2.5 as a function of reduced density p* and mole fraction x . n the mixture was expected to demix for p* > 0.66. At x n Constant volume Gibbs = 0.8, ensemble simulations were performed for an equal diameter mixture of 500 particles with total mole fraction x = 0.8 and total reduced density p* = 0.7. Starting at p* = \/2JS n 1.6, the reduced dipole moment was increased until the demixing phase transition occurred. A t p* = 1.8, the system was clearly mixed and at p* = 2.0, it was clearly demixed. A t p* = 1.9, the system showed large fluctuations from a mixed phase. We take this to be the critical reduced dipole moment. P* T h e transition occurred at a reduced pressure of « 3.7. T h e critical reduced dipole moment occurs at a critical mole fraction. Forstmann report the mole fraction at their critical point 0.8. Chen and Simulation results confirm that the mixture does not demix symmetrically (i.e., critical point at mole fraction x n = 0.5). T h e simulation results are consistent with a mole fraction Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres x n ft 0.8 at the critical reduced dipole moment. to the theoretical predictions. 127 Our results are qualitatively similar We started with a mixture which was predicted to be unstable. However, we found that the reduced dipole moment had to be increased from 1.6 to 1.8 (cooling the mixture) before demixing occured. Thus the theoretical predictions underestimate the critical reduced dipole moment. Chen et al. [31] have determined the stability of the mixture at p* = 0.8 as a function of the reduced dipole moment and mole fraction x . n becomes unstable for p* > 1.5 at x = 0.8. n We repeated our constant volume Gibbs ensemble simulations at total reduced density p* = 0.8. dipole moment p* = 1.8 at P* ft 5.8. T h e y report that the mixture We found the critical reduced Again the theoretical predictions underestimate the critical reduced dipole moment. We have attempted to find the reduced pressure at which the mixture becomes unstable at p* = 2.0. Again using a total of 500 particles with total mole fraction x n = 0.8, we increased the total reduced density until the system demixed. We found the reduced pressure to be P* ft 2.4 at p* = 2.0. T h e critical reduced dipole moments are tabulated as a function of reduced pressure in Table 5.4. Table 5.4: T h e critical reduced dipole moment as a function of reduced pressure for equal diameter mixture of neutral and dipolar hard spheres, mixture just before the demixing transition. p* Pi 1.8 P* 0.8 3.7 1.9 0.7 2.4 2.0 0.6 5.9 p* is the reduced density of the Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 5.5 128 The Critical Dipole Moment as a Function of Diameter of the Neutral Hard Spheres In the previous section we only considered the critical reduced dipole moment, but not the entire coexistence curve. We now describe our results for the coexistence curves. Gibbs ensemble simulations were performed at P* = 2.0 for a total of 500 particles. We performed a series of simulations at p* = 2.1,2.2,2.3, and 2.4. T h e results are tabulated in Table 5.5 and coexistence curves are plotted in F i g . 5.1. Although p d reSt agree within the uncertainty in each box, the uncertainty is quite large. This is due to the difficulty of exchanging a dipolar hard sphere between the two boxes in the Gibbs ensemble simulation. Table 5.5: Coexistence results for an equal diameter mixture of hard-sphere dipoles and neutral hard spheres at P* = 2.0. and fi d T h e first box is rich in dipolar hard spheres. p ,n res are the residual chemical potentials of the neutral and dipolar hard spheres, reSt respectively. P (N) pres d/kT pres n/kT </>*> I II I II I II I II 2.1 0.662(6) 0.546(7) 3.83(11) 3.84(11) -1.05(30) -1.25(45) 259(3) 241(3) 2.2 0.685(3) 0.546(6) 3.86(8) 3.86(9) -0.92(39) -1.36(10) 245(3) 255(3) 3.86(14) -2.15(70) -1.94(44) 240(6) 260(6) 3.81(9) -1.70(76) — 1.90(13) 229(4) 271(4) 2.3 0.697(9) 0.546(8) 3.86(15) 2.4 0.719(6) 0.543(5) 3.81(10) T h e mixture demixes into a mixture of spheres and dipoles (box / ) and a very dilute mixture of dipolar hard spheres in neutral hard spheres (box II). B o t h the pressure and chemical potential may be calculated for pure hard spheres using the Carnahan-Starling equation of state [see Equation (2.6.2)] [67]. For sufficiently large AT, an approximate Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 129 2.6 2.4 2.2 0 Figure 5.1: P* = 2.0. 0.2 0.4 0.6 0.8 1 Coexistence curves for mixtures of neutral and dipolar hard spheres at T h e solid, dotted, dashed, and long dashed curves are for reduced diame- ters of 1.0,0.9, 0.8, and 0.7 for the neutral hard spheres. T h e horizontal bars on the data points are the error bars giving the uncertainty in the coexisting mole fractions. Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 130 expression for the chemical potential for the neutral hard spheres in the mixture is p = kT\n(x p* n T h e excess chemical potential p ex n ex is a monotonically increasing function of reduced den- sity. T h e chemical potential of the neutral hard spheres is related to the probability of inserting an additional neutral hard sphere. From the perspective of the additional hard sphere, the dipoles are simply other hard spheres. T h e dipoles can only affect the neutral sphere chemical potential due to their effect on the excluded volume due to dipole-dipole interactions. Equation (5.5.1) becomes exact (in the thermodynamic limit) as x n —> 1. Knowing the reduced density p* = 0.546 of the neutral-rich phase, the Carnahan-Starling equation gives the coexisting reduced pressure P* = 2.02 and coexisting neutral hardsphere residual chemical potential p, ,n res the simulation values in Table 5.5. = 3.84. These values compare very well with Since the simulations were performed at reduced pressure P* = 2.0, the Carnahan-Starling equation could equivalently be used to predict the neutral-rich density. Confident that we were able to map out the coexistence curve the simulations were repeated with smaller neutral hard-sphere diameters. Still working at P* = have constructed the coexistence curves at a* = 0.9,0.8, and 0.7. became easier to perform for the small, neutral hard spheres. 2.0, we Particle exchanges However, we quickly ran into convergence problems for two reasons. As the diameter difference between the components became larger, it was harder to perform identity exchanges. Also, as the critical reduced dipole moment increased, convergence to equilibrium became Simulations at a* = 0.7 and p* = slower. 2.5 were at the limits of our ability to generate converged results. For a Gibbs ensemble simulation at p* = 2.4, we show in F i g . 5.2 a sample configuration from each of the two Gibbs boxes. T h e mixture with cr* = 0.8 was demixed n Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres at this reduced dipole moment. 131 We have chosen a configuration showing two dipoles (a short chain) in the neutral-rich phase. However, during much of the simulation, the neutral-rich phase was empty of dipoles. T h e coexistence curves are shown in F i g . 5.1. F r o m the coexistence curves, we find that the neutral-rich phase remains a very dilute mixture of dipolar hard spheres in neutral hard spheres as <r* decreases. However, the mole fraction of neutral hard spheres increases in the other phase as a* decreases at fixed p*. T h e critical dipole moment p* as a function of cr* is given in Table 5.6. These values were estimated from the coexistence curves in F i g . 5.1. T h e reduced dipole moment was varied in 0.1 increments and the critical values were determined by noting when the mixture demixed. A very bold linear extrapolation suggests p* « 3 at vanishing diameter of the neutral hard spheres. W h i l e this linear extrapolation cannot be taken too seriously, the high reduced dipole moment indicates that other transitions such as orientational ordering or freezing may occur before a gas-isotropic liquid transition [34]. Table 5.6: T h e critical reduced dipole moment as a function of the diameter of the neutral hard spheres for a mixture of neutral and dipolar hard spheres. * 1.0 M 2.0 0.9 2.1 0.8 2.2 0.7 2.3 c It would be very difficult, if p* c 3, to search for the gas-isotropic liquid transition for dipolar hard spheres by computer simulation, van Leeuwen and Smit have also searched for this transition using a different limiting procedure [30]. T h e y varied a parameter in the Hamiltonian, which smoothly transformed the system from a Stockmayer fluid to dipolar soft spheres. T h e extrapolated critical reduced dipole moment for dipolar soft 132 Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres spheres is close to our value. chaining of the dipoles. 133 However, the gas-isotropic liquid was pre-empted by a T h e y argue that a m i n i m u m amount of dispersion energy is required to observe a gas-isotropic liquid transition. Recall that no such transition was found for the hard-core liquid crystal model. However, this liquid crystal model provides a much simpler model to study. There are no long-range forces in the liquid crystal model, whereas the long-range dipolar forces likely give rise to chaining structures found in dipolar fluids. Also, the relevant temperatures are much more amenable to computer simulation in the liquid crystal case. In contrast, the demixing results have only shown that the gas-isotropic liquid transition for dipolar hard spheres, if it exists, occurs at a very low temperature. Chapter 6 Conclusions We have presented a model for an E R fluid in which the ions were treated explicitly. It was shown that the model could reproduce much of the behavior of real E R fluids. In particular, the electrorheological effect, which is the chaining of E R particles in an electric field, was reproduced. Our model differs from previous E R fluid models, which were based on point polarizable E R particles. In these models, it was assumed that chaining occurred at a field at which the net polarization varied linearly with the field. Our model showed that chaining occurred far from this linear response regime. It was also shown that non-spherical E R particles enhance the E R effect due to allowing for a larger dipole moment than a spherical E R particle of equal volume. In the prolate case, the symmetry axis essentially tracks the induced dipole which points along the electric field. Thus we can form an orientationally ordered system using an electric field even for slightly prolate particles. Without a field, prolate ellipsoids will only form a nematic for large (b/a > 3) elongations. Oblate ellipsoids behave quite differently since their symmetry axes align perpendicular to the field to maximize the dipole moment. It was also shown that it is possible to form a biaxial system from oblate ellipsoids for small length to breadth ratios (b/a < 1/3) at sufficiently high density and field. T h e spherical and ellipsoidal E R models are quite simple in that they involve only hard-core repulsions and ion-ion interactions. Simulations are limited by the number of E R particles and number of ions per E R particle which are computationally manageable as a simulation involves a total of (1 + Ni )N on 134 particles. We note that a model with Chapter 6. Conclusions 135 soft-core repulsions would be more amenable to molecular dynamics simulation methods. W i t h soft-core E R particles and M D simulation, we could get true dynamical information. Furthermore, the M D simulations could be performed under shear (non-equilibrium M D ) , in close analogy with the experimental set-up for measuring the apparent shear viscosity. We have performed a thorough investigation of some simple liquid crystal models using computer simulation. For Lennard-Jones particles modified by an anisotropic interaction [given in Equation (4.2.2)] with anisotropy parameter A = 0.15, simulations along the P* = 0.1 isobar. we have performed N P T Luckhurst and Romano have also studied this system at the same reduced pressure using N P T Monte Carlo simulations [78]. report an IN transition at T* ft 0.89. They We also find a transition to a nematic upon cooling the isotropic liquid, albeit at a lower reduced temperature of T* ft 0.75. However, it was found that the fee lattice melted at a higher reduced temperature, directly to an isotropic phase. T h e solid may have remained metastable beyond the true melting temperature, although the melting temperature from simulation showed little N dependence. Thus we cannot conclude from the simulations whether there is a thermodynamically stable nematic phase along the P* = 0.1 isobar at A = 0.15. Gibbs ensemble results for the gas-isotropic coexistence at A = 0.15 show that the coexistence is predominantly determined by the Lennard-Jones interaction. was true when the anisotropy parameter was increased to A = 0.3. T h e same We have used the Gibbs method to determine the gas-nematic coexistence and the IN coexistence at A = 0.3. Due to the stronger anisotropic interaction, the IN transition was clearly first order. Using a novel method for biasing particle insertions, we demonstrated that the Gibbs method could be used to determine IN coexistence. T h e main idea is that the orientation of the particle for a trial insertion is selected from the orientational distribution function of the receiving box. This method is of general applicability to systems with orientational degrees of freedom. It should work provided the density of each phase is not so high that Chapter 6. Conclusions the particle exchange rate is too low. 136 T h e density difference between the two phases must also be sufficiently large for the biasing method to work. T h e biasing method does increase the probability of an accepted exchange, but its main advantage is in not perturbing the orientational order of the nematic phase. Computationally, it only involves collecting the orientational distribution function of each box throughout the simulation, as well as using some subroutines to sample from these distribution functions. We have also investigated hard-sphere particles modified by the same anisotropic interaction. T h e anisotropy parameter was fixed at A = 0.15. Without the Lennard-Jones interaction, we were unable to find a gas-isotropic liquid coexistence. We suspect that no such coexistence exists. T h e stabilizing energy due to the condensation of the gas to an isotropic liquid phase does not compensate for the loss in entropy, thus the transition does not occur. However, we have determined an isotropic fluid-nematic coexistence using the Gibbs method. T h e conventional Gibbs method and the Gibbs method with rotationally biased insertions were compared using this system. Using the biasing method, we were able to perform simulations at temperatures much closer to the critical temperature than with the conventional Gibbs method. Otherwise, both methods gave identical results. T h e biasing method was applied only to one class of simple liquid crystal models. It would be interesting to use the Gibbs method with biased particle insertions to simulate the IN coexistence for hard ellipsoids of revolution, Gay-Berne particles [81], or particles interacting via a generalized Gaussian overlap potential [92]. In principle, we see no reason why this would not be successful, although the density difference between the two phases for each particular model may cause difficulties. It would also be interesting to revisit the theoretical results for the liquid crystal models now that we have a much clearer idea of the phase diagram. T h e Gibbs method was used to show that the predicted demixing transition for neutral and dipolar hard spheres does indeed occur. Our results are in qualitative agreement with Chapter 6. Conclusions 137 the integral equation theory predictions, although the theory underestimates the critical temperature. We have determined the range of temperatures and pressures for which the demixing phase transition occurs. Integral equation theory also predicts a gas-isotropic liquid transition for dipolar hard spheres, although computer simulations have failed to locate such a transition. In the limit of vanishing neutral hard spheres, a demixing transition resembles a gas-isotropic liquid transition for dipolar hard spheres. Thus we investigated the dependence of the critical temperature for the demixing transition as a function of the diameter of the neutral hard spheres at fixed pressure. A n extrapolation of the critical temperature to zero diameter gives an estimate of the gas-isotropic transition temperature, if such a transition exists. Although we were only able to construct the coexistence curve for a few diameters of the neutral hard spheres, an extrapolation to zero diameter indicates that the critical temperature would occur at quite a low temperature. This low temperature implies that it is difficult to find the gas-isotropic liquid transition via computer simulation due to convergence problems. We have noted that the hard-core liquid crystal model provides a more clear-cut example of a system without a distinct isotropic liquid phase. It would be interesting to investigate the demixing transition for a mixture of hard sphere ions and neutral hard spheres. T h e reason is that, unlike dipolar hard spheres, hard sphere ions are known to undergo a condensation [93,94]. Such work is underway by others [89]. One can also approach the gas-isotropic liquid phase transition using other limiting procedures. A n example is to have a dipole in the center of two hard spheres separated by a distance Ijcr. T h e two hard spheres may penetrate one another, but not other pairs. As / —> 0, one recovers the dipolar hard sphere case. Preliminary results show that this system has a gas-isotropic liquid transition down to I/a = 0.1 [89]. Bibliography [1] C . S. Powell, Sci. A m e r . A p r i l , (1994). [2] S. W . Depp and W . E . Howard, Sci. A m e r . M a r c h , (1993). [3] D . L . Klass and T . W . Martinek, J . A p p l . Phys. 38, 67 (1967). [4] H . Conrad, Y . Chen, and A . F . Sprecher, Int. J . M o d . Phys. B 6, 2575 (1992). [5] D . R. G a m o t a and F . E . Filisko, Int. J . M o d . Phys. B 6, 2595 (1992). [6] W . M . Winslow, J . A p p l . Phys. 2 0 , 1137 (1949). [7] H . Block and J . P. Kelley, J . Phys. D : A p p l . Phys. 2 1 , 1661 (1988). [8] T . C . Halsey and W . Toor, Phys. Rev. Lett. 65, 2820 (1990). [9] R . Tao and J . M . Sun, Phys. Rev. Lett. 67, 398 (1991). [10] R. Tao and J . M . Sun, Phys. Rev. A 44, R6181 (1991). [11] R . T . Bonnecaze and J . F . Brady, J . C h e m . Phys. 96, 2183 (1992). [12] M . W . Evans and D . M . Heyes, Phys. C h e m . 95, 5287 (1991). [13] D . J . Klingenberg, F . V a n Swol, and C . F . Zukoski, J . C h e m . Phys. 94, 6160 (1991) [14] D . J . Klingenberg, F . V a n Swol, and C . F . Zukoski, J . C h e m . Phys. 94, 6170 (1991) [15] Y . Chen, A . F . Sprecher, and H . Conrad, J . A p p l . Phys. 70, 6796 (1991). [16] T . C . Halsey, J . E . M a r t i n , and D . Adolf, Phys. Rev. Lett. 68 1519 (1992). [17] J . E . M a r t i n and R. A . Anderson, J . C h e m . Phys. 104, 4814 (1996). [18] D . Wei and G . N . Patey, Phys. Rev. Lett. 68, 2043 (1992). [19] D . Wei and G . N . Patey, Phys. Rev. E 47, 2954 (1993). [20] R . Eppenga and D . Frenkel, Molec. Phys. 52, 1303 (1984). [21] A . Stroobants, H . N . W . Lekkerkerker, and D . Frenkel, Phys. Rev. A 36, 2929 (1987) 138 Bibliography 139 [22] D . Frenkel, Molec. Phys. 60, 1 (1987). [23] M . P. Allen, L i q . Crystals 8, 499 (1990). [24] G . J . Zarragoicoechea, D . Levesque, and J . J . Weis, Molec. Phys. 75, 989 (1992). [25] J . A . C . Veerman and D . Frenkel, Phys. Rev. A 45, 5632 (1992). [26] J . J . Weis, D . Levesque, and G . J . Zarragoicoechea, Phys. Rev. Lett. 69, 913 (1992). [27] P. J . Collings, Liquid Crystals (Princeton University, Princeton, 1990). [28] S. Chandrasekhar, Liquid Crystals, 2nd Edition (Cambridge University, New York, 1992). [29] J . - M . Caillol, J . C h e m . Phys. 98, 9835 (1993). [30] M . E . van Leeuwen and B . Smit, Phys. Rev. Lett. 71, 3991 (1993). [31] X . S. Chen, M . Kasch, and F . Forstmann, Phys. Rev. Lett. 67, 2674 (1991). [32] X . S. Chen and F . Forstmann, Molec. Phys, 76, 1203 (1992). [33] M . S. Wertheim, J . C h e m . Phys. 55, 4291 (1971). [34] M . Kasch and F . Forstmann, J . C h e m . Phys. 99, 3037 (1993). [35] J . - P . Hansen and I. R. M c D o n a l d , Theory of Simple Liquids (Academic, San Diego, 1986) . [36] Simulation of Liquids and Solids, edited by B . Ciccotti, D . Frenkel, and I. R . M c Donald (North-Holland, Amsterdam, 1987). [37] M . P. Allen and D . J . Tildesley, Computer Simulation of Liquids (Oxford, New York, 1987) . [38] Computer Simulation in Materials Science, edited by M . Meyer and V . Pontikis (Kluwer, Dordrecht, 1991). [39] T . Cagin and B . M . Pettitt, M o i . Phys. 98A, 169 (1991). [40] C . L o and B . Palmer, J . C h e m . Phys. 102, 925 (1995). [41] D . J . Evans and G . P. Morriss, J . C h e m . Phys. 77, 63 (1983). [42] D . J . Evans and G . P. Morriss, Phys. Lett. 98A, 433 (1983). Bibliography 140 [43] D . J . Evans and G . P. Morriss, Statistical Mechanics of Nonequilibrium Liquids (Academic, San Diego, 1990). [44] S. L . Carnie and G . M . T o m e , A d v . C h e m . Phys. 56, 141 (1984). [45] H . M e t i u , Y . - T . L u , and Z . Zhang, Science 255, 1088 (1992). [46] J . S. Rowlinson and B . W i d o m , Molecular Theory of the Capillarity (Clarendon, Oxford, 1982). [47] S. Consta and R . K a p r a l , J . C h e m . Phys. 104, 4581 (1996). [48] D . W . Brenner, Phys. Rev. B 42, 9458 (1990). [49] K . Binder, in Computational Modeling of Polymers, edited by J . Bicerano (Marcel Dekker, New York, 1989). [50] M . P. Allen, M o l . Phys. 40, 1073 (1980). [51] J . P. Valleau, Computer Simulation in Materials Science, edited by M . Meyer and V . Pontikis (Kluwer, Dordrecht, 1991), Chap. 1.4. [52] J . P. Valleau and S. G . Whittington, in Statistical Mechanics, Part A, edited by B . J . Berne (Plenum, New York, 1977), Chap. 4. [53] J . M . Hammersley and D . C . Handscomb, Monte Carlo Methods, (Spottiswoode, Ballantyne & C o , London, 1965).' [54] N . Metropolis, A . W . Rosenbluth, M . N . Rosenbluth, A . H . Teller, and E . Teller, J . C h e m . Phys. 21, 1087 (1953). [55] D . Frenkel, Computer Simulation in Materials Science, edited by M . Meyer and V . Pontikis (Kluwer, Dordrecht, 1991), p. 110. [56] A . Z . Panagiotopoulos, Molec. Phys. 61, 813 (1987). [57] A . Z . Panagiotopoulos, N . Quirke, M . Stapleton, and D . J . Tildesley, Molec. Phys. 63, 527 (1988). [58] B . Smit, P h . De Smedt, and D . Frenkel, Molec. Phys. 68, 931 (1989). [59] B . Smit and D . Frenkel, Molec. Phys. 68, 951 (1989). [60] A . Z . Panagiotopoulos, Int. J . Thermophys. 10, 447 (1989). [61] L . F . R u l l , G . Jackson, and B . Smit, M o l . Phys. 85, 435 (1995). Bibliography 141 F . H . Stillinger and T . A . Weber, Phys. Rev. B 31, 5262 (1985). P. P. Ewald, A n n . Phys. 64, 253 (1921). S. W . De Leeuw, J . W . Perram, and E . R. Smith, Proc. R . Soc. L o n d . A 373, 27 (1980); 373, 57 (1980). B . W i d o m , J . C h e m . Phys. 39, 2808 (1963). J . W . Perram, M . S. Werthiem, J . L . Lebowitz, and G . 0 . Williams , C h e m . Phys. Lett. 105, 277 (1984). N . F . Carnahan and K . E . Starling, J . C h e m . Phys. 51, 635 (1969). L . C . Davis, Phys. Rev. A 46, R719 (1992). T . Chen, R . N . Zitter, and R . Tao, Phys. Rev. Lett. 68, 2555 (1992). R. Tao, Int. J . M o d . Phys. B 6, 2635 (1992). B . C . X u and K . C . Hass, J . C h e m . Phys. 98, 2258 (1992). H . Zhang and M . W i d o m , Phys. Rev. E 51, 2099 (1995). L . Eyges, The Classical Electromagnetic Field (Dover, New York, 1972). G . J . Zarragoicoechea, D . Levesque, and J . J . Weis, Molec. Phys. 78, 1475 (1993). P. A . Lebwohl and G . Lasher, Phys. Rev. A 6, 426 (1972); 7, 2222 (1973). C . Zannoni, The Molecular Physics of Liquid Crystals, edited by G . R . Luckhurst and G . W . Gray (Academic, London, 1979), Chap. 9. W . Maier and A . Saupe, Z . Naturf. A 13, 564 (1958); 14, 882 (1959); 15, 287 (1960). G . R . Luckhurst and S. Romano, Proc. R. Soc. A 373, 111 (1980). A . Perera, P. G . Kusalik, and G . N . Patey, Molec. Phys. 60, 77 (1987). A . Perera, G . N . Patey, and J . J . Weis, J . C h e m . Phys. 89, 6941 (1988). E . D e Miguel, L . F . R u l l , M . K . Chalam, and K . E . Gubbins, Molec. Phys. 71, 1223 (1990). D . Frenkel, in Phase Transitions in Liquid Crystals, edited by S. Martellucci and A . N . Chester (Plenun, New York, 1992), Chap. 5. P. H . Fries and G . N . Patey, J . C h e m . Phys. 82, 429 (1985). Bibliography 142 [84] E . De Miguel, L . F . R u l l , M . K . C h a l a m , K . E . Gubbins, and F . V a n Swol, M o l . Phys. 72, 593 (1991). [85] F . R . Cracknell, D . Nicholson, and N . G . Parsonage, Molec. Phys. 71, 931 (1990). [86] P. G . de Gennes, The Physics of Liquid Crystals (Oxford University, New York, 1974). [87] R . D . Mountain and A . H . Harvey, J . C h e m . Phys. 9 4 , 2238 (1991). [88] J . G . A m a r , Molec. Phys. 67, 739 (1989). [89] J . Shelley and G . N . Patey, unpublished. [90] P. G . Kusalik, J . C h e m . Phys. 9 3 , 3520 (1990). [91] D . Levesque and J . J . Weis, Phys. Rev. E 4 9 , 5131 (1994). [92] G . A y t o n and G . N . Patey, J . C h e m . Phys. 102, 9040 (1995). [93] G . Orkoulas and A . Z . Panagiotopoulos, J . C h e m . Phys. 101, 1452 (1994). [94] J . M . Caillol, J . C h e m . Phys. 100, 2161 (1994). [95] H . Goldstein, Classical Mechanics, 2nd Edition (Addison-Wesley, D o n Mills, 1980). Appendix A Rotation of a Particle Particles with orientational degrees of freedom must be randomly rotated as well as translated. We now describe how such rotations are performed. uniaxial particle is described by a unit vector it. T h e orientation of a To perform a rotation, the rotation matrix A, which rotates a unit vector along the z axis to the current unit vector tt, is constructed [95]. T h e unit vector ii is given by ^ cos 6 sin u sin <j> cos 9 ^ sin 9 (A.l) 9 where 9, (f> are the usual spherical coordinate system angles. T h e matrix A can be constructed by first rotating z through an angle 9 about the y axis. T h i s is performed by the matrix ( cos 9 0 0 1 0 0 cos A, = V — sin 9 sin6> ^ (A.2) 9 T h e resulting unit vector is then rotated through an angle (f> about the z axis. T h i s is performed by the matrix AA = cos (j) sin (j) — 0 sin cos 0 143 (j) (f) 0 0 1 \ (A.3) Appendix A. Rotation of a Particle B o t h Ag and A^ 144 may be constructed by considering the effect each must have on the column unit vectors x, y, z. T h e matrix A is then the product of A$ and ^ cos (j) cos A = AA 6 6 9 — sin < ^ > cos # cos <f> sin cos <f> sin <> / sin — sin ^ A unit vector it sin <f> 0 cos 9 Ag, 9 ^ (A.4) 6 is then generated in a spherical cap about the z axis on the unit sphere. z This unit vector is given by the angles 9 Z cf) z where Ri and R 2 = 7rPdA0 (A.5a) = (A.5b) 2TTR , 2 ave pseudo-random numbers uniformly selected from the unit interval and A9 € [0,1] determines the m a x i m u m angular step size. T h e unit vector is then ^ cos S z sin 9 ^ Z sin <j> sin 9 z cos T h e trial unit vector u Z (A.6) 9 Z is it = Ait, (A.7) Appendix B C h o o s i n g a N e w Ion P o s i t i o n Inside an E l l i p s o i d Confining spherical ions to the inside surface of an ellipsoid of revolution poses some difficulties due simply to geometrical considerations. Some of these difficulties are ad- dressed in this Appendix. A n ion must fit entirely within the ellipsoid. T h i s is a restriction on the size of the ion diameter <7,- . Consider a slice through the x = 0 plane of an ON ellipsoid. This gives the equation of an ellipse y 2 z — H 2 b 2 a' (B.l) = 1• For an oblate ellipsoid (b < a), this is an ellipse with foci at ±\J a — b . T h e m a x i m u m 2 2 radius of the ion is contained between one focus and the nearest point of the ellipse. It follows that b\ 2 ion,max - (6/a)" 1/3 1- (B.2a) - w For a prolate ellipsoid (b > a) the result is ion,max = (Va)" (B.2b) 1 7 3 T h i s restriction on ion size only becomes a concern for significantly non-spherical particles. For example, for an oblate ellipsoid with b/a = 1/3, the m a x i m u m value is <r* « 0.082. im E v e n if the ion will fit inside the ellipsoid, too large a value of <7, „ should be avoided 0 for the following reason. T o place an ion in an ellipsoid, a point is randomly chosen on the ellipsoid surface. This point is then the contact point of the ion with the inner surface 145 Appendix B. Choosing a New Ion Position Inside an Ellipsoid of the ellipsoid. 146 T h e probability of choosing a point on the ellipsoid surface is equal to the probability of choosing any other point on the surface. are determined by uniformly sampling the ellipsoid surface. Thus the ion positions Strictly speaking, it is the surface on which the ion centers of mass are confined that should be uniformly sampled. T h i s is a surface within the ellipsoid whose normal is everywhere c - / 2 from the ellipsoid t o n surface. T h i s inner surface is not an ellipsoid. Rather, it has a more complicated defining equation. If c,- on is too large, sampling from the ellipsoid to give the contact point will bias the ion center of mass towards the thin ends of the ellipsoid. a* expected to be a problem in practice for the small values of on T h i s biasing is not used in our study. A confined ion must be moved such that, in the absence of the B o l t z m a n n factor, all positions are equally probable. Such a condition satisfies the requirement of microscopic reversibility. Thus, we now describe how to generate a random position (x, y, z) uniformly chosen from the ellipsoid surface 4+4+S=i GT For an oblate ellipsoid GT A iate ob For a prolate ellipsoid 0^ (b < a) with eccentricity e = 1 is T Jo /2 sin = 47ra = 27ra + 7 r 6 e 2 (R3) 2 2 — (b/a) , the area of the ellipsoid 6(1 - e s i n 2 2 9) d8 1/2 2 - l n ^ . 1 —e (B.4) (b > a) with eccentricity e = 1 — (a/b) , the area of the ellipsoid is 2 /•7r/2 Aprdate = 4nab Jo = 2 7 r a + 27ra6 „ cos 6(1 — e sin 8)^ d8 2 , arcsin e 2 e 2 . (B.5 . To choose a random point on an ellipsoid, first the z coordinate is chosen. B y symmetry, only the case z > 0 needs to be considered. We define p(z) p(z)dz = (probability that a random point through Appendix B. Choosing a New Ion Position Inside an Ellipsoid 147 on the ellipsoid has a z component between z and z + dz) oc area of ellipsoid between z and z + dz . (B.6) T h e normalization condition / Jo (B.7) p(z)dz requires 1 p(z)dz = — (area of ellipsoid between z and z + dz) (B.8) Ai/2 where A / is half the ellipsoid area. We associate a random number R uniformly chosen T 2 from the unit interval with the probability (probability that z is between 0 and z) = R = — f A1/2 / 2 Je sin 0(1 - e sin Bfl de 2 2 2 (B.9) Integrating gives i Oyjl-e R = 7ror cos 0\fl n 2 o o sin 6 + 2 n I e , / e - In cos 6 + Vl - e sin 0 2 2 x (B.lOa) A 1/2 where ^ = (1 - • (B.lOb) A similar treatment of the prolate case gives R irab Ai/2 sin 0\/l — e smO + - arcsin(e sin 9) 2 This gives R as a function f(z/b) of zjb. function f (R) _1 We must invert this function to get z/b as a of R. This was done by tabulating R for a series of values of z/b and performing linear interpolation between points. ellipsoid, two random numbers £i and £ 2 ) To get a random point (x,y,z) on an where £ € [—1,1], are generated. T h e random Appendix B. Choosing a New Ion Position Inside an Ellipsoid 148 point is then given by z = 9 = frsignfo)/- ^!), 2TT|6| , / x = = (B.12b) zV 2 /2 (l-^-J / y (B.12a) 1 a cos 9 , (B.12c) zV 2 1-77 /2 a s ' m d - ( B 1 2 d ) Often, it is more convenient to move an ion within a region centered about the ion instead of to a point anywhere on the ellipsoid surface. To satisfy microscopic reversibility in this case requires that the area of the region is independent of the ion position. Thus it is necessary to generate a random position (x', y', z') uniformly chosen from a neighborhood about the ion on the ellipsoid surface. Given a point (x,y,z) on the ellipsoid, we calculate 8 = arctan (j-^j R = sign , 0e[O,27r] • (B.13a) (B.13b) We then choose two random numbers £i and £2 and calculate R' = R + ^Ar. , 9' = 9 + n£ Ar 2 where A r , (x',y',z'). o n R' ion ion €[-1,1] , € [0,1] is the m a x i m u m step size. (B.14a) (B.14b) R' and 9' are then used to generate T h e new point lies within a neighborhood of (x,y,z). B y construction, an interval AR taken from [—1,1] corresponds to a strip around the ellipsoid of the same area as an interval taken from any other range between [—1,1]. This implies that microscopic reversibility is satisfied, as required. Appendix C R o t a t i o n of Charge C o o r d i n a t e s T h e algorithm described for rotating a particle assumed a uniaxial particle. A n E R particle has hard-sphere ions confined to the E R particle's inner surface. One must be careful in rotating an E R particle as the ions break the uniaxial symmetry. T h i s is true for both a spherical particle, which contains an induced dipole, and an ellipsoidal particle, which has both an induced dipole and a symmetry axis. A spherical E R particle is rotated about its dipole. A n ellipsoidal E R particle is rotated about its dipole and symmetry axis. T h e particular axis about which to rotate is randomly chosen for each trial rotation. In this Appendix, we describe how the ions are rotated along with the E R particle. T h e ion position r ; o n is measured with respect to the center of the E R particle. T h e E R particle is rotated about its center. We have the current unit vector w, which may be determined in terms of the angles 6, (j>. We rotate u to lie along the z axis. Recall that to rotate the unit vector z along the z axis to u, we constructed a matrix A = A^Ag. A - 1 . Thus we want the inverse matrix T h i s is readily constructed from the inverse matrices cos (j) — sin V sin <f> \0 cj) cos <f> 0 o 149 Appendix C. Rotation of Charge Coordinates 150 and I cos An = 1 One can readily check that 9 0 - sin 1 0 0 V A^A^ 1 sin 9 0 • (C2) 9 cos I and AftAa 9 ^ = I. 1 Thus the inverse matrix is given by A- 1 One may note that A^ 1 = A; A; 1 1 . (C.3) is the transpose of A^ and similarly for If we only use the matrix A , Ag . 1 the E R particle ends up being rotated through an - 1 angle — <p in the z = 0 plane. To remove this rotation about the symmetry axis or induced dipole, we rotate the E R particle in its current position along the z axis through an angle <p about the z axis to cancel the undesired rotation. This is performed by the matrix A^. T h e entire operation from rotating the particle from its current position u to the z axis is performed by z = A+A^A; * . 1 (C.4) A s part of the rotation move, we may desire to rotate the particle through an angle tp £ [—7r, 7r] about the symmetry axis or induced dipole. We randomly choose an angle ip uniformly chosen from the interval [—irAip, irAip) where m a x i m u m angle. G [0,1] determines the This is set in combination with the m a x i m u m rotation to give an acceptance rate of roughly a half for trial rotations. T h i s rotation is performed by the matrix ^ cosip — smip sin xp cos ip 0 0 0 1 B ,i, = V 0 ^ (C.5) T h i s is just the matrix A^ except that the angle <p is replaced with ip. T h e operation is B,i,z (C.6) Appendix C. Rotation of Charge Coordinates 151 Although this does not appear to be a very exciting operation, recall that it is the ions that we are really interested in rotating. We then rotate the E R particle to its new trial position, given by trial angles 9 , <fi . Again, we remove the unwanted rotation. T h e entire operation of taking the particle along the z axis to the trial position it is performed by it = ApAe'A^z . (C.7) Collecting all these results, the operation for rotating it to it is given by it' = ApAg,A^}B^A^Ag 1 A^ 1 it . (C.8) This operation also rotates the ions from their original positions to their new trial positions.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A Monte Carlo study of fluids with orientational degrees...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A Monte Carlo study of fluids with orientational degrees of freedom Blair, Mark James 1996
pdf
Page Metadata
Item Metadata
Title | A Monte Carlo study of fluids with orientational degrees of freedom |
Creator |
Blair, Mark James |
Date Issued | 1996 |
Description | We have modeled an electrorheological (ER) fluid as hard-sphere particles, each with
smaller hard-sphere ions constrained to roll on the sphere's inside surface. When the
model E R fluid is placed in an electric field, each particle becomes polarized (due to rearrangement
of the ions confined within the particle) with a dipole moment depending on
both the field and interactions with its neighbors. Using NVT Monte Carlo simulations,
we have shown that our model can display chain formation as seen in real ER fluids.
Chaining occurred at a field where the net moment no longer varied linearly with the
field. This model was extended to include particles shaped as ellipsoids of revolution. In
the prolate case, slightly non-spherical particles were readily ordered by the field. In the
oblate case, the induced dipole is roughly perpendicular to the symmetry axis. Oblate
particles may then form a biaxial phase in an applied field.
NPT and Gibbs ensemble Monte Carlo simulations were performed for spherical particles
modified by an anisotropic potential of the form —4A£( |
Extent | 9290760 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-03-19 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0059608 |
URI | http://hdl.handle.net/2429/6260 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1996-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1996-147282.pdf [ 8.86MB ]
- Metadata
- JSON: 831-1.0059608.json
- JSON-LD: 831-1.0059608-ld.json
- RDF/XML (Pretty): 831-1.0059608-rdf.xml
- RDF/JSON: 831-1.0059608-rdf.json
- Turtle: 831-1.0059608-turtle.txt
- N-Triples: 831-1.0059608-rdf-ntriples.txt
- Original Record: 831-1.0059608-source.json
- Full Text
- 831-1.0059608-fulltext.txt
- Citation
- 831-1.0059608.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0059608/manifest