The Local Coulombic Monte Carlo Algorithm Applications to the Electric Double Layer by David Thompson B . S c . The University of Waterloo, 2004 A THESIS S U B M I T T E D IN P A R T I A L F U L F I 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 OF Master of Science in The Faculty of Graduate Studies (Physics) The University Of British Columbia August, 2007 © David Thompson 2007 Abstract A reformulation of the Coulomb problem, using a local Coulomb algorithm based on auxiliary fields, has been extended to slab and quasi-2D geome-tries. It has been implemented using Metropolis Monte Carlo and Gaussian charge interpolation functions. We have established the accuracy of the al-gorithm by generating effective pair potentials. Using this implementation, r the Gouy-Chapman problem was numerically resolved for constant potential slab boundaries. In the low coupling limit, we find excellent aggreement with analytic solutions. In the high coupling regime, we find agreement with the analytic theory in the limit of large wall separation. Using the contact value theorem, we calculate the pressure experienced by like-charged equipotential walls. The parameter space we consider pertains to many interesting bio-materials ranging from monovalent biomembranes to spermidine D N A . The numerical results show attractions mediated by counter-ions between the like-charged equipotential slab boundaries. We also extend the implementa-tion to allow for inhomogeneous dielectric backgrounds. The effect of a thin adsorbed layer of solvent is considered for an electrolyte system bounded by isolated electrodes. We show that a reduction in the dielectric value of this adsorbed layer results in a depletion of ions near the electrodes, even though the electrodes carry zero total charge. The applications considered show the versatility and accuracy of our implementation. Table of Contents Abstract ii Table of Contents i ii List of Tables vi List of Figures vii List of Programs ix Acknowledgements x Dedication . xi I Introduction 1 1 Importance of electrostatics and computational hurdles . 2 1.1 Periodic boundary conditions for bulk and slab systems . . . 4 1.2 An important model: The electric double layer 6 1.3 Limitations of existing methods . .' 7 1.4 History of the auxiliary-field technique 8 1.5 Outline of thesis . . . . . . 8 2 The Ewald summation method 9 2.1 The Ewald sum for 3D periodic boundary conditions . . . . 9 2.2 The Ewald sum for quasi-2D geometries . . 12 3 Electrostatics and constrained statistical mechanics .... 14 3.1 Homogeneous media 14 3.2 Inhomogeneous media . 17 3.3 Inhomogeneous media with metallic regions 18 iii Table of Contents 4 Analytic theory of the Gouy-Chapman problem 22 4.1 Electric double layer Hamiltonian 22 4.2 The weak coupling regime 24 4.2.1 One charged plate 30 4.2.2 Two charged plates 31 - 4.3 The strong coupling regime 32 4.3.1 Single charged plane 34 4.3.2 Two charged planes 34 4.4 The contact value theorem 35 4.5 Extensions 35 II A local Coulombic Monte Carlo algorithm 36 5 Numerical implementation of the local Coulomb algorithm 37 5.1 Metropolis Monte Carlo method 38 5.2 Discretization for bulk geometry 39 ' 5.2.1 Integrating transverse degrees of freedom 41 5.2.2 Integrating charge degrees of freedom 42 5.2.3 Surrounding medium 45 5.3 Discretization for slab geometries ' 47 5.3.1 Initialization 47 5.3.2 Surface charge fluctuations 48 5.3.3 Surface charge exchange 49 5.4 Quasi-2D geometry 49 5.5 Algorithm for a Monte Carlo sweep 50 5.6 Complexity 52 5.7 Sources of error 52 5.7.1 Green's function discretization 52 5.7.2 Aliasing 53 5.7.3 Cut off and charge neutrality violations 53 5.7.4 Short-range corrections 53 6 Analysis of accuracy 54 6.1 OK calibration of bulk geometries 55 6.2 OK calibration of quasi-2D geometry 59 7 An Application: The Gouy-Chapman problem 62 7.1 Weak coupling (Z -C 1) 64 7.2 Strong coupling (H » 1) '• • • 65 Table of Contents 7.3 Opt imiza t ion 67 7.4 Finite-size artifacts 68 7.5 Pressure between equipotential walls • • • • 70 8 The electric double layer: Inclusion of an inhomogeneous dielectric background 75 8.1 Implementation 75 8.2 Effect of mesh spacing 77 8.3 T h e Stern layer 80 9 Conclusion 83 Bibliography 84 List of Tables 5.1 Computational cost of each Monte Carlo update 52 6.1 Optimal parameters for bases of 5 3, 7 3 ; 9 3 ; and l l 3 58 • 7.1 Parameters used for E = 0.1 simulations 64 7.2 Parameters used for E = 100.0 simulations 66 7.3 Parameters used in testing finite-size artifacts 69 7.4 Parameters for pressure phase diagram simulations. 71 v i List of Figures 1.1 Periodic boundary condition for a rectangular simulation cell 4 1.2 Model systems with 2D symmetry 5 1.3 The method of images for two planar dielectric interfaces . . 6 1.4 Schematic of the electric double layer 6 2.1 Common bulk systems 10 2.2 Screening Gaussian functions for the Ewald sum 11 2.3 Common summation geometries for the Ewald method . . . . 11 2.4 Ewald method unit cell for the slab geometry .• . 13 3.1 Inhomogeneous ensemble with metallic regions 19 5.1 Schematic of the simulation mesh 40 5.2 Diagram of Link and Vertex classes 41 5.3 Plaquette definition 41 5.4 Plaquette update 42 5.5 Density fluctuation scheme 43 5.6 Possible charge density update 43 5.7 Different orientations of a Hamiltonian path 44 5.8 Initialization of the slab geometry 48 5.9 Implementation of surface charge fluctuations '. . 49 5.10 Implementation of surface charge exchange 50 6.1 Logarithmic inversion of g(r — r') 55 6.2 OK pair potential accuracy for bulk geometry 57 6.3 OK maximum pair potential error for bulk geometry 58 6.4 Green's function error for eq. (5.15) 59 6.5 OK pair potential accuracy for quasi^2D geometry 61 7.1 Implementation of conducting charged walls 63 7.2 Implementation of equipotential walls 63 7.3 Density for weak coupling simulations 65 List of Figures 7.4 Density for strong coupling simulations 66 7.5 Plaquette optimization 68 7.6 Finite-size artifacts 69 7.7 Pressure vs. d for equipotential planar walls 73 7.8 Phase diagram for pressure between equipotential walls, . . . 74 8.1 Implementation of the inhomogeneous electric double layer . 76 8.2 Pair potential for various mesh scales 78 8.3 Charge density for various mesh schemes 80 8.4 Charge density for Stern layer simulations 82 \ • •• • List of Programs 5.1 C++ code for generating an xy-oriented Hamiltonian path . . 46 5.2 C++ code for a MC sweep 51 ix Acknowledgements I would like to acknowledge my supervisor. Dr. Jorg Rottler. for constant guidance throughout this work. This work was supported in part by the National Science and Engineering Research Council of Canada (NSERC). It has been enabled by the use of WestGrid computing resources, which are funded in part by the Canada Foundation for Innovation, Alberta Innova-tion and Science, B C Advanced Education, and the participating research institutions. WestGrid equipment is provided by I B M , Hewlett Packard and SGI. X With love to Maria Part I Introduction i Chapter 1 Importance of electrostatics and computational hurdles Electrostatics plays a dominant role in many soft matter systems. For example, the stability of colloidal emulsions relies on Coulombic repulsion between colloidal particles[l, 2]. Many classes of materials such as polyelec-trolytes (charged polymer solutions) can only be understood through the inclusion of electrostatic interactions[3]. Often, the balance of these inter-actions is subtle, leading to startling behaviour. In biomaterials, this has led to many open questions concerning interactions between cellular mem-branes, transport of ions across cell walls, or the bundling of D N A within the cell's nucleus. These systems show astonishing phenomena ranging from at-tractions between like charged surfaces to the stable coiling of 10 1 0 electron charge within a diameter of 10/im[l, 2, 4. 5]. In biophysics, the systems studied are often complicated as well as macro-scopically large with scales of /urn to mm. Analytic theories are generally limited, making computer simulation the method of choice. Biophysics sim-ulations often include tens of thousands of charged particles using molec-ular dynamics (MD) or Monte Carlo (MC). Unlike short-range potentials, which model chemical bonds or Van Der Waals forces, Coulombic interac-tions cannot be truncated [6]. Thus a significant bulk of simulation time is spent calculating these interactions. This has made the aim of reducing the computational complexity of calculating the electrostatic energy of great importance. A naive sum over N charged particles would suggest the cost of calculating the Coulomb energy to be 0(N2). When periodic boundary conditions are applied this becomes even worse, as each charge must also interact with an infinite number of image charges. Some algorithms achieve scaling of O(N) but impose many restrictions on the systems they can reproduce. -These methods can typically be classified as one of the following techniques: 2 Chapter 1. Importance of electrostatics and computational hurdles summation techniques The Coulomb interaction is broken up into long-range terms and terms which may be truncated; an example is the Ewald sum. Often considered the standard in biophysics simulations, the Ewald sum makes use of the identity 1/r = erf(0r)/r + [1 — eri(0r)]/r. Terms composed of J~erf(^T") rapidly decay in real space while those made up of erf (/3r)/V quickly converge in Fourier-space[6]. The Ewald sum is often preferred due to its ease of implementation 3 and moderate efficiency, 0(N2). It is limited by its inability to handle large numbers of charges (> 103) or to deal with inhomogeneous di-electric functions e(r). It will be more thoroughly discussed in chapter 2. gr id /mesh techniques The electrostatic potential <p(r) is solved for by first discretizing the simulation space and then solving the general-ized Poisson equation V-e(r)Vc/>(r) = —p(r) on the grid, with p(r) the charge density and e(r) the dielectric function. Examples of these are particle-mesh Ewald and multigrid methods. For particle-mesh Ewald, the point charges are replaced by discrete charge distributions using interpolation functions. The charge density is Fourier-transformed and the electrostatic potential is solved for in Fourier-space. As with the Ewald sum, since the electrostatic interactions are resolved using Fourier decomposition the algorithm is unable to deal with inhomo-geneous dielectric functions e(r). The particle-mesh Ewald method may achieve scaling of 0(N log N) through the use of Fast-Fourier Transforms (FFTs)[6, 7, 8]. Multigrid approaches represents a gen-eral class of schemes used to solve elliptic partial differential equations on a grid. These methods are highly efficient and may achieve O(N) scaling, with a volume dependent cofactor.. They rely on the observa-tion that when solving elliptic partial differential equations on a grid, using relaxation methods, the residual error in the solution is most readily reduced when the grid spacing is similar to the -wavelength of the error. Thus by adaptively coarsening or refining the grid spacing, the error is efficiently reduced[9, 10]. By solving for the electrostatic potential in real space Multigrid offers greater flexibility than other mesh or summation techniques. constrained auxiliary-field techniques An effective Coulombic interac-tion is mediated between charged particles by the introduction of a constrained electric field. This algorithm is the subject of this thesis; we will discuss the theoretical foundation of this method in chapter 3. 3 Chapter 1.. Importance of electrostatics and computational hurdles 1.1 Periodic boundary conditions for bulk and slab systems Rather than simulating a macroscopic number of charged particles (~ 10 2 3), periodic boundary conditions (PBCs) are applied in biophysics simu-lations. By using PBCs , each charged particle is surrounded by an infinite number of chafges[2]. In this manner, the system's characteristics may be deduced from only hundreds to hundreds of thousands of charges. The ap-plication of periodic boundary conditions may be understood as an infinite replication of the simulation cell as shown in figure 1.1. o o o o o . o o o o o R " \ o ^ o o o o v. o O Figure 1.1: Periodic Boundary Condition for a Rectangular Simulation Cell. R, a direct lattice vector of the simulated lattice. When the characteristics of interest occur in the bulk of a material, P B C s are applied to all three dimensions. For behaviour near an interface, P B C s are only applied to two of the three dimensions. Examples of in-terfacial systems include electrolyte solutions between charged glass plates, lipid-bilayers, macro-ion membrane and membrane-membrane interactions in cellular biology, as well as interacting Wigner crystals in condensed mat-ter physics [1, 2, 4. 5, 11, 12]. Many simple quasi-2D models exist, a few of which are shown in figure 1.2. 4 Chapter 1. Importance of electrostatics and computational hurdles . • - * --~ Figure 1.2: Several simple model systems with 2D symmetry. A crude model of neutral lipid bilayers. point charges confined to a plane with a negative surface charge (top left). Membrane-membrane interactions in-volving counter ions, modeled by monovalent negative point charges con-fined between positively charged plates (top right). Macro-ion membrane interactions, modeled by a diffuse layer of monovalent point charges and a charged plane (bottom left). Polyelectrolytes between glass plates, modeled by charged strings constrained between charged planes (bottom right). For quasi-2D systems, the third dimension is either infinite or finite. When the third dimension is finite the simulation is said to have slab ge-ometry. Many algorithms have been tailored specifically for either bulk or quasi-2D geometries. Examples of 3D periodic methods include the Ewald sum and particle-mesh Ewald. For quasi-2D simulations, many summa-tion techniques have been developed including the Lekner-Sperb, Haut mart-Klein, and the two-dimension Ewald sum [13, 14, 15, 16, 17]. To simulate dielectric or metallic interfaces, boundary matrix or image charge methods must be be used in addition to the quasi-2D schemes. Boundary matrix methods discretize Maxwell's equations at an interface and solve for the in-duced electrostatic charge[18]. The alternative to boundary matrix methods is to ensure Maxwell's equations are satisfied at the interface(s) by introduc-ing image charges[19]. For a single interface this is often sufficient as only a finite number of image charges are required. However, for two or more the number of required image charges, even for a single charged particle. 5 Chapter 1. Importance of electrostatics and computational hurdles becomes infinite as is shown in figure 1.3. For s implic i ty many simulations util ize the method of images but truncate the images at a finite order, lead-ing to reduced accuracy and inferior scaling. Figure 1.3: M e t h o d of images for two planar symmetric dielectric interface 1.2 A n important model: The electric double layer A n important interfacial model is that of the electric double layer ( E D L ) . It consists of a charged wal l coupled to a diffuse region of mixed valency cation and anion species, as is shown in figure 1.4. Often included is a plane of nearest approach to the interface, which is used to model a th in layer of monovalent adsorbed anions. A common variant is interacting electric double layers, also shown in figure 1.4. + i e . mm o + l e mm i ° » o o ° ms 6 8 +2e +2e Figure 1.4: Schematic of the electric double layer (left) and interacting electric double layers (right). The E D L model has been the subject of much interest as it may be used on such varied systems as interfacial colloid emulsions and l ip id bilayers. A n important problem, known as the Gouy-Chapman problem[20], is the Chapter 1. Importance of electrostatics and computational hurdles determination of the ion density profiles. Analytic expressions are known in asymptotic limits, but a general theory remains unresolved. Computer simulation is a reliable means of predicting behaviour in the intermediate regimes as well as testing the validity of expansions and approximations. 1.3 Limitations of existing methods Realistic simulations remain limited by many factors. Most of the algo-rithms which provide suitable scaling cannot be applied directly to systems with slab geometry, or require additional computational complexity to re-solve interfacial boundaries. Methods reliant on Fourier transforms assume homogeneous dielectric functions. For many systems the effect of inhomoge-neous dielectrics is a reality and cannot be neglected as it leads to substantial modifications. For example, in the E D L , theory does not offer quantitative agreement with experiment, even in the known asymptotic limits, without the addition of a Stern layer[18]. The Stern layer is a thin layer of low di-electric coefficient that is believed to be caused by a strong polarization of the solvent near the dielectric interface. Contemporary simulations often incorporate hundreds of thousands of charged particles. To achieve simulations of this size, parallel computing-is necessary. Thus another scaling issue arises. Even though many meth-ods achieve impressive scaling with particle number, their efficiency peaks quickly as more processors are used. This is often due to the use of F F T s , which provide poor parallel performance due to the need to share the charge density across all processors at each time-step. This is indicative of a more severe limitation, non-locality. Most algorithms require a global solution of the potential after each update to the charge density. Thus global communi-cation between all processors is required following each update. An efficient parallel method would minimize this communication and ensure that only neighbouring nodes share information. Thus it is desirable to be able to compute the potential using only local updates. For Monte Carlo the issue of non-locality is even more problematic. For molecular dynamics global optimization remains reasonable, as the forces on all particles may be calculated simultaneously and the time stepped forward to achieve a new configuration in a single optimization. In contrast, in Monte Carlo a new arrangement requires N calculations of the electrostatic 7 Chapter 1. Importance of electrostatics and computational hurdles energy. This suggests that even the best electrostatic algorithms cost 0(N2) for Monte Carlo. 1 . 4 H i s t o r y o f t h e a u x i l i a r y - f i e l d t e c h n i q u e The auxiliary-field technique was proposed in 2002 by A . C . Maggs and collaborators[21]. It is a reformulation of the Coulomb.problem in terms of a constrained vector field, requiring local updates only. When implemented for Monte Carlo, it has a complexity O(N). It has been implemented for bulk and quasi-2D geometry lattice gases[21, 22, 23, 24]. As the formulation is entirely in real space, it is capable of treating inhomogeneous dielectrics. The behaviour of a lattice gas of mobile dielectrics has been studied in [24]. Extensions of the algorithm to simulate particles moving in a contin-uum have been implemented for bulk geometries[23, 25, 26]. In addition, a molecular dynamics method utilizing constrained auxiliary fields has also been developed[27] that also has a scaling O(N). 1 . 5 O u t l i n e o f t h e s i s For this thesis, we have developed an implementation of the auxiliary-field technique for application to the E D L . In chapters 2, 3, and 4 a review is given of the Ewald sum, the foundation of the auxiliary-field method, and the asymptotic limits of the Gouy-Chapman problem. A description of the implementation is given in chapter 5 and numerical accuracy is established in chapter 6. Our applications of the algorithm to the Gouy-Chapman problem and interacting E D L s with inhomogeneous dielectrics is presented in chapters 7 and 8. Throughout this thesis the electrostatic units used follow the SI convention. For conciseness, we present numerical results in Heaviside units with ks = eo = e = 1. 8 Chapter 2 The Ewald summation method In this chapter, we review the Ewald sum for systems with 2D and 3D periodicity as discussed in [6, 7. 8, 12, 13, 14, 16, 28, 29]. 2.1 T h e E w a l d s u m for 3 D p e r i o d i c b o u n d a r y c o n d i t i o n s Consider N charged point particles with charge {qj, q2,qN} at positions {rj, r - 2 , r / v } within a homogeneous dielectric medium of volume V and dielectric coefficient e. The electrostatic energy is 47T6 ^ r; - r; v ' i<j 11 Infinite replication of the simulation cell onto a three dimensional lattice defined by Bravais lattice vectors R, results in a sum which quickly converges with N. The energy of this replicated system is u = 1 L y y ' T [ — ( 2 . 2 ) with the prime indicating that for R = 0 the term with i = j is omitted. Unfortunately this sum is not absolutely convergent. In fact, for ^2,qi ^ 0 it does not converge at all. Fortunately for charge neutral .systems the electrostatic energy density does not diverge as long as a minimal length scale is imposed. However, the series is still not absolutely convergent, rather it converges conditionally. The value it converges to depends greatly upon the order in" which the terms are summed. This occurs because eq. (2.2) may be used to describe many different systems. In general, the application of 3D PBCs is intended to mimic a block of material, containing mobile 9 Chapter 2. The Ewald summation method charges embedded in a dielectric medium of dielectric constant es, such as shown in figure 2.1. Figure 2.1: Common boundary conditions used in conjunction with periodic boundary conditions for bulk systems. Vacuum (top, es = 1), tinfoil (middle, es = oo), and general dielectric (bottom). The electrostatic potential in the block is subject to boundary conditions at the interface. Therefore the energy contains a contribution from the interaction between the block's dipole moment and the induced charge at the boundary. Thus an absolutely convergent series, such as the Ewald sum, must include a term dependent on the simulation cell's dipole moment M , the summation geometry (hence the asymptotic shape of the block), and the permittivity of the surrounding medium, es/eo. The Ewald sum for 3D periodic systems is given by UEwald = TI direct + Urecip, + U BC (2-3) Uarect - ^ ^ W 4 v r e | | R | | (2-4) i.j R J * J l | k | | ^ 0 ' • V «, UBC = J(M,P,es) • (2.6) with 3 a real positive number, P the boundary geometry, and = — [12]. It may be understood as a factorization of the Coulomb problem into point charges screened by Gaussians, and interacting Gaussian distributions as shown in figure 2.2. 10 Chapter 2. The Ewald summation method Figure 2.2: The Ewald sum may be interpreted as the addition and sub-traction of Gaussian screening functions. The energy is resolved by the separation of terms as shown. The first term, 1)'direct, is the energy of the iV point charges, screened by Gaussians of width er == The second term, Urecip., contains a sum over all reciprocal lattice vectors k. The reciprocal lattice vectors are defined through the condition k-R = 27rn for arbitrary integer n, as in [30]. It may be seen as the energy of N interacting Gaussians of charge qi centered at the positions r^. Particle-mesh Ewald methods make use of this understanding by solving Poisson equation. eV2</>(r) = — p(r), on a grid for the N Gaus-sians, thus resolving Urecip.. The last term, J ( M , P , es) is the contribution to the energy due to the induced boundary charge. Common choices for es are shown in figure 2.1. For most simulations, spherical or plane-wise sum-mation (shown in figure 2.1) is used. For these cases J ( M , P , es) is given by ,7(M. P = spherical, es) J(M.P = slab,es) M 2 2Ve(es + l) 2V( (2.7) (2.8) " .... .... . . . .... Figure 2.3: Spherical (left) and plane-wise (right) summation geometries. 11 Chapter 2. The Ewald summation method 2.2 T h e E w a l d sum for quasi-2D geometries If we apply P B C s to only the x and y directions, the Coulomb energy for our N charged point particles is U 1 4TT£ N E E ' i<j Rj_ r, + Rj (2.9) with Rj^ a 2D direct lattice vector. For a rectangular lattice, has the form (nxLx.nyLy,{)) with Lx and Ly the linear dimensions of the unit cell. Upon use of similar identities to those used in the original Ewald derivation this may be expressed as El ±11 U E w a l , w = - ^ ^ m j 1 y <Mj cos(k • Tij) M l|k||#0 + e k Z i i erfc erfc I Bzij + 4e k 20 -0: + 20 1 N 7A E QiQj h3 Zij eric(0Zij) + 0 ^ -P2zl 0 N 1=1 Qi (2.10) with A the area of the 2D unit cell and Z y = z. Zj[12]. This equation is usually used only to establish the accuracy of other quasi-2D algorithms, as it is generally too expensive to calculate. For slab geometry systems, where the charges are confined to a finite region of width L , 3D periodic Ewald methods may be adapted to accurately approximate eq. (2.10). If a tetragonal-unit cell of volume NL3 and plane-wise summation is used, the energy of eq. (2.3) approaches that in eq. (2.10) for large N . The unit cell for this approximation is shown in figure 2.4. 12 Chapter 2. The Ewald summation method N x L L • i L Figure 2.4: Unit cell for the 3D-periodic Ewald method applied to the slab geometry Many simulations rely on this observation, as high quality 3D periodic Ewald methods have been developed. The leading order corrections have been established in [28, 29]. When the potential is resolved in this manner, it vanishes at infinite distance along the z direction. As was discussed in chapter 1, simulations wishing to incorporate dielectric or metallic interfaces must use boundary matrix or image charge methods in conjunction with this approximation. 13 Chapter 3 Electrostatics and constrained statistical mechanics In this chapter, we review a formulation of the partition function for a set of charges at positions r, with an electric field constrained to obey Gauss's law alone. We. show that this partition function is equivalent in many situations to the electrostatic partition function for the same system. This material has been developed by A . C . Maggs and collaborators in a series of papers [21, 23, 24, 25, 27, 31]. 3.1 Homogeneous media Consider a configuration of charges {gj,r-j} in a homogeneous dielectric medium. The electrostatic potential, d>(r), obeys Poisson's equation: V2r/>(r) = -p(r)/e and has energy U = / d J r ^(V</>)2 (3.1) (3.2) with p(r) = J2aifi(r ~ ri)- These electrostatics relations may be derived by minimizing an energy functional of a vector field E(r). Consider the energy functional . ^[E(r),0(r)] = y d 3 r We find the minima of T through e E 2 — - 0 { e V - E - p ( r ) (3.3) bT d 3 r eE • <5E - e0V-(5E d 3 r e V - E - p(r) =' / d 3re E + Vqb 5E d 3 r eV-E - p(r) (3.4) 14 Chapter 3. Electrostatics and constrained statistical mechanics The minima,(Ep, ( j ) p ) , of F\ E, (fi ] results in the equations of motion , Ep(r) =-V^,(r) (3.5) •V-Ep(r) = p(r)/e , • / (3.6) • (3.7) Upon substitution into T\ E, (fi ] we find the correct electrostatic energy T[ Ep, (fip] = J d ^ - ^ f ^ ( 3 . 8 ) This derivation serves as the motivation for a new formulation of electro-statics. Rather than minimize !?-"[ E, (fi } to solve Poisson's equation it is advantageous to use it as the actual energy density. Consider the canonical ensemble {g^r^} along with a field E(r). Instead of using the electrostatic energy density, we use T\ E. ]. The partition function is. Z N = T V l X p / d V V : D E e x p ( - / d 3 r ^ ) ^ V - e E W - ) <3'9) with A t the thermal wavelength of the charges. In addition to using T\ E, cfi ]. we have imposed the constraint that E(r) satisfy Gauss's law. With this contraint, E(r) may be decomposed into the minimum, E p , plus the curl of an arbitrary vector field Q(r). Clearly, this satisfies Gauss's law as V-E(r) = V-^Ep(r) + VxQ(r) = p(r)/e + V-VxQ(r) . = p(v)/e 15 Chapter 3. Electrostatics and constrained statistical mechanics (3.11) We may rewrite J d 3 r E 2 as j d 3 r E 2 = J d 3 r(E p + V x Q ) 2 = J d 3r(-Vr> p + V x Q ) 2 = y " d 3 r ^ ( V ( / ) p ) 2 - 2 V 0 p - V x Q + (VxQ) = j d 3 r (V0 p ) 2 + j d 3 r ( V x Q ) 2 + 2 jd3vebpV-VxQ - 2 j 4>pVxQ • dS = j dh{V<Pp)'2 + Jd3r{VxQ)2 We have used that the term § </>pVxQ • dS vanishes for periodic or constant potential boundary conditions. Substituting this into Zyv we arrive at the defining relation [21] • ZN = /TLW / D V V P E 5( EV-E(R) - ) ' I" f 3 6(E 2 + (VxQ) 2 ) ] x e x pri — — = ' f d 3 r " e * J - [ d 3 r ^ \ <3"12) x | P ( V x Q ) ( y ( V - V x Q ) e x p | - y d 3 r ^ ^ | = Zcoulombl { r i } ] x const. Thus aside from a constant, the partition function for the constrained field E(r). with energy functional T. is equal to a partition function constructed using Coulomb interactions. As the statistical weights of each configuration are independent of this constant, the two partition functions are equivalent, leading to the same expectation values for operators of {r^ }. The advantage of this formulation is that the energy density is locally defined in terms of D(r). The field D(r) in turn is only subject to the local constraint of Gauss's law. The evaluation of this partition function may be efficiently performed via Monte Carlo integration. We have developed an implementation, which is described in chapter 5. 16 Chapter 3. Electrostatics and constrained statistical mechanics 3.2 Inhomogeneous media The equivalence in eq. (3.12) is not limited to homogeneous dielectric materials[24]. For systems with a dielectric function e(r, {r-j}) we may ex-press F[ E , ( j f> ] and Gauss's law in terms of the displacement field D(r) = e(r. {r;})E(r). The constrained partition function for the inhomogeneous ensemble is Z n = iV!AF / d 3 r / V p D 6 { V ' D ( r ) _ p { r ) x exp<( — / d 3 r V2 (3.13) 2kBTe(r,{rl}) The general solution to V -D(r) = p(r), is D(r) = -e (r^(r ) )V^(r) + V x Q ( r ) (3.14) V-e(r,p(r))V0 p(r) = -p(r) (3.15) (3.16) With this decomposition, J d 3 r D 2 / 2 e can be rewritten as , 3 „ ° 2 _ / - ^ ( z i V ^ + V x Q ) 2 2e ' ^ ^ v ? ^ ^ , ( V x Q ) 5 2e = /d 3 r^M!+ / d 3 r ( V x Q ) 2 (3.17) 2 / 2e + / d V p V - V x Q j ^ V x Q • dS = | d : , r j ( V i ^ + | ( | 3 r ( V x Q ) ! Again we have utilized that the cross term dissappears for constant or pe-riodic boundary conditions. The partition function may be separated as before leading to N\XfN J x / D(VxQ)<) (V-VxQ) e J W B T ZCoulombl {ri} ] x ZFIUCI.I } (3.18) 17 Chapter 3. Electrostatics and constrained statistical mechanics For an e(r, {r;}) independent of the particle coordinates {r,}, ZpiuctX ] becomes constant and the equivalence is established. Throughout this thesis work we only consider e(r) independent of the charge-density. In general, for e(r, {rj}) dependent on the particle coordinates, it appears that this formu-lation leads to additional contributions to the free energy. By considering a pair of mobile dielectric inclusions in a vacuum, reference [24] was able to show that Z p i u c U leads to an interaction between the inclusions of 2>kBTv\8e\V20~e2 (47reo) 5 T/-/ \ o a y - i . u j u c i i ^ u c 2 , . V ( r ) = 7 7 — ^ (3-19) with v\ and V2 the volume of the inhomogeneities 5e\ and 5e-\. For a weakly dipolar material of volume vi and dipole moment pi , e^ /eo takes on the form 3fcjBT(47reo)2t'i substituting this form into the dielectric interaction we arrive at the Keesom potential which is the interaction potential for thermally agitated dipoles. Thus the additional terms in the free energy may be naturally interpreted as the interactions between fluctuating Langevin dipoles. Even though these in-teractions do not appear in the electrostatic solution, they are a desired addition to finite temperature systems. 3.3 Inhomogeneous media with metallic regions This formulation may even be extended to handle inhomogeneous media with metallic inclusions[31]. Consider an ensemble of charges in an inhomo-geneous dielectric background, with embedded conductors defined by sur-faces Si held at external potentials $ 2 . An example is shown in figure 3.1. 18 Chapter 3. Electrostatics and constrained statistical mechanics Figure 3.1: Inhomogeneous system with metallic regions. Before we may formulate the partition function in terms of a constrained field D(r) , we must determine what the appropriate free energy in electro-statics should be for Dirichlet boundary conditions, 4>(r) = $\ext\r). r G S,. So consider first the electrostatic energy, t/[E], of a system without metallic regions. A change of variables in terms of the potential requires that we determine the conjugate variable to <fr(r). Consider the variation in energy U[cj) + 8(p] = j d 3 r 3 e[V(4> + <5</>)]5 = U[(p}+ j d:5reV<P • VS4> + 0( Scj)2 ) = U\(p] + j d:irV-(-eV(50)(/> + 0{ Sep2 ) = U[(p} + j d : 5r Sp 0 + O( 5<j>2 ) (3.22) This shows that the conjugate variable to the potential is the charge density. Performing a Legendre transform to achieve the correct electrostatic free energy for our ensemble, including metallic interfaces, we see that f , 3 D 2 LEGENDRE f , 3 D 2 Jf jAexi) , c /o o • J ^ J d r _ 27 ' / d r 97 ~*r*'% 1 ' ^ where cr,(r) is the surface charge density of the ith equipotential surface. It is important to note that D is subject to Gauss's law throughout the 10 Chapter 3. Electrostatics and constrained statistical mechanics (3.26) medium and at the conducting boundaries. So that V-D(r) = p(r) r e V (3.24) D(r) • n?(r) = -a(r) r € Sj (3.25) for fij(r)' an outward facing unit normal to the surfaces, Sj. Following our experience in homogeneous media, we suspect that we can correctly derive the electrostatic equations of motion using the energy functional F[ D , o, <b, 3] = J d 3 r ( ^ - <f>{V D - p^j -z2<fs d S< (^*! 6 X t ) - S(fii • D + a*)) We confirm this through extremization, using 5T = j d 3 r ( ^jiP _ ^(y.<JD)^ +J^f ^ D • dS, -W d S i ^ - S ) ^ - f d 3 r(V-D - p)5<p^ y i dSi{hi • D + a^SE J i JSi = f d 3r(D/e + V0) • r5D + ^ / (H - )^<JD • dS* - w d S i ^ - s ) ^ - y d 3 r(V-D - p)8<p + j> dSi(nj • D + a ^ S Thus the extrema of T is defined by (3.27) -e ( r )V^r ) = D(r) r e V (3.28) V-D(r) = p(r) reV (3.29) M r ) - D ( r ) = ^(r) r e dSi (3.30) S(r) = </>(r) r € dS* (3.31) 2(r) = d>(-')(r) r e c\Sl (3.32) from which we arrive at the electrostatic equations of motion. The Legendre transform has achieved the desired effect of introducing Dirichlet boundary conditions. 20 Chapter 3. Electrostatics and constrained statistical mechanics Now consider the partition function with Dirichlet boundary conditions r(Dir.) _ 1 f fTjd3^ xVB 5( V--D-p)]j(vol 5( ni-D + Oi^ ( 3 . 3 3 ) > ^ - / < ^ ? £ ^ ) Let (fP be the electrostatic solution. We define erf (r) through the relation <7f(r) = -ni(r)-e(r)V<^(r), r e dS, ; ( 3 . 3 4 ) Substituting D = —eVcjP + V x Q and cr, = o\ + Oi into Zjy we get x ^ P ( V x Q ) <$( V-V.xQ) J J ^ P ^ S( hi • V x Q + a, )^ ( 3" 3 5) 2 /• ^ ^ ( e z t ) V ./ 2kBTe ^Js< kBT = ZCDoullmb\ 1 X ZF°iuct\ i r i ) \ As before, for a dielectric background independent of the particle coordinates (Dir ) we find that ZN is equivalent to the partition functions for an ensemble of particles interacting via Coulomb forces. However, the equivalent potential produced is now subject to Dirichlet boundary conditions along Sj. In the general case of mobile dielectrics, the term Z^^J has yet to be proven to result in interacting thermally agitated dipoles. However, the form of (Dir ) Zpiuct suggests equivalence. For the purpose of this work, proving the latter is unimportant as we are only concerned with systems with dielectric functions e(r) independent of the ion coordinates {ri}. This formulation is quite interesting as it presents an alternative to the use of either matrix boundary or image charges methods for dealing with metal-lic boundaries in the Coulomb problem. In addition, no approximations have been made to arrive at the statistical equivalence of these ensembles. 21 Chapter 4 Analytic theory of the Gouy-Chapman problem This chapter presents a review of field theory methods applied to the monovalent homogeneous EDL as developed by R. Netz and collaborators in [1. 2, 4. 5]. The asymptotic counter-ion density expressions reviewed are used to test our algorithm when applied to the slab geometry in chapter 7. 4 .1 Electric double layer Hamiltonian Consider N charges of charge —e confined to a region defined by the re-stricting function f2(r), interacting only via Coulomb interactions. In addi-tion to the charges we allow a surface charge, <r(r), to occupy the border of Q(r) . To ensure the system is charge neutral we have the condition -Ne + e J c\sr a(r) = 0 • (4.1) For example, for interacting EDLs, 0,(r) takes on the form = ° ^ i L ' (4.2, I 0, otherwise and a(r) = os5(z) + as5{Lz - z) (4.3) For simplicity we consider the system in a vacuum. The internal energy includes contributions from charge-charge, charge-wall, and wall-wall inter-22 Chapter 4. Analytic theory of the Gouy-Chapman problem actions. Thus the Hamiltonian is * - 2 N o N e 2mi • • + £ m r n ~ Id3™^ E - A — r IT (4-4) ^ 47re0 r 4 - Tj\\ J *r< 4TT€0 \\n - r + \jd3rdV;Vf(r,)„ - /d*rp(r)Mr) 2 7 4vreo 11^-^11 J We have introduced an external field h{r), coupled to the density operator p(r) = ~2~2f S(r — ri), to allow calculations of density distributions via { P { V ) ) = ^HrY ft=0 (4-5) Following reference [4] we do not couple h(r) to the charged surfaces. The reason for neglecting this term remains unknown. This system exhibits two classical length scales. The first, the Bjerrum length £ B . is defined as the length at which kBT is equal to the interaction between two unit charges. It is l B = S ^ T < 4 ' 6 ) The second, the Gouy-Chapman length p, is the distance at which a counter-ion has an electrostatic interaction with the boundary equal in strength to kBT. It is p = kBl — 1 6 G s (4.7) 2iroslB with us the average charge density along the boundary of $7(r). The coupling parameter, S, is the ratio of these lengths for the counter-ions. It is s = 7 ( 4- 8 ) 23 Chapter 4: Analytic theory of the Gouy-Chapman problem We now rewrite Ti as H kRT N - E 2m,iksT - e 2 y J d3rdVp(rMr-r>(r') + + e£B j d 3rd 3r'a(r)t;(r - r')p(r') - y y d3rd3rV(r)w(r - r')a(r') d 3r/5(r)/i(r) e2£BNv{0) (4.9) / d 3 r d 3 r ' 2 + y d 3r/i(r)p(r) + with v(r) = l /?\ When 7i is rescaled in terms of p and £ B it has a depen-dence on E only, as we will show in the following sections. The term with v(0) is included to remove self-interaction. 4.2 The weak coupling regime Weak coupling between the charged interface and the particles occurs when the Gouy-Chapman length is much greater than the Bjerrum length. In this case EE <C 1. We will apply this approximation following a reformu-lation of the partition function for our system. ep(r) — cr(r)j v(r — r') j^ ep(r') — er(r') e2£BNv(0) The partition function for our system in the canonical ensemble is Z N 1 m J J d 3 p ? ; d 3 r l f i ( r l ) ) exp We will show that this is equal to Z N d 3 r ? : ft(r;)\ [VA A? exp -W'[{ r i} ;A(r)] (4.10) (4.11) 24 Chapter 4. Analytic theory of the Gouy-Chapman problem with and -W'[{r,-},A(r)] = / d ' r d V A t r K V - O A l r ' ) + ~e j d 3rA(r)o-(r) - i J d3rA(r)/5(r) (4.12) f j 3 . w , Ne2£Bv(0) + J d 3r/i(r)p(r) + Z , = y P A e x p | - ^ 2 y d ^ d V A W u - ^ r - r O A C r O J (4.13) As before \ t is the thermal wavelength achieved through integration of the kinetic degrees of freedom. The transformation to a functional of the scalar field A(r) is known as the Hubbard-Stratonovich transform. After showing equivalence of the partition function in eq. (4.10) and eq. (4.11), we will transform to the grand-canonical ensemble. We will show that the free energy in this ensemble scales as E E - 1 . Thus a meanTfield approximation is appropriate for small H. We perform the saddle-point analysis and obtain a Poisson-Boltzmann equation for A(r), finding that A(r) is related to the electrostatic potential by a purely imaginary cofactor. We will utilize this relation by considering A(r) purely imaginary. Finally, we use the coupling to the external field to determine (p) for this approximation. Let ZP = YUJ P A e x p ( - W ' [ { r « } , A(r)]) (4.14) We rewrite Zp in terms of Fourier components using A(r) = i ^ A k e i k r A k = j d 3 r A ( r ) e - i k r (4.15) S(r) = i J2 e l k r *k,o = / d 3 r e " ' k r (4.16) V u and the relations A(r) = -A*( r ) A k = -A*_k (4.17) p(r) = p*(r) Pu = P-k (4.18) a(r) = a*(r) ak = a*_k (4.19) v(r) = v*(r) vk = v*_k (4.20) 25 Chapter 4. Analytic theory of the Gouy-Chapman problem From which we find d 3 r d V A ( r ) t ) - 1 ( r - r ' ) A ( r ' ) = - | d W A ' ^ - ^ r - r')A(r') = " E A k A k ' T 7 2 / d 3 r d V e - * V k ' r ' « - 1 ( r - r ' ) k,k' •' = - £ A £ A k , - ^ Id3Re"i(k-k')-R [dtAe-W+V-Wv-'iA) k,k' •' J • A = r - r ' k •/ = - ^ E i A k i V (4.21) and | d » r A ( p ) ( ^ - p ( r ) ) = ~ i 7 2 E A k ( v - ^ ) (4.22) So that we may rewrite W using H ' [ { r J ; A ( r ) ] - / d\h(r)p(r 1 Ne2£Bv(0) 91/2 £ 2 V 2 / r l f c { ^ l2 - 2 A B . A ; , k ( 2 i - p k k 1 91/2 £ 2 y 2 e 2 ^ k Ak - 'iez£BVk[ — - Pk e e2£BVk [ — - pu (4.23) 26 Chapter 4. Analytic theory of the Gouy-Chapman problem where we have used the identities in eq. (4.17) to (4.20). The last term is vke-iB\o-k „ 1 2 V2 Pk B 1 d 3rdV ep(r) — cr(r) v(r — r') ep(r') — cr(r') which follows in a similar fashion to eq. (4.21). So that ^=i/n(^)(n<»w)<v wf{r-" x exp \2V2 4 ^ e2eB\ A k - ie2eBvk{ — - p k (4.24) (4.25) Performing the change of variables 2 / A k - > A k + iq £Bvk p k e We find that Z P = /(II d 3 P i A ? e x P kBT Thus we arrive at eq. (4.11) (4.26) (4.27) ^ / ( n « ) / § M - * ' ^ H 27 Chapter 4. Analytic theory of the Gouy-Chapman problem Transforming to the grand-canonical ensemble results in the partition func-tion Zx = Y e»NZN = £ \$ZN N-Q N=0 = T J V A e X P { ~ 2 £ ^ j d 3 r d 3 r ' A ( r ) u _ 1 ( r " r')A(r') + %~ Jd3rA c^j 0 0 l f A d 3 / ] N\ x E 7 Y ? { A O J ~n(r)exp(h,(r) - zA(r) + -e2eBv(0)^ = ~klPAeXP{"27^2 / d ' r d V A W u - ^ r - r O A C r ' ) + ^ J d 3 r A a + A y d 3 r f i ( r )e / l ( r ) - i A ( r ) | ^ / »Aexp •£•1/ | - H [ A ] } (4.28) with Ao = eM the fugacity. and A = Aoee iBV(°)'2/\f.. We may simplify H using -y _ 1(r) = -|^ r<5(r) to get H = j d 3 rd 3 r 'A(r)f- 1 (r - r')A(r') - j d 3 r ^ A a + XQe'1'1^ , = " 2 / ^ 2 / d 3 r d V A ( r ) ^ ( J ( r - r')A(r') - j d 3 r ^ A a + A P ^ " ^ = " 8 ^ ? / d3'A(r)V»A(r) - /d3r(jA a + A Q e ^ ) d 3 r (VA) 2 - / d 3 r f - A a + XQe'1-' UfBe2 It is now advantageous to rescale H in terms of p, £B, and EL. Let (4.29) r = r/p ' (4.30) a(r) = pa(v)/as (4.31) A = • (4.32) fa/A = W\ (4-33) .28 Chapter 4. Analytic theory of the Gouy-Chapman problem Then d 3 r(VA(r)) 2 d 3 r : d 3r(VA(f)) d3fA(f)a(f) d 3 f A ^ ( f ) e h ( f ) - l A ( ^ a(r)A(r) 1 d 3rAfle' l ( r>- i A ( r ) 2nE 1 2 ^ E so that H[A(r)] . / i ( f ) - iA(f) (4.34) (4.35) (4.36) (VA(f)) 2 - zA(f)a(f) - XQ{r)e (4.37) =: ^W[A(f)] This suggests that in the limit E <C 1, 2^ is dominated by the saddle-point. The saddle-point is defined through-5H SA(T) = 0 H[A + 6A,h = 0] = « [A] + i / d ' r = .K[AJ + ± / d ' f - V A • V<5A - w i A + ?Afie-iA<5A + Q{5A2 (4.38) (4.39) 1 V 2 A -io + tXQe -iA 8A + 0(5A2) From which we get the saddle-point equation V 2 A s p + 2io = 2iAfie~ i A*p (4.40) This approximation has a very real'interpretation. The particle density for Z\ is (P(r)> = 5 ln 2 A 8h{r) h=0 i-Afi(r.) | p A e - A e x p ( - ^ l ) ( 4 . 4 1 ; \n(r)(e-iA) \n(r)e-iA"> 29 Chapter 4. Analytic theory of the Gouy-Chapman problem In terms of the rescaled ion the saddle-point equation now reads V 2 A s p + 2ib - 2ip . Let (j)sp = iAsp/2 =4>V 2^ s p = -p + a (4.42) From which we recognize a rescaled Poisson's equation for the electrostatic potential. Thus the use of the saddle-point approximation is equivalent to a mean-field theory where the ion-density decays with the potential. This equation may be derived very simply from a few assumption and is known as the Poisson-Boltzmann equation [32, 33, 34]. From electrostatics the potential </>(r) is related to the charge density through: In thermal equilibrium, the ion density adjusts according to the local value of the electrostatic potential. Using a Boltzmann weight for each ion, we find that the ion density is given by Thus the electrostatic potential is resolved through the Poisson-Boltzmann equation Eq. (4.45) results in eq. (4.42) for a homogeneous dielectric background when rescaled. The reason for following the field theoretic derivation to achieve this result, is to illuminate the fact that the Poisson-Boltzmann equation is in fact approximate; It is only applicable in the regime where S << 1. Deviations from the Poisson-Boltzmann equation are a result, of irregular ion densities, interactions between ions other than electrostatic interactions, and charge-charge correlations. 4.2.1 One charged plate For one charged plate, located at z = 0, we have V-e(r)V0(r) = ep(r) (4.43) (4.44) (4.45) (.4.46) 30 Chapter 4. Analytic theory of the Gouy-Chapman problem and cr(r) = 6(z). A n ansatz for iAsp is z A s p = 21n^l + A 1 / 2 / ) (4.47) which is trivially seen to be a solution of eq. (4.42). The rescaled density is p(~z)=(l + \1/2z) (4.48) We determine A through the normalization condition J dzp(z) = 1 - (4.49) from which we determine A = 1. Thus in the weak-coupling limit, the one-wall density is 4.2.2 Two charged plates For two charged plates, separated by a distance d in the z direction, we have ^ _ f l , 0 < z < d ^ )0, otherwise and a(f) = 6(z) + 6{d-z) (4.52) An ansatz for iAsp in eq. (4.42) is iAsp{z) = 21ncos(A 1/ 2[i - d/2]) (4.53) so that the density is given by ^ % o S ^ U ] ) (4-54) We find the value of A through the normalization condition Azp'z) = 2 (4.55) 31 Chapter 4. Analytic theory of the Gouy-Chapman problem which may be simplified to the transcendental equation A]/2tan(V/2J/2^ = 1 (4.56) Thus in the weak-coupling limit, the two-wall density is W - ^ U - « % > ] ) • (4-57) with A determined via eq. (4.56). 4.3 The strong coupling regime The ions are strongly bound to the charged perimeter when the Bjerrum length is-much larger than the Gouy-Chapman length. In this case, we apply perturbation theory in S _ 1 to determine the ion-density. We begin by rescaling the E D L Hamiltonian H - f- Pi c ^ W r r') I ^NvjO) kBT ^2m%kBT 2 ^ 1 j + 2 i i.j N •,. N + eiB J2 J d3ra(v)v(r - rf) + £ h(rz) i i J d3rd3rV(r>(r - r')a(r') 2 (4.58) 2niikBT 2 4 i,3 N + ^J2 j d3fff(f)t»(f - ri) + J > ( f i ) . i • i - j d:iM3r'd{r)v(r - f )a(r') 32 Chapter 4. Analytic theory of the Gouy-Chapman problem The grand-canonical partition function is then o o - ^ ) 7 n ( - ^ » ) x e x p j - ^ - ^ y d 3fd 3f / <r(f) 'u(f - f ' )cr(f ' ) | f ~ N 1 N f N ] = E ( ^ ) 7 n JV=0 N 7 J i v 7 x exp j - ^ y d 3fd 3f 'o'(f)u(f - f / ) a ( f ' ) | x e x p | - - ^ \ ( f i - r i ) + — J2 J d 3f5-(f)v(f - fi) + ^ / i ( r i ) j (4.59) The ion-density to leading order in A / E is 5\nZx\ (P(r)) = ft=0 <5fc(r) ^ Jd3h ^ ( r i ^ r - r O e x p ^ ^ d 3 f ' a ( r > ( r ' - f A ' / 1 3.Q(f) e x p ^ y d 3 f f f ( r>(r ' - f)^ (4.60) so that rescaled ion-density in the strong-coupling limit is p{r) ^XQ(r)exp(^ J d3r'o{r')v(?-r)^ (4.61) Thus the strong coupling partition function is the partition function for a single ion subject to c(r). 33 Chapter 4. Analytic theory of the Gouy-Chapman problem 4.3.1 Single charged plane For a single charged plane located at z = 0 we have fi(f) = ( 1 , ° ~ Z (4.62) ' [0, 0 < z ' and a(f) = 5{z) (4.63) The ion-density as given by eq. (4.61) takes the form ~p{~z) « n(z)e~2 (4.64) We have used that A = 1 from the normalization condition j dzp(z) = 1 (4.65) so that (p(z)) » 2 7 ^ ^ exp(-z/p) • (4.66) 4.3.2 Two charged planes We use the same Q(r) and <r(r) as in section 4.2.2. The integral 1 2TT So that ^ J dH'd{i')v{? - f) = 0 (4.67) ~p{~z) w 2fi(5)/d (4.68) Thus . (p{z)) « 4TrfBa2Q(z)p/d (4.69) Reference [4] shows that this approximation is valid only for small dfEL. For sufficiently large ci we expect the charged walls to decouple. Thus we expect to see two exponential profiles as described by the single wall theory. Unfortunately no analytic theory describes the cross-over in this behaviour. 34 Chapter 4. Analytic theory of the Gouy-Chapman problem 4.4 The contact value theorem The reason so much attention has been given to calculating the ion-density profile is due to the contact value theorem. This theorem relates the pres-sure, P. between two charged plates to the ion density at the interface. It is a general sum rule that be derived from statistical mechanical arguments. In terms of the rescaled ion-density it may be expressed as For weak coupling it has been shown that A > 1, and so the interaction between the plates is always repulsive. From eq. (4.68), we see that strong-coupling limit exhibits an attractive regime for d > 2. The charged plates are like charged, thus a simple analysis excluding counter-ions would suggest that the charged walls should repulse each other. Amazingly, the interaction between the walls mediated by the counter-ions can lead to attraction and even stable configurations. In principle this binding is similar to covalent bonds between like charged ions; however, the details differ significantly. 4.5 Extensions The analytic theories remains limited. The models assume homogeneous media as well as constant charge boundary conditions. In many cases these are not the most appropriate assumptions. For example, the hydrophilic surface of phospholipid bilayers may in many situation be treated as an equipotential [35]. As we have already discussed, the inclusion of dielectric inhomogeneities is also desired. In chapter 7, we confirm the results in equations (4.54) and (4.64). We then utilize the contact value theorem to predict the attractive regime for two equipotential planar walls. In chapter 8, we go beyond homogeneous dielectric media and consider the effects of an adsorbed solvent layer. (4.70) 35 Part II A local Coulombic Monte Carlo algorithm 36 Chapter 5 Numerical implementation of the local Coulomb algorithm c In this chapter, we introduce the numerical implementation developed for this work. Our simulation method is suitable for charges moving in a continuum of inhomogeneous dielectric media for bulk, slab or quasi-2D geometries. This is achieved through Monte Carlo integration of expectation values of operators A, using the constrained partition function z " = I ( n fi<r')d3r')i'D ' I ™ - * ) In chapter 3, we established the equivalence of this partition function with one constructed using Coulomb interactions. For constant charge boundary conditions we use the Hamiltonian W[{rl},D(r),{crl}] = | d 3 r ^ ! _ , (5.2) while for constant potential (Dirichlet) boundary conditions we use the Hamiltonian W[{rJ ,D(r) ,{aap"-) = / d h ^ ) ~ j • ^ When no charged surfaces are present, as in bulk simulations, the integration over J Dai is neglected. The integration is formulated in terms of four Monte Carlo updates • a particle move • a circular update of D(r) • a fluctuation in the local charge density of a surface 37 Chapter 5. Numerical implementation of the local Coulomb algorithm • an exchange of surface charge between metallic boundaries The particle, charge fluctuation, and exchange moves are all coupled to D(r) through the conditions n,; • D — —oi and V D = p. We begin by reviewing the Metropolis Monte Carlo method before getting into the details of the discretization scheme. 5.1 Metropolis Monte Carlo method Consider a system with coordinates described by a Hamiltonian Tt{qi). The thermodynamic expectation value of an operator A(qi) is given by n { A ) = hj (Ud9*)A(*)e~^ (5-4) i where Z is the appropriate partition function. Rosenbluth utilized that the expectation value calculated via a Markov chain asymptotically approaches (.4) [36]. A Markov chain is a string of configurations related by random updates, subject to the condition of detailed balance[37]: P(n)wn.m = P(m) (5.5) Here, P{n) is the probability of the system being in a configuration denoted by the index n and wn_m is the transition rate of the Markov chain from the state n to m. For the specified problem, the probability is given by with E(qlM) the energy of the configuration n. Metropolis et al. showed that detailed balance is achieved if the probability of accepting a random transition from the state {qi.n} to a new configuration {gj>m} is given by P<«*->(n-m) = 6 X PV " U (5.7) with AE the change of energy E(q^m) — E(qi_n). Thus the Metropolis Monte Carlo method may be summarizes as follows: 1. start with an initial random configuration {g^o} 2. perform an update to the configuration qjQ —» qj Q + 5 38 Chapter 5. Numerical implementation of the local Coulomb algorithm 3. if AE is less than 0, accept the update. Otherwise accept only if e-AE/kBi y r _ f o r r g (Q_ 2) a uniform random variable. 4. after attempting to update all degrees of freedom, </j, calculate the value of A. Use this instantaneous value to update a running average A. 5. repeat until A has sufficiently converged. The convergence rate of A may be optimized by ensuring that the acceptance rates are around 50%. In addition, the system exhibits a relaxation from its initial random distribution to a thermal configuration. It should be noted that the calculation of A should not proceed until after this relaxation. 5.2 Discretization for bulk geometry The simulation volume, V, is discretized, with charge density (pn) and dielectric function (en) defined on the vertices of a cubic mesh. In addition, a displacement field ( D n ) is defined on the links connecting the vertices. The grid voxels have a volume a 3 , so that vertices and links are located at n = a(i,j, k), vertices (5-8) n'x = a(i + 1/2, j, k), x-oriented field links (5-9) n'y = a(i,j + 1/2, k), y-oriented field links (5.10) n'z = a(i,j,k + 1/2), z-oriented field links (5-11) with i.j, k G Z. All displacement field variables ( D n ) take on the orientation of the link they sit on. For instance, a field defined on an x-oriented link would point in the positive x-directions. The particles' positions {r^} are not limited to vertices, rather they are allowed to continuously vary throughout V. subject to the restricting function O(r). The charge density on the vertices is determined via interpolation. This setup is illustrated in figure 5.1. 39 Chapter 5. Numerical implementation of the local Coulomb algorithm -t • - < i ) n ' ( 4 » 1 Figure 5.1: 2D sketch of the simulation mesh. Charge density and dielectric function sit on the vertices. The field variables are defined on the links. The shaded portions indicate the regions over which particles are interpolated. We have implemented V as an orthorhombic cell with linear dimensions aLx, aLy. and aLz. The imposition of PBCs on the x-direction is made by connecting the links at x=a(Lx — 1 + 1/2) to the vertices at x=0. Simi-lar conditions are applied to the y and z directions for full 3D periodicity. Gauss's law takes on the discrete form Pn (5-12) i { D r i * - o r • * + D " + ! 9 - o r 1 ' + o r 1 8 - Dr1*} = where we have used a seven-point divergence operator on D . A simple construction of the Hamiltonian in eq. (5.2) is n = Y. E ^ ( e V X D f ) 2 (5.13) n' e=x.y.z with e-*(e, n') = - + (5.14) 2 \ C n ' + | e e n ' - f e / For homogeneous media we have made use of an improved Hamiltonian proposed by A . C . Maggs[26]: n' e—x,y.z ^ ' This form leads to artifacts in the interaction potential which decay as l/h5 as discused in section 5.7.1. The scheme in eq. (5.13) has similar artifacts decaying as \/hA. Unfortunately, eq. (5.15) does not generalize in a simple way to inhomogeneous media. 40 Chapter 5. Numerical implementation of the local Coulomb algorithm The vertices and links are implemented in a similar fashion to linked-lists. Each vertex contains an e - 1 , a p, six pointers to connecting links, and six pointers to the neighbouring vertices. Each link contains a field D, an orientation, two pointers to the connecting vertices and two pointers to the nearest similarly oriented links. This is indicated in the diagram figure 5.2. Vertex - i p,e *Link[6] * Ver tex [6] Link D orientation *Link[2] * Ver tex [2] Figure 5.2: Diagram of the Link and Vertex classes. 5.2.1 Integrating transverse degrees of freedom As has been shown in chapter 3, the transverse modes of D(r) must be integrated to achieve statistical equivalence to the Coulomb problem. Con-sider a mesh, { D N , p n}, which satisfies Gauss's law. It is.advantageous to partition the field.links into plaquettes as shown in figure 5.3. 41 Chapter 5. Numerical implementation of the local Coulomb algorithm If we apply a circular update to the field links of a particular plaquette as shown in figure 5.4, we see that the field is modified without introducing any violations of Gauss's law. The flux into each vertex remains unchanged. 4 D3 3 4 D3-5 3 D4 D 2 — • D4-5 D2 + 5 1 D, 2 1 D,+5 2 Figure 5.4: Circular update to the displacement field which leaves the flux into each vertex unchanged. Thus we may integrate over the transverse modes of D n ' by introducing the following Monte Carlo update: 1. pick a random plaquette. 2. pick a random modification 5 € {—5mlax\Smax^) to apply to the pla-quette 3. accept the modification according to the Metropolis criterion The value of Smax^ is an optimization parameter which in general will depend on the temperature. The implementation of the plaquettes that we have developed is such that each one contains four pointers to link objects. For energy described by eq. (5.13) this leads to four modified terms while for eq. (5.15) twelve terms are involved. Thus this update is O( l ) . 5.2.2 Integrating charge degrees of freedom Consider first a modification to our mesh as shown in figure 5.5. The charge enclosed by a cube centered on vertex 1, of volume a 3 , is increased by as8. Thus §j D • dS] must correspondingly increase by a25. This is satisfied if D D + a6. -42 Chapter 5. Numerical implementation of the local Coulomb algorithm Pi D P 2 p.+8 D + a 5 p2-6 Figure 5.5: Simple density fluctuation satisfying Gauss's law Figure 5.6 illustrates a more complex charge update where density mod-ifications are made to n vertices. Let {5pi}, i 6 1,2, ...,n, be the change in the charge density of n connected vertices. If we increment the field con-necting nodes 1 and 2 by an amount abp\ then Gauss's law is satisfied at vertex 1. We may now proceed through the n vertices incrementing the field between i and i + 1 by ±a 2\I}=i °~Pji where the sign depends on the relative orientation of the link and a directed line between i and i + 1: As we pro-ceed along this path, Gauss's law is subsequently satisfied at each vertex. If we further require that the $Pi = 0, then by the time we reach the n t h vertex. Gauss's law is satisfied at all of the vertices. This may be seen by noting that had we extended the path by a single vertex, the increment to the link between n and n + 1 would be a 5pi = 0. ' 8p a (Sp, + 5p-,-t-...+ 8 p 1 5 ) 16 a8p l : Sp, 8p, 5p, a (Sp, + Sp,) Figure 5.6: A possible charge density update. 43 Chapter 5. Numerical implementation of the local Coulomb algorithm This algorithm for making modifications to the charge density allows us to formulate a way of moving charges in the continuum. We begin by inter-polating the charges to grid points n via Gaussian interpolation functions The sum over n is cut off at a value of ||r — [n]|| > Rc where [n] is the vertex closest to r, and Rc is a cut off radius measured in units of o. The linear dimension of the interpolation base is then p = 2Rc/o~ + 1 nodes. Thus each charge touches p3 vertices. For now lets assume that following our initial interpolation of all of the charges, we have made a global update such that the field satisfies Gauss's law. We may move a charge q from a position at r to r + Sr by interpolating —q at the position r and +q at the position r + 5r. This move introduces 0(p3) charge density updates {5pi}. Since the total charge interpolated is zero, we have £^ 5pi — 0. Thus we may use the above algorithm to re-impose Gauss's law. We only require a connected path between the 5pi. Figure 5.7 shows three such paths for the same region, these paths are known as Hamiltonian paths. A Hamiltonian path touches a set of nodes and has minimal length. The generalization to any orthorhombic region should be-obvious. Figure 5.7: Three different orientations of equivalent Hamiltonian paths. The C + _ l~ code for generating a Hamiltonian path which traverses pri-marily along xy-oriented planes is shown in program 5.1. The Hamiltonian path is stored in an array Move[ ]. It touches a total of nPts vertices and has dimension dimJ by dim_j in the xy-plane. For our algorithm, we choose the Hamiltonian path which touches primarily links oriented along the direction of the particle step. For example, for a particle move in the x-direction an xy-oriented Hamiltonian path is chosen. Thus we integrate the charge degrees of freedom with the following Monte Carlo move: (5.16) 44 Chapter 5. Numerical implementation of the local Coulomb algorithm 1. pick a random charge, i 6 [0,N — 1]. 2. pick a random orientation, e= x,y, or z, for the move and displacement " t V "max i °?nax J 3. interpolate a charge — ^ at and (p at r^ + 5e 4. perform a Hamiltonian path update for the charge density modifica-tion. To ensure we do not bias the field D n , we choose the orientation of the Hamiltonian to coincide with the orientation of the particle move. 5. accept the update according to the Metropolis criterion As with transverse updates the value S^ax^ may be tuned to optimize con-vergence. As the area affected includes 0(p3) vertices, this move has a complexity of 0(p3). In the coded implementation, we set an initial max-imum boundary on the size of 5maxm this w a y w e c a n pre-calculate all required Hamiltonian paths. All that remains is a means to satisfy Gauss's Law initially. We start with a zero initialized field. A Hamiltonian path update spanning the entire volume V results in a field configuration satisfying Gauss's law. There is one provision involved, the system's total charge must be zero. 5.2.3 Surrounding medium The algorithm give rise to some peculiar boundary conditions[25]. The dvnamics of the displacement field is such that it obeys for Q a rotational vector field (related very closely to the plaquettes) and J a current, such as described in the Hamiltonian path algorithm. When integrated, this may be used to define a dipole moment d. 1 /* d D(i) = - - y dtd 3 rJ = - - (5.18) d is related to the simulation cell's dipole moment M by d = M + gR, for R a Bravais lattice vector. There arises from d a contribution to the electrostatic energy of r f>2 d 2 45 Chapter 5. Numerical implementation of the local Coulomb algorithm Program 5.1 C++ code for generating an xy-oriented Hamiltonian path i = j = k. = iO = JO = 0; for(idx=0; idx < nPts; idx++){ i f ( i<dim_i-l && i>0 ){ i f ( i 0 < i){ iO = i ; i++; Move[idx] = INC_X; }/* i f 1.1 * / i f ( i 0 > i){ iO = i ; i — ; Move [idx] = DEC_X; }/* i f 1.2 */ }/* i f 1 * / else{ if(i0==i && i == dim_i - 1) { i0=i; i — ; Move[idx] = DEC_X; continue; >/* i f 2.1 * / if(i0==i && i == 0){ i0=i; i++; Move [idx] = INC_X; continue; }/* i f 2.2 */ iO = i ; i f ( j<dim_j-l && j>0 M i f ( j 0 < j){ j o = j ; j++; Move [idx] = INC_Y; }/* i f 2.3.1 * / if(JO > j){ JO = j:; j — ; Move [idx] = DEC_Y; }/* i f 2.3 .2 *./ }/* i f 2.3 */ else{ if(j0==j && j == dim_j-l){ j0=j; j — ; , Move [idx] = DEC_Y; continue; }/* i f 2.4.1 * / if(j0==j && j == ' O H . j0=j; Move[idx] = INC_Y; continue; }/* i f 2.4 .2 * / if(j0!=j){ •j0=j; k++; Move [idx] = INC_Z; }/* i f 2.4.3 */ }/* else 2.4 */ 37* else 2 * / }/* for idx */ 46 Chapter 5. Numerical implementation of the local Coulomb algorithm In reference [25] this extra term has been shown to result in a contribution to the free energy of f M i 1 / 1 / 3 ^ p Af=\w V (5.20) 10, otherwise If we refer back to eq. (2.7), we see that for V 1 / 3 -c ts the algorithm is equivalent to a Ewald sum with spherical summation geometry surrounded by a dielectric medium of e.s = 0. For V 1/-* » is, the boundary conditions correspond to a conducting surrounding medium, es = oo. 5.3 Discretization for slab geometries The discretization of slab geometry systems follows very closely to that of bulk geometry. We take the x-direction as finite and place charged surfaces on the vertices at the x = 0 and x = aLx planes. We use the restricting function n ( r ) = /1' • ap/2 < x < aLx - ap/2 ' 10, otherwise to ensure that the charges do not interpolate past the surfaces. The energy for constant-potential boundaries takes on the additional term a ^ M ^ o - S g ^ o O (5-22) By adding an additional layer of links to the x = aLx plane we may break the periodicity, in the x-direction. We require that these links do not contribute to the energy and do not- allow plaquette updates that would affect this condition. This is achieved by refusing any plaquettes which are xy or zx oriented at the x = aLx plane. Plaquette and particle updates then follow as before. 5.3.1 Initialization The slab geometry also requires that hi • D = — cr,, so that Hamiltonian paths no longer provide a means for initializing the grid. An alternative is to first migrate all charge violations away from the charged surfaces so that n, • D = —ai. We then follow with a Hamiltonian path which moves the accumulated violations along each vertex. This is shown in figure 5.8. As before, so long as the total charge on all vertices is zero, the mesh will satisfy Gauss's law following this update. 47 Chapter 5. Numerical implementation of the local Coulomb algorithm Figure 5.8: Initialization of the slab geometry. 5.3.2 Surface charge fluctuations To perform the integration over cr(r) we must introduce a Monte Carlo up-date which allows for a perturbation to the local surface charge density. This move must adhere to the constraint, • D = — <TJ, which when discretized becomes ID"+f* = - p n l x = 0 (5.23) l ^ x ' ^ = Pn\x=aLx (5.24) Consider the fluctuations shown in figure 5.9. These moves introduce the desired perturbation while maintaining Gauss's law and the boundary con-straints. Thus we may perform the integration over a(r) with J odSi held fixed with the following update: 1. pick a random surface i , i 6 0, aLx. 2. pick a random vertex j belonging to the surface, j G [0, LyLz] 3. pick a random direction, e= y, or z, for the charge migration and . displacement 5. e (-S^^^sl^1^) 4. perform the update according to figure 5.9 accepting according to the Metropolis criterion 48 Chapter 5. Numerical implementation of the local Coulomb algorithm - 5 a5 - 5 a5 a8 I 5 8 a5 -a5 Figure 5.9: Left, surface charge fluctuation on x = 0 plane. Right, fluctua-tion for x = aLx plane. The addition of local surface charge fluctuations does not allow for the equilibration of separate equipotentials. To achieve this, we introduce one more Monte Carlo update as suggested by Levrel [31]. This move allows charge density to flow from one equipotential to another. For the slab ge-ometry we choose the flow to be implemented in terms of x-directed lines from the surface at x = 0 to the surface at x = aLx. This is indicated in figure 5.10. The complexity of this move is 0(LX). 5.4 Quasi-2D geometry We achieve the general methods of quasi-2D geometry in the same manner as Particle-Mesh Ewald. We utilize an elongated tetragonal simulation cell and impose PBCs. We use the restricting function with n typically 3,4, or 5, and aLx = NL in figure 2.4. We also utilize an equivalence between our method and the Par tide-Mesh Ewald. Since both solve for the energy of interacting Gaussian clouds, the two methods should give identical results for low temperature up to a dipole term. Thus we may use the same corrections used for the Particle-Mesh Ewald when applied to quasi-2D geometries. Therefore, we correct for the summation order using the dipole term in eq. (2.8). 5.3.3 Surface charge exchange 0 < x < aLx/n otherwise (5.25) 49 Chapter 5. Numerical implementation of the local Coulomb algorithm Figure 5.10: Update exchanging of surface charge density. It is worth noting that this method exhibits one limitation. Since the sum of the charges on all vertices must be zero in order to satisfy Gauss's law, any terms involving charge in the Hamiltonian must be implemented explicitly on the mesh. In particular one cannot simulate a charged wall and counter-ions without including the charged wall on the mesh. 5.5 Algor i thm for a Monte Carlo sweep Program 5.2 shows C + + code for performing a Monte Carlo sweep. A n attempt to move each charge once is made in a random order. The proba-bility members m_dfProb of the Monte Carlo moves denoted m_mcinfo_Field, m_mcinfo_Particle, m_mcinfo_Surface_Fluct, m_mcinfo_Background_Field, and m_mcinfo_Surface_Exchange determine the relative ratio of plaquette, parti-cle, surface fluctuation, background field, and surface exchange updates per sweep. In this manner we can tune these proportions to achieve greater effi-ciency. The calls to member functions vTryJParticle, vTry JBackground_Field, vTry_Surface_Fluct, vTry_Field, vTry_SurfaceJExchange perform the Monte Carlo updates already described. 50 Chapter 5. Numerical implementation of the local Coulomb algorithm Program 5.2 C++ code for a MC sweep /* used to determine when to choose each type of MC move */ st a t i c REALTYPE df F i e l d = m_mcinfo_Field.m_dfProb, dfParticleMove = d f F i e l d + m_mcinfo_Particle.m_dfProb, dfBackgroundFieldMove = dfParticleMove + m_mcinfo_Background_Field.m_dfProb, dfSfluct = dfBackgroundFieldMove + m_mcinfo_Surface_Fluct.m_dfProb, •dfSSexchange = dfSfluct + m_mcinfo_Surface_Exchange.m_dfProb; REALTYPE r; /* used to determine when to choose which MC move */ int *grn_PtlList, nLeft = m_nParticles, nRP; grn_PtlList = new int[m_nParticles]; f o r ( i n t i = 0; i<m_nParticles; i++) grn_PtlList [i]=i; while(nLeft){ r = RANDOM; i f ( r <= dfField) vTry_Field(); else i f ( r > df F i e l d && r<=dfParticleMove){ /* t h i s code ensures each p a r t i c l e i s chosen once only */ nRP = int( RANDOM * ((REALTYPE)nLeft) ); vTry_Particle(grn_PtlList[nRP]); grn_PtlList[nRP]=grn_PtlList[—nLeft] ; }/* i f Par t i c l e move */ else i f ( r > dfParticleMove && r <=dfBackgroundFieldMove) vTry_Background_Field(); else i f ( r > dfBackgroundFieldMove && r<= dfSfluct) vTry_Surface_Fluct(); else i f ( r > dfSfluct && r<= dfSSexchange) vTry_Surface_Exchange(); }/* while */ vRecord_Sweep(); 51 Chapter 5. Numerical implementation of the local Coulomb algorithm 5.6 Complexity The complexity of each individual update has been established and is given in table 5.1. When performing a Monte Carlo sweep we update O(N) charges. 0(V) plaquettes and 0(LyLz) surface vertices. Thus the total cost of the algorithm is 0(p3N) + 0{V) + 0{LyLz) + 0(LxLyLz) = 0(p3N) + 0(V). For constant density simulations we have that TV oc V. Thus the total cost is 0{p3N) + 0{N). Table 5.1: Computational cost of each Monte Carlo update. 5.7 Sources of error This discretization is subject to many sources of inaccuracy. The sources and means through which we have reduced them are listed below. 5.7.1 Green's function discretization The electrostatic problem is often formulated in terms of a Green's func-tion. For homogeneous media G(r, r') = 4 7 r £ | | r _ r / | | so that the Fourier transform is G(k) = e~^k~2. The discretization of the energy density given in eq. (5.13) gives Plaquette Charge Surface Fluctuation Surface Exchange O(i) 0(P3) 0(1) Q(LX) (5.26) (5.27) e=x.y.z while that for eq. (5.15) leads to e=x.y.z (5.28) 52 Chapter 5. Numerical implementation of the local Coulomb algorithm When Taylor expanded, it may be seen that these inverse Green's functions exhibit errors which scale with k4/k4 and k4 respectively[26]. Thus we expect an error which decays as l / / i 3 for eq. (5.13) and 1/h5 for eq. (5.15). 5.7.2 Al i a s ing The interpolation of a continuous charge distribution onto a finite mesh leads to an important artifact in meshed simulations known as aliasing. There arises a periodic one-body potential, with periodicity of the mesh spacing. This leads to the charges being artificially trapped between vertices and introduces an unnatural dependence of the density on mesh spacing. Aliasing may be reduced through the use of large interpolation functions. Thus for Gaussians we expect this error to decrease as we increase the range of the Gaussian cr. Increasing a comes at a cost as the larger interpolating functions imply a cut in efficiency. 5.7.3 C u t off and charge neutral i ty violations A sum over eq. (5.16) for Gaussian interpolation does not add to one when a cut off is used. Thus one.expects the numerical sum of the charge density on the vertices to be non-zero. This suggests that violations of Gauss's law occur on the order of the total charge interpolated. Since violations of Gauss's law are not acceptable, we must rescale the Gaussians such that the interpolated charge is exact to within floating-point accuracy. This actually increases another source of error. Since the Gaussians do not vanish entirely at the cut off, the self-energy exhibits discrete jumps with periodicity of the grid. The size of these jumps is proportional to the value of the Gaussian at the cut off. By rescaling our distribution we have increased these errors. One means of decreasing these error is to use Gaussian which decay quickly. Thus one expects these discrete jumps to decrease with decreasing a and increasing cut off radius. 5.7.4 Short-range corrections By replacing point-charges by charge-distributions we have introduced yet one more source of inaccuracy. This substitution significantly modifies the short-range interactions. Fortunately we have made use of Gaussian inter-polation functions, so the energy of our simulation exhibits the same form as that in eq. (2.5) for the Ewald reciprocal sum. Thus for homogeneous me-dia, we may reduce these corrections by including the short-range correction form, eq. (2.4) for the Ewald sum. 53 Chapter 6 Analysis of accuracy We have rigorously tested the accuracy of the algorithm for bulk geome-try in terms of the Gaussian parameters a and p. This has been achieved through analysis of the electrostatic Green's function in real space. There are two main means through which the Green's function may be resolved. They are • Using that the pair correlation function g(r — r ' ) « exp{—G(r — r')/kBT} for dilute gases. We may determine the Green's function through logarithmic inversion of g(r — r ' ) . • measuring the minimized energy as a function of particle separation for a two particle system.. The minimization of the energy is performed through transverse updates only. This amounts to quenching the field temperature to OK. For all numerical results we express our data in reduced Heaviside units, where kB = f-o — e = 1. As an example, the Coulomb interaction is these units for two electrons is (47rr) _ 1 . For the finite temperature method, we have employed a cubic unit cell of volume (20a)3.. The boundary term in eq. (5.20) vanishes for lB <C 20a. By choosing a temperature of kBT = ^—- this is satisfied. Thus this method produces an energy equivalent to the Ewald reciprocal sum with tinfoil boundary conditions (e.s = oo). Figure 6.1 shows this for a two-particle system. Superposed on the data are the corresponding Ewald sums. These have been determined via numerical evaluation of eq. (2.5). The evaluated Ewald sums are well converged. Clearly, the statistics from the finite temperature method are too coarse for use in calibrating the discretization error. Each curve is the result of 30,000,000 M C sweeps taking almost 30 C P U hrs. As improving the statistics was not a reasonable option, we chose to perform the calibration using OK simulations. 54 Chapter 6. Analysis' of accuracy r [a] Figure 6.1: Interaction potential for kgT — The symbols are {cr, p 3} = {0.50a, 53} (•), {0.80a, 73} (o), {1.00a, 73} (A), and {1.00a, 73} (V). i.e. Different sized bases and ranges. The solid lines indicate corresponding Ewald sums. The offsets applied to each curve have been included for clarity. 6.1 OK calibration of bulk geometries The OK field method may be described in the following pseudo-code 0) set charge 1 = q:(0,0,0) and charge 2 = -q:(0,0,0) 1) i n i t i a l i z e the mesh with 1 and 2 2) pick a random plaquette 3) increment the plaquette by a random amount between -d and d 4) accept the increment i n 3) i f i t lowers the systems energy 5) repeat 3->4 for each plaquette keeping track of the acceptance rate and the accumulated change i n energy 6i) i f the acceptance rate i s > 60% mult iply d by 1.5 6 i i ) i f the acceptance rate i s < 40% divide d by 1.5 7i) i f the energy = 0.0 then we're done for th i s configuration 7 i i ) i f the | accumulated energy / system energy | < 1E-8 and the acceptance rate i s between 40% and 60% then we are 55 Chapter 6. Analysis of accuracy done for this configuration 7 i i i ) otherwise repeat 2->7ii for at least 10,000 sweeps 8) restore d to i t s o r i g i n a l value prior to the energy minimization 9) output the minimized energy for the given separation 10) step charge 1 and charge 2 away from each other along a particular direction 11) repeat 2->10 u n t i l charge 1 and charge 2 are a distance of 1/2 the linear dimension of the simulation c e l l away from each other t The scaling value of 1.5 for was determined through experimentation. We have conducted numerous OK simulations for cubic cells of size Lx = 20. As discussed, the boundary conditions for these simulations agree with the reciprocal part of the Ewald sum with es = 0. This is because in principle IB is infinite for T = OK. These simulations made use of the improved energy discretization in eq. (5.15). The particles were set to have a charge of e and — e and stepped away from each other along the lattice diagonal (1.1.1)/\/3. Figure 6.2 shows the difference between our generated energy and that of the Ewald methods for varying (a) and constant (b) interpolation base. The data is shown as a continuous line due to the fine resolution achieved. For (<r,p3) equal to (0.65a. 5 3) and (0.75a. 73) we see that the energy clearly exhibits a periodicity of \/3a. This is a result of aliasing due to the use of small a. The induced periodicity from aliasing is significantly reduced for a > 0.90a. For a > 0.90a. the error becomes dominated by short-range corrections. A third source of inaccuracy also becomes apparent for these choices of a. The discontinuities in the energy have periodicity \/3a and hence are periodic with the mesh. These jumps are a result of the residual non-vanishing charge density interpolated at the cut off radius of the Gaussians. By decreasing a/p, the cut off error is reduced. This is observed in the reversal of accuracy for (0.95a, 7 3) and (1.00a, 7 3). The balance between error induced from the use of cut offs and aliasing is clearly non-trivial. To determine the appropriate balance between these errors we have per-formed a full study of the parameter space for a and pA. The results are shown in figure 6.3. Each point represents the difference between the max-imum and minimum error for the energy. For example, the (0.65a, 5 3) 56 Chapter 6. Analysis of accuracy 0 I CO - 1 CO I 3: LU 0 -3 0 10 [a] Figure 6.2: Absolute error of the discretized energy for several parameter combinations, (a), the curves shown are Xcp3) = (0.65a;53), (0.75a, 73), (0.90a, 73), and (1.00a, 93) ordered from worst to best respectively, (b), the curves shown are (o,p3) = (0.80a, 73), (1.00a, 73), (0.95a, 73), and (0.90a, 73) ordered from least to most accurate. 57 Chapter 6. Analysis of accuracy curve in figure 6.2 a) has a maximum error of 8.0 x 10 4 a~ 1 and a min-imum error of —1.6 x 1 0 _ 3 a _ 1 . Thus it is represented by a point with max|(5c7| = 2.4 x 1 0 _ 3 a _ 1 . Table 6.1 gives the optimal value of a for inter-polation bases of 5 3 ; 7 3 . 9 3 and l l 3 . For most practical simulations the use of (0.90a, 7 3) is suggested, as larger base sizes prove expensive. a [a] 0(SU) [a"1] 5 3 0.75 IO" 3 7 3 • 0.90 10~ 4 9 3 1.15 IO" 5 l l 3 1.35 IO" 5 Table 6.1: Optimal parameters for bases of-5 , 7 , 9 , and 11 0.5 CO 3 to X CO E 0.5 0 1 r x 10^ + 5x5x5 + + i + t - + + 0.6 0.8 x 1 d"4 + 9x9x9 + + + + + + 1.2 0.5 0 x10'4 11x11x11 4 + + 1 + * + . 1.4 1 a [a] 1.2 1.4 Figure 6.3: Maximum error in the pair potential for bases of 5 3 to l l " 58 Chapter 6. Analysis of accuracy For comparison, we have also performed simulations using the discretized Hamiltonian in eq. (5.13). We have utilized the optimal parameter set (0.90a, 73) so as minimize all other sources of error. This is shown in figure 6.4. Clearly, the use of eq. (5.13) introduces a large correction at short ranges. This amounts to an increase in max|c>£/| by an order of magnitude. This is an undesired result but unavoidable for inhomogeneous media. 0 2 4 6 8 10 r ' [a] Figure 6.4: Pair potential accuracy for o = 0,90a and p3 = 7 3 using eq. (5.15) (dashed) and eq. (5.13) (solid). 6.2 OK calibration of quasi-2D geometry As described in section 5.4 we have employed a tetragonal unit cell to simulate quasi-2D systems. The unit cell is elongated in the x-direction and has dimension aLx = anL, aLy = a,L and aLz = ah. We use L = 15 and aspect ratios of n = 1,3,5 We have performed OK simulations for this configuration, as in the previous section. The boundary term, which is equivalent to a Ewald dipole term with es = 0, has been analytically subtracted and the plane-wise summation term in eq. (2.8) has been added. Thus our method achieves approximate equality to the 2D-periodic Ewald eq. (2.10). For comparison, we have evaluated eq. (2.10) numerically. Figure 6.5 shows a comparison of our quasi-2D pair potential to that of eq. (2.10). We displace the particles along the aperiodic direction in a), as well a periodic direction in b). Clearly for n = 1, the pair potential exhibits unacceptable inaccuracies in the aperiodic direction. For ratio n = 59 Chapter 6. Analysis of accuracy 3, we find that the pair potential is reasonably accurate in both the periodic and aperiodic directions. We would expect that as with 3D-periodic Ewald methods the accuracy should increase with the aspect ratio, n. However, for n larger than 3 the accuracy begins to decrease. This could be caused by a number of contributing factors: • The empty space is not being effectively updated. The use of a global field update should resolve this problem but its implementation in-creases the scaling of the algorithm and reduces the efficacy of the method. • Additional corrections to the 3D-periodic Ewald method need to be included. These corrections are described in [28, 29] but have not been implemented. The accuracy achieved is acceptable for most simulations but the quasi-2D method exhibits a significant efficiency problem. Buffering the system with empty media comes at a heavy price. For example in figure 6.5 b), the n = 1 simulation took 0.6 hours, the n = 3 simulation took 6.7 hours, and for n = 5 the simulation took 13.3 hours. This is due to the significant overhead of minimizing the field in empty space. A global update of the field may resolve this efficiency problem, however, we have not implemented one at this time. 60 Chapter 6. Analysis of accuracy CO 0 Q T3 CO LU Figure 6.5: Pair potential accuracy for quasi-2D geometry, o = 0.90a and p3 = 73. Charges separated along the aperiodic direction, a), n=l(dotted), 3(solid),5(dashed). In a) the n = 3 and 5 fall on top of each other. Charges separated along the periodic, b), n=l(dotted), 3(solid), 5(dashed). 61 Chapter 7 A n Applicat ion: The Gouy-Chapman problem We now apply the algorithm to the Gouy-Chapman problem as a means of demonstrating the algorithm's accuracy and versatility. We compare our results to the analytic solutions derived in chapter 4. Consider, as in chapter 4, a system with N point charges per volume V, confined between charged planar walls. The charges have charge — q and the walls have an average charge density as. Electroneutrality then gives qN + as Y i dS, = 0 • (7.1) i=l , 2 For a tetragonal unit cell with charged walls at x — 0 and x = d, this takes the form -qN + 2 C T L V = 0 (7.2) with aL the linear dimension of one of the square faces. By applying PBCs to the y and z directions, we achieve efficient convergence to the thermody-namic properties of a thin film. The coupling parameter, Gouy-Chapman and Bjerrum lengths for this system are d = d/n (7.3) „2 4 " ^ ( 7 4 ) In 1 _ a2L2 OSZB ~ * PB ^ NntB2 H a2L2 We have simulated this setup using the method described in chapter 5. The improved energy discretization in eq. (5.15) is used along with the Ewald short-range correction term (2.4). Since we cannot allow the charge to interpolate on to the walls, we extend the simulation cell by ap/2 at 62 Chapter 7. An Application: The Gouy-Chapman problem x = 0 and x = d. We then implement equipotential boundary conditions at x — —cip/2 and x = ap/2 + d. as indicated in figures 7.1 and 7.2. We thus have that aLx = d + ap. This is approximately equivalent to extending the model system by ap, but still restricting the particles to the region 0 < x < d. As —ap/2 < x < 0 and d < x < d + ap/2) remain empty, V 20(r) = 0 in these regions. Thus the potential at x = 0 and x = d is modified by at most a linear term. For all of the simulations in this chapter we use Gaussian interpolation functions with o = 0.90a and base 7 3. d 1 1 d Figure 7.1: Implementation of conducting charged walls. Figure 7.2: Implementation of equipotential walls. 63 Chapter 7. An Application: The Gouy-Chapman problem 7.1 Weak coupling (S < 1) Table 7.1 shows the simulations and corresponding Gouy-Chapman pa-rameters used in the charge density plot in figure 7.3. We have utilized constant-potential boundary conditions with zero potential at the charged walls. Since the charged walls are both at zero potential we expect a mini-mal effect from surface exchanges. The plates thus carry an average charge equivalent to the constant charge boundary conditions. Thus, we expect agreement with eq. (4.54). The rescaled fugacity, A, was resolved numer-ically for each system using eq. (4.56) in Maple. For clarity the data has been truncated prematurely for d = 7.6, 10.6 and 12.6. The corresponding weak-coupling expressions, PPB-, are superposed on each curve. The density plots are a result of 100,000 Monte Carlo sweeps. For the first 100 sweeps the plaquette, surface charge, charge exchange and particle step maximums are modified in a similar fashion to steps 6i and 6ii in the OK method. In this manner we achieve efficient convergence. In addition we begin accumulating the density after this optimization period. Figure 7.3, clearly shows that we have achieved excellent agreement in this limit. N L kBT [qh-'a-1] d[a] A time (hrs) 184 38 1.59 x 10" 1 4.6 0.23638 19.3 184 38 1.59 x 10" 1 7.6 0.10862 29.1 184 . 38 1.59 x 10" 1 10.6 0.06256 36.3 184 38 1.59 x 10" 1 13.6 0.04069 73.1 Table 7.1: Parameters used for H = 0.1 simulations. 64 Chapter 7. An Application: The Gouy-Chapman problem 1.2 0 I . . . . - 1 _ 1 1 2 ' 3 4 x [ j i ] Figure 7.3: Simulated ion density in the weak coupling limit (EE = 0.1) for d=4.6(D), 7.6(o), 10.6(A), and 13.6(V). The corresponding analytic solutions, PPB- obtained from eq. (4.54) are shown as lines. The inset shows the difference between the simulated results and the analytic solutions. 7.2 Strong coupling (S 3> 1) Unfortunately, the,two-wall strong coupling solution in eq. (4.68) results in a homogeneous density and is only valid for small d. For this regime we do not expect agreement since we have implemented the charge surfaces not at x = 0 or x = d but at x = —ap/2 and x = d + ap/2. In the small d regime we therefore expect artifacts from this approximations. Thus we have focused on the large d limit. In this limit, the two walls decouple and so the expected charge density is the simple exponential given in eq. (4.64). This is shown as the solid line in figure 7.4. We have simulated numerous systems varying d; the counter-ion density is shown in figure 7.4. The parameters used are given in table 7.2. We use the same method as for the weak coupling simulations, but perform 400,000 Monte Carlo sweeps to achieve better resolution near the charged interface. The exception is 65 Chapter 7. An Application: The Gouy-Chapman problem d= 28.75, for which we only performed 200,000 Monte Carlo updates due to the required simulation time. In addition, we have utilized an optimization described in the next section to limit the duration of the simulations . Rather than perform a complete sweep of all plaquettes we only update 10% per sweep. N L d P [ « ] kBT [ g 2 e " 1 a - 1 ] time (hrs) 199 50 115.00 0.20 3.9790 x IO" 3 33.3 199 100 57.50 0.40 1.9890 x IO" 3 56.3 106 146 28.75 0.80 9.9470 x 10~ 4 90.0 Table 7.2: Parameters used for E — 100.0 simulations. 1.8 0 i . . . 1 0 0.5 1 1.5 2 x [Hi Figure 7.4: Simulated ion density in the strong coupling limit (E = 100.0) for £2=115 (V) , 57.5 (o), and 28.75 (•) . The solid line is the solution in eq. (4.64)'for a single charged wall. The simulation parameters used are given in table 7.2. The agreement with eq. (4.64) becomes increasingly better as the coupling between the walls is decreased. 66 Chapter 7. An Application: The Gouy-Chapman problem Figure 7.4 clearly shows that as d is increased, the agreement with the exponential correspondingly becomes better. In-the extreme case of d=115, we show excellent agreement with eq. (4.64). Thus we have verified that for S > 1, we achieve agreement with the expected Gouy-Chapman solution. 7.3 Optimization The long durations in table 7.1 occur because we chose to update every plaquette per Monte Carlo sweep. The simulation times in table 7.2 are only comparable to those in table 7.1 because rather than updating each plaque-tte we have only updated 10%. For example, in the strong coupling \i = 0.8a simulation, the volume contains an extraordinary 1,939,756 plaquettes. In comparison, the largest weak coupling simulation contains only 326,344 pla-quettes. By updating only approximately 200,000 plaquettes per sweep for the strong coupling simulation, we have reduced the duration to be compa-rable to the weak coupling. As with other mesh reliant methods, the volume dependent scaling cofactor can be quite significant. To reduce this cofac-tor we have conducted simulations with different percentages of the total plaquettes updated per sweep. The results are shown in figure 7.5 with an identical setup to that of the largest weak coupling simulation. Figure 7.5, shows no difference between the 10% and 90% simulations! This is not an unexpected result because the surface fluctuations, exchanges, and particles moves are all coupled to the field variables. Thus a significant proportion of the integration over field degrees of freedom comes from these updates. For example in the weak coupling d — 13.6 simulation, for each Monte Carlo sweep we perform 184 particles updates, 2 • 382 surface fluctuations, and 382 surface exchanges. This results in 184 • 7 3 + 3 • 2 • 382 + (68 + 7)382 = 260, 076 field updates or 80% of all field variables. For the \.i = 0.8a strong cou-pling simulation, the same calculation yields 41%. Thus a reduction in the proportion of plaquettes updated is clearly acceptable. 67 Chapter 7. An Application: The Gouy-Chapman problem „ 0.8 CO w CM Figure 7.5: Weak coupling density for d=13.6. The proportion of total pla-quettes updated per sweep have been varied. Shown are 0.1%(x), 10%(D), 30%(o), 50%(A), 70%(V), and 90%(0). The solid line shows the corre-sponding analytic solution, pps, obtained from eq. (4.54). The inset shows the difference between the simulated results and the analytic solution. 7.4 Finite-size artifacts To ensure that our simulations have indeed converged to the thermody-namic limit, a series of strong coupling simulations were performed. These simulations all have p = 0.8a, d = 28.75 and 3 = 100. However, the num-ber of charges and area of the charged walls has been decreased. Thus any difference can be attributed to finite-size artifacts or statistics. The exact parameters are given in table 7.3. The method is the same as in the'previous strong coupling section. Figure 7.6, shows the percent difference between the N = 106 charge density and that of N = 21, 43, 58, 74, and 92. The percent difference clearly levels off to about 2% for N=92. The ex-pected order of the statistical error is around 1% to 2%. To understand this, 68 Chapter 7. An Application: The Gouy-Chapman problem consider the density plot in figure 7.4. For x > 1.5/z the density is approxi-mately 0.15 2TT£B°'2- For the curves in figure 7.6 we binned approximately 16,000 configurations. Using a bin size of 0.05//, we expect a total count of 0.15 • 0.05 • 16000 • N = 1207V used to generate each point for x > 1.5/x. The expected statistical percent error is thus roug hly ,/2/1207V « 1.3% for N=92. Since the total percent error is on this order, any finite-size effects are thus negligible. TV" L kBT { qh-'a-1 } time (hrs) 21 65 9.9470 x 10" 4 8.9 43 93 9.9470 x IO" 4 19.3 58 108 9.9470 x IO" 4 25.9 74 . 122 9.9470 x 10~ 4 35.1 92 136 9.9470 x IO" 4 40.8 106 146 9.9470 x IO" 4 90.0 Table 7.3: Parameters used in testing finite-size artifacts. For JV > 92, finite-size artifacts become negligible. , O 0.5 1 1.5 2 • * [ H ] Figure 7.6: Percent difference between the N = 106 density and 7V= 21 (x) , 43(D), 58(o), 74(A), and 92(V). 69 Chapter 7. An Application: The Gouy-Chapman problem 7.5 Pressure between equipotential walls Having established the accuracy of our implementation, it is interesting to see whether equipotential boundary conditions yield an attractive regime between grounded walls. We have conducted such a study using the pa-rameters in Table 7.4. It should be noted that the particle to Monte Carlo sweep ratio has been kept approximately constant to keep the resolution unchanged. Recall the contact value theorem P = /5(d) - 1 (7.7) 2ir£Ba2s Using the symmetry of the setup we have that /5(d) = p(0). To evaluate the pressure we generate rescaled density plots as in the previous sections and linearly extrapolate the first two points to x = 0, from which we subtract one. Figure 7.7 shows the resulting pressure for the simulation parameters in table 7.4. We have focused on 6 < S < 70 and 2 < d < 16 as this has been shown to be an area of particular interest for homogeneous charged walls[2]. This data has been gridded and paletted to generate the phase diagram in figure 7.8. All regions in figures 7.7 and 7.8 with P < 0 correspond to an attraction between the like-charged walls. Not surprisingly for large sepa-ration and small coupling the pressure between the interfaces approaches zero. Amazingly the system is stable or even attractive for small d in the case of large E. This exemplifies the importance of the counter-ions in this phenomena. 70 Chapter 7. An Application: The Gouy-Chapman problem 77 p [a] N L kBT [ g V ' a " 1 ] d Sweeps (103) 6 1.0 138 51 1.32629 x 1 0 " 0 2 2 210 6 1.0 138 51 1.32629 x I O " 0 2 3 210 6 1.0 138 51 1.32629 x I O " 0 2 4 210 6 1.0 138 51 1.32629 x I O " 0 2 6 210 6 1.0 138 51 1.32629 x I O " 0 2 8 210 6 1.0 138 51 1.32629 x I O " 0 2 11 210 6 1.0 138 51 1.32629 x I O - 0 2 16 210 10 1.0 156 70 7.95775 x I O " 0 3 2 190 10 1.0 156 70 7.95775 x I O " 0 3 3 190 10 1.0 156 70 7.95775 x I O " 0 3 4 190 10 1.0 156 70 7.95775 x I O " 0 3 6 190 10 1.0 156 70 7.95775 x 10~ 0 3 8 190 10 1.0 156 70 7.95775 x 10~ 0 3 11 190 10 1.0 156 70 7.95775 x I O " 0 3 16 190 16 1.0 106 73 4.97359 x I O " 0 3 2 270 16 1.0 106 73 4.97359 x I O - 0 3 3 270 16 1.0 106 73 4.97359 x I O " 0 3 4 270 16 1.0 106 73 4.97359 x I O - 0 3 6 270 16 1.0 106 73 4.97359 x I O " 0 3 8 270 16 1.0 106 73 4.97359 x I O " 0 3 11 270 16 1.0 106 73 4.97359 x I O " 0 3 16 270 26 1.0 97 89 3.06067 x 10~ 0 3 2 300 26 1.0 97 89 3.06067 x 10~ 0 3 3 300 26 1.0 97 89 . 3.06067 x I O - 0 3 4 300 26 1.0 97 89 3.06067 x I O " 0 3 6 300 26 1.0 97 89 3.06067 x I O " 0 3 8 300 26 1.0 97 89 3.06067 x I O " 0 3 11 300 26 1.0 97 89 3.06067 x I O " 0 3 16 300 43 0.5 133 67 3.70128 x I O " 0 3 4 220 43 0.5 133 67 3.70128 x I O " 0 3 8 220 43 0.5 133 67 3.70128 x I O " 0 3 12 220 43 0.5 133 67 3.70128 x I O " 0 3 16 220 70 0.5 144 89 2.27364 x 1 0 - ° 3 4 200 70 0.5 144 89 2.27364 x I O - 0 3 8 200 70 0.5 144 89 2.27364 x I O " 0 3 12 200 70 0.5 144 89 2.27364 x I O " 0 3 16 200 Table 7.4: Parameters for pressure phase diagram simulations. Chapter 7. An Application: The Gouy-Chapman problem Reference [1] gives numerous examples of p and S for charged biomate-rials. For charged membranes with divalent counter-ions one finds p = 1.1 Angstroms and 5 = 24.8. Figure 7.8 suggest an attractive regime for a pair of these membranes separated by a distance of approximately 3 Angstroms to 5 Angstroms. In comparison charged membranes with monovalent counter-ions have 5 = 3.1 and so we do not expect an attraction at all. Membranes .with trivalent counter-ions have E! = 83.7 and p = 0.8 Angstroms. Thus we expect attractions to be more predominant for biomembranes with trivalent counter-ions. Spermidine D N A has similar properties to a biomembrane with trivalent .counter-ions. It has S = 75.6 and p = 0.8 Angstroms. In-terestingly, it also has a radius of curvature that is 16Ap. This suggests that a local planar approximation is not too bad for spermidine D N A . Thus our simulations would suggest that the attraction mediated by counter-ions provides a mechanism for the bundling of this form of D N A . This study is an important unique result. As far as we are aware, it is the only study conducted on the Gouy-Chapman problem where constant-potential boundary conditions have been fully implemented. It clearly shows that attractions betweens like charged equipotentials are possible. For charged biomembranes with trivalent counter-ions as well as spermidine D N A it sug-gests a significant attractive regime. In particular,, it presents a possible mechanism for the binding of like-charged biomembranes and the bundling of spermidine D N A . The model we have considered neglects many details; however, the results of this study highlight the strength and versatility of our implementation. 72 Chapter 7. An Application: The Gouy-Chapman problem 0.4 -0.2 1 " V 1 i , 1 ; X § V X • G v. X • OEM G O DM> G $ . G A A 0 • A 0 1 . 0 2 6 10 14 Figure 7.7: Pressure between the like-charged equipotential walls as a func-tion of d for for H = 6 ( V ) , 10 (x) , 16 (•) , 26 (o), 43 (A) , and 70 (0). The parameters used are given in table 7.4. For the parameters shown no ana-lytic theory is applicable. This is a unique numerical study of the two-wall system with equipotential boundary conditions. 73 Chapter 7. An Application: The Gouy-Chapman problem 0.2 CM to 0 o 2 -0.2 16 16 32 64 Figure 7.8: Attractive and repulsive regime between like charged equipo-tential planar walls for 6 < s < 70 and 2 < d < 16. 71 Chapter 8 The electric double layer: Inclusion of an inhomogeneous dielectric background We now turn our attention to the effect of inhomogeneities in the dielec-tric background for interacting electric double layers. We study a model system of equal number of monovalent cation and anion species in solution. The ions are confined between two planar constant-charge conducting walls with zero total charge, separated by a distance d. Rather than being point particles, the ions are modeled as purely repulsive spheres of diameter D. The solvent is a neutral polarizable medium; it is modelled using a dielec-tric background. The solvent may be considered homogeneous for the most part with dielectric value €2, except near the conducting surface. Along the conducting boundaries, a layer of adsorbed solvent forms what is known as a Stern layer. The ions may not penetrate the Stern layer, which is usually less polarizable than the bulk solvent. We model the system's Stern layers using slabs of dielectric value ej as shown in figure 8.1. The thickness of the Stern layers are equal to the diameter of the neutral solvent molecules, which we have taken to be D for simplicity. This means the ions have a distance of closest approach to the walls of 3/2 D. This model may be understood as an electrolyte between zero charged conducting plates. It is an extension of the model studied by Henderson et al. in reference [18] where one planar wall was considered. 8.1 Implementat ion Figure 8.1 illustrates our method of simulating the two-wall model. In particular, the figure shows a 2D sketch of a possible setup for mesh spacing a = 1/2 D. As shown, the nodes are offset by a half lattice spacing from the 75 Chapter 8. The electric double layer: Inclusion of an inhomogeneous dielectric background E l ! 0 ® £ 2 © © © j i © ; ! 0 © ! © © © © ; h - ^ d ~ 1 V ( — ) J Figure 8.1: Illustration of the simulation setup. The bulk has a dielectric value of i%. The light grey regions have a dielectric value of e\. The dark grey indicates zero constant-charge conducting walls. The anions (©) and cations (©') are bound between the dotted lines. As indicated, we extend the simulation volume by half a mesh spacing on either side. reference frame so that the first plane of nodes sits at $ = —a/2 and the last at x = d + a/2. This offset is simply a matter of computational convenience. In general, we consider the regions —a/2 < x < 0 and d < x < d+ a/2 to be conducting but zero charged. The regions 0 < x < D and d — D < x < d are used to model the system's Stern layers and thus have a corresponding dielectric value of t\. We simulate this by implementing constant-charge conducting surfaces at x = —a/2 and x = d + a/2, as described for the slab geometry in chapter 5. In addition, the nodes with a/2 < x < D and d — D < x < d — a/2 take on the dielectric value t%. All other nodes have a dielectric value of eg. The simulation volume is thus (d + a) x aLy x aLz. We use a restricting function f l , V2D<x<d-ZI2D I 0. otherwise 76 Chapter 8. The electric double layer: Inclusion of an inhomogeneous dielectric background to ensure the ions remain between the Stern layers. For the cations and anions we use Gaussian interpolating functions with a = 0.9a and p 3 = 73. The spherical repulsive interactions are implemented using the zero-shifted cut off Leonard-Jones pair potential: Vu(r, M(?)'M?)6}+<> 0 < r < 2 i D ( g 2 ) [0, otherwise where we have used c — (4TVD) . This potential applies a smooth repulsion starting at r = 2&D and becoming dominant at 7- = D. For the simulations performed, we use €2 = 80eo and a temperature given by kBT = - i - A fn ^ = 4.288 x l O " 4 - ^ - ' (8.3) 2.32 47r80e0£> e0D v ; as in [18]. For an ion diameter of D = 3 Angstroms this corresponds to a temperature of 300 K . We use a simulation cell with d = 25D and aLy = aLz = 9D, containing 64 anions and 64 cations. This corresponds to an electrolyte solution of Molarity 1.94M. The Molarity of the solution is the number of Moles of cation or anion species per liter. For an aqueous solution, we take e\ = 6eo as in [18]. For comparison we also consider systems with ej = 40eo; 80eo and 120eo- Determining an appropriate mesh spacing is a non-trivial problem as it affects both accuracy and efficiency. We will address this issue in the following section. 8.2 Effect of mesh spacing Since our model exhibits spatial dielectric inhomogeneities, we cannot use the short range Ewald correction in eq. (2.4) or the improved energy discretization in eq. (5.15). Thus we expect the interactions mediated by the field to be incorrect at short ranges. By fortuitous choice of the mesh spacing a with respect to the ion diameter D , we may ensure that this incorrect potential is never sampled. This is possible only if any corrections vanish at a distance less than the ion diameter. By choosing an appropriate mesh spacing a, we ensure that the regime where any short range corrections need to be applied is dominated by the repulsive interaction in eq. (8.2). Thus the choice of D/a is extremely important, determining the extent that the incorrect short range potential affects the simulation. 77 Chapter 8. The electric double layer: Inclusion of an inhomogeneous dielectric background To determine an appropriate value of D/a, we have performed two par-ticle simulations with the field temperature quenched, as in chapter 6. The simulations use a cubic volume of (5.D) 3, bulk geometry and a homoge-neous dielectric background of eo- As before the Gaussian charge clouds are stepped away from each other along the diagonal (1, 1, l) /- \ /3. Recall that the expected pair potential is given by eq. (2.5). In addition, we expect a spherical dipole term given by eq. (2.7), with es = 0. Figure 8.2 shows simulated pair potentials for a = 1/2 D , 1/3 D, 1/4 D, and 1/5 D . 0 ^ -0.1 I Q 3 -0.2 -0.3 1 1 1 1 1 1 0 0.5 1 1.5 2 2.5 r [D] Figure 8.2: Bulk pair potential generated for a simulation volume of (5D)3 using the OK method described in chapter 6. The number of nodes inter-polated over has been kept constant at 7 3 . The Gaussian range varies with the mesh spacing, having a value of 0.90a. The curves shown are a = 1/2 D (dashed-dots), 1/3 D (dots), 1/4 D (bold-dashes) and 1/5 D (light-dashes). The grey region indicates where the spherical repulsive potential in eq. (8.2) is significant. The energy is discretized via eq. (5.13). We have kept the interpolation base fixed at p3 = 7 nodes and used or = 0.45 D, 0.3 D, 0.225 D, and 0.18 D corresponding to 0.90a. The discrete 78 Chapter 8. The electric double layer: Inclusion of an inhomogeneous dielectric background energy is calculated using eq. (5.13). The grey area indicates the region where r is less than or equal to D; this is the portion of the electrostatic pair potential that will be ignored when the repulsive interaction in eq.(8.2) is turned on. For a = D/4 the pair potential has converged beyond r = D to that of a = D /5 . This suggests that by using a discretization with a = D/4, we may attain sufficient accuracy at scales greater than D. To confirm this choice of Dfa we have performed multiple simulations of a 1.94M electrolyte, bound by conducting plates as described in the previous section. For simplicity, these simulations use 6] = £2 = 80eo- Our expecta-tion for the density of this system is a flat interior density with accretion near the conducting boundaries. The accretion comes from a decreased energy cost for having strong displacement fields near the charged bound-ary. Figure 8.4 shows the cation/anion density as a function of the distance from one of the conducting walls. These simulations represent a preliminary study. Clearly, the statistics are far too coarse. Close to the charged wall, the density appears to have converged; however, in the interior the density shows strong fluctuations. These fluctuations are most likely due to a de-creased ion diffusivity as the mesh spacing is decreased. For a = 1/2 D, 1/3 D and 1/4 D we have performed 500,000 M C sweeps. For a = 1/5 D we have performed 1,200,000 M C sweeps taking 107.5 C P U hrs. Thus achieving excellent convergence for a — 1/5 D remains a significant technical hurdle. However, these simulations do suggest that a mesh spacing of a = 1/4 D may be acceptable. 79 Chapter 8. The electric double layer: Inclusion of an inhomogeneous dielectric background 4.2 CO I Q — 3.6 <£' 3.0 2.4 x 10" A A ^ A A A A £ 8 0 B ^ O o o O ^ 0 ° ° ° o O O o o o @ § @ @ @ @ A A 4 6 8 x [ D ] 10 12 Figure 8v3: The cation/anion density for systems with homogeneous dielec-tric background, e\ = €2 = 80 eo, and planar conducting walls. Simulations with different mesh spacing are shown: a = 1/2 D (•), 1/3 D (o), 1/4 D (A) and 1/5 D (0). The plate separation is 25 D . The volume of the simulation cell is 25.2573 x 979 x 9D, containing 64 cations and 64 anions. 8.3 The Stern layer Using a mesh spacing of a = 1/4 73, we have performed a study of the affect of the dielectric value of the Stern layer on the 1.94M electrolyte system. We use the parameters and implementation described in section 8.1 and perform 1,200,000 MC sweeps. The cation/anion densities for these simulations are shown in figure 8.4. These results parallel the results of the single wall system studied in [18]. The ion density is strongly modified near the boundary by the dielectric value of the Stern layer. For a weakly polarized Stern layer (ej — 6eo) there is a strong depletion of the ion density close to the walls. In contrast, the e\ = 80eo and ej = 120eo simulations show a strong ion accumulation. The density is nearly flat for e] = 40er> These densities may be understood qualitatively in terms of image charges. The 80 Chapter 8. The electric double layer: Inclusion of an inhomogeneous dielectric background induced image charge at a dielectric interface for an ion has the same sign as the ion if the dielectric value decreases and the opposite sign if it increases. If we focus on a single wall, we see that each ion induces a charge at the ej-62 interface and at the ei-wall interface. For ej < €2, there is a repulsive interaction with the Stern layer interface and an attractive interaction with the wall. Thus one expects a depletion near the wall depending on the ratio ei/e2- For ei > 62, each ion exhibits two attractive interactions with the Stern layer and conducting wall interfaces. Thus it is not surprising that we see an accumulation of ions near the conducting wall. Alternatively, our algorithm gives a very nice qualitative explanation of. these results. We may think of each ion as a source of displacement field. The energy density D(r)/2e(r) suggests that the field energy associated with an ion is reduced in a region of high dielectric and increased in a low dielectric region. Thus we expect an accumulation of ions near the conducting walls. This affect is counteracted by the Stern layer if ei < £2-Though these simulations show reasonable agreement with the work in [18], it still remains a significant technical challenge to reduce the statistical errors. The message, however, is clear; the Stern layer for water results in a .depletion of the ion densities near conducting walls. We leave the challenge of optimizing convergence to future work. 81 Chapter 8. The electric double layer: Inclusion of an inhomogeneous dielectric background 3.6 [ q j ^ ^ Q ^ o x p a o o o c G o c o c c 0 ; ) 2.4 „ 3.6 CO Q Q. 2.4 3.6 2.4 3.6 2.4 x 10" o d) G G O x 10" C) v v V . , -,/ V V v V V V V v X 10 -2 b) ^ A A A ^ A ^ A A A A A A A A A A A A A . A x 10 -2 a) • • • • • • • • • • • • 2 4 6 8 10 12 x [D] Figure 8.4: The charge density for systems with Stern layers of thickness D = 4a and dielectric coefficient e\ — 6eo (a). 40eo (b). 80eo (c), and 120eo (d). The temperature used is kBT = 4.288 x l C T 4 e 2 e 0 1D~1 and e2 = 80e0-The plate separation is 25D. The volume of the simulation cell is 25.2573 x 9D x 973. containing 64 cations and 64 anions. 82 Chapter 9 Conclusion A reformulation of the Coulomb problem using auxiliary fields has been extended for the first time to slab and quasi-2D geometries. It has been im-plemented using Metropolis Monte Carlo and Gaussian interpolation func-tions. We have established the accuracy of the algorithm by generating effective pair potentials in chapter 6. For an interpolation base of 7 3 nodes the optimal Gaussian range is 0.90a and has an absolute error in the pair potential less than 2.4 x 10~ 4 a _ 1 . Using this implementation, the Gouy-Chapman problem was numerically resolved for constant-potential slab boundaries. In the low coupling limit we have found that the counter-ion density for constant-potential bound-ary conditions shows excellent agreement with the analytic Gouy-Chapman solution. In the high coupling regime we found agreement with the ana-lytic theory in the limit of large wall separation. Using the contact value theorem, we have conducted a thorough study of the phase behaviour of the pressure experienced by the equipotential walls in terms of the coupling constant £ and the rescaled wall separation d. The parameter space we have considered includes 6 < 5 < 70 and 2 < d < 16, pertaining to many interesting biomaterials ranging from monovalent biomembranes to spermi-dine DNA. The numerical results show that like-charged constant-potential slab boundaries may attract each other for E! > 20. This is a first result for constant-potential slab boundaries in the Gouy-Chapman problem. Finally, we have extended the implementation further to allow inhomoge-neous dielectric backgrounds. In chapter 8, we have simulated a reasonably realistic electrolyte system bounded by isolated conducting electrodes. The intent of the simulations was to investigate the role of the Stern layer, a thin layer of adsorbed solvent molecules, as well as to show the versatility of the algorithm. We have shown that a reduced dielectric value in the Stern layer, from the bulk dielectric value, leads to a depletion of ions near the electrodes. This study is a preliminary study so we do not expect full quantitative agreement with experiments. 83 Bibliography [1] A . Naji, S. Jungblut, G. Moreira, and R. Netz, Physica A 352, 131 (2005). [2] A . Moreira and R. Netz, in Novel Methods in Soft Matter Simulations (2004), pp. 245-278. [3] C. Holm and K . Kremer, A r X i v Condensed Matter e-prints (2002), cond-mat/0203541. [4] R. Netz, The European Physical Journal E 352, 557 (2005). [5] A . Moreira and R. Netz, Physical Review Letters 87, 078301 (2001), [6] C. Sagui and T. Darden, in Annual Review of Biophysics and Biomblec-ular Structure- (1999), pp. 155-179. [7] D. York and W . Yang, Journal of Chemical Physics 101, 3298 (1994). [8] Y . Shan, J . Klepeis, M . Eastwood, R. Dror, and D. Shaw, Journal of Chemical Physics 122, 054101 (2005). [9] W . Press, S. Teukolsky, W . Vetterling, and B . Flannery, Numerical Recipes in C (Cambridge University Press, New York, 1992), 2nd ed. [10] C. Sagui and T. Darden, Journal of Chemical Physics 114, 6578 (2001). [11] A . Tikhonov, Journal of Chemical Physics 124, 164704 (2006). [12] I. Yeh and M . Berkowitz, Journal of Chemical Physics 111, 3155 (1999). [13] S. De Leeuw and J . Perram, Molecular Physics 37, 1313 (1979). [14] A . Grzybowski, E . Gwozdz, and A . Brodka, Physical Review B 61, 6706 (2000). [15] A . Grzybowski and A . Brodka, Molecular Physics 101, 1079 (2003). [16] M . Mazars, Molecular Physics 103, 1241 (2005). 84 Bibliography [17] M . Mazars, Molecular Physics 103, 675 (2005). [18] D. Henderson, D. Gillespie, T. Nagy, and D. Boda, Molecular Physics 103, 2851 (2005). [19] J. Jackson, Classical Electrodynamics, (John Wiley and Sons, Inc., 1999), 3rd ed. [20] D. Chapman, Philos. Mag. 25, 475 (1913). [21] A. Maggs and V. Rossetto, Physical Review Letters 88, 196402 (2002). [22] A. Maggs, Journal of Chemical Physics 117, 1975 (2004). [23] L. Level, F. Alet, J. Rottler, and A. Maggs, Pramana - Journal of Physics 64, 1001 (2005), [24] A. Maggs, Journal of Chemical Physics 120, 3108 (2004). [25] J. Rottler and A. Maggs, Journal of Chemical Physics 120, 3119 (2004). [26] A. Maggs, Physical Review E 72, 040201 (2005). [27] J. Rottler and A. Maggs, Physical Review Letters 93, 170201 (2004). [28] A. Arnold, J. Joannis, and C. Holm, Journal of Chemical Physics .117, 2496 (2002). [29] A. Arnold, J. Joannis, and C. Holm, Journal of Chemical Physics 117, 2503 (2002). [30] N. Ashcroft and N. Mermin, Solid State Physics (W. B. Saunders Com-pany, 1976). [31] L. Levrel, Ph.D. thesis, L'Universite Paris 6 (2006). [32] J. Israelachvili, Intermolecular and Surface Forces (Academic Press, London, 1990). [33] G. Moy, B. Corry, S. Kuyucak, and S. Chung, Biophysical Journal 78, 2349 (2000). [34] B. Corry, S. Kuyucak, and S. Chung, Biophysical Journal 84, 3594 (2003). [35] R. Peitzsh and M . Eisenberg, Biophysical Journal 68, 729 (1995). 85 Bibliography [36] N . Metropolis, A. Rosenbluth, M . Rosenbluth, and A. Teller, Journal of Chemical Physics 21, 1087 (1953). [37] L. Reichl, A Modern Course in Statistical Physics (John Wiley and Sons, Inc., Toronto, 1998), 2nd ed. 86
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The local Coulombic Monte Carlo algorithm : applications...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The local Coulombic Monte Carlo algorithm : applications to the electric double layer Thompson, David 2007
pdf
Page Metadata
Item Metadata
Title | The local Coulombic Monte Carlo algorithm : applications to the electric double layer |
Creator |
Thompson, David |
Publisher | University of British Columbia |
Date Issued | 2007 |
Description | A reformulation of the Coulomb problem, using a local Coulomb algorithm based on auxiliary fields, has been extended to slab and quasi-2D geometries. It has been implemented using Metropolis Monte Carlo and Gaussian charge interpolation functions. We have established the accuracy of the algorithm by generating effective pair potentials. Using this implementation, the Gouy-Chapman problem was numerically resolved for constant potential slab boundaries. In the low coupling limit, we find excellent aggreement with analytic solutions. In the high coupling regime, we find agreement with the analytic theory in the limit of large wall separation. Using the contact value theorem, we calculate the pressure experienced by like-charged equipotential walls. The parameter space we consider pertains to many interesting biomaterials ranging from monovalent biomembranes to spermidine DNA. The numerical results show attractions mediated by counter-ions between the like-charged equipotential slab boundaries. We also extend the implementation to allow for inhomogeneous dielectric backgrounds. The effect of a thin adsorbed layer of solvent is considered for an electrolyte system bounded by isolated electrodes. We show that a reduction in the dielectric value of this adsorbed layer results in a depletion of ions near the electrodes, even though the electrodes carry zero total charge. The applications considered show the versatility and accuracy of our implementation. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-03-11 |
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.0084913 |
URI | http://hdl.handle.net/2429/32355 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2007-0615.pdf [ 6.14MB ]
- Metadata
- JSON: 831-1.0084913.json
- JSON-LD: 831-1.0084913-ld.json
- RDF/XML (Pretty): 831-1.0084913-rdf.xml
- RDF/JSON: 831-1.0084913-rdf.json
- Turtle: 831-1.0084913-turtle.txt
- N-Triples: 831-1.0084913-rdf-ntriples.txt
- Original Record: 831-1.0084913-source.json
- Full Text
- 831-1.0084913-fulltext.txt
- Citation
- 831-1.0084913.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-0084913/manifest