UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Nonlocal dielectric response of water Badali, Matthew 2012

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2012_fall_badali_matthew.pdf [ 600.56kB ]
JSON: 24-1.0073102.json
JSON-LD: 24-1.0073102-ld.json
RDF/XML (Pretty): 24-1.0073102-rdf.xml
RDF/JSON: 24-1.0073102-rdf.json
Turtle: 24-1.0073102-turtle.txt
N-Triples: 24-1.0073102-rdf-ntriples.txt
Original Record: 24-1.0073102-source.json
Full Text

Full Text

Nonlocal Dielectric Response of Water by Matthew Badali B.Sc., The University of Waterloo, 2010 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Physics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2012 c© Matthew Badali 2012 Abstract Water, possessing a permanent dipolar moment, generates a peculiar po- tential of mean force between solvated ions, the investigation of which is the goal of this thesis. Other authors have observed overscreening between two charges in water, leading to attractive regimes in otherwise repulsive potentials, and vice versa. In this thesis, atomistic molecular dynamics sim- ulations of Extended Simple Point Charge (SPC/E) water are used to gen- erate a nonlocal dielectric function that matches with other computational and experimental results from literature. Nonlocality means the consider- ation of the orientation of neighbouring water molecules when finding the displacement field caused by an electric field. Specifically, instead of the displacement field at a given point depending only on the electric field at that point, with nonlocality the displacement field involves an integral of the electric field over all space. To compare to this nonlocal dielectric function derived from simulations of only water, simulations of water with ions are performed. Simulations of two ions in water demonstrate the existence of nonlocality, also in agreement with current scientific understanding. Fur- ther, simulations of three ions show that the effective mean force in more complicated systems is not simply a summation of pair-wise forces but in- stead must account for multi-body effects. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . v 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Local Dielectrics . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Nonlocal Dielectrics . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.5 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Simulations with Ions . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Water Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 12 3.1 Linear Response Theory . . . . . . . . . . . . . . . . . . . . 12 3.2 Bulk Dielectric Constant . . . . . . . . . . . . . . . . . . . . 17 3.3 Deriving the Potential . . . . . . . . . . . . . . . . . . . . . . 18 3.4 Charge Density versus Dipole Correlation . . . . . . . . . . . 22 3.5 Necessity of Nonlocality . . . . . . . . . . . . . . . . . . . . . 22 3.6 Other Considerations . . . . . . . . . . . . . . . . . . . . . . 26 3.7 Three-Body Effects . . . . . . . . . . . . . . . . . . . . . . . 27 4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 iii List of Figures 2.1 Equilateral triangular arrangement of three ions showing the constraint forces output for data analysis. . . . . . . . . . . . 9 2.2 SPC/E water model with three point charges fixed relative to each other and a central LJ potential. . . . . . . . . . . . . . 11 3.1 Dielectric response function χ(k). The overscreening region, defined as χ > 1, is quite large, spanning from k< ≈ 0.3/Å to k> ≈ 16.4/Å. . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Dielectric function (k). The poles are at k< ≈ 0.3/Å and k> ≈ 16.4/Å. . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3 Electric potential of a point charge in nonlocal dielectric wa- ter, in units of kCal/mol per elementary charge. . . . . . . . 20 3.4 Screening factor found from figure 3.3, over whole domain. . 21 3.5 Effective potential energy between two ions solvated in water. 23 3.6 Radial distribution function (RDF) of oxygen atoms from wa- ter molecules about a central ion. . . . . . . . . . . . . . . . 25 3.7 Effective potential energy between charged or neutral parti- cles in water. . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.8 Effective force per ion between three ions, in the case of at- tractive interactions. . . . . . . . . . . . . . . . . . . . . . . 28 3.9 Effective force per ion between three ions, in the case of re- pulsive interactions. . . . . . . . . . . . . . . . . . . . . . . . 28 iv Acknowledgements Thanks to Dr Jörg Rottler and Dr Mona Berciu. v Chapter 1 Introduction 1.1 Motivation When dealing with electrostatics, the first approach of every physicist is to use Maxwell’s equations. Most of the interesting physical situations, how- ever, occur in a medium. In such cases there are the so-called ‘macroscopic’ Maxwell’s equations, and an effective Coulomb potential. The basic assumption behind Maxwell’s macroscopic equations is to treat a dielectric material as infinitely divisible, ie. as a homogeneous continuum. This inevitably fails, as dielectrics are made of matter and ultimately have a discrete microscopic length scale. It is not a significant issue in most solid dielectrics, because the material polarizes with electron density rearrange- ment practically continuously. One important example where it becomes a problem is in the case of water, which has a permanent dipole. In the presence of even a slight electric field a water molecule can re-orient itself, creating a large local electric displacement field. Water is the medium of primary interest in this thesis. This is in part because of its interesting and frankly weird properties, including its hydro- gen bonding, and consequent high surface tension and ordered clustering of liquid water. The lone pair electrons and greater electron affinity of the oxygen atom in H2O give it its permanent dipole moment. Water is also chosen because of its ubiquity. It is used heavily in industrial processes. One could not imagine a biological process occurring except in the presence of water. When doing atomistic computer simulations of biological systems, for instance protein folding problems, the simulation box length tends to scale with the length of the molecule of interest. However, the number of water molecules required to fill the box then goes as the volume, or the cube of the length, and quickly becomes prohibitively computationally expensive. The earliest simulations would simply consider the protein in a vacuum. A better way of dealing with all of these water molecules is to include them implicitly, as factors affecting how the relevant parts of the system interact. It is progress toward the creation of a better implicit water model that is 1 1.2. Local Dielectrics the main motivation for this project. The overall goal is to find a dielectric function that better represents the response of water to solvated charges, so that when given an arrangement of ions and partial charges in solution the electric fields can be accurately predicted, and thus the system can be evolved. Such a dielectric function would remove the need for explicit treatment of the water molecules. Of course, when two charged moieties are close together, on the or- der of Angstroms, each does not see the medium independent of the other, but rather the water orientation depends on the arrangement of the nearby molecules. A typical treatment of dielectric response assumes locality, lin- earity, and slow variation of the electric field in space. The purpose of this thesis is to investigate the first assumption, and to investigate methods of encapsulating nonlocality. Water, with its sub-nanometer-scale dipoles and hydrogen-bonding network, demonstrates such nonlocality. 1.2 Local Dielectrics We start by writing the familiar 1st Maxwell’s equation, Gauss’s Law[1]: o~∇ · ~E = ρ. (1.2.1) ~E is the electric field, ρ is the total charge density, and o is the electric constant, also called the permittivity of free space. The simplest approach to apply this equation in a medium is to pro- ceed from the so-called ‘macroscopic’ Maxwell’s equations. Equation 1.2.1 becomes ~∇ · ~D = ρf , (1.2.2) where ~D is the electric displacement field and ρf is the free charge density. By defining the polarization density field, ~P , and the bound charge density, ρb, which are related by the equation ~∇ · ~P = −ρb, (1.2.3) we can connect the regular, or ‘microscopic’, Maxwell’s equations with their ‘macroscopic’ analogues by ~D = o ~E + ~P (1.2.4) and ρ = ρb + ρf . (1.2.5) 2 1.2. Local Dielectrics The negative sign in equation 1.2.3 arises because the medium polarizes to oppose the applied electric field. One could instead define the induced charge density ρi such that ρi = −ρb. (1.2.6) Under the assumption of locality of the dielectric response we can also connect ~D and ~E by ~D(~r) = o(~r) ~E(~r), (1.2.7) with (~r) being the local dielectric function, assumed to depend only on the magnitude of ~E at that point. We also introduce the response function χ such that in a local dielectric we can write ~D(~r)χ(~r) = ~P (~r). (1.2.8) Finally we relate these two with χ(~r) = 1− 1 (~r) , or (~r) = 1 1− χ(~r) . (1.2.9) For reasons of stability,  cannot be between zero and one.[2, 3] In- tuitively, this is a sensible criterion: if an electric potential from a point charge is typically written as φ = 1 q 4pior , then with 0 ≤  < 1 outside a volume the effective charge would appear greater than q; this argument can be repeated indefinitely and thus the system is not stable. What is meant by (~r) is in fact the longitudinal component of the αβ(~r) tensor, and the same goes for χ. Because there are no magnetic fields in- volved in the systems investigated in this report, the displacement field should on average follow the direction of the electric field, thus we only concern ourselves with the longitudinal component (~r). More generally equation 1.2.7 can be written as Dα(~r) = ∑ β oαβ(~r)Eβ(~r). (1.2.10) Even more generally we can relax locality to write the equation Dα(~r) = ∑ β o ∫ d~r′αβ(~r, ~r′)Eβ(~r′). (1.2.11) This will be dealt with in the next section. 3 1.3. Nonlocal Dielectrics 1.3 Nonlocal Dielectrics Nonlocality means that at a given point the displacement field depends not only on the electric field at that point but also on the electric field elsewhere. In water this is because the orientation of the dipolar molecule depends very strongly on the orientation of its neighbours. This complicates the situa- tion, so to keep things reasonably tractable we start with equation 1.2.11, take only the longitudinal component of αβ, and assume homogeneity and translational invariance: ~D(~r) = o ∫ d~r′(~r − ~r′) ~E(~r′). (1.3.1) After a Fourier transform this can be written as ~D(~k) = o(~k) ~E(~k). (1.3.2) An analogous nonlocal χ also becomes local in Fourier space: ~D(~k)χ(~k) = ~P (~k). (1.3.3) To simplify matters even further we start by treating isotropic systems. Thus initially we will be looking for χ(k) and (k). These two functions remain related by equation 1.2.9: χ(k) = 1− 1 (k) , or (k) = 1 1− χ(k) . As stated before,  cannot be between zero and one. However, it can be less than zero. This phenomenon is called overscreening, wherein a medium overreacts to an applied electric field. This in turn can lead to oscillatory structure in the effective potential and even to attractive interactions in what would be exclusively repulsive ones in vacuum, or vice versa. At low k the dielectric function approaches its bulk value, a positive finite number, and at high k corresponding to the most minute spacial fluctuations  approaches the vacuum value of 1. Consequently if there is going to be an overscreening regime there must also be at least two poles or discontinuities in  as it jumps from greater than 1 to negative and back. This is why χ looks nicer: it starts from less than 1 at k = 0, goes greater than 1 in the overscreened region, and decays to 0 at long k. For stability reasons it can never be negative. 4 1.4. Previous Work 1.4 Previous Work There is a large body of literature on aqueous ion solvation.[4–6] H2O forms solvation shells around the ion with the negative oxygen pointing toward a positive ion or a positive hydrogen pointing toward a negative ion. The easiest way to treat the effective potential of a charged particle in water is to include a dielectric constant, but other authors [7, 8] have shown that this is insufficient on the scale of Angstroms, and I will show the same in this thesis. A nonlocal dielectric function is required. Kornyshev and his collaborators [3, 7, 9, 10] have done some work on de- riving a homogeneous, translationally invariant, and spherically-symmetric nonlocal dielectric function from pure water simulations. Using the Bopp- Janscó-Heinzinger (BJH) [11] water model they found overscreening in the intermediate k regime. The BJH model is an uncommon choice, in par- ticular because it uses partial charges two thirds to three quarters of more popular water models like TIP4P and SPC.[12] From their computational results they derived a dielectric response function that matches closely with the experimental work of Soper.[13, 14] Although Soper went on to say that his data could be off by a factor of 2 or more [15] the qualitative similarities still hold. The dielectric function presented in this thesis is a better match to the updated experimental results. For the investigation of a more general nonlocal dielectric response Hess et al. is the best example from the literature.[8, 16] Using a variety of sim- ple water and ion models they found the effective potential of mean force between two ions. This provided clear evidence of the need for nonlocality when considering water as it responds to electric fields; similar results are presented in this thesis. The authors explored the idea of multi-body effects being important, by increased ion concentration and observing a decrease in the bulk dielectric constant. They did not, however, explicitly look at sim- ulations of three ions, nor the potential of mean force over all wavevectors. 1.5 Thesis Outline The aim of this thesis is to provide a description of the dielectric response of water. Chapter 2 describes the numerical methods applied to get the results listed in Chapter 3. First, in section 3.1, we consider the simplest case, of a homogeneous, translationally invariant medium. Assuming a linear response to an applied electric field, a formulation is derived wherein the response and dielectric 5 1.5. Thesis Outline functions are extracted from a simulation of pure water. This (k) is used to find the effective electric potential generated by a point charge in the medium. Overscreening is observed, and there are large attractive regimes in the repulsive potential. Also in the pure water system, the bulk value of the dielectric is calculated and compared to literature values. Finally, the results of simulations with multiple ions are used to inves- tigate a more accurate effective potential energy between solvated charged particles. Having ions in these simulations breaks both the homogeneity and translational invariance of the system. Equation 1.2.11 is the most ap- propriate dielectric function in such a situation, but instead of extracting (~r, ~r′) the effective electric potential is found directly. These two ion results serve as a check; the linear response theory re- produces oscillatory structure in the effective potential between ions, yet remains a poor match when trying to emulate the real interaction between two ions in solution. What is more, simulations with three ions show that the two ion effective potential does not obey the superposition principle, and is therefore insufficient to describe the physics of charged particles interacting in water. 6 Chapter 2 Methodology 2.1 Molecular Dynamics Molecular dynamics (MD) simulations are carried out using the Large Atomic- Molecular Massively Parallel Simulator (LAMMPS).[19] Most are done on Westgrid, a partner member of Compute Canada.1 MD simulations, and indeed any in silico research, has been gaining popularity and utility along with the computational power made increasingly available over time. MD in particular is useful in that it helps find the trajectory of particles in complex systems when an analytic solution would be hard or impossible to find. A system of particles, in this case atoms alone or in molecules, is initialized in phase space. The chosen model dictates the interaction potential between pairs or groups of particles. By taking the derivative of the potential energy to find the force, Newton’s equations of motion can be solved numerically. The system is advanced for a short time to a new point in phase space, and the process is repeated. A common numerical integrator, used for the purposes of this thesis, is the velocity Verlet algorithm. For velocity Verlet the atomic positions are updated based on their current velocities and accelerations (derived from ~F = −∇U = m~a), their accelerations are updated based on the new positions, and the velocity is updated based on the new and current accelerations, before repeating. Specifically, with a potential U(~x) depending only on particle positions, the dynamic variables are updated as ~x[t+ ∆t] = ~x[t] + ∆t~v[t] + ∆t2~a[t]/2, ~a[t+ ∆t] = − 1 m ∇U (~x[t+ ∆t]) , ~v[t+ ∆t] = ~v[t] + ∆t (~a[t] + ~a[t+ ∆t]) /2. One underlying assumption of these MD simulations, to make it ap- plicable to thermodynamic properties like the dielectric response function, is that of ergodicity. The ergodic hypothesis states that an average of a 1http://www.westgrid.ca 7 2.2. Simulations with Ions thermodynamic quantity over time in a system’s trajectory is equivalent to an ensemble average in that system. The rationale behind this assumption is that each microstate is equally likely and accessible, so an average over time gives greater weighting to regions in phase space where there are more microstates. In a sufficiently long simulation the system visits all relevant regions of phase space, and will spend more time in the regions with higher microstate density. Unless otherwise stated, the parameters outlined in this section of the thesis are applied to all simulations. Simulations are run with periodic boundary conditions in a cubic box of length L = 62.14Å. Eight thou- sand water molecules are placed in the box giving a density of 0.03334 molecules/Å3 = 0.9974 g/cm3. The water molecules start on a lattice and are allowed to equilibrate for 10 ps before data collection commences. Sim- ulations are run for 10 ns with a time step of 2 fs. They evolve in an NVT fashion using a Nosé-Hoover thermostat held at 300.0 K. That is to say, the number of particles N, the volume V, and the temperature T are all kept constant over the course of the simulation. The thermostat does not keep the system at exactly 300K, but rather couples the system to a set of extra degrees of freedom. The extra degrees of freedom follow a logarithmic potential, which leads to a thermodynamic friction term in their equations of motion, which acts to return the entire system to the parameter T (taken here to be 300K). See reference [20] for more details. 2.2 Simulations with Ions For runs with two ions or more the distance between them is held fixed to better sample phase space. This is done using the SHAKE algorithm.[21] For two particles evolving together with the SHAKE algorithm, in each integration step the forces are calculated on each particle using the positions of all the particles, as described in the previous section. But before updating the accelerations, a constraint force is added along the axis between the two atoms, one that exactly cancels out the net force in that direction. In this way the two ions do not move relative to each other, but remain a fixed distance apart, floating like a dumbbell through the simulation. The algorithm is applied similarly to fix triangular arrangements of atoms, for instance in the rigid water molecules used in this thesis. In simulations with two ions this added force is recorded and averaged to find the effective mean force between the ions. The other simulation details are the same as those described above, except that L = 62.14Å and 8 2.2. Simulations with Ions N = 1000. The fixed distance ranges between 2 and 10 Å, at an interval of 0.2 Å. For three ion simulations the parameters are the same and the ions are arranged on the vertices of an equilateral triangle. Again the fixed distances differ by 0.2 Å and this time range between 2 and 8 Å. As for the force being output for analysis, this is done by selecting one atom to be ‘central’. The other two atoms are chosen such that they are symmetric, and so the constraint force between each of them and the central atom is on average the same. Each constraint force is found also with the SHAKE algorithm in a manner similar to that described above. The two constraint forces are averaged, together and over time, to report the effective mean force between the central ion and one other ion in this three ion arrangement. The force to be shown in the results section is then Fconstraint = (Fc1 + Fc2) /2 averaged over the simulation.  'central' ion + - - F c1 F c2 Figure 2.1: Equilateral triangular arrangement of three ions showing the constraint forces output for data analysis. Ions are given a mass of 22.990 or 35.453 atomic mass units (a.m.u.), the mass of sodium or chlorine, with oxygen and hydrogen being 15.9994 and 9 2.3. Water Model 1.00794 a.m.u. respectively. No difference is observed if only the ion mass is changed. The position of each atom is saved every 1000 or 5000 steps of their 10 ns trajectory, to be used for data analysis. 2.3 Water Model The water model chosen is called the extended Simple Point Charge (SPC/E) model.[22, 23] It is selected for its simplicity and frequency of use while still reproducing some global qualities of water such as pressure, bulk dielec- tric constant, and diffusion constant.[23] The SPC/E model is a rigid water molecule with the OH bonds fixed at 1Å and the HOH angle fixed at 109.47◦ via the SHAKE algorithm. The hydrogen atoms are represented as point charges with qH = +0.4238e, and the oxygen atoms are point charges with qO = -0.8476e, where e is the charge on an electron. The static dipole moment of one molecule is 0.489 e Å. Also centred on the oxygen is a Lennard-Jones (LJ) potential with parameters σ = 3.1656 Å and ε = 0.1554 kcal/mol as used in the equation ULJ(r) = 4ε [(σ r )12 − (σ r )6] (2.3.1) for a distance r between two oxygen atoms. See figure 2.2 for reference. The 1/r12 serves as a soft core repulsion between molecules, ensuring that they cannot overlap, and the −1/r6 is akin to the Van der Waals interaction. These same Lennard-Jones parameters are used between ions and the oxygen of water molecules. The ions themselves however have no LJ poten- tial between them. Real ions have repulsive cores so this is perhaps unreal- istic, but it allows for a clearer investigation of the electrostatic component of the interaction at small distances. 10 2.3. Water Model  q O q H q H 1.0Å 1.58Å 109.47º Figure 2.2: SPC/E water model with three point charges fixed relative to each other and a central LJ potential. 11 Chapter 3 Results and Discussion 3.1 Linear Response Theory By looking at the fluctuations in charge density the response to a test point charge can be calculated. It is possible to find the dielectric response of a homogeneous, translationally invariant medium solely from simulations of pure water. This section provides the details of how this can be done, then shows the result extracted from SPC/E simulations. We would like to calculate χ(k), the wavevector-dependent dielectric response of water, and use this to find (k), the dielectric function, using a series of configurations of a distribution of charges (ie. water molecules in a box). For each molecule, and in fact for each atom a, we know ~ra(t) over the course of the simulation. From Maxwell’s equations in a material with nonlocal effects we have equation 1.3.1 ~D(~r) = ∫ d~r′(~r − ~r′)o ~E(~r′) from which a Fourier transform gives ~D(~k) = o(~k) ~E(~k) Similarly equation 1.3.3 tells us that ~P (k) = χ(k) ~D(k) (3.1.1) with k = |~k| since the situation is also spherically symmetric. We start by defining an electric potential φ(~r) such that ~D(~r) = −o∇φ(~r). (3.1.2) We see that in Fourier space it is proportional to ~D(~k): ~D(~k) = −oi~k φ(~k). (3.1.3) 12 3.1. Linear Response Theory This deals with the right-hand side of 3.1.1. For the left-hand side we Fourier transform a combination of equations 1.2.3 and 1.2.6 to get i~k · ~P (~k) = ρi(~k). (3.1.4) We need to find a way to link ρi to the electric potential. It should be noted that the above equations are meant to refer to an equilibrium situation. We will link ρi to φ using simulation data, and so we will need to take a thermodynamic average of ρi. To make this link, start by defining H, the Hamiltonian of the system: H = Ho +H ′. (3.1.5) Ho represents the water-water interactions as well as translation and any other irrelevant physics. H ′ is the perturbation to the system, H ′ = ∫ d~r′ρ(~r′)φ(~r′), (3.1.6) which can be split using equations 1.2.5 and 1.2.6 into H ′ = ∫ d~r′ρf (~r′)φ(~r′)− ∫ d~r′ρi(~r′)φ(~r′). (3.1.7) In the case of simulations of only water, the water molecules form the in- duced charge and there is no free charge. Thus H ′, the perturbation in the system, has the linear response: H ′ = − ∫ d~r′ρi(~r′)φ(~r′). (3.1.8) This is where the thermodynamic average comes in. The induced charge density ρi showing up in equation 3.1.4, unlike in the Hamiltonian, is in fact a charge density ensemble average: ρi(~r) = 〈ρi(~r,~s)〉 = 1 Z ∫ d~s e−β(Ho(~s)+H ′(~s))ρi(~r,~s). (3.1.9) This is a phase space average with d~s being the infinitesimal volume element. The charge density of course depends on the positions of each atom, which along with the momenta is exactly the location ~s in phase space. β = 1kBT and Z is the partition function The brackets 〈A〉 represent an average of A over the perturbed system, and just as Ho is the Hamiltonian of the pure 13 3.1. Linear Response Theory water system without any external fields applied, 〈A〉o is the average of A over the unperturbed system. Now we make two approximations. The first is that the perturbation is small, resulting in little change of the partition function: Z ≈ Zo. This is correct to first order in H ′. The second is that in the Taylor expansion of e−βH′ we are only concerned up to the linear case. This allows us to write ρi(~r) ≈ 1 Zo ∫ d~s e−βHo(~s)(1− βH ′(~s))ρi(~r,~s) = 〈ρi(~r)〉o − β〈H ′ρi(~r)〉o = 〈ρi(~r)〉o + β ∫ d~r′φ(~r′)〈ρi(~r′)ρi(~r)〉o. (3.1.10) In equation 3.1.10 the first term averages out to zero since without a pertur- bation the system is unpolarized and charge neutral. The second term, on which we will focus more, is calculated from molecular dynamics simulations of a box of water molecules in equilibrium. By assuming translational invariance in the unperturbed, pure water system we can write 〈ρi(~r′)ρi(~r)〉o = 〈ρi(~R)ρi(~0)〉o (3.1.11) = 1 V ∫ d~q (2pi)3 ei~q·~R〈ρi(~q)ρi(−~q)〉o (3.1.12) where ~R is defined by ~R = ~r − ~r′. The V is a volume factor picked up by integrating over the position variable freed by translational invariance. Substituting equation 3.1.12 into a Fourier transformed equation 3.1.10 we find ρi(~k) = β ∫ d~re−i~k·~r ∫ d~r′φ(~r′) 1 V ∫ d~q (2pi)3 ei~q·(~r−~r ′)〈ρi(~q)ρi(−~q)〉o = β V ∫ d~q (2pi)3 ∫ d~re−i~r·(~k−~q) ∫ d~r′φ(~r′)e−i~q·~r ′〈ρi(~q)ρi(−~q)〉o = β V ∫ d~r′φ(~r′)e−i~k·~r ′〈ρi(~k)ρi(−~k)〉o = β V φ(~k)〈ρi(~k)ρi(−~k)〉o. (3.1.13) Returning to equation 3.1.1 equipped with equations 3.1.3, 3.1.4 and 14 3.1. Linear Response Theory 3.1.13 we can now write Pα(~k) = ρi(~k) i kα (3.1.14) = β i kαV φ(~k)〈ρi(~k)ρi(−~k)〉o (3.1.15) and χ(~k)Dα(~k) = χ(~k) ( −oi kα φ(~k) ) (3.1.16) = χ(~k) ( −oi kα φ(~k) ) (3.1.17) so β i kαV φ(~k)〈ρi(~k)ρi(−~k)〉o = χ(~k) ( −oi kα φ(~k) ) , (3.1.18) from which it is clear that χ(~k) = β ok2V 〈ρi(~k)ρi(−~k)〉o. (3.1.19) This is consistent with the definition in reference [9]. Using relation 1.2.9 the dielectric function turns out to be (~k) = ( 1− 1 kBT V 〈ρi(~k)ρi(−~k)〉o ok2 )−1 . (3.1.20) Thus the fluctuations in ρi lead to structure in the dielectric response. In a simulation of pure water ρi is calculated by ρi(~r) = 1 V ∑ a qaδ(~r − ~ra), (3.1.21) a sum of δ-functions at atomic positions ~ra weighted by atomic charge. Just as before, V is the volume of the periodic box. In Fourier space this conveniently becomes ρi(~k) = 1 V ∑ a qae −i~k·~ra . (3.1.22) The charge density correlation is calculated to find the response function χ. Since the system of pure water is rotationally invariant, quantities like the charge density or the response function depend only on |~k|, the magnitude of the wavevector. 15 3.1. Linear Response Theory  0  5  10  15  20  25  30  35  40  45  0  2  4  6  8  10  12  14  16  18 ch i (u nit les s) k (1/A) chi1 Figure 3.1: Dielectric response function χ(k). The overscreening region, defined as χ > 1, is quite large, spanning from k< ≈ 0.3/Å to k> ≈ 16.4/Å. Following this prescription, figure 3.1 shows the graph of the response function against wavevector obtained from MD simulations as described in Chapter 2. Both the qualitative shape and the order of magnitude of figure 3.1 are comparable to computational and experimental results in the literature.[7, 14] The calculated response function crosses χ = 1 around k< ≈ 0.3Å−1 and k> ≈ 16.4Å−1. This gives the bounds of the overscreened portion of the dielectric function. The response peaks around k∗ ≈ 3Å−1. In real space this is on the order of 2pi/k∗ ≈ 2Å. It is not clear at this time why this number comes up, as it is smaller than both the molecular diameter of water and the oscillations in the radial distribution function (see Section 3.5). In reference [7], Bopp et al. see a smaller overscreening region, between 1 and 12.5Å−1. Related to this thesis’ results, their χ reaches a maximum of 25, near 3Å−1. This is expected, as their model of water is similar to the SPC/E model but with the partial charges about three quarters of the oxygen and hydrogen charges used in this report.[11] With the perhaps naive assumption that the charge density follows a similar form regardless of the model, then one would expect the χ presented in this thesis to reach 16 3.2. Bulk Dielectric Constant a maximum of 25 times four thirds squared, i.e. 44, and at the same wavevector, which is indeed the case. Bopp et al. show a good match with experimental results [14], though later results from the same lab indicate that the correct data are approximately twice as large.[15] The regime of overscreening is perhaps more obvious in (k). In figure 3.2, which comes from the same simulation as figure 3.1, this is the region between the discontinuities, where (k) < 0. -40 -20  0  20  40  60  80  100  0  2  4  6  8  10  12  14  16  18 ep silo n ( un itle ss) k (1/A) dielectric function Figure 3.2: Dielectric function (k). The poles are at k< ≈ 0.3/Å and k> ≈ 16.4/Å. Equation 1.2.9 is used to create figure 3.2 except for the k = 0 value, which is calculated in the next section and added by hand. At this point the equation fails, and the bulk value of the dielectric must be calculated in a separate manner. There are only two values of (k) before the first pole, so extrapolating from these two alone is suspect, but nevertheless a linear extrapolation gives a bulk dielectric constant around 60, which is not unreasonable for water. 3.2 Bulk Dielectric Constant The SPC/E model has been validated as a reasonable emulation of water, re- producing such bulk quantities as the diffusion constant and latent heat.[23] 17 3.3. Deriving the Potential In order to evaluate the use of the model for electrostatic purposes, it is important to verify that at least the bulk dielectric constant matches with experimental or other simulated values. This same simulation of pure water is used, only the bulk dielectric value is derived from the polarization field, as calculated from the atomic positions and charges. The fluctuation of the net polarization about the mean polarization is used in the equation (k → 0) = bulk = 1 + 〈 ( ~M − 〈 ~M〉o )2〉o 3okBTV (3.2.1) to find the bulk dielectric constant of the simulated SPC/E water.[24, 25] Here ~M is the bulk value of ~P (k), found either as ~P (k → 0) or 1V ∫ d~r ~P (~r). And just as the charge density is a sum over charges, we can write ~M = 1 V ∑ j ~µj , (3.2.2) a sum over each molecular dipole ~µj = ∑ bj qbj~rbj , where bj indexes the atoms in the jth molecule. Literature values from MD simulations vary wildly, but generally place the value of the water bulk dielectric constant near 80 for systems around 300 K.[26] Using equation 3.2.1 the bulk dielectric constant is calculated to be (k = 0) = 75.92. This compares favourably with the experimental literature value of 74.76.[27] 3.3 Deriving the Potential Ultimately the results of this thesis will be used to find interactions between charged particles in solution. Therefore it is useful also to find the electric potential VE(r) at a distance r from a point charge in a nonlocal dielectric. We expect it to look like a modified Coulomb potential. We desire to find a VE such that ~E(~r) = −∇VE(~r) (3.3.1) meaning that ~E(~k) = −i~k VE(~k). (3.3.2) 18 3.3. Deriving the Potential Recalling equation 1.2.2 and applying a Fourier transform gives ρf (~k) = i~k · ~D(~k) = i~k · o(~k) ~E(~k) = i~k · o(~k) ( −i~kVE(~k) ) . This in turn gives VE(~k) = ρf (~k) o(~k)k2 (3.3.3) or VE(~r) = ∫ d~k (2pi)3 e−i~k·~r ρf (~k) o(~k)k2 (3.3.4) The electric potential can thus be calculated from any given free charge density ρf (~r) now that (k) is known. Taking equation 3.3.4 and assuming spherical symmetry we can write VE(r) = ∫ d3k (2pi)3 e−ik·rVE(k) = 1 (2pi)3 ∫ ∞ 0 k2dk ∫ 1 −1 d cos θ ∫ 2pi 0 dϕe−ikr cos θVE(k) = 1 (2pi)2 ∫ ∞ 0 k2dk ∫ 1 −1 d cos θ e−ikr cos θVE(k) = 2 (2pi)2 ∫ ∞ 0 k2dk VE(k) sin(kr) kr = 1 2pi2or ∫ ∞ 0 dk ρf (k) sin(kr) (k)k . (3.3.5) Consider the simplest distribution for free charge, that of a lone charge at the origin. For a charge q at the origin we get ρf (r) = qδ(r), which implies ρf (k) = q. Then VE(r) = q 2pi2or ∫ ∞ 0 dk sin(kr) k 1 (k) . (3.3.6) It is worth noting that ∫∞ 0 dk sin(kr) k = pi 2 sgn(r) so for (k) = 1 we recover the usual Coulomb potential. For computational ease equation 1.2.9 is used to rewrite equation 3.3.6 as VE(r) = q 2pi2or ( pi 2 − ∫ ∞ 0 dk sin(kr) k χ(k) ) . (3.3.7) 19 3.3. Deriving the Potential By inputting the dielectric response function χ(k) obtained from simulation we find the effective electric potential generated by a point charge, shown in figure 3.3. -200 -150 -100 -50  0  50  100  150  4  6  8  10  12  14  16 ele ctr ic p ote nti al (kC al/ mo l)/e r (A) effective potential Figure 3.3: Electric potential of a point charge in nonlocal dielectric water, in units of kCal/mol per elementary charge. Figure 3.4 allows for a comparison with a local dielectric. If we define a screening factor S(r) such that S(r) = bulk4pior q VE(r) (3.3.8) then a local medium with only a bulk dielectric constant would show a constant plot at S(r) = 1. Instead we see the screening factor vary orders of magnitude above and below, only settling toward the bulk value of 1 at long distances. When its magnitude is above bulk it indicates that the potential due to a point charge is even greater than the case of the same charge in a vacuum. Intuition states that the water acts to screen the charge; it is indeed overscreening if the medium ends up reinforcing the field to be greater than in the vacuum case. As seen in figure 3.3 there is indeed structure in the effective poten- tial. The observed oscillatory behaviour is consistent with an overscreening medium. Half of a period looks to be about 1Å, so the full period is on the 20 3.3. Deriving the Potential -1000 -800 -600 -400 -200  0  200  400  600  0  2  4  6  8  10  12  14  16 scr ee nin g f act or (un itle ss) r (A) screening factor1 Figure 3.4: Screening factor found from figure 3.3, over whole domain. order of r∗ ≈ 2Å. This corresponds to a wavevector of k∗ = 2pir∗ ≈ 3Å−1. Of course this is the wavevector expected to be most manifest, as it is the strongest contribution of χ(k). All of the wavevectors in the overscreening regime contribute, but the wavevector of maximal response contributes the most. Overscreening leads to oscillatory structure in the effective potential. The oscillations, if large enough, lead to the unintuitive result that there are some distances that the electric potential caused by a positive point charge is actually attractive to a positive test charge. It is obvious that if such a phenomenon really occurs it would lend an important correction to any MD simulations of implicitly solvated charges. There is structure in the effective electric potential. At long distances the limiting behaviour returns to a regular 1/r screened Coulomb potential. The structure in the effective potential, however, is present at all shorter distances, even going to a radius smaller than the size of a water molecule. This is because a point charge is used in the preceding derivation. A point charge is an unphysical representation of a charged particle in solution as it allows the solvent molecules to get arbitrarily close to the ion, and therefore the structure is present arbitrarily close to r = 0. 21 3.4. Charge Density versus Dipole Correlation 3.4 Charge Density versus Dipole Correlation It is the permanent dipole of water that lends itself to overscreening. There- fore one might conclude that a sensible approach to calculating a response function would include dipoles. Indeed it is possible to go through the derivation shown in the linear response theory section using the polariza- tion density field rather than charge density. The linear response of equation 3.1.8 could instead written as H ′ = ∫ d~r′ ~P (~r′) · ~D(~r′) and the derivation fol- lows analogously, only using ~P and ~D in place of ρi and φ. Equation 3.1.19 becomes χ(~k) = 1 o kBT V 〈~P (~k) · ~P (−~k)〉o. (3.4.1) However, this fails to accurately predict even the qualitative results of the dielectric function, and instead stays in an overscreened region in the k →∞ limit. Computationally, ~P (~r) is calculated as ~P (~r) = ∑ j ~µjδ(~rj)/V , a sum over each molecular dipole ~µj = ∑ b qb~rb. Placing the molecular dipole µj at only one position rj makes it a point dipole and neglects the finite structure of a water molecule. Consider equation 3.1.22. By rewriting the sum over atoms as a sum over molecules and atoms within the molecule we can write ρi(~k) = 1 V ∑ j,b qb exp{−i~k · (~rj + δ~rj,b)} ≈ 1 V ∑ j,b qb exp{−i~k · ~rj} ∑ n ( −i~k · δ~rj,b )n n! ≈M‖ +Q‖ +O‖ + . . . , which is a multipole expansion, with the longitudinal component of the dipole being only the first term in the expansion. It is therefore wrong to only include the correlations of the dipole moments when calculating (~k); the charge density is needed to properly account for all higher moment correlations that occur in the system. Reference [9] contains a further in depth discussion of this matter. 3.5 Necessity of Nonlocality We now consider a more general scenario, for a system which is neither homogeneous nor translationally invariant. To investigate this, simulations are run using a single pair of ions solvated in SPC/E water. The paired 22 3.5. Necessity of Nonlocality combinations of positive-positive, negative-negative, and positive-negative ions are kept at fixed distances between 2 and 10 Angstroms. We expect to find structure on the order of Angstroms in the effective Coulomb potential between two ions in water, based on the previous results and the typical size of a solvation shell. The force, Fconstraint, required to keep the ions fixed relative to each other is recorded and averaged over the run. It is integrated as U(r) = ∫ r rmax { 〈Fconstraint〉+ 2kBT x } dx+ q1q2 4piobulkrmax , r < rmax (3.5.1) to find the potential of mean force U(r) between the two ions.[8] The in- tegral runs up to rmax = 10Å at which point the curve begins to looks like the regular Coulomb energy with the dielectric constant of bulk water, (k = 0) = 75.92. The second term in the integral accounts for entropic contribution from more phase space being sampled by the farther separated ions. Figure 3.5 shows the results of the simulations. -4 -3 -2 -1  0  1  2  3  4  2  3  4  5  6  7  8  9  10 po ten tia l (k Ca l/m ol) ion separation (Angstrom) effective -- potentialeffective ++ potentialunmodified repulsive potentialeffective +- potentialunmodified attractive potential Figure 3.5: Effective potential energy between two ions solvated in water. It should be noted that in simulations with one positive and one nega- tive ion the constraint forces on the two were equal and opposite. In gen- eral, the above results compare favourably with simulations from another water model.[8] Both these results and the literature show the oscillatory 23 3.5. Necessity of Nonlocality behaviour. Furthermore, both the magnitude of the oscillations and the period of oscillation are comparable between the literature and this thesis. Importantly, there is an observed difference between the positive-positive and negative-negative pairings. What makes the results of figure 3.5 most interesting is the oscillatory structure of the effective potential. Just as with the application of linear response theory, these ‘true’ computational results show that the medium tends to overscreen the external charges. Furthermore, there are again the unintuitive attractive portions of the potential energy, even between like- charged ions. However, this is where the similarities end. While the linear response theory results of section 3.3 contains the same essential qualities as the computational result shown in figure 3.5, it fails in a few key regards. As previously discussed, it is unreasonable to treat an ion as a point charge; real ions possess a volume and displace the medium and other ions within that volume. While the oscillations are still observed, the number and positions of, for example, the attractive portions matches poorly with the two ion simulations. What is more, the magnitude of the linear response results is two orders of magnitude greater than these ‘real’ results. Finally, the results of section 3.3 predict that there should be the same potential for both positive-positive and negative-negative, and that the positive-negative combination should also have the same magnitude. It is obvious that the rapidly varying displacement field generated by an ion in water needs special treatment and must be represented by a nonlocal dielectric function. Even in the relatively simple case of two ions in solution, the effect of inhomogeneity cannot be neglected. Therefore a general (~r, ~r′) is needed to possibly describe the interaction of solvated charges in this dielectric medium. From these results, it is also clear that any single (~r, ~r′) will not account for all the details in the interaction. Notice the asymmetry between the positive-positive and negative-negative interactions. This discrepancy arises because water is not a point dipole; the positive partial charge in a water molecule is divided between the two hydrogen atoms, which allows for higher order multipole interactions between molecules. The effect this has on how water structures itself around ions can be unintuitive. Figure 3.6 is generated from simulations with only one ion per box, and shows the radial distribution of the oxygen atoms from water around the ion. The radial distribution function (RDF) of oxygen about the ion gives a measure of the amount of water at a given radial distance away. The RDF, or pair correlation function, gives the probability of finding a particle (in this case, oxygen) at a given distance from the reference particle (in this 24 3.5. Necessity of Nonlocality  0  1  2  3  4  5  6  7  8  2  3  4  5  6  7  8  9  10 g(r ) distance (Angstrom) RDF about positive ionRDF about negative ion Figure 3.6: Radial distribution function (RDF) of oxygen atoms from water molecules about a central ion. case, the ion), normalized by the probability of finding at that distance a particle of an ideal gas of the same density. Peaks in the RDF indicate shells of high density of oxygen about the ion. In both positive and negative ion cases there are multiple solvation shells, i.e. layers of water at various distances from the solvated ion. For a positive ion, the negative end of the water dipole, that is to say the oxygen atom, points towards the charge, as one would expect. But for a solvated negative ion, one of the two partially positive hydrogen atoms per molecule are drawn toward the charge, forming a hydrogen bond, and bringing the whole water molecule, including the oxygen, closer to the negative ion than in the positive case, such that the first peak in the ion-oxygen distribution function occurs closer. This phenomenon has been observed in other water models and is not unfamiliar.[4, 6] To obtain a complete description of how a medium responds to solvated charges, the dielectric function must account for the sign of charges causing the dielectric response. As an example, consider the longitudinal compo- nents of equation 1.2.11, ~D(~r) = o ∫ d~r′(~r, ~r′) ~E(~r′), with an electric field caused by a single solvated ion. The difference between a positive and neg- ative ion is solely a negative sign in ~E, which carries through the integral. But it is obvious from the results in this section that the field in water is 25 3.6. Other Considerations different by more than just a sign. Any progress from this point should be mindful of this detail. 3.6 Other Considerations For the most part, the results shown in figure 3.5 are due to water-water interactions of the solvation shells about ions, and the ions themselves at- tracting or repelling the oxygen and hydrogen atoms that make up the water molecules. But electrostatics of the ions are not the only factors. Figure 3.7 shows data from simulations of one or two ‘neutral’ ions. Simulations are run with the same parameters as those with two ions; the only difference is that one or both of the ‘ions’ in these simulations have zero charge. They still interact with the oxygen of water via a Lennard-Jones interaction, and they do not interact with each other.     	     




 Figure 3.7: Effective potential energy between charged or neutral particles in water. In such cases there are still clusters of water solvating the charged or neutral ‘ions’, which interact in their non-trivial ways. Consider the neutral- neutral curve in figure 3.7. Placing one neutral ‘ion’ in solution disrupts the hydrogen-bonding network of pure water, at least within the radius of exclusion that comes from the LJ potential of the ‘ion’. With two such ‘ions’ 26 3.7. Three-Body Effects far away from each other, this would be twice the energy cost. But as the two ‘ions’ move closer to each other, especially as their exclusion volumes begin to overlap, the pair together take up less volume, and break the hydrogen- bonding network less. Thus it is energetically favourable for the ‘ions’ to be closer, and indeed it is seen that the potential curves downward at smaller distances. On the other hand, a charged ion in solution is energetically favourable; while the hydrogen-bonding network of water still gets disrupted, this is because the water instead orients itself with respect to the ion. This is es- pecially true for negative ions, with which the water forms hydrogen-bonds. Therefore one would imagine that putting a neutral ‘ion’ close to one of these solvated ions would disrupt the water shell and therefore be energeti- cally unfavourable, especially for a negative ion. This is exactly what is seen at smaller r in figure 3.7. There is also slight oscillatory behaviour seen in all of the potentials, but this probably has to do with overlap of solvation shells, and either way is small enough as to not explain the oscillations seen in figure 3.3. 3.7 Three-Body Effects How well do the results of Section 3.5 fare in more complicated systems? Specifically, the question is whether the two-ion potential described is su- perimposable, or are there multi-body effects when there are more than two ions. Is the potential in an arbitrary situation simply a sum of pairwise interactions, or are there higher order terms? Similar to the case with two ions, simulations are run with three ions, held fixed relative to each other as the vertices of an equilateral triangle. Figures 3.8 and 3.9 show how some particular triangular arrangements com- pare to the two ion forces. One ion of the triangle is chosen as the ‘central’ atom, and the constraint force of the other two symmetric ions is averaged. The difference between the last two lines in figure 3.8, the green and purple ones, is that the first (green) shows the average force on the central negative ion from each of two positive ions, and the second (purple) shows the reverse (as in figure 2.1). If the potential of mean force between solvated ions could be decomposed as a sum of pairwise interactions, figures 3.8 and 3.9 would show the same two ion and three ion plots. It is evident from the presented data that elec- trostatics from ions in water are not merely superimposable. The solvation shells must overlap in a complicating fashion. But even this is not consistent 27 3.7. Three-Body Effects -4 -3 -2 -1  0  1  2  2  3  4  5  6  7  8  9  10 for ce (kC al/ mo l-A ng str om ) ion separation (Angstrom) effective +- force from two ionsunmodified attractive forceeffective +- force from three ionseffective -+ force from three ions Figure 3.8: Effective force per ion between three ions, in the case of attractive interactions. -2 -1  0  1  2  3  4  2  3  4  5  6  7  8  9  10 for ce (kC al/ mo l-A ng str om ) ion separation (Angstrom) effective ++ force from two ionseffective -- force from two ionsunmodified repulsive forceeffective ++ force from three ionseffective -- force from three ions Figure 3.9: Effective force per ion between three ions, in the case of repulsive interactions. 28 3.7. Three-Body Effects in all cases. In figure 3.8 it is seen that the force felt between solvated ions of opposite charge is exaggerated in the presence of a third ion. It still shows the same oscillatory behaviour with repulsive regimes but the magnitude of effective force is much greater. Similarly in figure 3.9, the force between two positive ions is compounded. This contrasts with the negative ion case, where the interaction between two negative ions is hardly perturbed when a third negative ion is present. Whatever the reason, the effective interaction between two ions in solu- tion is confounded when a third ion is present. While this does not discount the results presented in figure 3.5, it does indicate that those results are in- sufficient to describe anything more than the simplest interactions. Pairwise interactions are not enough to describe the physics of a complicated system. 29 Chapter 4 Conclusion By assuming homogeneity and translational invariance, a linear response theory applied to a simulation of pure SPC/E water uses charge density correlations to produce a dielectric function (k), seen in figure 3.2. This dielectric function has the correct limits at bulk and infinite k and shows an overscreening regime between k< ≈ 0.3Å−1 and k< ≈ 16.4Å−1. Experiment and theory literature support such a dielectric function.[7, 14, 15] From (k) the effective electric potential due to a point charge in solution is calculated and presented in figure 3.3. There are counter-intuitive regions of attraction in otherwise repulsive potentials, and oscillatory structure is seen. This is a good starting point but does not accurately reflect the field generated by solvated charges. In particular, the magnitude of the calculated potential from a point charge is greater than the charge would generate in a vacuum. Atomistic simulations of two ions in solution give a more realistic poten- tial energy (figure 3.5). This too shows oscillatory behaviour modifying a Coulomb interaction in a uniform polarizable medium. There is an asymme- try between two negative ions and two positive ions interacting, suggesting that even the charge must be considered when modeling the potential be- tween two ions. Of course, the level of symmetry of the pure water case has also been reduces, such that a homogeneous, translationally invariant dielectric function is insufficient to model the effective interactions. Simulations of three ions in solution show that two ion potentials cannot simply be added together (figures 3.9 and 3.8). This failure of superposition of potentials is a crucial result. Complicated systems cannot be decomposed to pairwise interactions, and any future work must take this into account. The next step would be to move beyond the assumption of translational invariance. In the situation with even one ion present it acts as an inho- mogeneity. Then instead of (~r − ~r′) one would have to use the generalized (~r, ~r′), or more likely (~r,~k′). In this case one of the variables, say ~r, would be used to represent where in space one is regarding relative to the ion, and would depend only on the distance from the ion, since spherical symmetry is retained in a single ion system. The second variable over which one would 30 Chapter 4. Conclusion integrate (or multiply, if the Fourier transform has been applied to make the relation local in reciprocal space) would then depend on the angle away from the ion-~r axis. Introducing another ion breaks the symmetry further, such that the full (~r,~k′) depending on two three-dimensional vectors would need to be used. Still, the approach used here would work almost as well in these more complicated scenarios. One thing is clear: for simulations involving charged particles interacting on the lengthscale of Angstroms in water, a nonlocal dielectric response is necessary to account for the interesting physics. Due to its highly polar na- ture water cannot be treated by conventional means to simulate electrostatic interactions. 31 Bibliography [1] J. D. Jackson, Classical Electrodynamics. New Jersey, USA: John Wiley & Sons, Inc., 3 ed., 1991. [2] O. Dolgov, D. Kirzhnits, and E. Maksimov, “On an admissible sign of the static dielectric function of matter,” Rev. Mod. Phys., vol. 53, no. 1, pp. 81–94, 1981. [3] A. Kornyshev, S. Leikin, and G. Sutmann, “”Overscreening” in a polar liquid as a result of coupling between polarization and density fluctua- tions,” Electrochimica Acta, vol. 42, no. 5, pp. 849–865, 1997. [4] L. X. Dang, J. E. Rice, J. Caldwell, and P. A. Kollman, “Ion solvation in polarizable water: molecular dynamics simulations,” J. Am. Chem. Soc., vol. 113, pp. 2481–2486, 1991. [5] S. Koneshan, J. C. Rasaiah, R. Lynden-Bell, and S. Lee, “Solvent struc- ture, dynamics, and ion mobility in aqueous solutions at 25c,” J. Phys. Chem., vol. 102, pp. 4193–4204, 1998. [6] J. C. Rasaiah and R. Lynden-Bell, “Computer simulation studies of the structure and dynamics of ions and non-polar solutes in water,” Phil. Trans. R. Soc. Lond. A, vol. 359, pp. 1545–1574, 2001. [7] P. A. Bopp, A. A. Kornyshev, and G. Sutmann, “Static nonlocal dielec- tric function of liquid water,” Phys. Rev. Lett., vol. 76, no. 8, pp. 1280– 1283, 1996. [8] B. Hess, C. Holm, and N. van der Vegt, “Modeling multibody effects in ionic solutions with a concentration dependent dielectric permittivity,” Phys. Rev. Lett., vol. 96, p. 147801, 2006. [9] P. A. Bopp, A. A. Kornyshev, and G. Sutmann, “Frequency and wave- vector dependent dielectric function of water: Collective modes and relaxation spectra,” J. Chem. Phys., vol. 109, pp. 1939–1958, 1998. 32 Bibliography [10] A. A. Kornyshev and G. Sutmann, “Nonlocal nonlinear static dielec- tric response of polar liquids,” Journal of Electroanalytical Chemistry, vol. 450, pp. 143–156, 1998. [11] P. Bopp, G. Jancso, and K. Heinzinger, “An improved potential for non- rigid water molecules in the liquid phase,” Chem. Phys. Lett., vol. 98, pp. 129–133, 1983. [12] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, “Comparison of simple potential functions for simulating liquid,” J. Chem. Phys., vol. 79, p. 926, 1983. [13] A. Soper, C. Andreani, and M. Nardone, “Reconstruction of the ori- entational pair-correlation function from neutron-diffraction data: The case of liquid hydrogen iodide,” Phys. Rev. E, vol. 47, no. 4. [14] A. Soper, “Orientational correlation function for molecular liquids: The case of liquid water,” J. Chem. Phys., vol. 101, no. 8, pp. 6888–6901, 1994. [15] A. Soper, “The radial distribution functions of water and ice from 220 to 673 k and at pressures up to 400 mpa,” Chem. Phys., vol. 258, pp. 121–137, 2000. [16] B. Hess, C. Holm, and N. van der Vegt, “Osmotic coefficients of atom- istic nacl (aq) force fields,” J. Chem. Phys., vol. 124, p. 164509, 2006. [17] J. Rottler and B. Krayenhoff, “Numerical studies of nonlocal electro- static effects on the sub-nanoscale,” J. Phys.: Condens. Matter, vol. 21, p. 255901, 2009. [18] J. Rottler and A. Maggs, “Long-ranged electrostatics from local algo- rithms,” Soft Matter, vol. 7, pp. 3260–3267, 2011. [19] S. Plimpton, “Fast parallel algorithms for short-range molecular dy- namics,” J. Comp. Phys., vol. 117, pp. 1–19, 1995. [20] W. G. Hoover, “Canonical dynamics: Equilibrium phase-space distri- butions,” Phys. Rev. A, vol. 31, no. 3, pp. 1695–1697, 1985. [21] J.-P. Ryckaert, G. Ciccotti, and H. J. Berendsen, “Numerical integra- tion of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes,” J. Comp. Phys., vol. 23, pp. 327–341, 1977. 33 Bibliography [22] F. Stillinger and A. Rahman, “Revised central force potentials for wa- ter,” J. Chem. Phys., vol. 68, no. 2, pp. 666–670, 1978. [23] H. Berendsen, J. Grigera, and T. Straatsma, “The missing term in effective pair potentials,” J. Phys. Chem., vol. 91, pp. 6269–6271, 1987. [24] S. de Leeuw, J. Perram, and E. Smith, “Simulation of electrostatic systems in periodic boundary conditions. i. lattice sums and dielectric constants,” Proc. R. Soc. London, Ser. A, vol. 373, no. 1752, pp. 27–56, 1980. [25] M. Neumann, “Dipole moment fluctuation formulas in computer simu- lations of polar systems,” Molecular Physics, vol. 50, no. 4, pp. 841–858, 1983. [26] M. Rami Reddy and M. Berkowitz, “The dielectric constant of spc/e water,” Chem. Phys. Lett., vol. 155, no. 2, pp. 173–176, 1989. [27] C. Wohlfarth, CRC Handbook of Chemistry and Physics. Internet Ver- sion, 92 ed., 2012. 34


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items