A M O N T E C A R L O S T U D Y O F F L U I D S W I T H O R I E N T A T I O N A L D E G R E E S O F 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 T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S ( D E P A R T M E N T O F C H E M I S T R Y ) We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H 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 my 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) The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V 6 T 1Z1 Date: Abstract 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 rear-rangement 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 par-ticles modified by an anisotropic potential of the form —4A£(<r / r ) 6 P 2 ( cos7). 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. At 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. Our 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 vi i List of Figures vii i Table of Symbols x Acknowledgement xv 1 Introduction 1 2 Computer Simulation 8 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 Hard Spheres 40 2.7 Orientational Degrees of Freedom and Order Parameters 41 3 Mode l Electrorheological Fluid 46 3.1 Introduction 46 iv 3.2 The Model 49 3.3 Computational Details 51 3.4 Determination of Ordering and Structure 53 3.5 Model 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 Fluid of Oblate Ellipsoids as a Function of Density 82 3.7.4 Summary 85 4 Liquid Crystal Models 87 4.1 Introduction 87 4.2 Model 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 Summary 115 5 Mixtures of Neutral and Dipolar Hard Spheres 118 5.1 Introduction 118 5.2 Model and Computational Details 120 5.3 Computational Tests 123 5.4 The Critical Dipole Moment as a Function of Pressure 126 5.5 The Critical Dipole Moment as a Function of Diameter of the Neutral Hard Spheres 128 v 6 Conclusions 134 Bibliography 138 Appendices 143 A Rotation of a Particle 143 B Choosing a New Ion Position Inside an Ellipsoid 145 C Rotation of Charge Coordinates 149 v i Lis 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* o n = 0.05, and E* = 100 . . . 59 3.6 Results for b/a = 1, />* = 0.5, = 0.5, < n = 0.1, and E* = 100 60 3.7 Results for p* = 0.75, q* = 0.35, a*on = 0.05, and E* = 40 62 3.8 Results for 6/a = 1, p* = 0.75, q* = 0.35, and <r*on = 0.05 63 3.9 Results for b/a = 3, p* = 0.75, q* = 0.35, and a*on = 0.05 75 3.10 Results for p* = 0.75, q* - 0.35, a*on - 0.05, and E* = 40 77 3.11 Results for b/a = 0.333, q> = 0.35, a* o n = 0.05, and E* = 100 85 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 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 16 3.1 Sample configurations for N = 128 and N = 120with b/a = 1, p* = 0.5, q* = 0.5, and a*ion = 0.1 at E* = 100 61 3.2 The mean square displacement for b/a = 1, p* — 0.75, q* = 0.35, and < n = 0-05 64 3.3 (a) (/j.*) and (b) (Pj) versus E* for b/a = 1, p* = 0.75, q* = 0.35, and < n = 0-05 66 3.4 (a) gn(z) and (b) g(r) for b/a = 1, p* = 0.75, q* = 0.35, and cr*on = 0.05. . 68 3.5 Dipole moment distributions for b/a = 1, p* = 0.75, q* = 0.35, and C = 0-05 69 3.6 Sample configurations for b/a — 1, p* = 0.75, = 0.35, and a*on = 0.05 at (a) E* = 50 and (b) E* = 60 70 3.7 (P 2) versus E* for 6/a = 3, p* = 0.75, ? * = 0.35, and a*on = 0.05 72 3.8 0H(z; u) for 6/a = 3, p* = 0.75, = 0.35, and a*on = 0.05 73 3.9 Sample configurations for b/a = 3, p* = 0.75, 9* = 0.35, and a*on = 0.05 at (a) E* = 0, (b) E* = 10, and (c) E* = 100 74 3.10 l/D vs. £ * 2 for (a) b/a = 1 and (b) b/a = 3 at p* = 0.75, 4* = 0.35, and < n = 0.05 76 3.11 (a) p \ (b) (P 2 ) , and (c) (cos 7) versus b/a for / = 0.75, q* = 0.35, a*on = 0.05 at £ * = 40 78 3.12 g\\(z] u) for p* = 0.75, q* = 0.35, < o r i = 0.05 at £ * = 40 80 viii 3.13 Sample configuration for b/a = 0.333, p* = 0.75, q* = 0.35, and a*on = 0.05 at E* = 40 81 3.14 p/pmax versus E* for a single particle for q* = 0.35, and a*on = 0.05. . . . 83 3.15 (a) g\\(z;n) and (b) g\\(z;u) for b/a = 0.333, q* = 0.35, and a*on = 0.05 at E* = 100 84 4.1 (a) p* and (b) (P 2) versus T * at P* = 0.1 for soft model with A = 0.15. . 96 4.2 (a) g(r) and (b) g220{r) at P* = 0.1 for soft model with A = 0.15 98 4.3 Mean square displacements at P* = 0.1 for soft model with A = 0.15. . . 99 4.4 p* versus T* for N = (a) 32, (b) 64, (c) 108, and (d) 150 for soft model with A = 0.15 101 4.5 (P 2) versus T* for N = (a) 32, (b) 64, (c) 108, and (d) 150 for soft model 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. . 104 4.7 Gas-liquid coexistence curve for soft model with A = 0.15 108 4.8 (a) p* and (b) (P 2) versus T* for soft model with A = 0.3 110 4.9 Phase diagram for soft model with A = 0.3 I l l 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, xn = 0.8, P* = 2.0, and p* = 2.4 132 IX Table of Symbols 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(r) Radial distribution function. g220(r) Term in expansion of full angle dependent distribution func-9\\(z) tion. Distribution function measuring ordering alonj 1 induced 9\\(zlu) dipole of spherical E R particle. Distribution function measuring ordering along symmetry 9\\{z\n) axis ellipsoidal E R particle. Distribution function measuring ordering alonj I induced G dipole of ellipsoidal E R particle. Gibbs free energy fiN. h Planck's constant. k Boltzmann's constant. K Total kinetic energy. L Box length. m Particle mass. X M Net moment, or total polarization, of the system. n Integer component lattice vector. nmax Maximum length of vector n included in lattice vector sum. N Number of particles. Nion Number of ions confined in each E R particle. P Pressure. P1 Polarization order parameter. P2 Nematic order parameter. P2N Nematic order parameter estimated from largest eigenvalue of ordering matrix Q. P2(a;) Second order Legendre polynomial. q Ion charge. Q Classical partition function. Q Ordering matrix. Q22 Biaxial order parameter. r Particle position. rcut Spherical cut-off separation for pair potential. nj Pair separation (under minimum image convention). non 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.). u a(12) Anisotropic contribution to pair potential. UDD Dipole-dipole pair potential. x i UHS Hard sphere pair potential. U{ON Coulombic (ion-ion) pair potential. ULJ Lennard-Jones pair potential. U Total potential energy. UF Total field energy for E R fluid ( - M • E). Ui Total interaction energy for E R fluid (U — Us). Ucor 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. Zs Classical configuration integral in scaled coordinates. a Convergence parameter in Ewald sum. f3 1/kT 7 Shear strain. 7 Strain rate. 8 Shell size for calculating hard core contribution to the pres-sure. A r Maximum translational step size along each coordinate axis. A r i o n Maximum 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 Maximum volume change. AO Maximum angular step size. Aip Maximum 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. e0 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 Thermal wavelength. p Chemical potential. /j, Dipole (scalar quantity is the dipole moment). pex Excess chemical potential p — pid-pid Ideal chemical potential kT In pA3. pmax Maximum possible induced dipole in E R particle. pres Residual chemical potential p — kTln(A/a)3. nd Limiting distribution for Markov chain. p Number density N/V. a Particle diameter. Xlll Gion Diameter of ion confined in E R particle, r Shear stress. r 0 Yield stress. f2 Particle orientation either in spherical polar coordinates (j>) or Euler angles (9, <f), ip). (...) Indicates an ensemble average. * Superscript indicates a reduced quantity. x i v Acknowledgement I thank Professor Gren Patey for providing both guidance and freedom throughout my time at U B C . Many 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 Ayton, who has provided computer code and assistance with the word processing of this thesis. I must also thank Kurt 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 This thesis covers three loosely related projects. The first project was the modeling of, and investigation of the orientational ordering in, an electrorheological (ER) fluid. The second project was an investigation of the isotropic-nematic (IN) phase transition for a simple liquid crystal model. The 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 (MC) 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 orienta-tional 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. To 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. The 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. Of course, we cannot get any dynamical information from M C calculations. The truth of the matter is that liquid crystals and E R fluids are intrinsically inter-esting. 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. The 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. The 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. He took out two patents based on E R fluids before publishing a paper on the Induced Fibration of Suspensions in 1949 [6]. Winslow himself first realized the dependence of shear viscosity on the square of the electric field. With 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 minimum 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. From 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. The dipole-dipole 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. The 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 maximum 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. Non-spherical E R particles introduce a liquid crystal component to the E R fluid. The dependence of the orientation of the E R particle's induced dipole on particle shape was also studied. The 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. The 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 par-ticles, 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. Two 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. As we continue to heat the liquid crystal, the particles will start to spin in all directions. At 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. The simplest type of liquid crystal is a nematic, in which there is no long-range translational order but orientational order in one direction. The 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. The orientational order in the liquid layer need not be perpendicular to the plane. The 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). The 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 Chapter 1. Introduction 6 [18,19]. The 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 display-ing 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 rela-tively new Gibbs ensemble Monte Carlo method (Gibbs method). The Gibbs method is an extension of the conventional M C method but designed to maintain thermal, mechan-ical, 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. Much 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. This 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. The mixture undergoes a Chapter 1. Introduction 7 demixing phase transition, in which it separates into a neutral-rich phase and an isotropic, dipole-rich phase. The demixing transition provides another route by which to investigate the existence of a gas-isotropic phase for dipolar hard spheres. The 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. At 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. The 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 Computer 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 (MD) 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 Simulation of Liquids and Solids [36]. The 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 corre-lated) 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. The reason is as fol-lows. The number of particles in a system is determined by its chemical potential. When 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. The rapid growth in the calculating speed of computers encourages longer, larger, and more complicated sim-ulations. Developments in the methodology allows one to select the most convenient ensemble with which to work [36,41,42]. Wi th 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. Quantum mechanical M C and M D calcula-tions have also been done [36]. The interaction potentials between particles can be made realistic enough to both reproduce and predict experimental results [38,48]. Special tech-niques 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. Two 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 By restricting ourselves to M C simulations, we have excluded the study of the dynam-ical properties of the systems under consideration. This may seem to be an unnecessary sacrifice, considering that one can determine both the equilibrium and dynamical prop-erties 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 equilib-rium properties from such simulations, the investigation of equilibrium properties need not be so straightforward. Using the assumption that time averages are equal to ensem-ble 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 Stat is t ical Mechanics and the Monte Car lo M e t h o d Classical equilibrium statistical mechanics can be expressed entirely in terms of multi-dimensional integrals and ratios of multidimensional integrals [51,52]. The 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. The pseudo-random numbers are used to sample values in the space of the in-tegral (the region over which the integration is performed). The 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. The integral Chapter 2. Computer Simulation 11 is defined as Z(N,V,T) = J...jexp(-/3U(rN))drN, (2.1.1) where rN is an abbreviation for the positions 7*1, r 2 , ...f;v of the TV particles, /? = 1/fcT 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 configu-ration 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 com-plicated, 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 ( n s = f---fO(rN)eM-PU(rN))drN ( 2 2 ] { ' f---fexp(-f3U{rN))drN ' 1 ' ' ' 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. The 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. Each state in the Markov chain is a pos-sible configuration. The 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 7rd(rN) oc exp(-(3U(rN)) . (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 Boltzmann 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? What 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. The 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 maximum 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. The 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 ) en-semble. The 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 av-erage of the pressure equals the desired pressure P. In the N P T ensemble, the limiting distribution is ird{rN,V) oc exp [-B(U+ PV - NkTInV)] oc VNe-W+pv^ (2.1.4) Chapter 2. Computer Simulation 14 and the probability in going from one configuration to the next is min( l , exp(Aiu)), where and A V is a value randomly chosen from the interval [—AVmax, A V 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. The probability is set such that changes. AVmax is adjusted so that roughly half of the attempted volume changes are accepted. When 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. Of 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. The 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- The estimate of the mean is (2.1.5) a reasonable fraction (1/10 to 1/4) of the total number of steps are attempted volume (0) (2.1.6) The standard deviation of the mean is a. (2.1.7a) Chapter 2. Computer Simulation 15 where a0 = J^-rKO2) - (0?\ • (2.1.7b) V nb - 1 The standard deviation o~o is somewhat independent of rib, but the standard deviation of the mean am 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 am from our estimate roughly 65% of the time. We choose to report (0) ± 3crm. • We may be confident that the mean lies within 3<7m roughly 99% of the time. Reporting the larger error has two advantages. The first is that if values disagree by more than the reported error, we may be reasonably certain that the difference is meaningful. The second advantage is that in practice, all simulations were for rib — 10 blocks, so that ao ~ 3<rm. 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. The 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 The 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 Fig. 2.1. The 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. The P vs. T diagram is made up of a set of coexistence curves. Any 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. The density difference implies a discontinuity in the derivative of the Gibbs free energy G with respect to pressure, since 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. This 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. This requires the location of phase transitions. If an N V T simulation is started at a tempera-ture 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. This does not happen [55]. Instead, the system remains in a single phase. The 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. Recall that the method for performing N P T simulations was discussed in section 2.1. This corresponds more closely to a search through the P vs. T phase diagram, rather than through the T vs. p phase diagram. As the coexistence curve (2.2.1) 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. The system may remain mechanically stable, (dV/dP)x < 0 , beyond 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 simu-lations. 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 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 There is a simulation method which allows for the simulation of two phases in equi-librium with one another. Thus it is a wonderful way for mapping out phase diagrams. The three conditions for equilibrium between two phases are equal temperatures, equal pressures, and equal chemical potentials. The Gibbs ensemble Monte Carlo method is de-signed to establish these three conditions between the two simulation boxes. The Gibbs method, developed in 1987, is a relatively new simulation technique [56]. The 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 be-tween the boxes. The 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 (2.2.2) give results over a range of temperatures, densities, etc. from a single simulation [51]. Chapter 2. Computer Simulation 19 to a two phase system, one should eventually get a different phase is each of the two boxes. The phases will be in equilibrium with one another. The Gibbs method reduces to a simulation of two equivalent N V T systems for an initial density corresponding to a single phase system. The 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. The original version of the Gibbs method used macroscopic, thermodynamic argu-ments to derive the probability of accepting each type of move. The Gibbs method was put on a firmer footing using statistical mechanical arguments [57]. The partition func-tion, 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. They derive the acceptance probabilities starting from the parti-tion 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] where the superscripts / and / / label each of the two simulation boxes. The probability for the acceptance of any type of move is min( l , exp(Aiw)). For example, for a volume exp(zw) = exp In { WIN"] Chapter 2. Computer Simulation 20 exchange of AV between boxes / and 77, the expression for Aw is A „ = N * l n ( £ + £ V ) + N n l n ( Y ! L ^ _ ? A V . _ . ( 2 . 2 . 4 ) For particle exchanges (say for exchanging a particle from box / to box II), the expression for Aw is / Ni yll\ Aw = In [j^yr) ~ P&U1 - (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 nexch 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. They 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. The 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) is A w =ln (/vJTl) +ln (TVFTt) ~ ^ U * - B ~ • ( 2 - 2 " 6 ) One of the two boxes (/ or II) is chosen with equal probability. The 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. The important point is that the selection algorithm must satisfy the condition of microscopic reversibility [61]. The 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. The advantage is that one is free to choose the component for which to perform particle exchanges. The sensible choice is to pick the component for which particle exchanges are easiest. The 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. The answer is supplied by the Gibbs phase rule (which actually was, unlike the Gibbs ensemble, developed by Josiah Willard Gibbs). The phase rule says that the number of thermodynamic degrees of freedom / (independent intensive thermodynamic variables) is given by / = 2 + n c - n p , (2.2.7) where nc is the number of components and np 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 thermo-dynamic 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. Thus an N P T sim-ulation is perfectly valid, even for a pure system. 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 parti-cle 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. The benefit of exchanging neutral spheres only becomes better as the diameter of the neutral spheres shrinks. 2.3 Pair Potentials and Boundary Conditions Consider a system of TV interacting particles. The potential energy may be written as the sum over individual particles, pair interactions, triplet interactions etc. to give U(rN) = J2n(rl) + J2<rl,rJ)+ £ u(r%, r i ? r k ) . . . . (2.3.1) i i<j i<j<k 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 three-body 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. The first sum accounts for the effects of applied fields. The second sum is a sum of the pair potential for the interacting particles. For a translationally invariant system, the pair potential may be written as where = Tj — V{. The subscript ij is sometimes dropped for convenience where it is understood that we mean the pair separation. The 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 0 r < a UHs(r) = i (2.3.2) oo r > a , where a is the hard-sphere diameter, and the Lennard-Jones potential, ULj(r) = 4e CT\12 la (2.3.3) . r, where — e is the minimum 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". The pair potential may depend on the orientation of the particles. In this case, we use the notation u(12) where 12 = r, J?2 and fl represents the Euler angles describing the orientation of a particle. When 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. The simplest way remove surface effects is to impose minimum image (MI) boundary conditions. To impose M I boundary conditions, one calculates the pair separations using the MI 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 . The 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 1 / 3 is the 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 c u t , the correction to the total potential energy is [37] where p = N/V is the number density. The 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. The correction assumes that particles are uncorrelated beyond a separation of rcut. Thus, the correction is not valid for a solid phase or for rcut too small. O f course, in practice, rcut must not exceed half the box length. For some ensembles, the box size varies, so it is most convenient to always set rcut to half the box Chapter 2. Computer Simulation 25 length. In this case, rcut 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. The evaluation may be carried out as done for short-ranged potentials. For short-range potentials, minimum image (MI) boundary conditions are imposed. The problem with long-range potentials is that the minimum 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. The long-range potentials used in this thesis are the coulombic (ion-ion) potential uim(r) = ^ , (2.3.5) where q is the ion charge, and the dipolar (dipole-dipole) potential uDD(12) = ^ , (2.3.6) where / / is the dipole. The 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 0 , where eo is the permittivity of free space. The 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; The minimum image (MI) expres-sion 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. The 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) The 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. The lattice sum is a sum over integer component vectors n — ny, nz) , where nx, ny, nz = 0, ± 1 , ± 2 , . . . . (2.3.9) The lattice sum is only conditionally convergent, meaning that the result depends on the order of summation. The 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! To obtain a unique limit for a conditionally convergent sum, the method of summation (the order in which the sum is performed) must be defined. The total potential energy of the central cube may be written as a sum over approximately spherical shells of cubic images. The 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. The 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. The opposite limit is a surrounding vacuum (e' = 1), where the extra term makes its maximum contribution. Both the conducting and vacuum limits of P B C are used. Although avoiding the problems of minimum 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 ^ i j n \ iJ ' \ ^ e - v 2 n 2 / a 2 + « J2 5 — l5Zftexp(27rm-s,-)| 2 / rijto 7 r n i + / ( O y E f c - . - ] 2 } , (2-3.10) where positions and separations are in units of the box length (s,- = r , - / L and S y = rij/L). The complimentary error function is defined as erfc(x) = - = / e~* dt . (2.3.11) The 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) The potential energy has been expressed as the sum of two converging series. The 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 maximum allowed error) beyond half the box length. The maximum length nmax 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. nmax is taken sufficiently large Chapter 2. Computer Simulation 28 so that the truncation error for both sums match each other. As a final simplification, consider the sum over n. The terms have the symmetry property that the nth term equals the — nth 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. The final term in Equation (2.3.10) is proportional to the net moment, M = ( 2 - 3 - 1 4 ) t or total polarization, of the system. A computationally efficient form for the potential energy with periodic boundary conditions is Tit >\ 1 I eiic(asij) rtmax * e - 7 r 2 n 2 / a 2 + £ T , — | V] qi e x p ( 2 7 r m • S j ) a P-1 + f ( ^ M \ (2.3.15) Typically, a value of a ~ 5 is sufficient so that only the first term in the real space (complementary error function) sum need be included. Of course this implies including many more terms in the reciprocal space sum. A value of n2max = 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. The Ewald expression is readily generalized to a rectangular box. In this case, the shortest side is taken as the box length L. The box dimensions L x , L y , L z are given in Chapter 2. Computer Simulation 29 terms of L so that L\ = L{/L where i = x,y,z. At least one of L'X,L', L'z will be 1, and all three will be one for a cubic simulation cell. The 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 1 j erfc(o;3tj) h [i<j SiJ Y v m a x -k g — ir2v2/a2 + TT77 ——T~ I2>exp(27rii/-a,-) | 2 27T + m ^ M 2 . (2.3.17) A few comments are in order about the scaled coordinates. Both in computer simulation and theoretical derivations (see, for example, section 2.4), it is convenient to work with scaled coordinates s = r/V1/3. For the most part, we employ a cubic simulation cell, in which case, L = V 1 / 3 and the scaled coordinates s — rjL have the same meaning as before. A rectangular box has V 1 / 3 / L, where L is the shortest box length. For a rectangular box, we have chosen to use s = r/L, as opposed to s = r/V1/3. Of course, this makes no difference to the final results. The periodic summation may also be performed for the dipolar pair potential. The computationally efficient Ewald expression for dipoles (in a cubic box) is [37] L [i<j SiJ SU Umax * ^ e — K 2 n 2 / a 2 + J2 2 I - n) exp(27rm • s{) \2 n # o n i ^ E ^ + f ^ M 2 , (2.3.18a) Chapter 2. Computer Simulation 30 where B(s) = erfc(a.s) + 2as (2.3.18b) and C{s) = B(s) + 4(as) 3 — — 6 (2.3.18c) The remaining quantities all have the same meaning as in Equation (2.3.15) for ions. 2.4 Chemica l Potent ia l Recall that three conditions must be satisfied for equilibrium between two phases of a system. The three conditions are equal temperatures, equal pressures and equal chem-ical potentials (for each component in a mixture) in each phase. The 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. What 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. To prove this point, there is no better example than the chemical potential. The 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. The 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 The total energy is E(rN,pN) = K(PN) + U(rN), (2.4.2a) where K(PN) = 7^EP1 (2-4.2b) is the total kinetic energy and m is the particle mass. The integrations over the canonical momenta pN are readily performed to yield We define the thermal wavelength as and the classical configuration integral as Z(N, V,T) = J...J e-WUr" . (2.4.5) We may work with scaled coordinates s = V _ 1 / 3 r , such that dr = Vds, so that the limits of integration no longer depend on volume to obtain ^ ~ N\A3NZS yN = QidZa , (2.4.6) where Q u is the ideal gas partition function and Zs = J ...je~^sN^dsN . (2.4.7) We note that Zs still depends on the volume through the total potential energy. For an ideal gas, the potential energy U is zero, thus Zs = 1. Chapter 2. Computer Simulation 32 The appropriate thermodynamic potential for the canonical ensemble is the Helmholtz free energy A — E — TS. The connection to statistical mechanics is through A(N,V,T) = -kTlnQ . (2.4.8) From Equation (2.4.6), we may write the free energy as A = Ald + Aex , (2.4.9a) where — = N ln A3/V + ln NI . (2.4.9b) kT Using Stirling's approximation (TV! m N ln N — AC), we find - ^ - = l n / 9 A 3 - i . (2.4.10) NkT y y ! From 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 3 . (2.4.14) We now look at the excess contribution. Since , u v / N,T Chapter 2. Computer Simulation 33 converting to a partial derivative with respect to density gives, It follows that dp ( A^T")N,T {pkr) p' (2.4.17) As a function of density, the excess contribution to the free energy may be expressed as Aex{p) A(p) - Ald(p) NkT NkT + / ' - 1 I ^ • (2-4.18) Aex(po) j" (JP_ _ ) d£ NkT JP0 [p'kT ) p' In the limit po —* 0, Aex(p0) also goes to 0 and we assume the integral is well-behaved to 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. The chemical potential is = A{N + 1,V,T)-A(N,V,T) = _ j f c r i n % ± 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 The 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. The 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 right-hand-side is pid/kT. The 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 _ / 3 A L / + . Thus we get the expression lx = u i d - k T \ n ( e - ^ ) N , (2.4.24) due to Widom [65], for the chemical potential. This is the particle insertion method. We may define the residual chemical potential ures through [59] H = kTln(^j +ures , (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 pex) which must be equal for two phases in equilibrium. The particle diameter cr has been included to keep the arguments of the logarithms dimensionless. V* = V/cr3 is the 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+ is of course still the energy for inserting an (TV + l)th ghost particle. 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 fac-tor. The 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. The pressure may be broken into P = Pid + Pex, where the ideal pressure is PM = NkT/V and the excess pressure is ' - r (£L- (£L-The pair potential may include both a continuous part and a discontinuous, hard-core term. The excess pressure may be broken into Pex = Pex,hard + Pex,soft to account for these two contributions to the total pressure. We first consider Pex,soft, which is identified as ex,so ft The potential energy often has the form Pex,sofi = - [ % ) • (2-5-6) 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 = L3, one may write The quantity being summed no longer depends on the volume. At 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). The 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. The hard-core contribution to the pressure may be identified with the first term on the right-hand-side of Equation (2.5.5), however, this is of little use for explicitly evaluating Pex,hard-The hard-core contribution to the pressure is most easily calculated from the virial equa-tion (which in turn is derived from the thermodynamic equation of state) — = l - - £ - [ g ( r ) r ( ^ ) d V , (2.5.11) NkT 6kT J y y ' \dr) K ! where g(r) is the radial distribution function. The radial distribution function is defined by 9(r) > ( 2 - 5 - 1 2 ) where 6(r) is the Dirac delta function. The radial distribution function is the ratio of the angle-averaged local density about a particle to the bulk density. The 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 The derivative in the integral is awkward, so one typically converts it to a derivative of a step function, giving Pe*Mrd = Y^^I9{r) ePUHS {jre~PUHS)r*dr • (2-5-14) The derivative of the step function gives a delta function, thus PesMrd = YkTP29(°)°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 Ns which are within a shell of width ACT about each hard sphere, we get an estimate for g(cr). The estimate is g(a) « ° A 2 Ns . (2.5.16) The factor 2/(N — 1) accounts for the fact that we are using all N(N — l ) /2 separations of the N particle system. Ns is estimated using the contact function F where Ns = J2F(i,j). (2.5.17) The 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) The maximum value for F is then Fmax ~ 1 + — • (2-5.19) Chapter 2. Computer Simulation 39 In general, the contact function has the property that t < 1 if 1 and 2 overlap W W = 1 if 1 and 2 are in contact (2.5.20a) > 1 if 1 and 2 do not overlap The expression for the hard-sphere compressibility (PV/NkT) is PV NkT = 1 + 7^7 l i m -z-Y.F(hJ) ^ 3N A < T - O Ao- f£ y ' J ' (2.5.21) Some workers define 8 = 2Aa/a and sum over 1 < F < 1 + 8. The problem has been reduced from extrapolating g(r) to contact to choosing a reasonable value for Aa/a. A typical value is Aa/a 0.01. The contact function method for calculating the pressure is surprisingly insensitive to the value up to moderate reduced densities (up to p* ?s 0.5). At higher reduced density, one must reach a compromise between the number of counted separations and the uncertainty. The 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 parti-cles. In this thesis, we consider hard ellipsoids of revolution. The 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 P* = Pa3/e. With a hard-core interaction, one typically reports the reduced pressure P* — Pa3/kT. Note that the reduced pressure depends on a. Thus if we arbitrarily measured the particle diameter as 2a, the reduced pressure would change. The 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 P* z= —O* 2 ex,hard o > (2.5.22) 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 xA _|_ (-^L) xAxB 9AB{O-AB) + (—'] x2B gBB(o-B) V cr A / W A / where / 9 * = /0<r ,^ O ^ B = (cr^ + 0 \ B ) / 2 and £; = Ni/N (i = A, B) 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. Much 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. The virial (or pressure) equation is op oo ^ = l + £ i W T > " . (2.6.1) P n=l Carnahan and Starling observed a pattern in the virial coefficients for a hard-sphere system. They were then able to perform the summation, giving BP 1 + 7] + n2 - n3 P (i-v)3 ' where n — \p*• From Equation (2.4.19), one can show that Aex r/(4 - 3r/) (2.6.2) NkT~ (1-n)2 • (2-6-3) Chapter 2. Computer Simulation 41 To this, we add — 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 8TT2 for 3 angles . As was the case for the translational conjugate momenta, we can readily integrate over the rotational conjugate momenta. The classical partition function for a particle with two rotational degrees of freedom is «=m (^£)"^ /e-""'",^oW- (2-7-la> where Qrot = h2/87r2Ik . (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 r o<, labeled QA, 0 B , and 0 c . For the partition Chapter 2. Computer Simulation 42 function for particles with three angular degrees of freedom, one replaces (T/QROT)N in Equation (2.7.1a) with (2.7.2) \ e A e B & c J ' In either case, Widom's expression for the chemical potential may be written as p = Hid,r0t ~ kTl*(^ / e-^dsNdnN^j , (2.7.3) where the effects of the rotational momenta are absorbed into pid,rot- The 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) The 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. In general, the order parameter may not be so simple. 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>. (2-7-5) i i The quantity calculated in simulations is the ensemble average of P l 5 (Pi) = (M)lN(p) . (2.7.6) P i measures the degree to which the dipoles are pointing in the same direction. Essen-tially, 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 P2 order parameter is defined as ^ = ^ E ^ ( 3 c o s 2 ^ - l ) , (2.7.7) where #8 is the angle between the unit vector iii along the symmetry axis of the ith particle and the director d (defined below). In practice, however, the order parameter is calculated using the ordering matrix [20]. The ordering matrix is Q = ^ E ^ ( 3 A A - J ) . (2.7.8) It has three eigenvalues, in decreasing order, A + , A 0 , A _ . The 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 P2. In the isotropic phase, P2 is zero. Due to finite size effects, the estimate of P2 in the isotropic phase from computer simulation is not quite zero. Rather, the eigenvalues in the isotropic phase have the N dependence A + = 0 Bv) • A° = 0 (if) ' M D x- = 0 (TF) ' < 2 ' 7 ' 9 ) In the nematic phase, the eigenvalues have the N dependence K = P2 + 0(Ly A o = _ l p 2 + 0 ( - ^ ) , and A _ = - I f t + 0 ( - i = ) . (2.7.10) The best estimate of P2 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 0, as it is easiest to distinguish between an isotropic and weakly nematic 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]. The difference between the two Chapter 2. Computer Simulation 44 estimates is small, but large enough to be annoying for the purpose of quantitative comparisons. The notation P2N is used to indicate when we compare with published results using A + 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. In this case, we must measure the biaxial order parameter Ql2 [23]. In general, the orientation of a biaxial particle is described by the Euler angles 0,<f>,ip. The biaxial order parameter is then defined as [23] '(1 + cos 2 di) Ql\2 — cos 2<j>i cos 2ipi — cos 9{ sin 2(j){ sin 2ipi (2.7.11) 2 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, Y, Z (as measured in the lab-fixed frame). Consider the unit vectors iix, iiy, iiz pointing along the axes of the biaxial particle. We may then construct three ordering matrices Q X \ QYY, Q Z \ where QXX = ^ E ^ ( 3 « f « f - I) , (2-7-12) with similar expressions for QVV and QZZ. The eigenvectors corresponding to the largest eigenvalues of each of the three matrices then define X, Y, Z. Our 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. The unit vector ii = (ux, uy, uz) for the particle is then rotated into the 2 = 0 plane, giving a vector if = . 1 (UX,Uy,0) \lul + uy = « , < , 0 ) - (2.7.13) The vector uv is readily constructed by the requirement that ux • iiy = 0. The result is ii y = «,-<,o) Chapter 2. Computer Simulation 45 0) • (2.7.14) The matrices QXX, QYY and the corresponding eigenvectors X, Y ave readily constructed, from which one can calculate Q\2 using the expression Q\2 = | ( X • QXX • X + Y • QYY • Y - X • QYY • X - Y • QXX • Y) . (2.7.15) Chapter 3 M o d e l Electrorheological F l u i d 3.1 Int roduct ion Rheology relates the way a fluid deforms or flows under the various forces which may be applied to that fluid. Electrorheological (ER) fluids are those whose rheology varies dramatically with applied electric field. Although the field dependence of the fluid prop-erties 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. The 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. The dipole in each particle tends to align with the field, giving an orientationally ordered system. The dipole-dipole interactions between the particles further enhances the orien-tational order. The 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 Chapter 3. Model Electrorheological Fluid 47 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 Bingham plastic). The plastic is described by the equation T = T 0 + 77o7 , (3.1.1) for T greater than the yield stress T 0 , where 770 is the plastic viscosity. A n E R fluid is observed to behave as a Bingham plastic in an electric field such that T ( E , 7 ) = T 0 ( E ) + 77O7, (3.1.2) for shearing perpendicular to the field direction. It is typically found that 770 does not depend on field, TO(0) = 0 (Newtonian fluid in zero field) and T0(E) OC EN where n ^ 2 ( E R effect). The enhanced apparent viscosity is defined as 7 7 (E ,7 ) - 77(0,7) = ^ • (3.1.3) 7 Showing that this quantity behaves as E2 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]. They 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]. The liquid phase of an E R fluid has been studied within the mean spherical approximation, and the enhancement of the dipole moment due to many-body effects was investigated [71]. The 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 Widom have studied the force between particles for the magnetic analogue of an E R fluid [72]. They 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. When polarized by a uniform applied field, such particles will distort to a prolate ellipsoidal shape. The distortion is understood by the competition between surface tension and field energy. The 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. The net charge on each large sphere is zero but each sphere will have its own field-induced dipole moment. The 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. Both 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. The unit of length, a, is defined through a3 = 8a2b. Physically, this is equivalent to taking the unit of length as the diameter of a sphere of equal volume to the ellipsoid. The unit of energy, ekT, is the product of the dielectric constant e of the solvent, Boltzmann's constant k, and temperature T. The dielectric constant e of the solvent is taken to be unity. 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. Our 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 Nion 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 iV , - o n , the reduced density p* = Ncr3/V, the reduced ion charge q* = q/^4:Tre0ekTa, the reduced ion diameter o*on = <Tion/a, and the reduced electric field E* = ^4:Tre0ea3/kTE must be set. We wish to make a connection between the reduced parameters employed here and typical experimental values. Typical E R particle diameters range from lpm to 100/mi [11]. Since this covers quite a large range, for now we take a — dpm, where d is a variable parameter. At room temperature (T = 300K), the reduced fundamental charge is q* ~ 0 .236J - 1 / 2 . A typical electric field at which the E R effect is observed is l k V / m m . In reduced units, this field gives E* w 164d 3/ 2 . The reduced field is sensitive to the particular value of d and can range from E* « 164 for d = 1 to E* « 164,000 for J = 100. We typically report reduced particle energies, i.e., in units of ekT. The reduced per particle energy u* for an infinite chain of point dipoles is [10] u* « — 2Ap*2, where p* = pIV'4irtoekT'a3 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. The maximum 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 Umax ~ 2^ion(i* bY1/3 max (3.2.1) where the quantity in square brackets is /*. For a spherical particle with a*on <C 1, p*max as Nlonq*/2. 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. The dimensionless parameter p*2 has been used to characterize E R fluids [8,71]. For p*2 1, chaining behavior is expected. Unfortunately, with present day state-of-the-art computer workstations,, a realistic practical limit to values of / V , o n is about ten. To make chain formation energetically favorable, we should have q* > 2p*max/Nion. For 7V t o n = 8, this requires q* > 0.16. Consider splitting a pair of ions using an electric field. Assuming cr*on <C 1, the total gain in reduced energy is q*E* whereas the loss due to the lower coulombic interaction is 1*2/a*on- Thus the field required to break up the pair is roughly given by equating these two energies, which implies E* « q*/u*on. cr*on is a particularly awkward parameter to 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 in 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 are given in Appendix B . 3 . 3 C o m p u t a t i o n a l D e t a i l s Periodic boundary conditions were imposed on the simulation cell and the Ewald sum-mation method [63,64] was employed to calculate the total potential energy. Recall that the Ewald summation method was discussed in section 2.3. The total energy U was broken up into a self energy Us and an interaction energy Ui, as defined below, and the field energy UF- The total potential energy of the system of ions confined within the E R particles was calculated using the Ewald summation method. The self energy Us was then defined as k i<3 i,j€k Chapter 3. Model Electrorheological Fluid 52 where u,j = qiqj/rij is the usual coulombic pair potential. The summation index k labels an E R particle, and i and j label ions within the kth E R particle. We then use Ui — U — Us — UF to determine the interaction energy. The dipole of the kth E R particle is /•** = £ ?fr,-. (3.3.2) The field energy is then UF = - £ E-iik = -E • £ nk = -E • M . (3.3.3) k k Note that for a polarizable system, the net lowering in energy due to the field is —\E-M [73] (assuming the net moment M grows linearly with E). For a non-polarizable system, the result is — E • M. 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 • M and resulting in a net lowering in energy of approximately — \E • M. The 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 . To 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. The 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. On average, a M C sweep was composed of TV translations, 2JV rotations (TV rotations for spheres) and JViV,-o n ion rolls. For non-spherical particles, rotations were performed about both the symmetry axis and the induced dipole. Each time, the type of move to perform was randomly selected with the appropriate probability. The step sizes were adjusted to give an acceptance ratio of ~ 0.5 for each of the three types of moves. Run 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 De te rmina t ion of Order ing and Structure As the E R fluid becomes polarized, each E R particle acquires an induced dipole p. The vector sum of all TV dipoles gives the net moment M. In section 2.7, we defined the normalized polarization per particle as (Pi) = (M)/N(p), 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. The orientational order of the symmetry axes is described by the (P 2) order parameter [20]. Recall from section 2.7 that it is defined as (P 2) = (i(3cos 20-1)} , (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 2) is zero in an isotropic phase but increases towards unity for an orientationally ordered phase. The 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)P^a/2Y ' ( 3 A 2 a ) where zij = • r{j\ , r\3 = \ri3 - z^e{\ , (3.4.2b) 6(x) is the Dirac delta function, 9(x) is the Heaviside step function, and e; is the unit vector along some direction associated with the ith 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)\2) (where 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)j2) is proportional to the number of 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 M o d e l B u i l d i n g The present E R models are characterized by a large number of parameters. The 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). The results are weakly dependent on the ion diameter. However, the depen-dence is weak only when the ions within each E R particle are relatively mobile. When 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 </0 («*/> (Pi) uncertainty ± 0.010 ± 0.04 ± 0.03 ± 0.6 ± 0.02 0.01 0.639 0.88 -0.67 -31.0 0.97 0.05 0.604 0.85 -0.60 -29.2 0.97 0.1 0.558 0.67 -0.49 -26.9 0.96 A set of simulations was performed varying the number of ions per E R particle (Nion) while keeping the total charge (q*Nion) and the reduced density of the ions {Niona*3n) fixed. The choice of quantities to keep fixed is somewhat arbitrary. For example, instead of fixing the charge density, the ion size could be held constant. The 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. The results at E* = 10 and Chapter 3. Model Electrorheological Fluid 56 E* = 100 are shown in Tables 3.2 and Tables 3.3, respectively. At low field, the results depend quite strongly on the number of ions. At high field, the dependence is not as dramatic. 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. At 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 o n . A l l further simulations were performed with 8 ions per particle. Table 3.2: Results for 6/a = !,/>* = 0.5, an d E* = 10. q* Nlon a*on (p*) K) <«}> (Pi) uncertainty ± 0.010 ± 0.12 ± 0.11 ± 0.1 ± 0.02 1.0 4 0.126 1.075 0.5 8 0.1 0.829 0.25 16 0.0794 0.585 -4.22 -1.57 -0.71 -2.12 -0.85 -0.33 -9 .9 -7 .4 -5 .0 0.92 0.89 0.85 Table 3.3: Results for b/a = 1, p* = 0.5, and E* = 100 as a function of number of ions per E R particle, q* and a*on were adjusted so that charge density and packing fraction within E R particle remained fixed. a* N- a* 1 lyion uton ( « 5 > (Pi) uncertainty ± 0.001 ± 0.01 ± 0.3 ± 0.1 ± 0.001 1.0 4 0.126 1.696 6.30 -14.0 -168.8 0.996 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 57 3 . 6 Computa t iona l Tests We first checked how our results depend on the choice of boundary conditions (see Ta-ble 3.4). The 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 (vac-uum) and e' —> oo (a conductor). The difference is that in the vacuum case, there is a positive contribution to the energy proportional to M2. 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. Which value of d should be used? Ideally, a value of e' should be chosen consistent with the dielectric constant of the system. A typical E R fluid is composed of materials of low dielectric constant. 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. The 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. The results were then checked that both initial configurations lead to the same equilibrium Chapter 3. Model Electrorheological Fluid 58 Table 3.4: Results for b/a = 1, p* = 0.75, q* = 0.35, a*on = 0.05, and E* = 20. B C indicates the boundary condition used. For the Ewald method, listed are a, nmax and e'. B C (Pi) uncertainty ± 0.004 ± 0.03 ± 0.03 ± 0.1 ± 0.008 5,14,1 0.848 0.63 -0.23 -16.0 0.941 3,14,oo 0.900 0.98 -1.22 -17.1 0.950 5,14,oo 0.905 1.02 -1.42 -17.2 0.950 7,14,oo 0.906 1.03 -1.44 -17.2 0.951 5,25,oo 0.907 1.04 -1.44 -17.3 0.950 values. For certain sets of parameters, this was the case. At 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*on. For the smaller ion diameter, the ions form bound pairs within the E R particle at zero field. At 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. The dependence of the results on the shape of the simulation cell was investigated. From the length of the simulations, it was difficult to determine whether a chained system was in 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. The 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. The 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. The results are given in Table 3.5. The 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 = 0.05, and E* = 100 box ratio initial configuration </*•> < « 5 > (w/> (PI) uncertainty ± 0.001 ± 0.05 ± 0.17 ± 0.1 ± 0.001 1:1:1 random 1:1:1 lattice 6.53:3.27:4 lattice 1.242 1.242 1.248 5.37 5.38 5.71 -5.08 -123.3 0.993 -5.10 -123.3 0.993 -6.86 -123.9 0.993 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). The 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 60 results. The 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 Fig. 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*on = 0.1, and E* = 100. TV K> (Pi) uncertainty ± 0.001 ± 0.04 ± 0.1 ± 0.1 ± 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 / 3 and 2b* — (b/a)2/3. Some tests for N dependence were performed at the extreme limits of particle shape considered in this thesis. The 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. The 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*on = 0.1 at E* = 100. Field direction is indicated by arrow in middle of box. Chapter 3. Model Electrorheological Fluid 62 Table 3.7: Results for p* — 0.75, q* = 0.35, a*on = 0.05, and E* = 40 N b/a </**> («s> (u*i) <«*•> (Pi) (P*) (cos 7) uncertainty ± 0.003 ± 0.04 ± 0.11 ± 0.1 ± 0.003 ± 0.001 ± 0.0008 64 3 128 3 2.501 2.513 6.41 6.56 -3.68 -4.74 -99.1 -99.7 0.991 0.992 0.970 0.972 0.9985 0.9986 64 2 128 2 1.856 1.859 5.15 5.18 -4.02 -3.99 -73.4 -73.5 0.988 0.988 0.950 0.951 0.9948 0.9949 64 0.333 128 0.333 1.687 1.683 4.15 4.09 -3.43 -3.44 -66.6 -66.4 0.988 0.987 0.0348 0.0354 3.7 Results and Discussion The configurations in Fig . 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. Two points should be noted on this matter. The periodic boundary conditions imply that a chain stretching across the simulation cell represents an infinitely long chain. Also, while chains in 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. The behavior of spherical and ellipsoidal E R fluid models is now investigated more systematically. 3.7.1 E R Behavior as a Funct ion of Elec t r ic F i e ld 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. The system remains in a fluid phase for all but the highest field studied. At E* = 100, the system was in a chained phase. The motion of the chains was slow, and they appeared to wriggle about an equilibrium position. The mean square displacements, shown in Fig. 3.2, confirm this picture. The mean square displacement at E* = 100 appears to level off, indicating a solid phase. The 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*on = 0.05 as a function of E*. The number in brackets is the uncertainty (three standard deviations of the mean) in the final digit. E* </0 < « 5 > <«J> (Pi) uncertainty ± 0.002 ± 0.04 ± 0.13 ± 0.1 0 0.351 -2.00 -0.08 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) 100 1.242 5.37 -5.08 -123.3 0.9930(6) The average dipole moment (the magnitude of the dipole) and the order parameter (Pi) are plotted as a function of the reduced field in Fig. 3.3(a) and Fig. 3.3(b), respec-tively. The 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 64 MC sweeps Figure 3.2: The mean square displacement for b/a = 1, p* — 0.75, q* = 0.35, and a*on — 0.05. Solid, dotted, dashed, and long-dashed curves are for E* — 0,20,50, and 100. 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. At 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. The small number of ions within the E R particle will give it a dipole at any instant in time. At each instant, the dipole moment may be large, even though the vector sum of the dipole over time is zero. The 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 maximum 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 proportional to the field E), that chain formation occurs. From Fig. 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]. At high field, the total energy is dominated by the field energy with the self and interaction energies roughly cancelling each other out. The repulsion from an ion's like-charged 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. The 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 Fig . 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 Fig. 3.2. Chapter 3. Model Electrorheological Fluid 67 1 diameter and 2 diameters. g(r) is given in Fig. 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. The dipole moment distributions are plotted in Fig. 3.5. The dipole moment distribu-tions 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 chain-ing. The 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 Fig. 3.6. The 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. The 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. At 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 68 Figure 3.4: (a) g\\(z) and (b) g(r) and for same system as in Fig. 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 Fig. 3.2. Solid, dotted, dashed, and long-dashed curves are for E* — 0,20,50 and 100. Chapter 3. Model Electrorheological Fluid 7 0 (A) (B) Figure 3.6: Sample configurations for same system as in Fig. 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. The 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 Fig. 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 Fig. 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 Fig. 3.9). At E* — 0, the system is disordered. At a small reduced field of E* = 10, the system is strongly ordered, but not chained. At E* — 100, the system is chained. The 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 sim-ulations, 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 2) as a function of field for b/a = 3, p* = 0.75, q* = 0.35, and a*on = 0.05. Chapter 3. Model Electrorheological Fluid Figure 3.8: g\\(z\u) for same system as in Fig. 3.7. The solid, dotted, and dashed are for E* = 20, 50, and 100. Chapter 3. Model Electrorheological Fluid 74 Figure 3.9: Sample configurations for b/a = 3, p* - 0.75, q* = 0.35, and cr*on = 0.05 at (a) E* = 0, (b) E* = 10, and (c) E* = 100. Chapter 3. Model Electrorheological Fluid 75 Table 3.9: Results for b/a = 3, p* = 0.75, q* = 0.35, and cr*on = 0.05 as a function of reduced field. Results are for 128 particles in a 1:1:2 box. E* </0 < « 5 > (Pi) (P*) uncertainty ± 0.013 ± 0.04 ± 0.16 ± 0.1 0 0.420 -1.84 -0.11 0.10(2) 0.15(4) 5 0.838 -0.90 -0.51 -3.6 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||) . What we have measured is the trace of the diffusion tensor (D = | D | | + f - D j J -Initially, the curves appear quite linear. As chains are formed, the viscosity may grow faster than E2. 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 Behavior as a Funct ion of Par t ic le Shape The effect of E R particle shape on the electrorheological behavior was investigated. Re-sults as a function of b/a at E* = 40 are presented in Table 3.10. The field was strong enough to orient the E R particles, but not to cause chain formation. The mean square displacements indicated that each system was in a fluid phase under the chosen con-ditions. The dipole moment as a function of b/a is shown in Fig. 3.11(a). The shape Chapter 3. Model Electrorheological Fluid 76 2 h-1 / D 1 h 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ( a ) ; 1 1 1 1 1 1 1 1 1 1 1 I I I 1 I I I 0 20 40 60 80 100 E*2 / 100 3 1111111111111111111 o5 1 / D 0 20 40 60 80 100 E*2 / 100 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 fb\2/S Umax ~ 2Nionq* \a) ~ Nionq*b* (prelates) (3.7.1a) 1 / 6 \ ~ 1 / 3 fCax ~ 2 N i m q \ a ) ~Ni°n(L*a* (oblates). (3.7.1b) For prolates, the reduced dipole moment grows as (b/a)2/3 and for oblates it grows as (b/a)~ll3. This is due to the choice of length scale <7, where a is the diameter of a sphere of equal volume to the ellipsoids. When 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 («*/> <«F> (Pi) (cos 7) uncertainty ± 0.004 ± 0.05 ± 0.13 ± 0.2 ± 0.003 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 (P2) in Fig. 3.11(b), even slightly prolate particles were readily ordered by the field. The 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. The particles tend to maximize their dipole moment by orienting themselves with the longest Chapter 3. Model Electrorheological Fluid 78 3 I I I I I | I I I I | I I I I | I M 1 2.5 ' - (a) ' ' ' ' ' ' ' ' ' I ' I ' I I I T1 -3 0.2 F-_ l 1 1 1 1 1 1 1 11 1 | 1 l - l j l 1 • _l 1 1 1 1 1 1 1 —n — -1 1 1 1 1 1 1 1 "3 1 2 3 b / a Figure 3.11: (a) The 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 p* = 0.75, q* = 0.35, a*on = 0.05, and E* = 40. The dashed curve in (a) shows f^max' 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 Fig. 3.11(c). It is almost a step function. The 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 Fig. 3.12 for various particle shapes. Recall that the notation indicates that g\\(z) 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. The situation for oblate particles is shown by a sample configuration for 128 particles with elongation b/a = 0.333. Fig. 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. The dipole moment increased from the sphere to oblate ellipsoid to prolate ellipsoid. However, when scaled by the maximum dipole moment for each shape, the behavior of the dipole moment as a function of field was similar for each shape (see Fig. 3.14). The dipole moment shows a rapid increase over the same range of the electric field for each particle shape. The 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. The 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 Fig. 3.11. The 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 81 Figure 3.13: Sample configuration for b/a = 0.333, p* = 0.75, q* = 0.35, and a\on = 0.05 at E* = 40. Chapter 3. Model Electrorheological Fluid 82 the dipole moment. 3.7.3 Order ing of an E R F l u i d of Oblate El l ipsoids as a Funct ion of Dens i ty 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. The results for b/a = 0.333 as a function of density are presented in 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. The biaxial order parameter (Q\2) w a s calculated for the oblate ellipsoids. Recall that Q\2 was discussed in 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. The only feature of note in g\\(z; p) is that the peak becomes broadened as the density is increased. The 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. The 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. The 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. The mean square displacements indicate that Figure 3.14: Scaled dipole moment as a function of reduced electric field for a single par-ticle for q* = 0.35, and a*on = 0.05. Solid, dashed, and dot-dashed curves for b/a = 3,1, and 0.333. Chapter 3. Model Electrorheological Fluid 84 10 8 6 gn(z;M) A 4 2 0 0 0.5 1 1.5 2 z/cr gn(z;u) -" " l 1 i i i i i i i i i i i i i i i_ II II II l;i (b) j = III r i\ w \ I \ 11 11 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 0 0.5 1 1.5 2 z/a Figure 3.15: (a) g\\(z; p) and (b) g\\(z;u) for 6/a = 0.333, q* = 0.35, and a*on = 0.05 at E* = 100. Solid, dotted, and dashed curves are for p* = 0.75,0.9, and 1.05. 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*on = 0.05, and E* = 100 as a function of density. p* </**> («5> <«F> (Qh) (Pi) (cos 7) uncertainty ± 0.001 ± 0.04 ± 0.09 ± 0.1 ± 0.10 ± 0.0009 ± 0.0003 0.75 1.848 7.23 -6.96 -184.0 0.22 0.9958 0.0190 0.9 1.840 6.93 -6.55 -183.1 0.48 0.9952 0.0201 1.05 1.836 6.79 -6.53 -182.7 0.79 0.9953 0.0205 3.7.4 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. The interesting properties of an E R fluid are due to the chaining of the E R particles along the field direction. Our 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. The 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. The 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 L i q u i d Crys ta l Mode l s 4.1 In t roduct ion 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 where A is a constant and the angle dependence enters through P2(cos7), which is defined in section 4.2. This potential was chosen to allow comparisons between simulation results and predictions of Maier-Saupe theory [77]. Of course, one could not expect to simulate a true nematic liquid on a lattice. 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. The model is not the best one for a real liquid crystal. Nevertheless, it contains the key features required for an IN transition. Our interest in this model was to use the Gibbs ensemble Monte Carlo method to simulate the IN coexistence. The conventional [75,76] u(12) AP 2 (cos7) , 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. At 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 hys-teresis 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. To 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. The 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 gas-isotropic liquid transition is found for dipolar hard spheres [29,30]. The 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 Computa t iona l Detai ls We consider axially symmetric particles which interact via a pair potential of the form « ( 1 2 ) = u 0 ( r ) + ix0(12) , (4.2.1) where uo(r) is a spherically symmetric interaction, u a(12) is an anisotropic term and r is the distance between particles 1 and 2. We consider two forms for u0(r). For a soft-core model, we take u0(r) to be the usual Lennard-Jones potential ULJ{T). In addition, we also consider a hard-core model where u0(r) is the usual hard-sphere potential u#s(r). The anisotropic contribution to the potential is u „ ( 1 2 ) = -4Ae ( ^ ) 6 P 2 ( C O S 7 ) , (4.2.2) where P 2 ( c o s 7 ) is the usual second order Legendre polynomial [P2(cc) = (3x 2 — l)/2] and 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. The potential was spherically truncated for separations greater than half the box length. The usual corrections were applied to the Lennard-Jones contribution to the energy [37]. The correction for the Lennard-Jones contribution to the total potential en-ergy was given by Equation (2.3.4). We now discuss the correction due to the anisotropic term in the pair potential. The 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) is a function only of separation r for large r. It follows that g(l2) — g(r) for large r in the isotropic phase. Assuming g(r) ~ 1 beyond the cut-off separation rcut, the anisotropic contribution to the correction for the potential integrates to zero in the isotropic phase. The same argument cannot be made in the nematic phase. The projection g220(r) is one term in the general expansion of the pair distribution function [18,83]. In the nematic phase, the projection g220(r) must satisfy [18] 9220(r) - 5 (P 2 ) 2 , r - o o , (4.2.3) where (P2) is the equilibrium order parameter [20]. The correction for the anisotropic energy is then X(P2)2Ucor-6, 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 (P2) beforehand. In practice, this requires determining (P2) self-consistently. We have tested the effect of including the correction in an orientationally ordered phase. The difference made by including the correction should be most apparent in the anisotropic contribution to the energy for a small system. The 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 P2 order parameter which is zero in an isotropic phase Chapter 4. Liquid Crystal Models 91 and has a maximum 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 g220(r). g(r) is of little help in distinguishing between the isotropic and nematic phases. As previously noted, however, g220{r) decays to zero in the isotropic phase but has a finite value for large r in the nematic phase. Technically, large r should be less than the scale of spatial director fluctuations, although this is not a problem in 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 nexch particle exchange moves are also performed. nexch was adjusted so that between 1 and 3 percent of the particles were exchanged during each M C sweep. The 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]. Why 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. The problem is that even if one finds room to insert a particle in 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. The 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 at-tempted 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. The problem was that the insertion of randomly oriented particles into the nematic phase constantly pushed it into the isotropic phase. The 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 ro-tationally bias the inserted particles. The 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 co-existence curve where the density difference was sufficiently large that each box retained its identity long enough to collect useful statistical averages. At any particular moment, a phase with orientational degrees of freedom has a di-rector 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 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. At first, /(cos 9) is quite noisy, but it quickly settles down to a smooth curve. The orien-tation of the particle to be inserted is then sampled from the distribution function of the receiving box. The details on how this was performed are given below. Another technical detail is that a production run was typically started from an equilibrium configuration. The 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. They 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. This 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. They 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 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 (4.2.5) of 2. Chapter 4. Liquid Crystal Models 94 given by min( l , Pexch,bias) where _ / J(cosfl) *exch,bias — .;;/ n\±exch i I t . ^ . U I r J (cos6 i ) and Pexc/i is the expression for an unbiased exchange [57]. Recall from section 2.2 that Pexch = exp(Au;) where Aw is given by Equation (2.2.5). We must evaluate the two quantities fI(cos9) and fn(cosd). Evaluating fI(cos9) is straightforward. The 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 fI(cos9). Determining fn(cos9) is slightly more complicated. The orientation of the particle to be inserted into box II is sampled from the orientational distribution function fn(cos9) of box II. One method of performing this sampling is as follows. A random number R is uniformly selected from the unit interval. This value is equated with the area under the curve fn(cos9) from cos0 to 1, R= C f11 (cos 9')d cos 9'. (4.2.7) ./cos 8 The 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. In other words, since f11 (cos 9) is accumulated throughout the simulation, the only unknown in Equa-tion (4.2.7) is cos 9. For example, if the receiving box happened to be an isotropic phase, then f11 (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 pTes from the trial particle insertions during the Gibbs simulation. Since the particle insertions are rotationally biased, this bias must be removed to recover the correct estimate of the residual chemical potential. The equation with this correction is pres = -kT\n^{N^i)fexp-^u+^ , (4.2.8) Chapter 4. Liquid Crystal Models 95 where AU+ = U(N + 1) — U(N) 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 Resul ts and Discussion 4 . 3 . 1 Resul ts 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). The reduced pressure P* = P a 3 / e was set at 0.1. As 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. The system solidifies into an fee lattice. There are certain numbers of particles (4i3, i = 1, 2, 3 . . .) for which the fee lattice exactly fits the cubic simulation cell. Of 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. Upon 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 jump in density signaling the transition to a nematic phase [see Fig. 4.1(a)]. Nevertheless, the order parameter (P2N) 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 Fig. 4.1(b)]. We continued to cool this nematic phase and there Chapter 4. Liquid Crystal Models 96 Figure 4.1: 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. The dashed squares are the M C results of Luckhurst and Romano [78]. Chapter 4. Liquid Crystal Models 97 was a transition to a solid at T* = 0.6. The system spontaneously solidified into a (slightly imperfect) fee lattice at this reduced temperature. The 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* ~ 0.85, and a freezing branch with a freezing from a nematic liquid at T* ~ 0.6. The 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. Along the freezing branch, there is an IN transition. 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 g220(r) for typical isotropic, nematic and solid phases in Fig. 4.2. We also include a plot of the mean square displacement versus the number of M C sweeps (see Fig. 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. The 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 g220(r) for the same system as in figure 4.1. The solid, dotted, and dashed curves are for T* = 0.5,0.7, and 0.825. Chapter 4. Liquid Crystal Models 99 0.4 0.3 -0.2 0.1 0 1 1 1 1 1 1 1 1 1 1 I <[r(s)-r(0)]2 / / 1 | 1 1 1 | 1 > / a 2 1 1 — / / / / / / / 1 / / 1 f / - 1 ..-**" — ~ 1 :/ — / / V — I I 1 I 1 1 1 1 1 I 1 I 1 1 1 I I 1 0 200 400 600 800 1000 MC sweeps Figure 4.3: The mean square displacement as a function of the number of M C sweeps for the same system as in figure 4.1. The 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 100 this model at P* = 0.1 and A = 0.15. Their results are compared with ours in Fig . 4.1. Although their densities correspond to our IN 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. The embedded unit vectors were then aligned in a single direction and the anisotropic potential was turned on. The system was then heated in stages at constant pressure. The 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. The reduced density and (P2) 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 Fig. 4.4(c)]. Upon cooling the isotropic phase, the system shows some kind of transition at T* = 0.75 since it Chapter 4. Liquid Crystal Models 101 1.1 1 0.9 0.8 0.7 111111111111111111111111 ' i i i i I i i i i I i i i i I i i i i I m f 0.5 0.6 0.7 0.8 0.9 1 r 1.1 1 F-0.9 F -0.8 h 111111111111111111111111 (c) 3 0 7 hi 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 0.5 0.6 0.7 0.8 0.9 1 1.1 1 = P - 0.9 I-111111111111111111111111 (b) j 0.8 r 0 7 "l I I I I I I I I I I I I I I I I I I I I I M 0.5 0.6 0.7 0.8 0.9 1 r 0.7 111111111111111111111111 (d) j " i i i i I i i i i I i i i i I i i i i I i i i r 0.5 0.6 0.7 0.8 0.9 1 r 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 102 0.5 0.6 0.7 0.8 0.9 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. This 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 Fig. 4.4(a)]. The 32 particle system never formed a nematic phase. Any 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. This 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 Fig. 4.4(b)]. The liquid phase at T* = 0.75 has a small order parameter of (P2) ~ 0.28 [see Fig. 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 nematic-solid 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. The 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 de-pendence. In Fig. 4.6, both u* and u*a are plotted as a function of 1/N. The energies for the isotropic phase at T* = 0.9 show no noticeable N dependence. As expected, Chapter 4. Liquid Crystal Models 104 10 8 -<u> 6 4 -2 -0 T—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—r L B - B B - - ^ " B B IB--IS— BM-Bt—->•—H---P—+—!—>—+•••!—j-B I I 0 0.01 0.02 0.03 1/N Figure 4.6: 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. At 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. At T* = 0.7, there is a significant difference between the small and large N results. At larger N, the system is somewhat glassy and cannot freeze onto the lattice. At 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. The 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. Upon 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 jump in density. It is known that the IN 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. The discontinuity in density for the IN transition at A = 0.15 could not be resolved by the simulation results. At 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. For 32 particles, we found no nematic phase whatsoever. 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. The 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 TjN = 0.89 ± 0 . 0 1 and report that Maier-Saupe theory predicts TfN = 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. The 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]. Density functional theory predicts a transition to a nematic phase at TfN = 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%. The structure of the higher density phase, whether it is a nematic or solid, is an input for the density functional calculations. The 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. The 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 TjN = 0.89 is consistent with the integral equation theory result of T§L = 0.786 and is in reasonable agreement with the density functional theory result of TfN = 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. The 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. This 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. MS is Maier-Saupe theory [78], L R is Luckhurst and Romano [78], SL is stability limit [79], and D F T is density functional theory [80]. The Monte Carlo results (LR) were at reduced pressure P* = 0.1. The number in brackets is the uncertainty (three standard deviations of the mean) in the final digit. The theoretical predictions (SL and D F T ) were for fixed isotropic reduced density p*IN = 0.79. MS L R SL D F T 1IN 0.811 0.89(1) 0.786 0.875 P*N 0.790(2) 1.04 Pi 0.780(3) 0.79 0.79 We now discuss the Gibbs ensemble simulation results. At A = 0.15, we have deter-mined 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. The larger number of particles was used close to the critical point. The gas-liquid coex-istence curve was constructed and is shown in Fig. 4.7. The critical point values were determined for the gas-isotropic transition at A = 0.15. 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. At the low-est reduced temperature of T* = 0.7, the order parameter in the coexisting liquid was (P2) ~ 0.6. The 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 1.4 1.2 - o ho T* 1 0.8 0.6 i i i i i i I i i i I i i i • H i i i I i i l i i i 0 0.2 0.4 0.6 0.8 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 109 4.3.2 Results at A = 0.3 At 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. The system was started on an fee lattice and the temperature was slowly increased. The solid phase remained stable up to a reduced temperature of T* = 1.18. The 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. The isotropic phase was stable up to the gas-isotropic 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. The N P T results are shown in Fig. 4.8. To 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 Fig. 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 Fig. 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 1.2 1.4 r 1 0.8 0.6 <P»>n „ 0.4 0.2 0 0.8 1 1.2 1.4 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. At 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 I I I I I I I I I I I l - l 0 0.2 0.4 0.6 0.8 1 * P 1.1 1.2 1.3 r 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. The slightly larger system size of 600 particles helped to decrease the frequency of this occurrence. The 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. From Fig . 4.9(a) and Fig. 4.9(b), one can readily identify the triple point and two critical points for this system. 4 . 3 . 3 Results for the Hard-core 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. The coexistence results are plotted in Fig. 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. The 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 0.4 0.6 0.8 1.2 1 J I M | 1 1 1 1 | I I I I . y ( b ) ] P* 0.8 i x ~ 1 J / I -0.6 / — 0.4 " l 1 1 1 / l , , , , 1 , , , " 0.25 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. The dashed line in (b) gives the points at which a 108 particle system melts from a fee lattice. This 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. To test for this possibility, an fee lattice of 108 particles was heated under constant pressure until it melted. By repeating the process at various pressures, we generated the dashed curve in Fig. 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. prea is the residual chemical potential. (/) is the value of the orientation distribution function /(cos 0) extrapolated to cos $ — 1. The final column is the percentage of the total number of particles exchanged per M C sweep. p* Pi -es 1 kT I 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.286 0.53(2) 0.43(27) 0.22(5) 0.21(10) ( A ) % exchanged 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 21.1(10) 1.8(1) 0.030(9) 0.841(7) 1.9 The integral equation calculations were performed for this hard-core model at p* = 0.8 at A = 0.15 [79]. The temperature was lowered until the isotropic system became ori-entationally unstable. This gives the stability limit temperature for the IN transition. Chapter 4. Liquid Crystal Models 115 The reduced stability limit temperature using the R H N C closure was TgL = 0.85. Recall that integral equation theory with the R H N C closure predicted T§L = 0.786 for the soft-core model at p* = 0.79. This value for T£L 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. The integral equa-tion 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 tempera-tures well below the values used in the integral equation calculations. Our 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. The p vs. T diagrams for the soft-core model for A = 0.15 and A = 0.3 are shown in Figs. 4.7 and 4.9(a), respectively. The p vs. T diagram for the hard-core model is shown in Fig. 4.10(a). Both 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. Thus 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 116 case as well, is an IN coexistence. This coexistence appears as a small hump in the coexisting liquid curve [right-hand-side of Fig. 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. The phase diagram for the hard-core model is quite different from that for the soft-core model. Replacing the Lennard-Jones particles by hard spheres, with no anistropic potential, results in the disappearance of the gas-liquid coexistence. In other words, hard spheres only have a fluid and solid phase. 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. As 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]. The 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. The 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. The reason that a gas-isotropic liquid coexistence does not exist for the hard-core model is as follows. The gas phase has essentially zero configurational energy. The 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 Neu t ra l and Dipo la r H a r d Spheres 5.1 In t roduct ion The 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]. The fluid may undergo a gas-liquid transition as well as a demixing transition. The Lennard-Jones mixtures of Panagiotopoulos et al. provided examples of both types of transition. They 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. The 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. At 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 non-additive 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 . The demixing phase transition for a mixture of neutral (n) and dipolar (d) hard spheres has been investigated using integral equation theory [31,32]. The 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. The transition may occur at a temperature too low to get reliable simulation results. The condensation of hard-sphere ions is not much affected by the presence of neutral hard spheres [89]. The 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 demix-ing 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. The 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 re-duced critical dipole moment would increase as the diameter of the neutral hard spheres was decreased. The 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 gas-isotropic 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 and Computa t iona l Detai ls A mixture consists of Nd dipolar hard spheres of diameter ad and Nn neutral hard spheres of diameter an. The diameter of the dipolar hard spheres is taken as the unit of length ( a d = 1). The ratio of component diameters is a* = an/ad (reduced diameter of the neutral hard spheres). The total number of particles is N = Nd + Nn and the reduced density is p* = pa\, where p = N/V is the total number density. The mole fractions are Xd — Nd/N and xn = Nn/N = 1 — Xd- Each dipolar hard sphere has a reduced dipole moment p* = p/'y4nCotkl'ad. Equivalently, a reduced temperature may be defined 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 con-ditions were employed using the Ewald summation method [63,64]. The numerical param-eters used for the Ewald method were a convergence parameter a = 5 with a maximum lattice vector length of nmax = 14 and a surrounding continuum of dielectric constant e' = 1. 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. The Gibbs method involves two boxes, labeled I and II. One may specify the values of Nd, Nn, and the reduced total volume V* = V/ad for the two boxes. Instead of initially specifying the total volume (thus the total density), we may specify the reduced pressure P* = Pad/kT. Along with u* and cr*, this completely characterizes the system. Trial translations are performed for each particle within each box. The 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^ = p1/ and pln = pi1). 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. The 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. Numer-ically, this did not much affect the equilibrium values. The 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. The 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. The 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 123 5.3 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 Lennard-Jones mixtures. The Lennard-Jones mixture was composed of two components, labeled A and B. The 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 and EAB = \J^A^B- For a total mole fraction of XA,totai = 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. The 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. The standard mixing rules were used to determine the Lennard-Jones cross parameters (see text). The first row shows the published results [57]. The number in brackets is the uncertainty (three standard deviations of the mean) in the final digit(s). X A,total xA (P*. ) (N) I II I II I II 0.3 0.433 0.567 0.666(14) 0.69(1) 0.67(1) 0.70(3) 0.195(17) 0.207(5) 0.20(1) 0.22(2) 0.848(15) 0.836(9) 0.841(11) 0.836(6) 0.098(4) 0.099(1) 0.098(6) 0.100(3) 116(5) 484(5) 299(12) 301(12) 431(30) 169(30) As 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. The 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. The 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 dipole-rich 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. (x n) (AO / II I II I II 0.712(4) 0.685(7) 0.745(6) 0.683(3) 0.674(34) 0.394(28) 0.948(39) 0.968(14) 270(4) 230(4) 146(5) 354(5) The 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. The 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 xn (x (AO I 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 ( iV = 250), it was possible for the dipoles to form a chain spanning the entire simulation cell. The 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. Instead, one might see clusters, strands, or rings of dipoles. We used the largest system size that was computationally practical to avoid, as much as possible, any unphysical dipole structures. At high reduced dipole moment, the dipole-rich phase may become ferroelectric. Sim-ulations 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 The C r i t i c a l Dipo le M o m e n t as a Funct ion of Pressure 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 xn. At xn = 0.8, the mixture was expected to demix for p* > 0.66. Constant volume Gibbs ensemble simulations were performed for an equal diameter mixture of 500 particles with total mole fraction xn = 0.8 and total reduced density p* = 0.7. Starting at p* = \/2JS 1.6, the reduced dipole moment was increased until the demixing phase transition occurred. At p* = 1.8, the system was clearly mixed and at p* = 2.0, it was clearly demixed. At p* = 1.9, the system showed large fluctuations from a mixed phase. We take this to be the critical reduced dipole moment. The transition occurred at a reduced pressure of P* « 3.7. The critical reduced dipole moment occurs at a critical mole fraction. Chen and Forstmann report the mole fraction at their critical point 0.8. Simulation results confirm that the mixture does not demix symmetrically (i.e., critical point at mole fraction xn = 0.5). The simulation results are consistent with a mole fraction Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 127 xn ft 0.8 at the critical reduced dipole moment. Our results are qualitatively similar to the theoretical predictions. 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 xn. They report that the mixture becomes unstable for p* > 1.5 at xn = 0.8. We repeated our constant volume Gibbs ensemble simulations at total reduced density p* = 0.8. We found the critical reduced dipole moment p* = 1.8 at P* ft 5.8. Again the theoretical predictions underestimate the critical reduced dipole moment. We have attempted to find the reduced pressure at which the mixture becomes un-stable at p* = 2.0. Again using a total of 500 particles with total mole fraction xn = 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. The critical reduced dipole moments are tabulated as a function of reduced pressure in Table 5.4. Table 5.4: The critical reduced dipole moment as a function of reduced pressure for equal diameter mixture of neutral and dipolar hard spheres, p* is the reduced density of the mixture just before the demixing transition. p* Pi P* 5.9 1.8 0.8 3.7 1.9 0.7 2.4 2.0 0.6 Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 128 5 . 5 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. The results are tabulated in Table 5.5 and coexistence curves are plotted in Fig. 5.1. Although preStd 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. The first box is rich in dipolar hard spheres. pres,n and fireStd are the residual chemical potentials of the neutral and dipolar hard spheres, respectively. P </>*> pres n/kT pres d/kT (N) I II I II I II I II 2.1 2.2 2.3 2.4 0.662(6) 0.685(3) 0.697(9) 0.719(6) 0.546(7) 0.546(6) 0.546(8) 0.543(5) 3.83(11) 3.86(8) 3.86(15) 3.81(10) 3.84(11) 3.86(9) 3.86(14) 3.81(9) -1.05(30) -0.92(39) -2.15(70) -1.70(76) -1.25(45) -1.36(10) -1.94(44) — 1.90(13) 259(3) 245(3) 240(6) 229(4) 241(3) 255(3) 260(6) 271(4) The 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). Both 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 0.2 0.4 0.6 0.8 1 Figure 5.1: Coexistence curves for mixtures of neutral and dipolar hard spheres at P* = 2.0. The 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. The 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 The excess chemical potential pex is a monotonically increasing function of reduced den-sity. The 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. The 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 xn —> 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 hard-sphere residual chemical potential p,res,n = 3.84. These values compare very well with the simulation values in Table 5.5. 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* = 2.0, we have constructed the coexistence curves at a* = 0.9,0.8, and 0.7. Particle exchanges became easier to perform for the small, neutral hard spheres. 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 slower. Simulations at a* = 0.7 and p* = 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 Fig. 5.2 a sample config-uration from each of the two Gibbs boxes. The mixture with cr* = 0.8 was demixed pn = kT\n(xnp* ex n Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 131 at this reduced dipole moment. 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. The coexistence curves are shown in Fig. 5.1. From 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*. The critical dipole moment p* as a function of cr* is given in Table 5.6. These values were estimated from the coexistence curves in Fig. 5.1. The 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. While 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: The critical reduced dipole moment as a function of the diameter of the neutral hard spheres for a mixture of neutral and dipolar hard spheres. * M c 1.0 2.0 0.9 2.1 0.8 2.2 0.7 2.3 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]. They varied a parameter in the Hamiltonian, which smoothly transformed the system from a Stockmayer fluid to dipolar soft spheres. The extrapolated critical reduced dipole moment for dipolar soft 132 Chapter 5. Mixtures of Neutral and Dipolar Hard Spheres 133 spheres is close to our value. However, the gas-isotropic liquid was pre-empted by a chaining of the dipoles. They argue that a minimum 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. The 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 + Nion)N particles. We note that a model with 134 Chapter 6. Conclusions 135 soft-core repulsions would be more amenable to molecular dynamics simulation methods. Wi th 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, we have performed N P T simulations along the P* = 0.1 isobar. Luckhurst and Romano have also studied this system at the same reduced pressure using N P T Monte Carlo simulations [78]. They report an IN transition at T* ft 0.89. 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. The 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. The same was true when the anisotropy parameter was increased to A = 0.3. 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. The 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 136 the particle exchange rate is too low. The density difference between the two phases must also be sufficiently large for the biasing method to work. The 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. The 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. The 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. The 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. The 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. The 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 neu-tral 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. The 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. The 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]. Bibl iography [1] C . S. Powell, Sci. Amer. A p r i l , (1994). [2] S. W . Depp and W . E . Howard, Sci. Amer. M a r c h , (1993). [3] D . L . Klass and T . W . Martinek, J . Appl . Phys. 38, 67 (1967). [4] H . Conrad, Y . Chen, and A . F . Sprecher, Int. J . Mod. Phys. B 6, 2575 (1992). [5] D . R. Gamota and F . E . Filisko, Int. J . Mod. Phys. B 6, 2595 (1992). [6] W . M . Winslow, J . Appl . Phys. 2 0 , 1137 (1949). [7] H . Block and J . P. Kelley, J . Phys. D: Appl . 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 . Chem. Phys. 96, 2183 (1992). [12] M . W . Evans and D . M . Heyes, Phys. Chem. 95, 5287 (1991). [13] D. J . Klingenberg, F . Van Swol, and C . F . Zukoski, J . Chem. Phys. 94, 6160 (1991) [14] D . J . Klingenberg, F . Van Swol, and C . F . Zukoski, J . Chem. Phys. 94, 6170 (1991) [15] Y . Chen, A . F . Sprecher, and H . Conrad, J . Appl . Phys. 70, 6796 (1991). [16] T . C . Halsey, J . E . Martin, and D. Adolf, Phys. Rev. Lett. 68 1519 (1992). [17] J . E . Mart in and R. A . Anderson, J . Chem. 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 iq . 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 . Chem. 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 . Chem. Phys. 55, 4291 (1971). [34] M . Kasch and F . Forstmann, J . Chem. Phys. 99, 3037 (1993). [35] J.-P. Hansen and I. R. McDonald, 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, Moi . Phys. 98A, 169 (1991). [40] C . Lo and B . Palmer, J . Chem. Phys. 102, 925 (1995). [41] D . J . Evans and G . P. Morriss, J . Chem. 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 . Tome, Adv. Chem. Phys. 56, 141 (1984). [45] H . Metiu, Y . - T . L u , and Z. Zhang, Science 255, 1088 (1992). [46] J . S. Rowlinson and B . Widom, Molecular Theory of the Capillarity (Clarendon, Oxford, 1982). [47] S. Consta and R. Kapral , J . Chem. 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, Mol . 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 & Co, London, 1965).' [54] N . Metropolis, A . W . Rosenbluth, M . N . Rosenbluth, A . H . Teller, and E . Teller, J . Chem. 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, Ph . 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 . Rul l , G . Jackson, and B. Smit, Mol . Phys. 85, 435 (1995). Bibliography 141 F . H . Stillinger and T . A . Weber, Phys. Rev. B 31, 5262 (1985). P. P. Ewald, Ann. Phys. 64, 253 (1921). S. W . De Leeuw, J . W . Perram, and E . R. Smith, Proc. R. Soc. Lond. A 373, 27 (1980); 373, 57 (1980). B . Widom, J . Chem. Phys. 39, 2808 (1963). J . W . Perram, M . S. Werthiem, J . L . Lebowitz, and G . 0 . Williams , Chem. Phys. Lett. 105, 277 (1984). N . F . Carnahan and K . E . Starling, J . Chem. 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 . Mod. Phys. B 6, 2635 (1992). B . C . X u and K . C . Hass, J . Chem. Phys. 98, 2258 (1992). H . Zhang and M . Widom, 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 . Chem. Phys. 89, 6941 (1988). E . De Miguel, L . F . Rul 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 . Chem. Phys. 82, 429 (1985). Bibliography 142 [84] E . De Miguel, L . F . Rull , M . K . Chalam, K . E . Gubbins, and F . Van Swol, Mol . 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 . Chem. Phys. 9 4 , 2238 (1991). [88] J . G . Amar, Molec. Phys. 67, 739 (1989). [89] J . Shelley and G . N . Patey, unpublished. [90] P. G . Kusalik, J . Chem. Phys. 9 3 , 3520 (1990). [91] D . Levesque and J . J . Weis, Phys. Rev. E 4 9 , 5131 (1994). [92] G . Ayton and G . N . Patey, J . Chem. Phys. 102, 9040 (1995). [93] G . Orkoulas and A . Z. Panagiotopoulos, J . Chem. Phys. 101, 1452 (1994). [94] J . M . Caillol, J . Chem. Phys. 100, 2161 (1994). [95] H . Goldstein, Classical Mechanics, 2nd Edition (Addison-Wesley, Don 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. The orientation of a uniaxial particle is described by a unit vector it. 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]. The unit vector ii is given by ^ cos 6 sin 9 ^ u sin <j> sin 9 cos 9 (A . l ) where 9, (f> are the usual spherical coordinate system angles. The matrix A can be constructed by first rotating z through an angle 9 about the y axis. This is performed by the matrix ( cos 9 0 sin6> ^ A, = V (A.2) 0 1 0 — sin 9 0 cos 9 The resulting unit vector is then rotated through an angle (f> about the z axis. This is performed by the matrix AA = cos (j) — sin (j) 0 sin (j) cos (f) 0 0 0 1 \ (A.3) 143 Appendix A. Rotation of a Particle 144 Both Ag and A^ may be constructed by considering the effect each must have on the column unit vectors x, y, z. The matrix A is then the product of A$ and Ag, ^ cos (j) cos 9 — sin <f> cos <f> sin 9 ^ A = A6A6 sin <^> cos # cos <f> sin </> sin 9 — sin ^ 0 cos 6 (A.4) A unit vector itz is then generated in a spherical cap about the z axis on the unit sphere. This unit vector is given by the angles 9Z = 7rPdA0 cf)z = 2TTR2 , (A.5a) (A.5b) where Ri and R2 ave pseudo-random numbers uniformly selected from the unit interval and A9 € [0,1] determines the maximum angular step size. The unit vector is then ^ cos Sz sin 9Z ^ sin <j>z sin 9Z cos 9Z (A.6) The trial unit vector u is it = Ait, (A.7) A p p e n d i x B Choosing a N e w Ion Pos i t ion Inside an E l l i p so id 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. This is a restric-tion on the size of the ion diameter <7,-ON. Consider a slice through the x = 0 plane of an ellipsoid. This gives the equation of an ellipse 2 2 y z — H = 1 • a' b2 (B. l ) For an oblate ellipsoid (b < a), this is an ellipse with foci at ±\J a2 — b2. The maximum radius of the ion is contained between one focus and the nearest point of the ellipse. It follows that - (6 /a)" 1 / 3 ion,max b\2 1 - -w (B.2a) For a prolate ellipsoid (b > a) the result is ion,max = (Va)" 1 7 3 (B.2b) This restriction on ion size only becomes a concern for significantly non-spherical par-ticles. For example, for an oblate ellipsoid with b/a = 1/3, the maximum value is <r*im « 0.082. Even if the ion will fit inside the ellipsoid, too large a value of <7,0„ should be avoided for the following reason. To 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 146 of the ellipsoid. The probability of choosing a point on the ellipsoid surface is equal to the probability of choosing any other point on the surface. Thus the ion positions are determined by uniformly sampling the ellipsoid surface. Strictly speaking, it is the surface on which the ion centers of mass are confined that should be uniformly sampled. This is a surface within the ellipsoid whose normal is everywhere c t - o n / 2 from the ellipsoid surface. This 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. This biasing is not expected to be a problem in practice for the small values of a*on used in our study. A confined ion must be moved such that, in the absence of the Boltzmann 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- (R3) GT GT 0 ^ For an oblate ellipsoid (b < a) with eccentricity e = 1 — (b/a)2, the area of the ellipsoid is Aobiate = 4 7 r a 2 T / 2 sin 6(1 - e2 sin 2 9)1/2d8 Jo = 2 7 r a 2 + 7 r 6 2 - l n ^ . (B .4) e 1 — e For a prolate ellipsoid (b > a) with eccentricity e = 1 — (a/b)2, the area of the ellipsoid is /•7r/2 Aprdate = 4nab cos 6(1 — e2 sin 2 8)^ d8 Jo „ , arcsin e . = 2 7 r a 2 + 27ra6 . (B .5 e 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) through p(z)dz = (probability that a random point 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) The normalization condition / p(z)dz Jo requires 1 p(z)dz = — (area of ellipsoid between z and z + dz) Ai/2 (B.7) (B.8) where AT/2 is half the ellipsoid area. We associate a random number R uniformly chosen from the unit interval with the probability (probability that z is between 0 and z) = R = — f / 2 sin 0(1 - e2 sin 2 Bfl2de (B.9) A 1 / 2 Je Integrating gives R = 7ror A1/2 0\fl n i o o n I - e , / e cos cos Oyjl-e2 sin 2 6 + In 6 + Vl - e 2 s i n 2 0 x where ^ = (1 - • (B.lOa) (B.lOb) A similar treatment of the prolate case gives R irab Ai/2 sin 0\/l — e2smO + - arcsin(e sin 9) This gives R as a function f(z/b) of zjb. We must invert this function to get z/b as a function f_1(R) of R. This was done by tabulating R for a series of values of z/b and performing linear interpolation between points. To get a random point (x,y,z) on an ellipsoid, two random numbers £i and £ 2 ) where £ € [—1,1], are generated. The random Appendix B. Choosing a New Ion Position Inside an Ellipsoid 148 point is then given by z = frsignfo)/-1^!), (B.12a) 9 = 2TT|6| , (B.12b) / z2V/2 x = ( l - ^ - J a cos 9 , (B.12c) / z2V/2 y = 1 - 7 7 a s ' m d - ( B- 1 2 d) Often, it is more convenient to move an ion within a region centered about the ion in-stead 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 neigh-borhood about the ion on the ellipsoid surface. Given a point (x,y,z) on the ellipsoid, we calculate 8 = arctan (j-^j , 0 e [ O , 2 7 r ] (B.13a) R = sign • (B.13b) We then choose two random numbers £i and £2 and calculate R' = R + ^Ar.ion, R' € [ -1 ,1 ] (B.14a) 9' = 9 + n£2Arion , (B.14b) where A r , - o n € [0,1] is the maximum step size. R' and 9' are then used to generate (x',y',z'). The 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. A p p e n d i x C Rota t ion of Charge Coordinates The 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. This 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. The 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. The ion position r ; o n is measured with respect to the center of the E R particle. The 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. Thus we want the inverse matrix A - 1 . This is readily constructed from the inverse matrices cos (j) sin <f> 0 \ — sin cj) cos <f> 0 V o 149 Appendix C. Rotation of Charge Coordinates 150 and I cos 9 0 - sin 9 ^ An1 = V 0 1 0 • ( C 2 ) sin 9 0 cos 9 I and AftAa1 = I. Thus the inverse matrix is One can readily check that A^A^1 given by A-1 = A;1 A;1 . (C.3) One may note that A^1 is the transpose of A^ and similarly for Ag1. If we only use the matrix A - 1 , the E R particle ends up being rotated through an 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^. The 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) As 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 G [0,1] determines the maximum angle. This is set in combination with the maximum rotation to give an acceptance rate of roughly a half for trial rotations. This rotation is performed by the matrix ^ cosip — smip 0 ^ B ,i, = sin xp cos ip 0 0 0 1 This is just the matrix A^ except that the angle <p is replaced with ip. The operation is V (C.5) 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. The 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^Ag1 A^1 it . (C.8) This operation also rotates the ions from their original positions to their new trial posi-tions.
- 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